From: kanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d> Date: Sun, 28 Nov 2010 23:21:29 +0000 (+0000) Subject: cleanup X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=73c6e2399b2379e91c4eed666a8da616833a037a;p=dealii-svn.git cleanup git-svn-id: https://svn.dealii.org/trunk@22871 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/lac/precondition_block.h b/deal.II/include/deal.II/lac/precondition_block.h index 5d75ce1ae3..dfb125d5d2 100644 --- a/deal.II/include/deal.II/lac/precondition_block.h +++ b/deal.II/include/deal.II/lac/precondition_block.h @@ -126,6 +126,22 @@ class PreconditionBlock * are equal to save memory. */ bool same_diagonal; + /** + * Choose the inversion + * method for the blocks. + */ + typename PreconditionBlockBase<inverse_type>::Inversion inversion; + + /** + * The if #inversion is SVD, + * the threshold below which + * a singular value will be + * considered zero and thus + * not inverted. This + * parameter is used in the + * call to LAPACKFullMatrix::compute_inverse_svd(). + */ + double threshold; }; @@ -309,6 +325,9 @@ class PreconditionBlock unsigned int block_size () const; /** + * @deprecated Use size() + * instead. + * * The number of blocks of the * matrix. */ @@ -371,11 +390,6 @@ class PreconditionBlock */ double relaxation; - /** - * Number of blocks. - */ - unsigned int nblocks; - /** * The permutation vector */ @@ -547,53 +561,21 @@ class PreconditionBlockJacobi : public virtual Subscriptor, Accessor accessor; }; - /** - * Make initialization function - * publicly available. - */ PreconditionBlock<MATRIX, inverse_type>::initialize; - - /** - * Make function of base class public again. - */ PreconditionBlock<MATRIX, inverse_type>::clear; - - /** - * Make function of base class public again. - */ PreconditionBlock<MATRIX, inverse_type>::empty; - - /** - * Make function of base class public again. - */ PreconditionBlock<MATRIX, inverse_type>::el; - - /** - * Make function of base class public again. - */ PreconditionBlock<MATRIX, inverse_type>::set_same_diagonal; - - /** - * Make function public. - */ PreconditionBlock<MATRIX, inverse_type>::invert_diagblocks; - - /** - * Make function public. - */ PreconditionBlock<MATRIX, inverse_type>::block_size; - + PreconditionBlockBase<inverse_type>::size; + PreconditionBlockBase<inverse_type>::inverse; + PreconditionBlockBase<inverse_type>::inverse_householder; + PreconditionBlockBase<inverse_type>::inverse_svd; /** - * Make function public. + * @deprecated Use size() instead */ PreconditionBlock<MATRIX, inverse_type>::n_blocks; - /** - * Make function accessible to iterator. - */ - PreconditionBlock<MATRIX, inverse_type>::inverse; - /** - * Make function accessible. - */ PreconditionBlock<MATRIX, inverse_type>::set_permutation; /** @@ -749,42 +731,18 @@ class PreconditionBlockSOR : public virtual Subscriptor, /** * Make type publicly available. */ - using PreconditionBlock<MATRIX,inverse_type>::AdditionalData; - - /** - * Make initialization function - * publicly available. - */ - using PreconditionBlock<MATRIX, inverse_type>::initialize; - - /** - * Make function of base class public again. - */ - using PreconditionBlock<MATRIX, inverse_type>::clear; - - /** - * Make function of base class public again. - */ - using PreconditionBlock<MATRIX, inverse_type>::empty; - - /** - * Make function of base class public again. - */ - using PreconditionBlock<MATRIX, inverse_type>::el; - - /** - * Make function of base class public again. - */ - using PreconditionBlock<MATRIX, inverse_type>::set_same_diagonal; - - /** - * Make function of base class public again. - */ - using PreconditionBlock<MATRIX, inverse_type>::invert_diagblocks; - /** - * Make function accessible. - */ - using PreconditionBlock<MATRIX, inverse_type>::set_permutation; + PreconditionBlock<MATRIX,inverse_type>::AdditionalData; + PreconditionBlock<MATRIX, inverse_type>::initialize; + PreconditionBlock<MATRIX, inverse_type>::clear; + PreconditionBlock<MATRIX, inverse_type>::empty; + PreconditionBlockBase<inverse_type>::size; + PreconditionBlockBase<inverse_type>::inverse; + PreconditionBlockBase<inverse_type>::inverse_householder; + PreconditionBlockBase<inverse_type>::inverse_svd; + PreconditionBlock<MATRIX, inverse_type>::el; + PreconditionBlock<MATRIX, inverse_type>::set_same_diagonal; + PreconditionBlock<MATRIX, inverse_type>::invert_diagblocks; + PreconditionBlock<MATRIX, inverse_type>::set_permutation; /** * Execute block SOR @@ -957,46 +915,28 @@ class PreconditionBlockSSOR : public virtual Subscriptor, */ PreconditionBlockSSOR (); - /** - * Make type publicly available. - */ - using PreconditionBlockSOR<MATRIX,inverse_type>::AdditionalData; + // Keep AdditionalData accessible + PreconditionBlockSOR<MATRIX,inverse_type>::AdditionalData; + // The following are the + // functions of the base classes + // which we want to keep + // accessible. /** * Make initialization function * publicly available. */ - using PreconditionBlockSOR<MATRIX,inverse_type>::initialize; - - /** - * Make function of base class public again. - */ - using PreconditionBlockSOR<MATRIX,inverse_type>::clear; - - /** - * Make function of base class public again. - */ - using PreconditionBlockSOR<MATRIX,inverse_type>::set_permutation; - - /** - * Make function of base class public again. - */ - using PreconditionBlockSOR<MATRIX, inverse_type>::empty; - - /** - * Make function of base class public again. - */ - using PreconditionBlockSOR<MATRIX, inverse_type>::el; - - /** - * Make function of base class public again. - */ - using PreconditionBlockSOR<MATRIX,inverse_type>::set_same_diagonal; - - /** - * Make function of base class public again. - */ - using PreconditionBlockSOR<MATRIX,inverse_type>::invert_diagblocks; + PreconditionBlockSOR<MATRIX,inverse_type>::initialize; + PreconditionBlockSOR<MATRIX,inverse_type>::clear; + PreconditionBlockBase<inverse_type>::size; + PreconditionBlockBase<inverse_type>::inverse; + PreconditionBlockBase<inverse_type>::inverse_householder; + PreconditionBlockBase<inverse_type>::inverse_svd; + PreconditionBlockSOR<MATRIX,inverse_type>::set_permutation; + PreconditionBlockSOR<MATRIX, inverse_type>::empty; + PreconditionBlockSOR<MATRIX, inverse_type>::el; + PreconditionBlockSOR<MATRIX,inverse_type>::set_same_diagonal; + PreconditionBlockSOR<MATRIX,inverse_type>::invert_diagblocks; /** * Execute block SSOR @@ -1039,21 +979,6 @@ class PreconditionBlockSSOR : public virtual Subscriptor, #ifndef DOXYGEN -template<class MATRIX, typename inverse_type> -inline -PreconditionBlock<MATRIX, inverse_type>::AdditionalData:: -AdditionalData (const unsigned int block_size, - const double relaxation, - const bool invert_diagonal, - const bool same_diagonal) - : - relaxation (relaxation), - block_size(block_size), - invert_diagonal(invert_diagonal), - same_diagonal(same_diagonal) -{} - - template<class MATRIX, typename inverse_type> inline bool PreconditionBlock<MATRIX, inverse_type>::empty () const @@ -1068,7 +993,7 @@ template<class MATRIX, typename inverse_type> inline unsigned int PreconditionBlock<MATRIX, inverse_type>::n_blocks () const { - return nblocks; + return this->size(); } diff --git a/deal.II/include/deal.II/lac/precondition_block.templates.h b/deal.II/include/deal.II/lac/precondition_block.templates.h index 4e28fa9d23..67b0ab6087 100644 --- a/deal.II/include/deal.II/lac/precondition_block.templates.h +++ b/deal.II/include/deal.II/lac/precondition_block.templates.h @@ -26,6 +26,22 @@ DEAL_II_NAMESPACE_OPEN +template<class MATRIX, typename inverse_type> +PreconditionBlock<MATRIX, inverse_type>::AdditionalData:: +AdditionalData (const unsigned int block_size, + const double relaxation, + const bool invert_diagonal, + const bool same_diagonal) + : + relaxation (relaxation), + block_size(block_size), + invert_diagonal(invert_diagonal), + same_diagonal(same_diagonal), + inversion(PreconditionBlockBase<inverse_type>::gauss_jordan), + threshold(0.) +{} + + template <typename number> PreconditionBlockBase<number>::~PreconditionBlockBase () {} @@ -67,8 +83,9 @@ void PreconditionBlock<MATRIX,inverse_type>::initialize ( Assert (A->m()%bsize==0, ExcWrongBlockSize(bsize, A->m())); blocksize=bsize; relaxation = parameters.relaxation; - nblocks = A->m()/bsize; - this->reinit(nblocks, blocksize, parameters.same_diagonal); + const unsigned int nblocks = A->m()/bsize; + this->reinit(nblocks, blocksize, parameters.same_diagonal, + parameters.inversion); if (parameters.invert_diagonal) invert_diagblocks(); @@ -92,7 +109,7 @@ void PreconditionBlock<MATRIX,inverse_type>::initialize ( Assert (A->m()%bsize==0, ExcWrongBlockSize(bsize, A->m())); blocksize=bsize; relaxation = parameters.relaxation; - nblocks = A->m()/bsize; + const unsigned int nblocks = A->m()/bsize; this->reinit(nblocks, blocksize, parameters.same_diagonal); if (parameters.invert_diagonal) @@ -158,7 +175,7 @@ void PreconditionBlock<MATRIX,inverse_type>::invert_permuted_diagblocks( // of the unkowns. M_cell = 0; - for (unsigned int cell=0; cell<nblocks; ++cell) + for (unsigned int cell=0; cell<this->size(); ++cell) { const unsigned int cell_start = cell*blocksize; for (unsigned int row_cell=0; row_cell<blocksize; ++row_cell) @@ -193,10 +210,10 @@ void PreconditionBlock<MATRIX,inverse_type>::invert_permuted_diagblocks( case PreconditionBlockBase<inverse_type>::householder: this->inverse_householder(cell).initialize(M_cell); break; - case PreconditionBlockBase<inverse_type>::svd: - this->inverse_svd(cell) = M_cell; - this->inverse_svd(cell).compute_inverse_svd(0.); - break; + case PreconditionBlockBase<inverse_type>::svd: + this->inverse_svd(cell) = M_cell; + this->inverse_svd(cell).compute_inverse_svd(0.); + break; default: Assert(false, ExcNotImplemented()); @@ -241,7 +258,7 @@ void PreconditionBlock<MATRIX,inverse_type>::forward_step ( FullMatrix<number> M_cell(this->blocksize); // Loop over all blocks - for (unsigned int cell=0; cell < this->nblocks; ++cell) + for (unsigned int cell=0; cell < this->size(); ++cell) { const unsigned int permuted_block_start = permuted ? permutation[block_start] @@ -328,11 +345,11 @@ void PreconditionBlock<MATRIX,inverse_type>::backward_step ( // row, column are the global numbering // of the unkowns. unsigned int row, row_cell; - unsigned int block_end=this->blocksize * this->nblocks; + unsigned int block_end=this->blocksize * this->size(); number2 b_cell_row; FullMatrix<number> M_cell(this->blocksize); - for (unsigned int cell=this->nblocks; cell!=0 ;) + for (unsigned int cell=this->size(); cell!=0 ;) { --cell; const unsigned int block_start = block_end - this->blocksize; @@ -433,13 +450,27 @@ void PreconditionBlock<MATRIX,inverse_type>::invert_diagblocks() } if (this->store_diagonals()) this->diagonal(0) = M_cell; - this->inverse(0).invert(M_cell); + switch(this->inversion) + { + case PreconditionBlockBase<inverse_type>::gauss_jordan: + this->inverse(0).invert(M_cell); + break; + case PreconditionBlockBase<inverse_type>::householder: + this->inverse_householder(0).initialize(M_cell); + break; + case PreconditionBlockBase<inverse_type>::svd: + this->inverse_svd(0) = M_cell; + this->inverse_svd(0).compute_inverse_svd(1.e-12); + break; + default: + Assert(false, ExcNotImplemented()); + } } else { M_cell = 0; - for (unsigned int cell=0; cell<nblocks; ++cell) + for (unsigned int cell=0; cell<this->size(); ++cell) { const unsigned int cell_start = cell*blocksize; for (unsigned int row_cell=0; row_cell<blocksize; ++row_cell) @@ -462,7 +493,21 @@ void PreconditionBlock<MATRIX,inverse_type>::invert_diagblocks() if (this->store_diagonals()) this->diagonal(cell) = M_cell; - this->inverse(cell).invert(M_cell); + switch(this->inversion) + { + case PreconditionBlockBase<inverse_type>::gauss_jordan: + this->inverse(cell).invert(M_cell); + break; + case PreconditionBlockBase<inverse_type>::householder: + this->inverse_householder(cell).initialize(M_cell); + break; + case PreconditionBlockBase<inverse_type>::svd: + this->inverse_svd(cell) = M_cell; + this->inverse_svd(cell).compute_inverse_svd(1.e-12); + break; + default: + Assert(false, ExcNotImplemented()); + } } } this->inverses_computed(true); @@ -543,7 +588,7 @@ void PreconditionBlockJacobi<MATRIX,inverse_type> if (!this->inverses_ready()) { FullMatrix<number> M_cell(this->blocksize); - for (unsigned int cell=0; cell < this->nblocks; ++cell) + for (unsigned int cell=0; cell < this->size(); ++cell) { for (row=cell*this->blocksize, row_cell=0; row_cell<this->blocksize; @@ -569,7 +614,7 @@ void PreconditionBlockJacobi<MATRIX,inverse_type> } } else - for (unsigned int cell=0; cell < this->nblocks; ++cell) + for (unsigned int cell=0; cell < this->size(); ++cell) { for (row=cell*this->blocksize, row_cell=0; row_cell<this->blocksize; @@ -730,7 +775,7 @@ void PreconditionBlockSOR<MATRIX,inverse_type>::forward ( // inverses were not precomputed FullMatrix<number> M_cell(this->blocksize); - for (unsigned int cell=0; cell < this->nblocks; ++cell) + for (unsigned int cell=0; cell < this->size(); ++cell) { const unsigned int permuted_block_start = permuted ? this->permutation[block_start] @@ -833,11 +878,11 @@ void PreconditionBlockSOR<MATRIX,inverse_type>::backward ( // row, column are the global numbering // of the unkowns. unsigned int row, row_cell; - unsigned int block_end=this->blocksize * this->nblocks; + unsigned int block_end=this->blocksize * this->size(); number2 b_cell_row; FullMatrix<number> M_cell(this->blocksize); - for (unsigned int cell=this->nblocks; cell!=0 ;) + for (unsigned int cell=this->size(); cell!=0 ;) { --cell; const unsigned int block_start = block_end - this->blocksize; @@ -991,7 +1036,7 @@ void PreconditionBlockSSOR<MATRIX,inverse_type>::vmult (Vector<number2> &d const double scaling = (2.-this->relaxation)/this->relaxation; // Multiply with diagonal blocks - for (unsigned int cell=0; cell < this->nblocks; ++cell) + for (unsigned int cell=0; cell < this->size(); ++cell) { unsigned int row = cell*this->blocksize; @@ -1022,7 +1067,7 @@ void PreconditionBlockSSOR<MATRIX,inverse_type>::Tvmult (Vector<number2> & const double scaling = (2.-this->relaxation)/this->relaxation; // Multiply with diagonal blocks - for (unsigned int cell=0; cell < this->nblocks; ++cell) + for (unsigned int cell=0; cell < this->size(); ++cell) { unsigned int row = cell*this->blocksize;