Reference documentation for deal.II version 8.1.0
precondition_block.templates.h
1 // ---------------------------------------------------------------------
2 // @f$Id: precondition_block.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__precondition_block_templates_h
18 #define __deal2__precondition_block_templates_h
19 
20 
21 #include <deal.II/base/config.h>
23 #include <deal.II/base/logstream.h>
24 #include <deal.II/base/memory_consumption.h>
25 #include <deal.II/lac/householder.h>
26 #include <deal.II/lac/precondition_block.h>
27 #include <deal.II/lac/vector.h>
28 #include <deal.II/lac/full_matrix.h>
29 
31 
32 
33 template<class MATRIX, typename inverse_type>
35 AdditionalData (const size_type block_size,
36  const double relaxation,
37  const bool invert_diagonal,
38  const bool same_diagonal)
39  :
40  relaxation (relaxation),
41  block_size(block_size),
42  invert_diagonal(invert_diagonal),
43  same_diagonal(same_diagonal),
44  inversion(PreconditionBlockBase<inverse_type>::gauss_jordan),
45  threshold(0.)
46 {}
47 
48 
49 template <typename number>
51 {}
52 
53 
54 template <class MATRIX, typename inverse_type>
57  blocksize(0),
58  A(0, typeid(*this).name())
59 {}
60 
61 
62 template <class MATRIX, typename inverse_type>
64 {}
65 
66 
67 template <class MATRIX, typename inverse_type>
69 {
71  blocksize = 0;
72  A = 0;
73 }
74 
75 
76 template <class MATRIX, typename inverse_type>
78  const MATRIX &M,
79  const AdditionalData parameters)
80 {
81  const size_type bsize = parameters.block_size;
82 
83  clear();
84  Assert (M.m() == M.n(), ExcNotQuadratic());
85  A = &M;
86  Assert (bsize>0, ExcIndexRange(bsize, 1, M.m()));
87  Assert (A->m()%bsize==0, ExcWrongBlockSize(bsize, A->m()));
88  blocksize=bsize;
89  relaxation = parameters.relaxation;
90  const unsigned int nblocks = A->m()/bsize;
91  this->reinit(nblocks, blocksize, parameters.same_diagonal,
92  parameters.inversion);
93 
94  if (parameters.invert_diagonal)
95  {
96  if (permutation.size() == M.m())
97  invert_permuted_diagblocks(permutation, inverse_permutation);
98  else
99  invert_diagblocks();
100  }
101 }
102 
103 
104 template <class MATRIX, typename inverse_type>
106  const MATRIX &M,
107  const std::vector<size_type> &permutation,
108  const std::vector<size_type> &inverse_permutation,
109  const AdditionalData parameters)
110 {
111  set_permutation(permutation, inverse_permutation);
112  initialize(M, parameters);
113 }
114 
115 template <class MATRIX, typename inverse_type>
117  const std::vector<size_type> &permutation,
118  const std::vector<size_type> &inverse_permutation)
119 {
120  Assert (A!=0, ExcNotInitialized());
121  Assert (blocksize!=0, ExcNotInitialized());
122 
123  const MATRIX &M=*A;
124  Assert (this->inverses_ready()==0, ExcInverseMatricesAlreadyExist());
125  AssertDimension (permutation.size(), M.m());
126  AssertDimension (inverse_permutation.size(), M.m());
127 
128  FullMatrix<inverse_type> M_cell(blocksize);
129 
130  if (this->same_diagonal())
131  {
132  deallog << "PreconditionBlock uses only one diagonal block" << std::endl;
133 
134  for (size_type row_cell=0; row_cell<blocksize; ++row_cell)
135  {
136  typename MATRIX::const_iterator entry = M.begin(row_cell);
137  const typename MATRIX::const_iterator row_end = M.end(row_cell);
138  while (entry != row_end)
139  {
140  if (entry->column() < blocksize)
141  M_cell(row_cell, entry->column()) = entry->value();
142  ++entry;
143  }
144  }
145  if (this->store_diagonals())
146  this->diagonal(0) = M_cell;
147  switch (this->inversion)
148  {
150  this->inverse(0).invert(M_cell);
151  break;
153  this->inverse_householder(0).initialize(M_cell);
154  break;
156  this->inverse_svd(0) = M_cell;
157  this->inverse_svd(0).compute_inverse_svd(0.);
158  break;
159  default:
160  Assert(false, ExcNotImplemented());
161 
162  }
163  }
164  else
165  {
166  // cell_row, cell_column are the
167  // numbering of the blocks (cells).
168  // row_cell, column_cell are the local
169  // numbering of the unknowns in the
170  // blocks.
171  // row, column are the global numbering
172  // of the unkowns.
173  M_cell = 0;
174 
175  for (unsigned int cell=0; cell<this->size(); ++cell)
176  {
177  const size_type cell_start = cell*blocksize;
178  for (size_type row_cell=0; row_cell<blocksize; ++row_cell)
179  {
180  const size_type urow = row_cell + cell_start;
181 
182  const size_type row = permutation[urow];
183 
184  typename MATRIX::const_iterator entry = M.begin(row);
185  const typename MATRIX::const_iterator row_end = M.end(row);
186 
187  for (; entry != row_end; ++entry)
188  {
189  //if (entry->column()<cell_start)
190  if (inverse_permutation[entry->column()]<cell_start)
191  continue;
192 
193  const size_type column_cell = inverse_permutation[entry->column()]-cell_start;
194  if (column_cell >= blocksize)
195  continue;
196  M_cell(row_cell, column_cell) = entry->value();
197  }
198  }
199 
200  if (this->store_diagonals())
201  this->diagonal(cell) = M_cell;
202  switch (this->inversion)
203  {
205  this->inverse(cell).invert(M_cell);
206  break;
208  this->inverse_householder(cell).initialize(M_cell);
209  break;
211  this->inverse_svd(cell) = M_cell;
212  this->inverse_svd(cell).compute_inverse_svd(0.);
213  break;
214  default:
215  Assert(false, ExcNotImplemented());
216 
217  }
218  }
219  }
220  this->inverses_computed(true);
221 }
222 
223 
224 
225 template <class MATRIX, typename inverse_type>
226 template <typename number2>
228  Vector<number2> &dst,
229  const Vector<number2> &prev,
230  const Vector<number2> &src,
231  const bool transpose_diagonal) const
232 {
233  Assert (this->A!=0, ExcNotInitialized());
234 
235  const MATRIX &M=*this->A;
236 
237  if (permutation.size() != 0)
238  Assert (permutation.size() == M.m() || permutation.size() == this->size(),
239  ExcMessage("Permutation vector size must be equal to either the number of blocks or the dimension of the system"));
240 
241  const bool permuted = (permutation.size() == M.m());
242  const bool cell_permuted = (permutation.size() == this->size());
243 
244 // deallog << "Permutation " << permutation.size();
245 // if (permuted) deallog << " point";
246 // if (cell_permuted) deallog << " block";
247 // deallog << std::endl;
248 
249  Vector<number2> b_cell(this->blocksize), x_cell(this->blocksize);
250 
251  // cell_row, cell_column are the
252  // numbering of the blocks (cells).
253  // row_cell, column_cell are the local
254  // numbering of the unknowns in the
255  // blocks.
256  // row, column are the global numbering
257  // of the unkowns.
258  size_type row, row_cell;
259  number2 b_cell_row;
260  // The diagonal block if the
261  // inverses were not precomputed
262  FullMatrix<number> M_cell(this->blocksize);
263 
264  // Loop over all blocks
265  for (unsigned int rawcell=0; rawcell < this->size(); ++rawcell)
266  {
267  const unsigned int cell = cell_permuted ? permutation[rawcell] : rawcell;
268  const size_type block_start = cell*this->blocksize;
269  const size_type permuted_block_start = permuted
270  ? permutation[block_start]
271  : block_start;
272 
273 // deallog << std::endl << cell << '-' << block_start
274 // << '-' << permuted_block_start << (permuted ? 't' : 'f') << '\t';
275 
276  for (row = permuted_block_start, row_cell = 0;
277  row_cell < this->blocksize;
278  ++row_cell, ++row)
279  {
280 // deallog << ' ' << row;
281  const typename MATRIX::const_iterator row_end = M.end(row);
282  typename MATRIX::const_iterator entry = M.begin(row);
283 
284  b_cell_row=src(row);
285  for (; entry != row_end; ++entry)
286  {
287  const size_type column = entry->column();
288  const size_type inverse_permuted_column = permuted
289  ? inverse_permutation[column]
290  : column;
291  b_cell_row -= entry->value() * prev(column);
292 //TODO:[GK] Find out if this is really once column and once permuted
293  if (!this->inverses_ready()
294  && inverse_permuted_column >= block_start
295  && inverse_permuted_column < block_start + this->blocksize)
296  {
297  const size_type column_cell = column - block_start;
298  if (transpose_diagonal)
299  M_cell(column_cell, row_cell) = entry->value();
300  else
301  M_cell(row_cell, column_cell) = entry->value();
302  }
303  }
304  b_cell(row_cell)=b_cell_row;
305  }
306  if (this->inverses_ready())
307  {
308  if (transpose_diagonal)
309  this->inverse_Tvmult(cell, x_cell, b_cell);
310  else
311  this->inverse_vmult(cell, x_cell, b_cell);
312  }
313  else
314  {
315  Householder<number> house(M_cell);
316  house.least_squares(x_cell,b_cell);
317  }
318 
319  // distribute x_cell to dst
320  for (row=permuted_block_start, row_cell=0;
321  row_cell<this->blocksize;
322  ++row_cell, ++row)
323  dst(row) = prev(row) + this->relaxation*x_cell(row_cell);
324  }
325 }
326 
327 
328 template <class MATRIX, typename inverse_type>
329 template <typename number2>
331  Vector<number2> &dst,
332  const Vector<number2> &prev,
333  const Vector<number2> &src,
334  const bool transpose_diagonal) const
335 {
336  Assert (this->A!=0, ExcNotInitialized());
337 
338  const MATRIX &M=*this->A;
339 
340  if (permutation.size() != 0)
341  Assert (permutation.size() == M.m() || permutation.size() == this->size(),
342  ExcMessage("Permutation vector size must be equal to either the number of blocks or the dimension of the system"));
343 
344  const bool permuted = (permutation.size() == M.m());
345  const bool cell_permuted = (permutation.size() == this->size());
346 
347  Vector<number2> b_cell(this->blocksize), x_cell(this->blocksize);
348 
349  // cell_row, cell_column are the
350  // numbering of the blocks (cells).
351  // row_cell, column_cell are the local
352  // numbering of the unknowns in the
353  // blocks.
354  // row, column are the global numbering
355  // of the unkowns.
356  size_type row, row_cell;
357  number2 b_cell_row;
358 
359  FullMatrix<number> M_cell(this->blocksize);
360  for (unsigned int rawcell=this->size(); rawcell!=0 ;)
361  {
362  --rawcell;
363  const unsigned int cell = cell_permuted ? permutation[rawcell] : rawcell;
364  const size_type block_start = cell*this->blocksize;
365  const size_type block_end = block_start + this->blocksize;
366  const size_type permuted_block_start = permuted
367  ? permutation[block_start]
368  : block_start;
369  for (row = permuted_block_start, row_cell = 0;
370  row_cell<this->blocksize;
371  ++row_cell, ++row)
372  {
373  const typename MATRIX::const_iterator row_end = M.end(row);
374  typename MATRIX::const_iterator entry = M.begin(row);
375 
376  b_cell_row=src(row);
377  for (; entry != row_end; ++entry)
378  {
379  const size_type column = entry->column();
380  const size_type inverse_permuted_column = permuted
381  ? inverse_permutation[column]
382  : column;
383  b_cell_row -= entry->value() * prev(column);
384  if (!this->inverses_ready()
385  && inverse_permuted_column < block_end
386  && column >= block_start)
387  {
388  const size_type column_cell = column - block_start;
389  // We need the
390  // transpose of the
391  // diagonal block,
392  // so we switch row
393  // and column
394  // indices
395  if (transpose_diagonal)
396  M_cell(column_cell, row_cell) = entry->value();
397  else
398  M_cell(row_cell, column_cell) = entry->value();
399  }
400  }
401  b_cell(row_cell)=b_cell_row;
402  }
403  if (this->inverses_ready())
404  {
405  if (transpose_diagonal)
406  this->inverse_Tvmult(cell, x_cell, b_cell);
407  else
408  this->inverse_vmult(cell, x_cell, b_cell);
409  }
410  else
411  {
412  Householder<number> house(M_cell);
413  house.least_squares(x_cell,b_cell);
414  }
415 
416 
417  // distribute x_cell to dst
418  for (row=permuted_block_start, row_cell=0;
419  row_cell<this->blocksize;
420  ++row_cell, ++row)
421  dst(row) = prev(row) + this->relaxation*x_cell(row_cell);
422  }
423 }
424 
425 
426 template <class MATRIX, typename inverse_type>
429 {
430  return blocksize;
431 }
432 
433 
434 template <class MATRIX, typename inverse_type>
436 {
437  Assert (A!=0, ExcNotInitialized());
438  Assert (blocksize!=0, ExcNotInitialized());
439 
440  const MATRIX &M=*A;
441  Assert (this->inverses_ready()==0, ExcInverseMatricesAlreadyExist());
442 
443  FullMatrix<inverse_type> M_cell(blocksize);
444 
445  if (this->same_diagonal())
446  {
447  deallog << "PreconditionBlock uses only one diagonal block" << std::endl;
448  for (size_type row_cell=0; row_cell<blocksize; ++row_cell)
449  {
450  typename MATRIX::const_iterator entry = M.begin(row_cell);
451  const typename MATRIX::const_iterator row_end = M.end(row_cell);
452  while (entry != row_end)
453  {
454  if (entry->column() < blocksize)
455  M_cell(row_cell, entry->column()) = entry->value();
456  ++entry;
457  }
458  }
459  if (this->store_diagonals())
460  this->diagonal(0) = M_cell;
461  switch (this->inversion)
462  {
464  this->inverse(0).invert(M_cell);
465  break;
467  this->inverse_householder(0).initialize(M_cell);
468  break;
470  this->inverse_svd(0) = M_cell;
471  this->inverse_svd(0).compute_inverse_svd(1.e-12);
472  break;
473  default:
474  Assert(false, ExcNotImplemented());
475  }
476  }
477  else
478  {
479  M_cell = 0;
480 
481  for (unsigned int cell=0; cell<this->size(); ++cell)
482  {
483  const size_type cell_start = cell*blocksize;
484  for (size_type row_cell=0; row_cell<blocksize; ++row_cell)
485  {
486  const size_type row = row_cell + cell_start;
487  typename MATRIX::const_iterator entry = M.begin(row);
488  const typename MATRIX::const_iterator row_end = M.end(row);
489 
490  for (; entry != row_end; ++entry)
491  {
492  if (entry->column()<cell_start)
493  continue;
494 
495  const size_type column_cell = entry->column()-cell_start;
496  if (column_cell >= blocksize)
497  continue;
498  M_cell(row_cell, column_cell) = entry->value();
499  }
500  }
501 
502  if (this->store_diagonals())
503  this->diagonal(cell) = M_cell;
504  switch (this->inversion)
505  {
507  this->inverse(cell).invert(M_cell);
508  break;
510  this->inverse_householder(cell).initialize(M_cell);
511  break;
513  this->inverse_svd(cell) = M_cell;
514  this->inverse_svd(cell).compute_inverse_svd(1.e-12);
515  break;
516  default:
517  Assert(false, ExcNotImplemented());
518  }
519  }
520  }
521  this->inverses_computed(true);
522 }
523 
524 
525 
526 template <class MATRIX, typename inverse_type>
528  const std::vector<size_type> &p,
529  const std::vector<size_type> &i)
530 {
531  Assert (p.size() == i.size(), ExcDimensionMismatch(p.size(), i.size()));
532 
533  if (this->inverses_ready())
534  {
535  AssertDimension(p.size(), this->size());
536  }
537 
538  permutation.resize(p.size());
539  inverse_permutation.resize(p.size());
540  for (unsigned int k=0; k<p.size(); ++k)
541  {
542  permutation[k] = p[k];
543  inverse_permutation[k] = i[k];
544  }
545 }
546 
547 
548 template <class MATRIX, typename inverse_type>
549 std::size_t
551 {
552  return (sizeof(*this)
555 }
556 
557 
558 
559 
560 /*--------------------- PreconditionBlockJacobi -----------------------*/
561 
562 
563 template <class MATRIX, typename inverse_type>
564 template <typename number2>
567  const Vector<number2> &src,
568  bool adding) const
569 {
570  Assert(this->A!=0, ExcNotInitialized());
571 
572  const MATRIX &M=*this->A;
573 
574  Vector<number2> b_cell(this->blocksize), x_cell(this->blocksize);
575 
576  // cell_row, cell_column are the
577  // numbering of the blocks (cells).
578  // row_cell, column_cell are the local
579  // numbering of the unknowns in the
580  // blocks.
581  // row, column are the global numbering
582  // of the unkowns.
583  size_type row, row_cell, begin_diag_block=0;
584 
585  if (!this->inverses_ready())
586  {
587  FullMatrix<number> M_cell(this->blocksize);
588  for (unsigned int cell=0; cell < this->size(); ++cell)
589  {
590  for (row=cell*this->blocksize, row_cell=0;
591  row_cell<this->blocksize;
592  ++row_cell, ++row)
593  {
594  b_cell(row_cell)=src(row);
595  for (size_type column_cell=0, column=cell*this->blocksize;
596  column_cell<this->blocksize; ++column_cell, ++column)
597  M_cell(row_cell,column_cell)=M(row,column);
598  }
599  Householder<number> house(M_cell);
600  house.least_squares(x_cell,b_cell);
601  // distribute x_cell to dst
602  for (row=cell*this->blocksize, row_cell=0;
603  row_cell<this->blocksize;
604  ++row_cell, ++row)
605  if (adding)
606  dst(row)+=x_cell(row_cell);
607  else
608  dst(row)=x_cell(row_cell);
609 
610  begin_diag_block+=this->blocksize;
611  }
612  }
613  else
614  for (unsigned int cell=0; cell < this->size(); ++cell)
615  {
616  for (row=cell*this->blocksize, row_cell=0;
617  row_cell<this->blocksize;
618  ++row_cell, ++row)
619  {
620  b_cell(row_cell)=src(row);
621  }
622  this->inverse_vmult(cell, x_cell, b_cell);
623  // distribute x_cell to dst
624  for (row=cell*this->blocksize, row_cell=0;
625  row_cell<this->blocksize;
626  ++row_cell, ++row)
627  if (adding)
628  dst(row)+=x_cell(row_cell);
629  else
630  dst(row)=x_cell(row_cell);
631 
632  begin_diag_block+=this->blocksize;
633  }
634  dst *= this->relaxation;
635 }
636 
637 
638 template <class MATRIX, typename inverse_type>
639 template <typename number2>
642  const Vector<number2> &src) const
643 {
644  do_vmult(dst, src, false);
645 }
646 
647 
648 template <class MATRIX, typename inverse_type>
649 template <typename number2>
652  const Vector<number2> &src) const
653 {
654  do_vmult(dst, src, false);
655 }
656 
657 
658 template <class MATRIX, typename inverse_type>
659 template <typename number2>
662  const Vector<number2> &src) const
663 {
664  do_vmult(dst, src, true);
665 }
666 
667 
668 template <class MATRIX, typename inverse_type>
669 template <typename number2>
672  const Vector<number2> &src) const
673 {
674  do_vmult(dst, src, true);
675 }
676 
677 
678 template <class MATRIX, typename inverse_type>
679 template <typename number2>
682  const Vector<number2> &src) const
683 {
685  typename VectorMemory<Vector<number2> >::Pointer aux(mem);
686  aux->reinit(dst);
687 
688  this->forward_step(*aux, dst, src, false);
689  dst = *aux;
690 }
691 
692 
693 template <class MATRIX, typename inverse_type>
694 template <typename number2>
697  const Vector<number2> &src) const
698 {
700  typename VectorMemory<Vector<number2> >::Pointer aux(mem);
701  aux->reinit(dst);
702 
703  this->backward_step(*aux, dst, src, true);
704  dst = *aux;
705 }
706 
707 
708 
709 
710 /*--------------------- PreconditionBlockSOR -----------------------*/
711 
712 
713 template <class MATRIX, typename inverse_type>
716 
717 {}
718 
719 template <class MATRIX, typename inverse_type>
722 
723 {}
724 
725 template <class MATRIX, typename inverse_type>
726 template <typename number2>
728  Vector<number2> &dst,
729  const Vector<number2> &src,
730  const bool transpose_diagonal,
731  const bool) const
732 {
733  Assert (this->A!=0, ExcNotInitialized());
734 
735  const MATRIX &M=*this->A;
736  const bool permuted = (this->permutation.size() != 0);
737  if (permuted)
738  {
739  Assert (this->permutation.size() == M.m(), ExcDimensionMismatch(this->permutation.size(), M.m()));
740  }
741 
742  Vector<number2> b_cell(this->blocksize), x_cell(this->blocksize);
743 
744  // cell_row, cell_column are the
745  // numbering of the blocks (cells).
746  // row_cell, column_cell are the local
747  // numbering of the unknowns in the
748  // blocks.
749  // row, column are the global numbering
750  // of the unkowns.
751  size_type row, row_cell, block_start=0;
752  number2 b_cell_row;
753  // The diagonal block if the
754  // inverses were not precomputed
755  FullMatrix<number> M_cell(this->blocksize);
756 
757  for (unsigned int cell=0; cell < this->size(); ++cell)
758  {
759  const size_type permuted_block_start = permuted
760  ? this->permutation[block_start]
761  :block_start;
762 
763  for (row = permuted_block_start, row_cell = 0;
764  row_cell < this->blocksize;
765  ++row_cell, ++row)
766  {
767  const typename MATRIX::const_iterator row_end = M.end(row);
768  typename MATRIX::const_iterator entry = M.begin(row);
769 
770  b_cell_row=src(row);
771  for (; entry != row_end; ++entry)
772  {
773  const size_type column = entry->column();
774  const size_type inverse_permuted_column = permuted
775  ? this->inverse_permutation[column]
776  : column;
777 
778  if (inverse_permuted_column < block_start)
779  b_cell_row -= entry->value() * dst(column);
780  else if (!this->inverses_ready() && column < block_start + this->blocksize)
781  {
782  const size_type column_cell = column - block_start;
783  if (transpose_diagonal)
784  M_cell(column_cell, row_cell) = entry->value();
785  else
786  M_cell(row_cell, column_cell) = entry->value();
787  }
788  }
789  b_cell(row_cell)=b_cell_row;
790  }
791  if (this->inverses_ready())
792  {
793  if (transpose_diagonal)
794  this->inverse_Tvmult(cell, x_cell, b_cell);
795  else
796  this->inverse_vmult(cell, x_cell, b_cell);
797  }
798  else
799  {
800  Householder<number> house(M_cell);
801  house.least_squares(x_cell,b_cell);
802  }
803 
804  // distribute x_cell to dst
805  for (row=permuted_block_start, row_cell=0;
806  row_cell<this->blocksize;
807  ++row_cell, ++row)
808  dst(row)=this->relaxation*x_cell(row_cell);
809 
810  block_start+=this->blocksize;
811  }
812 }
813 
814 
815 template <class MATRIX, typename inverse_type>
816 template <typename number2>
818  Vector<number2> &dst,
819  const Vector<number2> &src,
820  const bool transpose_diagonal,
821  const bool) const
822 {
823  Assert (this->A!=0, ExcNotInitialized());
824 
825  const MATRIX &M=*this->A;
826  const bool permuted = (this->permutation.size() != 0);
827  if (permuted)
828  {
829  Assert (this->permutation.size() == M.m(), ExcDimensionMismatch(this->permutation.size(), M.m()));
830  }
831 
832  Vector<number2> b_cell(this->blocksize), x_cell(this->blocksize);
833 
834  // cell_row, cell_column are the
835  // numbering of the blocks (cells).
836  // row_cell, column_cell are the local
837  // numbering of the unknowns in the
838  // blocks.
839  // row, column are the global numbering
840  // of the unkowns.
841  size_type row, row_cell;
842  size_type block_end=this->blocksize * this->size();
843  number2 b_cell_row;
844 
845  FullMatrix<number> M_cell(this->blocksize);
846  for (unsigned int cell=this->size(); cell!=0 ;)
847  {
848  --cell;
849  const size_type block_start = block_end - this->blocksize;
850  // Collect upper triangle
851  const size_type permuted_block_start = (this->permutation.size() != 0)
852  ? this->permutation[block_start]
853  : block_start;
854  for (row = permuted_block_start, row_cell = 0;
855  row_cell<this->blocksize;
856  ++row_cell, ++row)
857  {
858  const typename MATRIX::const_iterator row_end = M.end(row);
859  typename MATRIX::const_iterator entry = M.begin(row);
860 
861  b_cell_row=src(row);
862  for (; entry != row_end; ++entry)
863  {
864  const size_type column = entry->column();
865  const size_type inverse_permuted_column = permuted
866  ? this->inverse_permutation[column]
867  : column;
868  if (inverse_permuted_column >= block_end)
869  b_cell_row -= entry->value() * dst(column);
870  else if (!this->inverses_ready() && column >= block_start)
871  {
872  const size_type column_cell = column - block_start;
873  // We need the
874  // transpose of the
875  // diagonal block,
876  // so we switch row
877  // and column
878  // indices
879  if (transpose_diagonal)
880  M_cell(column_cell, row_cell) = entry->value();
881  else
882  M_cell(row_cell, column_cell) = entry->value();
883  }
884  }
885  b_cell(row_cell)=b_cell_row;
886  }
887  if (this->inverses_ready())
888  {
889  if (transpose_diagonal)
890  this->inverse_Tvmult(cell, x_cell, b_cell);
891  else
892  this->inverse_vmult(cell, x_cell, b_cell);
893  }
894  else
895  {
896  Householder<number> house(M_cell);
897  house.least_squares(x_cell,b_cell);
898  }
899 
900 
901  // distribute x_cell to dst
902  for (row=permuted_block_start, row_cell=0;
903  row_cell<this->blocksize;
904  ++row_cell, ++row)
905  dst(row)=this->relaxation*x_cell(row_cell);
906  block_end = block_start;
907 
908  }
909 }
910 
911 
912 template <class MATRIX, typename inverse_type>
913 template <typename number2>
916  const Vector<number2> &src) const
917 {
918  forward(dst, src, false, false);
919 }
920 
921 
922 template <class MATRIX, typename inverse_type>
923 template <typename number2>
926  const Vector<number2> &src) const
927 {
928  forward(dst, src, false, true);
929 }
930 
931 
932 template <class MATRIX, typename inverse_type>
933 template <typename number2>
936  const Vector<number2> &src) const
937 {
938  backward(dst, src, true, false);
939 }
940 
941 
942 template <class MATRIX, typename inverse_type>
943 template <typename number2>
946  const Vector<number2> &src) const
947 {
948  backward(dst, src, true, true);
949 }
950 
951 
952 
953 template <class MATRIX, typename inverse_type>
954 template <typename number2>
957  const Vector<number2> &src) const
958 {
959  this->forward_step(dst, dst, src, false);
960 }
961 
962 
963 template <class MATRIX, typename inverse_type>
964 template <typename number2>
967  const Vector<number2> &src) const
968 {
969  this->backward_step(dst, dst, src, true);
970 }
971 
972 
973 
974 
975 //---------------------------------------------------------------------------
976 
977 
978 template <class MATRIX, typename inverse_type>
981 
982 {}
983 
984 
985 template <class MATRIX, typename inverse_type>
986 template <typename number2>
988  const Vector<number2> &src) const
989 {
990  Vector<number2> help;
991  help.reinit(dst);
992 
993  this->forward(help, src, false, false);
994 
995  Vector<inverse_type> cell_src(this->blocksize);
996  Vector<inverse_type> cell_dst(this->blocksize);
997  const double scaling = (2.-this->relaxation)/this->relaxation;
998 
999  // Multiply with diagonal blocks
1000  for (unsigned int cell=0; cell < this->size(); ++cell)
1001  {
1002  size_type row = cell*this->blocksize;
1003 
1004  for (size_type row_cell=0; row_cell<this->blocksize; ++row_cell)
1005  cell_src(row_cell)=help(row+row_cell);
1006 
1007  this->diagonal(cell).vmult(cell_dst, cell_src);
1008 
1009  for (size_type row_cell=0; row_cell<this->blocksize; ++row_cell)
1010  help(row+row_cell) = scaling * cell_dst(row_cell);
1011  }
1012 
1013  this->backward(dst, help, false, false);
1014 }
1015 
1016 template <class MATRIX, typename inverse_type>
1017 template <typename number2>
1019  const Vector<number2> &src) const
1020 {
1021  Vector<number2> help;
1022  help.reinit(dst);
1023 
1024  this->backward(help, src, true, false);
1025 
1026  Vector<inverse_type> cell_src(this->blocksize);
1027  Vector<inverse_type> cell_dst(this->blocksize);
1028  const double scaling = (2.-this->relaxation)/this->relaxation;
1029 
1030  // Multiply with diagonal blocks
1031  for (unsigned int cell=0; cell < this->size(); ++cell)
1032  {
1033  size_type row = cell*this->blocksize;
1034 
1035  for (size_type row_cell=0; row_cell<this->blocksize; ++row_cell)
1036  cell_src(row_cell)=help(row+row_cell);
1037 
1038  this->diagonal(cell).Tvmult(cell_dst, cell_src);
1039 
1040  for (size_type row_cell=0; row_cell<this->blocksize; ++row_cell)
1041  help(row+row_cell) = scaling * cell_dst(row_cell);
1042  }
1043 
1044  this->forward(dst, help, true, false);
1045 }
1046 
1047 
1048 template <class MATRIX, typename inverse_type>
1049 template <typename number2>
1052  const Vector<number2> &src) const
1053 {
1054  this->forward_step(dst, dst, src, false);
1055  this->backward_step(dst, dst, src, false);
1056 }
1057 
1058 
1059 template <class MATRIX, typename inverse_type>
1060 template <typename number2>
1063  const Vector<number2> &src) const
1064 {
1065  this->backward_step(dst, dst, src, true);
1066  this->forward_step(dst, dst, src, true);
1067 }
1068 
1069 
1070 
1071 DEAL_II_NAMESPACE_CLOSE
1072 
1073 #endif
void Tvmult(Vector< number2 > &, const Vector< number2 > &) const
void Tvmult(Vector< number2 > &, const Vector< number2 > &) const
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:858
void backward_step(Vector< number2 > &dst, const Vector< number2 > &prev, const Vector< number2 > &src, const bool transpose_diagonal) const
void step(Vector< number2 > &dst, const Vector< number2 > &rhs) const
::ExceptionBase & ExcMessage(std::string arg1)
void backward(Vector< number2 > &, const Vector< number2 > &, const bool transpose_diagonal, const bool adding) const
AdditionalData(const size_type block_size, const double relaxation=1., const bool invert_diagonal=true, const bool same_diagonal=false)
void vmult(Vector< number2 > &, const Vector< number2 > &) const
void step(Vector< number2 > &dst, const Vector< number2 > &rhs) const
void initialize(const MATRIX &A, const AdditionalData parameters)
void do_vmult(Vector< number2 > &, const Vector< number2 > &, bool adding) const
void Tstep(Vector< number2 > &dst, const Vector< number2 > &rhs) const
void forward_step(Vector< number2 > &dst, const Vector< number2 > &prev, const Vector< number2 > &src, const bool transpose_diagonal) const
std::size_t memory_consumption() const
void step(Vector< number2 > &dst, const Vector< number2 > &rhs) const
virtual void reinit(const size_type N, const bool fast=false)
#define Assert(cond, exc)
Definition: exceptions.h:299
void forward(Vector< number2 > &, const Vector< number2 > &, const bool transpose_diagonal, const bool adding) const
void Tvmult_add(Vector< number2 > &, const Vector< number2 > &) const
void invert_permuted_diagblocks(const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation)
::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
void set_permutation(const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation)
PreconditionBlockBase< inverse_type >::Inversion inversion
void vmult(Vector< number2 > &, const Vector< number2 > &) const
void Tstep(Vector< number2 > &dst, const Vector< number2 > &rhs) const
PreconditionBlock(bool store_diagonals=false)
::ExceptionBase & ExcNotImplemented()
::ExceptionBase & ExcNotInitialized()
void Tvmult(Vector< number2 > &, const Vector< number2 > &) const
void vmult_add(Vector< number2 > &, const Vector< number2 > &) const
void vmult_add(Vector< number2 > &, const Vector< number2 > &) const
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void vmult(Vector< number2 > &, const Vector< number2 > &) const
SmartPointer< const MATRIX, PreconditionBlock< MATRIX, inverse_type > > A
double least_squares(Vector< number2 > &dst, const Vector< number2 > &src) const
void Tstep(Vector< number2 > &dst, const Vector< number2 > &rhs) const
void Tvmult_add(Vector< number2 > &, const Vector< number2 > &) const
types::global_dof_index size_type