1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15
16 #ifndef dealii_precondition_block_h
17 #define dealii_precondition_block_h
18
19
20 #include <deal.II/base/config.h>
21
22 #include <deal.II/base/exceptions.h>
23 #include <deal.II/base/smartpointer.h>
24 #include <deal.II/base/subscriptor.h>
25
26 #include <deal.II/lac/precondition_block_base.h>
27
28 #include <vector>
29
30 DEAL_II_NAMESPACE_OPEN
31
32 /*! @addtogroup Preconditioners
33 *@{
34 */
35
36
37 /**
38 * Base class for actual block preconditioners. This class assumes the
39 * <tt>MatrixType</tt> consisting of invertible blocks of @p blocksize on the
40 * diagonal and provides the inversion of the diagonal blocks of the matrix.
41 * It is not necessary for this class that the matrix be block diagonal;
42 * rather, it applies to matrices of arbitrary structure with the minimal
43 * property of having invertible blocks on the diagonal. Still the matrix must
44 * have access to single matrix entries. Therefore, BlockMatrixArray and
45 * similar classes are not a possible matrix class template arguments.
46 *
47 * The block matrix structure used by this class is given, e.g., for the DG
48 * method for the transport equation. For a downstream numbering the matrices
49 * even have got a block lower left matrix structure, i.e. the matrices are
50 * empty above the diagonal blocks.
51 *
52 * @note This class is intended to be used for matrices whose structure is
53 * given by local contributions from disjoint cells, such as for DG methods.
54 * It is not intended for problems where the block structure results from
55 * different physical variables such as in the Stokes equations considered in
56 * step-22.
57 *
58 * For all matrices that are empty above and below the diagonal blocks (i.e.
59 * for all block diagonal matrices) the @p BlockJacobi preconditioner is a
60 * direct solver. For all matrices that are empty only above the diagonal
61 * blocks (e.g. the matrices one gets by the DG method with downstream
62 * numbering) @p BlockSOR is a direct solver.
63 *
64 * This first implementation of the @p PreconditionBlock assumes the matrix
65 * has blocks each of the same block size. Varying block sizes within the
66 * matrix must still be implemented if needed.
67 *
68 * The first template parameter denotes the type of number representation in
69 * the sparse matrix, the second denotes the type of number representation in
70 * which the inverted diagonal block matrices are stored within this class by
71 * <tt>invert_diagblocks()</tt>. If you don't want to use the block inversion
72 * as an exact solver, but rather as a preconditioner, you may probably want
73 * to store the inverted blocks with less accuracy than the original matrix;
74 * for example, <tt>number==double, inverse_type=float</tt> might be a viable
75 * choice.
76 *
77 * @see
78 * @ref GlossBlockLA "Block (linear algebra)"
79 */
80 template <typename MatrixType,
81 typename inverse_type = typename MatrixType::value_type>
82 class PreconditionBlock : public virtual Subscriptor,
83 protected PreconditionBlockBase<inverse_type>
84 {
85 private:
86 /**
87 * Define number type of matrix.
88 */
89 using number = typename MatrixType::value_type;
90
91 /**
92 * Value type for inverse matrices.
93 */
94 using value_type = inverse_type;
95
96 public:
97 /**
98 * Declare type for container size.
99 */
100 using size_type = types::global_dof_index;
101
102 /**
103 * Parameters for block preconditioners.
104 */
105 class AdditionalData
106 {
107 public:
108 /**
109 * Constructor. Block size must be given since there is no reasonable
110 * default parameter.
111 */
112 AdditionalData(const size_type block_size,
113 const double relaxation = 1.,
114 const bool invert_diagonal = true,
115 const bool same_diagonal = false);
116
117 /**
118 * Relaxation parameter.
119 */
120 double relaxation;
121
122 /**
123 * Block size.
124 */
125 size_type block_size;
126
127 /**
128 * Invert diagonal during initialization.
129 */
130 bool invert_diagonal;
131
132 /**
133 * Assume all diagonal blocks are equal to save memory.
134 */
135 bool same_diagonal;
136 /**
137 * Choose the inversion method for the blocks.
138 */
139 typename PreconditionBlockBase<inverse_type>::Inversion inversion;
140
141 /**
142 * The if #inversion is SVD, the threshold below which a singular value
143 * will be considered zero and thus not inverted. This parameter is used
144 * in the call to LAPACKFullMatrix::compute_inverse_svd().
145 */
146 double threshold;
147 };
148
149
150 /**
151 * Constructor.
152 */
153 PreconditionBlock(bool store_diagonals = false);
154
155 /**
156 * Destructor.
157 */
158 ~PreconditionBlock() override = default;
159
160 /**
161 * Initialize matrix and block size. We store the matrix and the block size
162 * in the preconditioner object. In a second step, the inverses of the
163 * diagonal blocks may be computed.
164 *
165 * Additionally, a relaxation parameter for derived classes may be provided.
166 */
167 void
168 initialize(const MatrixType &A, const AdditionalData parameters);
169
170 protected:
171 /**
172 * Initialize matrix and block size for permuted preconditioning.
173 * Additionally to the parameters of the other initialize() function, we
174 * hand over two index vectors with the permutation and its inverse. For the
175 * meaning of these vectors see PreconditionBlockSOR.
176 *
177 * In a second step, the inverses of the diagonal blocks may be computed.
178 * Make sure you use invert_permuted_diagblocks() to yield consistent data.
179 *
180 * Additionally, a relaxation parameter for derived classes may be provided.
181 */
182 void
183 initialize(const MatrixType & A,
184 const std::vector<size_type> &permutation,
185 const std::vector<size_type> &inverse_permutation,
186 const AdditionalData parameters);
187
188 /**
189 * Set either the permutation of rows or the permutation of blocks,
190 * depending on the size of the vector.
191 *
192 * If the size of the permutation vectors is equal to the dimension of the
193 * linear system, it is assumed that rows are permuted individually. In this
194 * case, set_permutation() must be called before initialize(), since the
195 * diagonal blocks are built from the permuted entries of the matrix.
196 *
197 * If the size of the permutation vector is not equal to the dimension of
198 * the system, the diagonal blocks are computed from the unpermuted entries.
199 * Instead, the relaxation methods step() and Tstep() are executed applying
200 * the blocks in the order given by the permutation vector. They will throw
201 * an exception if length of this vector is not equal to the number of
202 * blocks.
203 *
204 * @note Permutation of blocks can only be applied to the relaxation
205 * operators step() and Tstep(), not to the preconditioning operators
206 * vmult() and Tvmult().
207 *
208 * @note It is safe to call set_permutation() before initialize(), while the
209 * other order is only admissible for block permutation.
210 */
211 void
212 set_permutation(const std::vector<size_type> &permutation,
213 const std::vector<size_type> &inverse_permutation);
214
215 /**
216 * Replacement of invert_diagblocks() for permuted preconditioning.
217 */
218 void
219 invert_permuted_diagblocks();
220
221 public:
222 /**
223 * Deletes the inverse diagonal block matrices if existent, sets the
224 * blocksize to 0, hence leaves the class in the state that it had directly
225 * after calling the constructor.
226 */
227 void
228 clear();
229
230 /**
231 * Check whether the object is empty.
232 */
233 bool
234 empty() const;
235
236 /**
237 * Read-only access to entries. This function is only possible if the
238 * inverse diagonal blocks are stored.
239 */
240 value_type
241 el(size_type i, size_type j) const;
242
243 /**
244 * Stores the inverse of the diagonal blocks in @p inverse. This costs some
245 * additional memory - for DG methods about 1/3 (for double inverses) or 1/6
246 * (for float inverses) of that used for the matrix - but it makes the
247 * preconditioning much faster.
248 *
249 * It is not allowed to call this function twice (will produce an error)
250 * before a call of <tt>clear(...)</tt> because at the second time there
251 * already exist the inverse matrices.
252 *
253 * After this function is called, the lock on the matrix given through the
254 * @p use_matrix function is released, i.e. you may overwrite of delete it.
255 * You may want to do this in case you use this matrix to precondition
256 * another matrix.
257 */
258 void
259 invert_diagblocks();
260
261 /**
262 * Perform one block relaxation step in forward numbering.
263 *
264 * Depending on the arguments @p dst and @p pref, this performs an SOR step
265 * (both reference the same vector) of a Jacobi step (both different
266 * vectors). For the Jacobi step, the calling function must copy @p dst to
267 * @p pref after this.
268 *
269 * @note If a permutation is set, it is automatically honored by this
270 * function.
271 */
272 template <typename number2>
273 void
274 forward_step(Vector<number2> & dst,
275 const Vector<number2> &prev,
276 const Vector<number2> &src,
277 const bool transpose_diagonal) const;
278
279 /**
280 * Perform one block relaxation step in backward numbering.
281 *
282 * Depending on the arguments @p dst and @p pref, this performs an SOR step
283 * (both reference the same vector) of a Jacobi step (both different
284 * vectors). For the Jacobi step, the calling function must copy @p dst to
285 * @p pref after this.
286 *
287 * @note If a permutation is set, it is automatically honored by this
288 * function.
289 */
290 template <typename number2>
291 void
292 backward_step(Vector<number2> & dst,
293 const Vector<number2> &prev,
294 const Vector<number2> &src,
295 const bool transpose_diagonal) const;
296
297
298 /**
299 * Return the size of the blocks.
300 */
301 size_type
302 block_size() const;
303
304 /**
305 * Determine an estimate for the memory consumption (in bytes) of this
306 * object.
307 */
308 std::size_t
309 memory_consumption() const;
310
311 /**
312 * @addtogroup Exceptions
313 * @{
314 */
315
316 /**
317 * For non-overlapping block preconditioners, the block size must divide the
318 * matrix size. If not, this exception gets thrown.
319 */
320 DeclException2(ExcWrongBlockSize,
321 int,
322 int,
323 << "The blocksize " << arg1 << " and the size of the matrix "
324 << arg2 << " do not match.");
325
326 /**
327 * Exception
328 */
329 DeclException0(ExcInverseMatricesAlreadyExist);
330
331 //@}
332
333 protected:
334 /**
335 * Size of the blocks. Each diagonal block is assumed to be of the same
336 * size.
337 */
338 size_type blocksize;
339
340 /**
341 * Pointer to the matrix. Make sure that the matrix exists as long as this
342 * class needs it, i.e. until calling @p invert_diagblocks, or (if the
343 * inverse matrices should not be stored) until the last call of the
344 * preconditoining @p vmult function of the derived classes.
345 */
346 SmartPointer<const MatrixType, PreconditionBlock<MatrixType, inverse_type>> A;
347 /**
348 * Relaxation parameter to be used by derived classes.
349 */
350 double relaxation;
351
352 /**
353 * The permutation vector
354 */
355 std::vector<size_type> permutation;
356
357 /**
358 * The inverse permutation vector
359 */
360 std::vector<size_type> inverse_permutation;
361 };
362
363
364
365 /**
366 * Block Jacobi preconditioning. See PreconditionBlock for requirements on the
367 * matrix. This class satisfies the
368 * @ref ConceptRelaxationType "relaxation concept".
369 *
370 * @note Instantiations for this template are provided for <tt>@<float@> and
371 * @<double@></tt>; others can be generated in application programs (see the
372 * section on
373 * @ref Instantiations
374 * in the manual).
375 */
376 template <typename MatrixType,
377 typename inverse_type = typename MatrixType::value_type>
378 class PreconditionBlockJacobi
379 : public virtual Subscriptor,
380 private PreconditionBlock<MatrixType, inverse_type>
381 {
382 private:
383 /**
384 * Define number type of matrix.
385 */
386 using number = typename MatrixType::value_type;
387
388 public:
389 /**
390 * Declare type for container size.
391 */
392 using size_type = types::global_dof_index;
393
394 /**
395 * Standard-conforming iterator.
396 */
397 class const_iterator
398 {
399 private:
400 /**
401 * Accessor class for iterators
402 */
403 class Accessor
404 {
405 public:
406 /**
407 * Constructor. Since we use accessors only for read access, a const
408 * matrix pointer is sufficient.
409 */
410 Accessor(const PreconditionBlockJacobi<MatrixType, inverse_type> *matrix,
411 const size_type row);
412
413 /**
414 * Row number of the element represented by this object.
415 */
416 size_type
417 row() const;
418
419 /**
420 * Column number of the element represented by this object.
421 */
422 size_type
423 column() const;
424
425 /**
426 * Value of this matrix entry.
427 */
428 inverse_type
429 value() const;
430
431 protected:
432 /**
433 * The matrix accessed.
434 */
435 const PreconditionBlockJacobi<MatrixType, inverse_type> *matrix;
436
437 /**
438 * Save block size here for further reference.
439 */
440 size_type bs;
441
442 /**
443 * Current block number.
444 */
445 size_type a_block;
446
447 /**
448 * Iterator inside block.
449 */
450 typename FullMatrix<inverse_type>::const_iterator b_iterator;
451
452 /**
453 * End of current block.
454 */
455 typename FullMatrix<inverse_type>::const_iterator b_end;
456
457 // Make enclosing class a friend.
458 friend class const_iterator;
459 };
460
461 public:
462 /**
463 * Constructor.
464 */
465 const_iterator(
466 const PreconditionBlockJacobi<MatrixType, inverse_type> *matrix,
467 const size_type row);
468
469 /**
470 * Prefix increment.
471 */
472 const_iterator &
473 operator++();
474
475 /**
476 * Dereferencing operator.
477 */
478 const Accessor &operator*() const;
479
480 /**
481 * Dereferencing operator.
482 */
483 const Accessor *operator->() const;
484
485 /**
486 * Comparison. True, if both iterators point to the same matrix position.
487 */
488 bool
489 operator==(const const_iterator &) const;
490 /**
491 * Inverse of <tt>==</tt>.
492 */
493 bool
494 operator!=(const const_iterator &) const;
495
496 /**
497 * Comparison operator. Result is true if either the first row number is
498 * smaller or if the row numbers are equal and the first index is smaller.
499 */
500 bool
501 operator<(const const_iterator &) const;
502
503 private:
504 /**
505 * Store an object of the accessor class.
506 */
507 Accessor accessor;
508 };
509
510 /**
511 * import functions from private base class
512 */
513 using typename PreconditionBlock<MatrixType, inverse_type>::AdditionalData;
514 using PreconditionBlock<MatrixType, inverse_type>::initialize;
515 using PreconditionBlock<MatrixType, inverse_type>::clear;
516 using PreconditionBlock<MatrixType, inverse_type>::empty;
517 using PreconditionBlock<MatrixType, inverse_type>::el;
518 using PreconditionBlock<MatrixType, inverse_type>::invert_diagblocks;
519 using PreconditionBlock<MatrixType, inverse_type>::block_size;
520 using PreconditionBlockBase<inverse_type>::size;
521 using PreconditionBlockBase<inverse_type>::inverse;
522 using PreconditionBlockBase<inverse_type>::inverse_householder;
523 using PreconditionBlockBase<inverse_type>::inverse_svd;
524 using PreconditionBlockBase<inverse_type>::log_statistics;
525 using PreconditionBlock<MatrixType, inverse_type>::set_permutation;
526
527 /**
528 * Execute block Jacobi preconditioning.
529 *
530 * This function will automatically use the inverse matrices if they exist,
531 * if not then BlockJacobi will need much time inverting the diagonal block
532 * matrices in each preconditioning step.
533 */
534 template <typename number2>
535 void
536 vmult(Vector<number2> &, const Vector<number2> &) const;
537
538 /**
539 * Same as @p vmult, since Jacobi is symmetric.
540 */
541 template <typename number2>
542 void
543 Tvmult(Vector<number2> &, const Vector<number2> &) const;
544 /**
545 * Execute block Jacobi preconditioning, adding to @p dst.
546 *
547 * This function will automatically use the inverse matrices if they exist,
548 * if not then BlockJacobi will need much time inverting the diagonal block
549 * matrices in each preconditioning step.
550 */
551 template <typename number2>
552 void
553 vmult_add(Vector<number2> &, const Vector<number2> &) const;
554
555 /**
556 * Same as @p vmult_add, since Jacobi is symmetric.
557 */
558 template <typename number2>
559 void
560 Tvmult_add(Vector<number2> &, const Vector<number2> &) const;
561
562 /**
563 * Perform one step of the Jacobi iteration.
564 */
565 template <typename number2>
566 void
567 step(Vector<number2> &dst, const Vector<number2> &rhs) const;
568
569 /**
570 * Perform one step of the Jacobi iteration.
571 */
572 template <typename number2>
573 void
574 Tstep(Vector<number2> &dst, const Vector<number2> &rhs) const;
575
576 /**
577 * Iterator starting at the first entry.
578 */
579 const_iterator
580 begin() const;
581
582 /**
583 * Final iterator.
584 */
585 const_iterator
586 end() const;
587
588 /**
589 * Iterator starting at the first entry of row @p r.
590 */
591 const_iterator
592 begin(const size_type r) const;
593
594 /**
595 * Final iterator of row @p r.
596 */
597 const_iterator
598 end(const size_type r) const;
599
600
601 private:
602 /**
603 * Actual implementation of the preconditioner.
604 *
605 * Depending on @p adding, the result of preconditioning is added to the
606 * destination vector.
607 */
608 template <typename number2>
609 void
610 do_vmult(Vector<number2> &, const Vector<number2> &, bool adding) const;
611
612 friend class Accessor;
613 friend class const_iterator;
614 };
615
616
617
618 /**
619 * Block SOR preconditioning. This class satisfies the
620 * @ref ConceptRelaxationType "relaxation concept".
621 *
622 * The functions @p vmult and @p Tvmult execute a (transposed) block-SOR step,
623 * based on the blocks in PreconditionBlock. The elements outside the diagonal
624 * blocks may be distributed arbitrarily.
625 *
626 * See PreconditionBlock for requirements on the matrix. The blocks used in
627 * this class must be contiguous and non-overlapping. An overlapping Schwarz
628 * relaxation method can be found in RelaxationBlockSOR; that class does not
629 * offer preconditioning, though.
630 *
631 * <h3>Permutations</h3>
632 *
633 * Optionally, the entries of the source vector can be treated in the order of
634 * the indices in the permutation vector set by #set_permutation (or the
635 * opposite order for Tvmult()). The inverse permutation is used for storing
636 * elements back into this vector. This functionality is automatically enabled
637 * after a call to set_permutation() with vectors of nonzero size.
638 *
639 * @note The diagonal blocks, like the matrix, are not permuted! Therefore,
640 * the permutation vector can only swap whole blocks. It may not change the
641 * order inside blocks or swap single indices between blocks.
642 *
643 * <h3>Instantiations</h3>
644 *
645 * @note Instantiations for this template are provided for <tt>@<float@> and
646 * @<double@></tt>; others can be generated in application programs (see the
647 * section on
648 * @ref Instantiations
649 * in the manual).
650 */
651 template <typename MatrixType,
652 typename inverse_type = typename MatrixType::value_type>
653 class PreconditionBlockSOR
654 : public virtual Subscriptor,
655 protected PreconditionBlock<MatrixType, inverse_type>
656 {
657 public:
658 /**
659 * Declare type for container size.
660 */
661 using size_type = types::global_dof_index;
662
663 /**
664 * Default constructor.
665 */
666 PreconditionBlockSOR();
667
668 /**
669 * Define number type of matrix.
670 */
671 using number = typename MatrixType::value_type;
672
673 /**
674 * import types and functions from protected base class.
675 */
676 using typename PreconditionBlock<MatrixType, inverse_type>::AdditionalData;
677 using PreconditionBlock<MatrixType, inverse_type>::initialize;
678 using PreconditionBlock<MatrixType, inverse_type>::clear;
679 using PreconditionBlockBase<inverse_type>::size;
680 using PreconditionBlockBase<inverse_type>::inverse;
681 using PreconditionBlockBase<inverse_type>::inverse_householder;
682 using PreconditionBlockBase<inverse_type>::inverse_svd;
683 using PreconditionBlock<MatrixType, inverse_type>::invert_diagblocks;
684 using PreconditionBlock<MatrixType, inverse_type>::set_permutation;
685 using PreconditionBlockBase<inverse_type>::log_statistics;
686
687 /**
688 * Execute block SOR preconditioning.
689 *
690 * This function will automatically use the inverse matrices if they exist,
691 * if not then BlockSOR will waste much time inverting the diagonal block
692 * matrices in each preconditioning step.
693 *
694 * For matrices which are empty above the diagonal blocks BlockSOR is a
695 * direct solver.
696 */
697 template <typename number2>
698 void
699 vmult(Vector<number2> &, const Vector<number2> &) const;
700
701 /**
702 * Execute block SOR preconditioning.
703 *
704 * Warning: this function performs normal @p vmult without adding. The
705 * reason for its existence is that BlockMatrixArray requires the adding
706 * version by default. On the other hand, adding requires an additional
707 * auxiliary vector, which is not desirable.
708 *
709 * @see vmult
710 */
711 template <typename number2>
712 void
713 vmult_add(Vector<number2> &, const Vector<number2> &) const;
714
715 /**
716 * Backward application of vmult().
717 *
718 * In the current implementation, this is not the transpose of vmult(). It
719 * is a transposed Gauss-Seidel algorithm applied to the whole matrix, but
720 * the diagonal blocks being inverted are not transposed. Therefore, it is
721 * the transposed, if the diagonal blocks are symmetric.
722 */
723 template <typename number2>
724 void
725 Tvmult(Vector<number2> &, const Vector<number2> &) const;
726
727 /**
728 * Execute backward block SOR preconditioning.
729 *
730 * Warning: this function performs normal @p vmult without adding. The
731 * reason for its existence is that BlockMatrixArray requires the adding
732 * version by default. On the other hand, adding requires an additional
733 * auxiliary vector, which is not desirable.
734 *
735 * @see vmult
736 */
737 template <typename number2>
738 void
739 Tvmult_add(Vector<number2> &, const Vector<number2> &) const;
740
741 /**
742 * Perform one step of the SOR iteration.
743 */
744 template <typename number2>
745 void
746 step(Vector<number2> &dst, const Vector<number2> &rhs) const;
747
748 /**
749 * Perform one step of the transposed SOR iteration.
750 */
751 template <typename number2>
752 void
753 Tstep(Vector<number2> &dst, const Vector<number2> &rhs) const;
754
755 protected:
756 /**
757 * Constructor to be used by PreconditionBlockSSOR.
758 */
759 PreconditionBlockSOR(bool store);
760
761 /**
762 * Implementation of the forward substitution loop called by vmult() and
763 * vmult_add().
764 *
765 * If a #permutation is set by set_permutation(), it will automatically be
766 * honored by this function.
767 *
768 * The parameter @p adding does not have any function, yet.
769 */
770 template <typename number2>
771 void
772 forward(Vector<number2> &,
773 const Vector<number2> &,
774 const bool transpose_diagonal,
775 const bool adding) const;
776
777 /**
778 * Implementation of the backward substitution loop called by Tvmult() and
779 * Tvmult_add().
780 *
781 * If a #permutation is set by set_permutation(), it will automatically be
782 * honored by this function.
783 *
784 * The parameter @p adding does not have any function, yet.
785 */
786 template <typename number2>
787 void
788 backward(Vector<number2> &,
789 const Vector<number2> &,
790 const bool transpose_diagonal,
791 const bool adding) const;
792 };
793
794
795 /**
796 * Block SSOR preconditioning. This class satisfies the
797 * @ref ConceptRelaxationType "relaxation concept".
798 *
799 * The functions @p vmult and @p Tvmult execute a block-SSOR step, based on
800 * the implementation in PreconditionBlockSOR. This class requires storage of
801 * the diagonal blocks and their inverses.
802 *
803 * See PreconditionBlock for requirements on the matrix. The blocks used in
804 * this class must be contiguous and non-overlapping. An overlapping Schwarz
805 * relaxation method can be found in RelaxationBlockSSOR; that class does not
806 * offer preconditioning, though.
807 *
808 * @note Instantiations for this template are provided for <tt>@<float@> and
809 * @<double@></tt>; others can be generated in application programs (see the
810 * section on
811 * @ref Instantiations
812 * in the manual).
813 */
814 template <typename MatrixType,
815 typename inverse_type = typename MatrixType::value_type>
816 class PreconditionBlockSSOR
817 : public virtual Subscriptor,
818 private PreconditionBlockSOR<MatrixType, inverse_type>
819 {
820 public:
821 /**
822 * Declare type for container size.
823 */
824 using size_type = types::global_dof_index;
825
826 /**
827 * Define number type of matrix.
828 */
829 using number = typename MatrixType::value_type;
830
831 /**
832 * Constructor.
833 */
834 PreconditionBlockSSOR();
835
836 // Keep AdditionalData accessible
837 using typename PreconditionBlockSOR<MatrixType, inverse_type>::AdditionalData;
838
839 // The following are the
840 // functions of the base classes
841 // which we want to keep
842 // accessible.
843 /**
844 * Make initialization function publicly available.
845 */
846 using PreconditionBlockSOR<MatrixType, inverse_type>::initialize;
847 using PreconditionBlockSOR<MatrixType, inverse_type>::clear;
848 using PreconditionBlockBase<inverse_type>::size;
849 using PreconditionBlockBase<inverse_type>::inverse;
850 using PreconditionBlockBase<inverse_type>::inverse_householder;
851 using PreconditionBlockBase<inverse_type>::inverse_svd;
852 using PreconditionBlockBase<inverse_type>::log_statistics;
853 using PreconditionBlockSOR<MatrixType, inverse_type>::set_permutation;
854 using PreconditionBlockSOR<MatrixType, inverse_type>::empty;
855 using PreconditionBlockSOR<MatrixType, inverse_type>::el;
856 using PreconditionBlockSOR<MatrixType, inverse_type>::invert_diagblocks;
857
858 /**
859 * Execute block SSOR preconditioning.
860 *
861 * This function will automatically use the inverse matrices if they exist,
862 * if not then BlockSOR will waste much time inverting the diagonal block
863 * matrices in each preconditioning step.
864 */
865 template <typename number2>
866 void
867 vmult(Vector<number2> &, const Vector<number2> &) const;
868
869 /**
870 * Same as vmult()
871 */
872 template <typename number2>
873 void
874 Tvmult(Vector<number2> &, const Vector<number2> &) const;
875
876 /**
877 * Perform one step of the SOR iteration.
878 */
879 template <typename number2>
880 void
881 step(Vector<number2> &dst, const Vector<number2> &rhs) const;
882
883 /**
884 * Perform one step of the transposed SOR iteration.
885 */
886 template <typename number2>
887 void
888 Tstep(Vector<number2> &dst, const Vector<number2> &rhs) const;
889 };
890
891 /*@}*/
892 //---------------------------------------------------------------------------
893
894 #ifndef DOXYGEN
895
896 template <typename MatrixType, typename inverse_type>
897 inline bool
empty()898 PreconditionBlock<MatrixType, inverse_type>::empty() const
899 {
900 if (A == nullptr)
901 return true;
902 return A->empty();
903 }
904
905
906 template <typename MatrixType, typename inverse_type>
907 inline inverse_type
el(size_type i,size_type j)908 PreconditionBlock<MatrixType, inverse_type>::el(size_type i, size_type j) const
909 {
910 const size_type bs = blocksize;
911 const unsigned int nb = i / bs;
912
913 const FullMatrix<inverse_type> &B = this->inverse(nb);
914
915 const size_type ib = i % bs;
916 const size_type jb = j % bs;
917
918 if (jb + nb * bs != j)
919 {
920 return 0.;
921 }
922
923 return B(ib, jb);
924 }
925
926 //---------------------------------------------------------------------------
927
928 template <typename MatrixType, typename inverse_type>
929 inline PreconditionBlockJacobi<MatrixType, inverse_type>::const_iterator::
Accessor(const PreconditionBlockJacobi<MatrixType,inverse_type> * matrix,const size_type row)930 Accessor::Accessor(
931 const PreconditionBlockJacobi<MatrixType, inverse_type> *matrix,
932 const size_type row)
933 : matrix(matrix)
934 , bs(matrix->block_size())
935 , a_block(row / bs)
936 , b_iterator(&matrix->inverse(0), 0, 0)
937 , b_end(&matrix->inverse(0), 0, 0)
938 {
939 // This is the end accessor, which
940 // does not have a valid block.
941 if (a_block == matrix->size())
942 return;
943
944 const size_type r = row % bs;
945
946 b_iterator = matrix->inverse(a_block).begin(r);
947 b_end = matrix->inverse(a_block).end();
948
949 AssertIndexRange(a_block, matrix->size());
950 }
951
952
953 template <typename MatrixType, typename inverse_type>
954 inline typename PreconditionBlockJacobi<MatrixType, inverse_type>::size_type
955 PreconditionBlockJacobi<MatrixType,
row()956 inverse_type>::const_iterator::Accessor::row() const
957 {
958 Assert(a_block < matrix->size(), ExcIteratorPastEnd());
959
960 return bs * a_block + b_iterator->row();
961 }
962
963
964 template <typename MatrixType, typename inverse_type>
965 inline typename PreconditionBlockJacobi<MatrixType, inverse_type>::size_type
966 PreconditionBlockJacobi<MatrixType,
column()967 inverse_type>::const_iterator::Accessor::column() const
968 {
969 Assert(a_block < matrix->size(), ExcIteratorPastEnd());
970
971 return bs * a_block + b_iterator->column();
972 }
973
974
975 template <typename MatrixType, typename inverse_type>
976 inline inverse_type
977 PreconditionBlockJacobi<MatrixType,
value()978 inverse_type>::const_iterator::Accessor::value() const
979 {
980 Assert(a_block < matrix->size(), ExcIteratorPastEnd());
981
982 return b_iterator->value();
983 }
984
985
986 template <typename MatrixType, typename inverse_type>
987 inline PreconditionBlockJacobi<MatrixType, inverse_type>::const_iterator::
const_iterator(const PreconditionBlockJacobi<MatrixType,inverse_type> * matrix,const size_type row)988 const_iterator(
989 const PreconditionBlockJacobi<MatrixType, inverse_type> *matrix,
990 const size_type row)
991 : accessor(matrix, row)
992 {}
993
994
995 template <typename MatrixType, typename inverse_type>
996 inline
997 typename PreconditionBlockJacobi<MatrixType, inverse_type>::const_iterator &
998 PreconditionBlockJacobi<MatrixType, inverse_type>::const_iterator::
999 operator++()
1000 {
1001 Assert(*this != accessor.matrix->end(), ExcIteratorPastEnd());
1002
1003 ++accessor.b_iterator;
1004 if (accessor.b_iterator == accessor.b_end)
1005 {
1006 ++accessor.a_block;
1007
1008 if (accessor.a_block < accessor.matrix->size())
1009 {
1010 accessor.b_iterator =
1011 accessor.matrix->inverse(accessor.a_block).begin();
1012 accessor.b_end = accessor.matrix->inverse(accessor.a_block).end();
1013 }
1014 }
1015 return *this;
1016 }
1017
1018
1019 template <typename MatrixType, typename inverse_type>
1020 inline const typename PreconditionBlockJacobi<MatrixType, inverse_type>::
1021 const_iterator::Accessor &
1022 PreconditionBlockJacobi<MatrixType, inverse_type>::const_iterator::
1023 operator*() const
1024 {
1025 return accessor;
1026 }
1027
1028
1029 template <typename MatrixType, typename inverse_type>
1030 inline const typename PreconditionBlockJacobi<MatrixType, inverse_type>::
1031 const_iterator::Accessor *
1032 PreconditionBlockJacobi<MatrixType, inverse_type>::const_iterator::
1033 operator->() const
1034 {
1035 return &accessor;
1036 }
1037
1038
1039 template <typename MatrixType, typename inverse_type>
1040 inline bool
1041 PreconditionBlockJacobi<MatrixType, inverse_type>::const_iterator::
1042 operator==(const const_iterator &other) const
1043 {
1044 if (accessor.a_block == accessor.matrix->size() &&
1045 accessor.a_block == other.accessor.a_block)
1046 return true;
1047
1048 if (accessor.a_block != other.accessor.a_block)
1049 return false;
1050
1051 return (accessor.row() == other.accessor.row() &&
1052 accessor.column() == other.accessor.column());
1053 }
1054
1055
1056 template <typename MatrixType, typename inverse_type>
1057 inline bool
1058 PreconditionBlockJacobi<MatrixType, inverse_type>::const_iterator::
1059 operator!=(const const_iterator &other) const
1060 {
1061 return !(*this == other);
1062 }
1063
1064
1065 template <typename MatrixType, typename inverse_type>
1066 inline bool
1067 PreconditionBlockJacobi<MatrixType, inverse_type>::const_iterator::
1068 operator<(const const_iterator &other) const
1069 {
1070 return (accessor.row() < other.accessor.row() ||
1071 (accessor.row() == other.accessor.row() &&
1072 accessor.column() < other.accessor.column()));
1073 }
1074
1075
1076 template <typename MatrixType, typename inverse_type>
1077 inline
1078 typename PreconditionBlockJacobi<MatrixType, inverse_type>::const_iterator
begin()1079 PreconditionBlockJacobi<MatrixType, inverse_type>::begin() const
1080 {
1081 return const_iterator(this, 0);
1082 }
1083
1084
1085 template <typename MatrixType, typename inverse_type>
1086 inline
1087 typename PreconditionBlockJacobi<MatrixType, inverse_type>::const_iterator
end()1088 PreconditionBlockJacobi<MatrixType, inverse_type>::end() const
1089 {
1090 return const_iterator(this, this->size() * this->block_size());
1091 }
1092
1093
1094 template <typename MatrixType, typename inverse_type>
1095 inline
1096 typename PreconditionBlockJacobi<MatrixType, inverse_type>::const_iterator
begin(const size_type r)1097 PreconditionBlockJacobi<MatrixType, inverse_type>::begin(
1098 const size_type r) const
1099 {
1100 AssertIndexRange(r, this->A->m());
1101 return const_iterator(this, r);
1102 }
1103
1104
1105
1106 template <typename MatrixType, typename inverse_type>
1107 inline
1108 typename PreconditionBlockJacobi<MatrixType, inverse_type>::const_iterator
end(const size_type r)1109 PreconditionBlockJacobi<MatrixType, inverse_type>::end(
1110 const size_type r) const
1111 {
1112 AssertIndexRange(r, this->A->m());
1113 return const_iterator(this, r + 1);
1114 }
1115
1116 #endif // DOXYGEN
1117
1118 DEAL_II_NAMESPACE_CLOSE
1119
1120 #endif
1121