Reference documentation for deal.II version 8.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
sparse_vanka.templates.h
1 // ---------------------------------------------------------------------
2 // @f$Id: sparse_vanka.templates.h 30040 2013-07-18 17:06:48Z maier @f$
3 //
4 // Copyright (C) 1999 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__sparse_vanka_templates_h
18 #define __deal2__sparse_vanka_templates_h
19 
20 
21 #include <deal.II/base/memory_consumption.h>
22 #include <deal.II/base/thread_management.h>
23 #include <deal.II/lac/sparse_vanka.h>
24 #include <deal.II/lac/full_matrix.h>
25 #include <deal.II/lac/sparse_matrix.h>
26 #include <deal.II/lac/vector.h>
27 
28 #include <algorithm>
29 #include <map>
30 
32 
33 template<typename number>
35  const std::vector<bool> &selected,
36  const bool conserve_mem,
37  const unsigned int n_threads)
38  :
39  matrix (&M, typeid(*this).name()),
40  conserve_mem (conserve_mem),
41  selected (selected),
42  n_threads (n_threads),
43  inverses (M.m(), 0)
44 {
45  Assert (M.m() == M.n(), ExcNotQuadratic ());
46  Assert (M.m() == selected.size(), ExcDimensionMismatch(M.m(), selected.size()));
47 
48  if (conserve_mem == false)
50 }
51 
52 
53 template<typename number>
55 {
56  typename std::vector<SmartPointer<FullMatrix<float>,SparseVanka<number> > >::iterator i;
57  for (i=inverses.begin(); i!=inverses.end(); ++i)
58  {
59  FullMatrix<float> *p = *i;
60  *i = 0;
61  if (p != 0) delete p;
62  }
63 }
64 
65 
66 template <typename number>
67 void
69 {
70 #ifndef DEAL_II_WITH_THREADS
71  compute_inverses (0, matrix->m());
72 #else
73  const size_type n_inverses = std::count (selected.begin(),
74  selected.end(),
75  true);
76 
77  const size_type n_inverses_per_thread = std::max(n_inverses / n_threads,
78  static_cast<size_type> (1U));
79 
80  // set up start and end index
81  // for each of the
82  // threads. note that we have
83  // to work somewhat to get this
84  // appropriate, since the
85  // indices for which inverses
86  // have to be computed may not
87  // be evenly distributed in the
88  // vector. as an extreme
89  // example consider numbering
90  // of DoFs by component, then
91  // all indices for which we
92  // have to do work will be
93  // consecutive, with other
94  // consecutive regions where we
95  // do not have to do something
96  std::vector<std::pair<size_type, unsigned int> > blocking (n_threads);
97 
98  unsigned int c = 0;
99  unsigned int thread = 0;
100  blocking[0].first = 0;
101 
102  for (size_type i=0; (i<matrix->m()) && (thread+1<n_threads); ++i)
103  {
104  if (selected[i] == true)
105  ++c;
106  if (c == n_inverses_per_thread)
107  {
108  blocking[thread].second = i;
109  blocking[thread+1].first = i;
110  ++thread;
111 
112  c = 0;
113  };
114  };
115  blocking[n_threads-1].second = matrix->m();
116 
117  typedef void (SparseVanka<number>::*FunPtr)(const size_type,
118  const size_type);
119  const FunPtr fun_ptr = &SparseVanka<number>::compute_inverses;
120 
121  // Now spawn the threads
122  Threads::ThreadGroup<> threads;
123  for (unsigned int i=0; i<n_threads; ++i)
124  threads += Threads::new_thread (fun_ptr, *this,
125  blocking[i].first,
126  blocking[i].second);
127  threads.join_all ();
128 #endif
129 }
130 
131 
132 template <typename number>
133 void
135  const size_type end)
136 {
137  // set-up the vector that will be used
138  // by the functions which we call
139  // below.
140  std::vector<size_type> local_indices;
141 
142  // traverse all rows of the matrix
143  // which are selected
144  for (size_type row=begin; row<end; ++row)
145  if (selected[row] == true)
146  compute_inverse (row, local_indices);
147 }
148 
149 
150 
151 template <typename number>
152 void
154  std::vector<size_type> &local_indices)
155 {
156  // first define an alias to the sparsity
157  // pattern of the matrix, since this
158  // will be used quite often
159  const SparsityPattern &structure
160  = matrix->get_sparsity_pattern();
161 
162  const size_type row_length = structure.row_length(row);
163 
164  inverses[row] = new FullMatrix<float> (row_length, row_length);
165 
166  // collect the dofs that couple
167  // with @p row
168  local_indices.resize (row_length);
169  for (size_type i=0; i<row_length; ++i)
170  local_indices[i] = structure.column_number(row, i);
171 
172  // Build local matrix
173  inverses[row]->extract_submatrix_from (*matrix,
174  local_indices,
175  local_indices);
176 
177  // Compute inverse
178  inverses[row]->gauss_jordan();
179 }
180 
181 
182 template<typename number>
183 template<typename number2>
184 void
186  const Vector<number2> &src) const
187 {
188  // first set output vector to zero
189  dst = 0;
190  // then pass on to the function
191  // that actually does the work
192  apply_preconditioner (dst, src);
193 }
194 
195 
196 
197 template<typename number>
198 template<typename number2>
199 void
201  const Vector<number2> &src,
202  const std::vector<bool> *const dof_mask) const
203 {
204  Assert (dst.size() == src.size(),
205  ExcDimensionMismatch(dst.size(), src.size()));
206  Assert (dst.size() == matrix->m(),
207  ExcDimensionMismatch(dst.size(), src.size()));
208 
209  // first define an alias to the sparsity
210  // pattern of the matrix, since this
211  // will be used quite often
212  const SparsityPattern &structure
213  = matrix->get_sparsity_pattern();
214 
215 
216  // store whether we shall work on
217  // the whole matrix, or only on
218  // blocks. this variable is used to
219  // optimize access to vectors a
220  // little bit.
221  const bool range_is_restricted = (dof_mask != 0);
222 
223  // space to be used for local
224  // systems. allocate as much memory
225  // as is the maximum. this
226  // eliminates the need to
227  // re-allocate memory inside the
228  // loop.
229  FullMatrix<float> local_matrix (structure.max_entries_per_row(),
230  structure.max_entries_per_row());
231  Vector<float> b (structure.max_entries_per_row());
232  Vector<float> x (structure.max_entries_per_row());
233 
234  std::map<size_type, size_type> local_index;
235 
236  // traverse all rows of the matrix
237  // which are selected
238  const size_type n = matrix->m();
239  for (size_type row=0; row<n; ++row)
240  if ((selected[row] == true) &&
241  ((range_is_restricted == false) || ((*dof_mask)[row] == true)))
242  {
243  const size_type row_length = structure.row_length(row);
244 
245  // if we don't store the
246  // inverse matrices, then alias
247  // the entry in the global
248  // vector to the local matrix
249  // to be used
250  if (conserve_mem == true)
251  {
252  inverses[row] = &local_matrix;
253  inverses[row]->reinit (row_length, row_length);
254  };
255 
256  b.reinit (row_length);
257  x.reinit (row_length);
258  // mapping between:
259  // 1 column number of all
260  // entries in this row, and
261  // 2 the position within this
262  // row (as stored in the
263  // SparsityPattern object
264  //
265  // since we do not explicitly
266  // consider nonsysmmetric sparsity
267  // patterns, the first element
268  // of each entry simply denotes
269  // all degrees of freedom that
270  // couple with @p row.
271  local_index.clear ();
272  for (size_type i=0; i<row_length; ++i)
273  local_index.insert(std::pair<size_type, size_type>
274  (structure.column_number(row, i), i));
275 
276  // Build local matrix and rhs
277  for (std::map<size_type, size_type>::const_iterator is=local_index.begin();
278  is!=local_index.end(); ++is)
279  {
280  // irow loops over all DoFs that
281  // couple with the present DoF
282  const size_type irow = is->first;
283  // index of DoF irow in the matrix
284  // row corresponding to DoF @p row.
285  // runs between 0 and row_length
286  const size_type i = is->second;
287 
288  // copy rhs
289  b(i) = src(irow);
290 
291  // for all the DoFs that irow
292  // couples with
293  // number of DoFs coupling to
294  // irow (including irow itself)
295  for (typename SparseMatrix<number>::const_iterator p=matrix->begin(row);
296  p != matrix->end(row); ++p)
297  {
298  // find out whether this DoF
299  // (that couples with @p irow,
300  // which itself couples with
301  // @p row) also couples with
302  // @p row.
303  const std::map<size_type, size_type>::const_iterator js
304  = local_index.find(p->column());
305  // if not, then still use
306  // this dof to modify the rhs
307  //
308  // note that if so, we already
309  // have copied the entry above
310  if (js == local_index.end())
311  {
312  if (!range_is_restricted ||
313  ((*dof_mask)[p->column()] == true))
314  b(i) -= p->value() * dst(p->column());
315  }
316  else
317  // if so, then build the
318  // matrix out of it
319  if (conserve_mem == true)
320  (*inverses[row])(i,js->second) = p->value();
321  }
322  }
323 
324  // Compute new values
325  if (conserve_mem == true)
326  inverses[row]->gauss_jordan();
327 
328  // apply preconditioner
329  inverses[row]->vmult(x,b);
330 
331  // Distribute new values
332  for (std::map<size_type, size_type>::const_iterator is=local_index.begin();
333  is!=local_index.end(); ++is)
334  {
335  const size_type irow = is->first;
336  const size_type i = is->second;
337 
338  if (!range_is_restricted ||
339  ((*dof_mask)[irow] == true))
340  dst(irow) = x(i);
341  // do nothing if not in
342  // the range
343  }
344 
345  // if we don't store the
346  // inverses, then unalias the
347  // local matrix
348  if (conserve_mem == true)
349  inverses[row] = 0;
350  }
351 }
352 
353 
354 
355 template <typename number>
356 std::size_t
358 {
359  std::size_t mem = (sizeof(*this) +
361  for (size_type i=0; i<inverses.size(); ++i)
362  mem += MemoryConsumption::memory_consumption (*inverses[i]);
363 
364  return mem;
365 }
366 
367 
368 
369 
370 template <typename number>
372  const std::vector<bool> &selected,
373  const unsigned int n_blocks,
374  const BlockingStrategy blocking_strategy,
375  const bool conserve_memory,
376  const unsigned int n_threads)
377  :
378  SparseVanka<number> (M, selected, conserve_memory, n_threads),
379  n_blocks (n_blocks),
380  dof_masks (n_blocks,
381  std::vector<bool>(M.m(), false))
382 {
383  compute_dof_masks (M, selected, blocking_strategy);
384 }
385 
386 
387 template <typename number>
388 void
390  const std::vector<bool> &selected,
391  const BlockingStrategy blocking_strategy)
392 {
393  Assert (n_blocks > 0, ExcInternalError());
394 
395  const size_type n_inverses = std::count (selected.begin(),
396  selected.end(),
397  true);
398 
399  const size_type n_inverses_per_block = std::max(n_inverses / n_blocks,
400  static_cast<size_type> (1U));
401 
402  // precompute the splitting points
403  std::vector<std::pair<size_type, size_type> > intervals (n_blocks);
404 
405  // set up start and end index for
406  // each of the blocks. note that
407  // we have to work somewhat to get
408  // this appropriate, since the
409  // indices for which inverses have
410  // to be computed may not be evenly
411  // distributed in the vector. as an
412  // extreme example consider
413  // numbering of DoFs by component,
414  // then all indices for which we
415  // have to do work will be
416  // consecutive, with other
417  // consecutive regions where we do
418  // not have to do something
419  if (true)
420  {
421  unsigned int c = 0;
422  unsigned int block = 0;
423  intervals[0].first = 0;
424 
425  for (size_type i=0; (i<M.m()) && (block+1<n_blocks); ++i)
426  {
427  if (selected[i] == true)
428  ++c;
429  if (c == n_inverses_per_block)
430  {
431  intervals[block].second = i;
432  intervals[block+1].first = i;
433  ++block;
434 
435  c = 0;
436  };
437  };
438  intervals[n_blocks-1].second = M.m();
439  };
440 
441  // now transfer the knowledge on
442  // the splitting points into the
443  // vector<bool>s that the base
444  // class wants to see. the way how
445  // we do this depends on the
446  // requested blocking strategy
447  switch (blocking_strategy)
448  {
449  case index_intervals:
450  {
451  for (unsigned int block=0; block<n_blocks; ++block)
452  std::fill_n (dof_masks[block].begin()+intervals[block].first,
453  intervals[block].second - intervals[block].first,
454  true);
455  break;
456  };
457 
458  case adaptive:
459  {
460  // the splitting points for
461  // the DoF have been computed
462  // above already, but we will
463  // only use them to split the
464  // Lagrange dofs into
465  // blocks. splitting the
466  // remaining dofs will be
467  // done now.
468 
469  // first count how often the
470  // Lagrange dofs of each
471  // block access the different
472  // dofs
473  Table<2,size_type> access_count (n_blocks, M.m());
474 
475  // set-up the map that will
476  // be used to store the
477  // indices each Lagrange dof
478  // accesses
479  std::map<size_type, size_type> local_index;
480  const SparsityPattern &structure = M.get_sparsity_pattern();
481 
482  for (size_type row=0; row<M.m(); ++row)
483  if (selected[row] == true)
484  {
485  // first find out to
486  // which block the
487  // present row belongs
488  unsigned int block_number = 0;
489  while (row>=intervals[block_number].second)
490  ++block_number;
491  Assert (block_number < n_blocks, ExcInternalError());
492 
493  // now traverse the
494  // matrix structure to
495  // find out to which
496  // dofs number the
497  // present index wants
498  // to write
499  const size_type row_length = structure.row_length(row);
500  for (size_type i=0; i<row_length; ++i)
501  ++access_count[block_number][structure.column_number(row, i)];
502  };
503 
504  // now we know that block @p i
505  // wants to write to DoF @p j
506  // as often as
507  // <tt>access_count[i][j]</tt>
508  // times. Let @p j be allotted
509  // to the block which
510  // accesses it most often.
511  //
512  // if it is a Lagrange dof,
513  // the of course we leave it
514  // to the block we put it
515  // into in the first place
516  for (size_type row=0; row<M.m(); ++row)
517  if (selected[row] == true)
518  {
519  unsigned int block_number = 0;
520  while (row>=intervals[block_number].second)
521  ++block_number;
522 
523  dof_masks[block_number][row] = true;
524  }
525  else
526  {
527  // find out which block
528  // accesses this dof
529  // the most often
530  size_type max_accesses = 0;
531  unsigned int max_access_block = 0;
532  for (unsigned int block=0; block<n_blocks; ++block)
533  if (access_count[block][row] > max_accesses)
534  {
535  max_accesses = access_count[block][row];
536  max_access_block = block;
537  };
538  dof_masks[max_access_block][row] = true;
539  };
540 
541  break;
542  };
543 
544  default:
545  Assert (false, ExcInternalError());
546  };
547 }
548 
549 
550 
551 template <typename number>
552 template <typename number2>
554  const Vector<number2> &src) const
555 {
556  dst = 0;
557 
558  // if no blocking is required, pass
559  // down to the underlying class
560  if (n_blocks == 1)
561  this->apply_preconditioner (dst, src);
562  else
563  // otherwise: blocking requested
564  {
565 #ifdef DEAL_II_WITH_THREADS
566  // spawn threads. since
567  // some compilers have
568  // trouble finding out
569  // which 'encapsulate'
570  // function to take of all
571  // those possible ones if
572  // we simply drop in the
573  // address of an overloaded
574  // template member
575  // function, make it
576  // simpler for the compiler
577  // by giving it the correct
578  // type right away:
579  typedef void (SparseVanka<number>::*mem_fun_p)
580  (Vector<number2> &,
581  const Vector<number2> &,
582  const std::vector<bool> *const) const;
583  const mem_fun_p comp
584  = &SparseVanka<number>::template apply_preconditioner<number2>;
585  Threads::ThreadGroup<> threads;
586  for (unsigned int block=0; block<n_blocks; ++block)
587  threads += Threads::new_thread (comp,
588  *static_cast<const SparseVanka<number>*>(this),
589  dst, src,&dof_masks[block]);
590  threads.join_all ();
591 #else
592  for (unsigned int block=0; block<n_blocks; ++block)
593  this->apply_preconditioner (dst, src,
594  &dof_masks[block]);
595 #endif
596  }
597 }
598 
599 
600 
601 template <typename number>
602 std::size_t
604 {
605  std::size_t mem = SparseVanka<number>::memory_consumption();
606  for (size_type i=0; i<dof_masks.size(); ++i)
607  mem += MemoryConsumption::memory_consumption (dof_masks[i]);
608  return mem;
609 }
610 
611 
612 
613 DEAL_II_NAMESPACE_CLOSE
614 
615 #endif
void vmult(Vector< number2 > &dst, const Vector< number2 > &src) const
SparseVanka(const SparseMatrix< number > &M, const std::vector< bool > &selected, const bool conserve_memory=false, const unsigned int n_threads=multithread_info.n_threads())
size_type m() const
STL namespace.
Thread< RT > new_thread(const std_cxx1x::function< RT()> &function)
const_iterator begin() const
SparseBlockVanka(const SparseMatrix< number > &M, const std::vector< bool > &selected, const unsigned int n_blocks, const BlockingStrategy blocking_strategy, const bool conserve_memory=false, const unsigned int n_threads=multithread_info.n_threads())
#define Assert(cond, exc)
Definition: exceptions.h:299
std::size_t memory_consumption(const T &t)
size_type n() const
std::size_t memory_consumption() const
unsigned int row_length(const size_type row) const
size_type column_number(const size_type row, const unsigned int index) const
types::global_dof_index size_type
Definition: sparse_vanka.h:140
void compute_dof_masks(const SparseMatrix< number > &M, const std::vector< bool > &selected, const BlockingStrategy blocking_strategy)
const SparsityPattern & get_sparsity_pattern() const
std::size_t memory_consumption() const
const_iterator end() const
void compute_inverse(const size_type row, std::vector< size_type > &local_indices)
size_type max_entries_per_row() const
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Definition: table.h:33
::ExceptionBase & ExcInternalError()
void apply_preconditioner(Vector< number2 > &dst, const Vector< number2 > &src, const std::vector< bool > *const dof_mask=0) const
std::size_t size() const
void vmult(Vector< number2 > &dst, const Vector< number2 > &src) const