From: guido Date: Wed, 18 Sep 2002 14:58:47 +0000 (+0000) Subject: PreconditionBlock... operates on SparseMatrixEZ X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=7dd6692616b40d98156031e0c24583a312328252;p=dealii-svn.git PreconditionBlock... operates on SparseMatrixEZ git-svn-id: https://svn.dealii.org/trunk@6451 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/news/2002/c-3-4.html b/deal.II/doc/news/2002/c-3-4.html index c85c213373..51bbb5f414 100644 --- a/deal.II/doc/news/2002/c-3-4.html +++ b/deal.II/doc/news/2002/c-3-4.html @@ -234,6 +234,17 @@ contributor's names are abbreviated by WB (Wolfgang Bangerth), GK

lac

    +
  1. Changed: The classes PreconditionBlock, PreconditionBlockJacobi, PreconditionBlockSOR, and PreconditionBlockSSOR have changed their + template signature. The first template argument is now the matrix + type, not just a number type. +
    + (GK 2002/09/18) +

    +
  2. New: class BlockVector has a function collect_sizes(), very much as BlockSparsityPattern. This allows diff --git a/deal.II/lac/include/lac/precondition_block.h b/deal.II/lac/include/lac/precondition_block.h index c00ae7f564..b251d65737 100644 --- a/deal.II/lac/include/lac/precondition_block.h +++ b/deal.II/lac/include/lac/precondition_block.h @@ -75,9 +75,15 @@ template class Vector; * * @author Ralf Hartmann, Guido Kanschat, 1999, 2000 */ -template +template class PreconditionBlock : public virtual Subscriptor { + private: + /** + * Define number type of matrix. + */ + typedef typename MATRIX::value_type number; + public: /** * Constructor. @@ -102,7 +108,7 @@ class PreconditionBlock : public virtual Subscriptor * parameter for derived classes * may be provided. */ - void initialize (const SparseMatrix& A, + void initialize (const MATRIX& A, const unsigned int block_size, const double relaxation = 1.); @@ -243,7 +249,7 @@ class PreconditionBlock : public virtual Subscriptor * @p{vmult} function of the * derived classes. */ - SmartPointer > A; + SmartPointer A; /** * Relaxation parameter to be * used by derived classes. @@ -293,31 +299,37 @@ class PreconditionBlock : public virtual Subscriptor * See @ref{PreconditionBlock} for requirements on the matrix. * @author Ralf Hartmann, Guido Kanschat, 1999, 2000 */ -template +template class PreconditionBlockJacobi : public virtual Subscriptor, - private PreconditionBlock + private PreconditionBlock { + private: + /** + * Define number type of matrix. + */ + typedef typename MATRIX::value_type number; + public: /** * Make initialization function * publicly available. */ - PreconditionBlock::initialize; + PreconditionBlock::initialize; /** * Make function of base class public again. */ - PreconditionBlock::clear; + PreconditionBlock::clear; /** * Make function of base class public again. */ - PreconditionBlock::set_same_diagonal; + PreconditionBlock::set_same_diagonal; /** * Make function public. */ - PreconditionBlock::invert_diagblocks; + PreconditionBlock::invert_diagblocks; /** * Execute block Jacobi @@ -384,31 +396,37 @@ class PreconditionBlockJacobi : public virtual Subscriptor, * See @ref{PreconditionBlock} for requirements on the matrix. * @author Ralf Hartmann, Guido Kanschat, 1999, 2000 */ -template +template class PreconditionBlockSOR : public virtual Subscriptor, - protected PreconditionBlock + protected PreconditionBlock { + private: + /** + * Define number type of matrix. + */ + typedef typename MATRIX::value_type number; + public: /** * Make initialization function * publicly available. */ - PreconditionBlock::initialize; + PreconditionBlock::initialize; /** * Make function of base class public again. */ - PreconditionBlock::clear; + PreconditionBlock::clear; /** * Make function of base class public again. */ - PreconditionBlock::set_same_diagonal; + PreconditionBlock::set_same_diagonal; /** * Make function of base class public again. */ - PreconditionBlock::invert_diagblocks; + PreconditionBlock::invert_diagblocks; /** * Execute block SOR @@ -447,10 +465,16 @@ class PreconditionBlockSOR : public virtual Subscriptor, * See @ref{PreconditionBlock} for requirements on the matrix. * @author Ralf Hartmann, Guido Kanschat, 1999, 2000 */ -template +template class PreconditionBlockSSOR : public virtual Subscriptor, - private PreconditionBlockSOR + private PreconditionBlockSOR { + private: + /** + * Define number type of matrix. + */ + typedef typename MATRIX::value_type number; + public: /** * Constructor. @@ -460,22 +484,22 @@ class PreconditionBlockSSOR : public virtual Subscriptor, * Make initialization function * publicly available. */ - PreconditionBlockSOR::initialize; + PreconditionBlockSOR::initialize; /** * Make function of base class public again. */ - PreconditionBlockSOR::clear; + PreconditionBlockSOR::clear; /** * Make function of base class public again. */ - PreconditionBlockSOR::set_same_diagonal; + PreconditionBlockSOR::set_same_diagonal; /** * Make function of base class public again. */ - PreconditionBlockSOR::invert_diagblocks; + PreconditionBlockSOR::invert_diagblocks; /** * Execute block SSOR diff --git a/deal.II/lac/include/lac/precondition_block.templates.h b/deal.II/lac/include/lac/precondition_block.templates.h index 8ea6a1ea42..0def0191be 100644 --- a/deal.II/lac/include/lac/precondition_block.templates.h +++ b/deal.II/lac/include/lac/precondition_block.templates.h @@ -24,8 +24,8 @@ #include -template -PreconditionBlock::PreconditionBlock (): +template +PreconditionBlock::PreconditionBlock (): blocksize(0), A(0), store_diagonals(false), @@ -33,8 +33,8 @@ PreconditionBlock::PreconditionBlock (): {}; -template -PreconditionBlock::~PreconditionBlock () +template +PreconditionBlock::~PreconditionBlock () { if (var_inverse.size()!=0) var_inverse.erase(var_inverse.begin(), var_inverse.end()); @@ -43,8 +43,8 @@ PreconditionBlock::~PreconditionBlock () } -template -void PreconditionBlock::clear () +template +void PreconditionBlock::clear () { if (var_inverse.size()!=0) var_inverse.erase(var_inverse.begin(), var_inverse.end()); @@ -56,8 +56,8 @@ void PreconditionBlock::clear () } -template -void PreconditionBlock::initialize (const SparseMatrix &M, +template +void PreconditionBlock::initialize (const MATRIX &M, const unsigned int bsize, const double relax) { @@ -71,9 +71,9 @@ void PreconditionBlock::initialize (const SparseMatrix +template const FullMatrix& -PreconditionBlock::inverse(unsigned int i) const +PreconditionBlock::inverse(unsigned int i) const { if (same_diagonal) return var_inverse[0]; @@ -83,9 +83,9 @@ PreconditionBlock::inverse(unsigned int i) const } -template +template const FullMatrix& -PreconditionBlock::diagonal(unsigned int i) const +PreconditionBlock::diagonal(unsigned int i) const { Assert(store_diagonals, ExcDiagonalsNotStored()); @@ -97,37 +97,37 @@ PreconditionBlock::diagonal(unsigned int i) const } -template -unsigned int PreconditionBlock::block_size() const +template +unsigned int PreconditionBlock::block_size() const { return blocksize; } -template +template void -PreconditionBlock::set_same_diagonal() +PreconditionBlock::set_same_diagonal() { Assert(var_inverse.size()==0, ExcInverseMatricesAlreadyExist()); same_diagonal = true; } -template +template bool -PreconditionBlock::inverses_ready() const +PreconditionBlock::inverses_ready() const { return (var_inverse.size() != 0); } -template -void PreconditionBlock::invert_diagblocks() +template +void PreconditionBlock::invert_diagblocks() { Assert (A!=0, ExcNotInitialized()); Assert (blocksize!=0, ExcNotInitialized()); - const SparseMatrix &M=*A; + const MATRIX &M=*A; Assert (var_inverse.size()==0, ExcInverseMatricesAlreadyExist()); const unsigned int n_cells = M.m()/blocksize; @@ -147,11 +147,12 @@ void PreconditionBlock::invert_diagblocks() for (unsigned int row_cell=0; row_cell::const_iterator entry = M.begin(row_cell); - const SparseMatrix::const_iterator row_end = M.end(row_cell); - while(entry != row_end && entry->column() < blocksize) + typename MATRIX::const_iterator entry = M.begin(row_cell); + const typename MATRIX::const_iterator row_end = M.end(row_cell); + while(entry != row_end) { - M_cell(row_cell, entry->column()) = entry->value(); + if (entry->column() < blocksize) + M_cell(row_cell, entry->column()) = entry->value(); ++entry; } } @@ -202,18 +203,18 @@ void PreconditionBlock::invert_diagblocks() for (unsigned int row_cell=0; row_cell::const_iterator entry = M.begin(row); - const SparseMatrix::const_iterator row_end = M.end(row); + typename MATRIX::const_iterator entry = M.begin(row); + const typename MATRIX::const_iterator row_end = M.end(row); - while (entry != row_end && entry->column()column()column()-cell_start; if (column_cell >= blocksize) - break; + continue; M_cell(row_cell, column_cell) = entry->value(); - ++entry; } } @@ -226,9 +227,9 @@ void PreconditionBlock::invert_diagblocks() -template +template unsigned int -PreconditionBlock::memory_consumption () const +PreconditionBlock::memory_consumption () const { unsigned int mem = sizeof(*this); for (unsigned int i=0; i::memory_consumption () const /*--------------------- PreconditionBlockJacobi -----------------------*/ -template +template template -void PreconditionBlockJacobi +void PreconditionBlockJacobi ::do_vmult (Vector &dst, const Vector &src, bool adding) const @@ -267,10 +268,10 @@ void PreconditionBlockJacobi // parse error for gcc for the // exceptions that do not take // args... - typedef PreconditionBlock BaseClass; + typedef PreconditionBlock BaseClass; Assert(this->A!=0, ExcNotInitialized()); - const SparseMatrix &M=*this->A; + const MATRIX &M=*this->A; const unsigned int n_cells=M.m()/this->blocksize; Vector b_cell(this->blocksize), x_cell(this->blocksize); @@ -337,9 +338,9 @@ void PreconditionBlockJacobi } -template +template template -void PreconditionBlockJacobi +void PreconditionBlockJacobi ::vmult (Vector &dst, const Vector &src) const { @@ -347,9 +348,9 @@ void PreconditionBlockJacobi } -template +template template -void PreconditionBlockJacobi +void PreconditionBlockJacobi ::Tvmult (Vector &dst, const Vector &src) const { @@ -357,9 +358,9 @@ void PreconditionBlockJacobi } -template +template template -void PreconditionBlockJacobi +void PreconditionBlockJacobi ::vmult_add (Vector &dst, const Vector &src) const { @@ -367,9 +368,9 @@ void PreconditionBlockJacobi } -template +template template -void PreconditionBlockJacobi +void PreconditionBlockJacobi ::Tvmult_add (Vector &dst, const Vector &src) const { @@ -379,9 +380,9 @@ void PreconditionBlockJacobi /*--------------------- PreconditionBlockSOR -----------------------*/ -template +template template -void PreconditionBlockSOR::vmult (Vector &dst, +void PreconditionBlockSOR::vmult (Vector &dst, const Vector &src) const { // introduce the following typedef @@ -400,17 +401,13 @@ void PreconditionBlockSOR::vmult (Vector &ds // parse error for gcc for the // exceptions that do not take // args... - typedef PreconditionBlock BaseClass; + typedef PreconditionBlock BaseClass; Assert (this->A!=0, ExcNotInitialized()); - const SparseMatrix &M=*this->A; + const MATRIX &M=*this->A; const unsigned int n_cells=M.m()/this->blocksize; - const SparsityPattern &spars = M.get_sparsity_pattern(); - const unsigned int *rowstart = spars.get_rowstart_indices(); - const unsigned int *columns = spars.get_column_numbers(); - Vector b_cell(this->blocksize), x_cell(this->blocksize); // cell_row, cell_column are the @@ -420,7 +417,7 @@ void PreconditionBlockSOR::vmult (Vector &ds // blocks. // row, column are the global numbering // of the unkowns. - unsigned int row, column, row_cell, begin_diag_block=0; + unsigned int row, column, row_cell, block_start=0; number2 b_cell_row; if (!this->inverses_ready()) @@ -428,44 +425,55 @@ void PreconditionBlockSOR::vmult (Vector &ds FullMatrix M_cell(this->blocksize); for (unsigned int cell=0; cellblocksize, row_cell=0; - row_cellblocksize; + for (row = block_start, row_cell = 0; + row_cell < this->blocksize; ++row_cell, ++row) { + const typename MATRIX::const_iterator row_end = M.end(row); + typename MATRIX::const_iterator entry = M.begin(row); + b_cell_row=src(row); - for (unsigned int j=rowstart[row]; jcolumn(); + if (column < block_start) + b_cell_row -= entry->value() * dst(column); + else if (column < block_start + this->blocksize) + { + const unsigned int column_cell = column - block_start; + M_cell(row_cell, column_cell) = entry->value(); + } + } b_cell(row_cell)=b_cell_row; - for (unsigned int column_cell=0, column=cell*this->blocksize; - column_cellblocksize; - ++column_cell, ++column) - M_cell(row_cell,column_cell)=M(row,column); } M_cell.householder(b_cell); M_cell.backward(x_cell,b_cell); // distribute x_cell to dst - for (row=cell*this->blocksize, row_cell=0; + for (row=block_start, row_cell=0; row_cellblocksize; ++row_cell, ++row) dst(row)=this->relaxation*x_cell(row_cell); - begin_diag_block+=this->blocksize; + block_start+=this->blocksize; } } else for (unsigned int cell=0; cellblocksize, row_cell=0; + for (row = block_start, row_cell = 0; row_cellblocksize; ++row_cell, ++row) { + const typename MATRIX::const_iterator row_end = M.end(row); + typename MATRIX::const_iterator entry = M.begin(row); + b_cell_row=src(row); - for (unsigned int j=rowstart[row]; jcolumn(); + if (column < block_start) + b_cell_row -= entry->value() * dst(column); + } b_cell(row_cell)=b_cell_row; } this->inverse(cell).vmult(x_cell, b_cell); @@ -475,14 +483,14 @@ void PreconditionBlockSOR::vmult (Vector &ds ++row_cell, ++row) dst(row)=this->relaxation*x_cell(row_cell); - begin_diag_block+=this->blocksize; + block_start+=this->blocksize; } } -template +template template -void PreconditionBlockSOR::Tvmult (Vector &dst, +void PreconditionBlockSOR::Tvmult (Vector &dst, const Vector &src) const { // introduce the following typedef @@ -501,17 +509,13 @@ void PreconditionBlockSOR::Tvmult (Vector &d // parse error for gcc for the // exceptions that do not take // args... - typedef PreconditionBlock BaseClass; + typedef PreconditionBlock BaseClass; Assert (this->A!=0, ExcNotInitialized()); - const SparseMatrix &M=*this->A; + const MATRIX &M=*this->A; const unsigned int n_cells=M.m()/this->blocksize; - const SparsityPattern &spars = M.get_sparsity_pattern(); - const unsigned int *rowstart = spars.get_rowstart_indices(); - const unsigned int *columns = spars.get_column_numbers(); - Vector b_cell(this->blocksize), x_cell(this->blocksize); // cell_row, cell_column are the @@ -521,8 +525,8 @@ void PreconditionBlockSOR::Tvmult (Vector &d // blocks. // row, column are the global numbering // of the unkowns. - unsigned int row, column, row_cell; - unsigned int end_diag_block=this->blocksize *n_cells; + unsigned int row, column, row_cell, block_start = 0; + unsigned int block_end=this->blocksize *n_cells; number2 b_cell_row; if (!this->inverses_ready()) @@ -532,23 +536,26 @@ void PreconditionBlockSOR::Tvmult (Vector &d { unsigned int cell = (unsigned int) icell; // Collect upper triangle - for (row=cell*this->blocksize, row_cell=0; + for (row = block_start, row_cell = 0; row_cellblocksize; ++row_cell, ++row) { + const typename MATRIX::const_iterator row_end = M.end(row); + typename MATRIX::const_iterator entry = M.begin(row); + b_cell_row=src(row); - for (unsigned int j=rowstart[row]; j= end_diag_block) - b_cell_row -= M.global_entry(j) * dst(column); + const unsigned int column = entry->column(); + if (column >= block_end) + b_cell_row -= entry->value() * dst(column); + else if (column >= block_start) + { + const unsigned int column_cell = column - block_start; + M_cell(row_cell, column_cell) = entry->value(); + } } - b_cell(row_cell)=b_cell_row; - // Collect diagonal block - for (unsigned int column_cell=0, column=cell*this->blocksize; - column_cellblocksize; ++column_cell, ++column) - M_cell(row_cell,column_cell)=M(row,column); } M_cell.householder(b_cell); M_cell.backward(x_cell,b_cell); @@ -558,7 +565,8 @@ void PreconditionBlockSOR::Tvmult (Vector &d ++row_cell, ++row) dst(row)=this->relaxation*x_cell(row_cell); - end_diag_block-=this->blocksize; + block_end -= this->blocksize; + block_start += this->blocksize; } } else @@ -569,13 +577,16 @@ void PreconditionBlockSOR::Tvmult (Vector &d row_cellblocksize; ++row_cell, ++row) { + const typename MATRIX::const_iterator row_end = M.end(row); + typename MATRIX::const_iterator entry = M.begin(row); + b_cell_row=src(row); - for (unsigned int j=rowstart[row]; j= end_diag_block) - b_cell_row -= M.global_entry(j) * dst(column); - } + const unsigned int column = entry->column(); + if (column >= block_end) + b_cell_row -= entry->value() * dst(column); + } b_cell(row_cell)=b_cell_row; } this->inverse(cell).vmult(x_cell, b_cell); @@ -585,29 +596,30 @@ void PreconditionBlockSOR::Tvmult (Vector &d ++row_cell, ++row) dst(row)=this->relaxation*x_cell(row_cell); - end_diag_block-=this->blocksize; + block_end -= this->blocksize; + block_start += this->blocksize; } } //----------------------------------------------------------------------// -template -PreconditionBlockSSOR::PreconditionBlockSSOR () +template +PreconditionBlockSSOR::PreconditionBlockSSOR () { this->store_diagonals = 1; } -template +template template -void PreconditionBlockSSOR::vmult (Vector &dst, +void PreconditionBlockSSOR::vmult (Vector &dst, const Vector &src) const { Vector help; help.reinit(dst); - PreconditionBlockSOR::vmult(help, src); + PreconditionBlockSOR::vmult(help, src); Vector cell_src(this->blocksize); Vector cell_dst(this->blocksize); @@ -628,18 +640,18 @@ void PreconditionBlockSSOR::vmult (Vector &d help(row+row_cell)=this->relaxation*cell_dst(row_cell); } - PreconditionBlockSOR::Tvmult(dst, help); + PreconditionBlockSOR::Tvmult(dst, help); } -template +template template -void PreconditionBlockSSOR::Tvmult (Vector &dst, +void PreconditionBlockSSOR::Tvmult (Vector &dst, const Vector &src) const { Vector help; help.reinit(dst); - PreconditionBlockSOR::Tvmult(help, src); + PreconditionBlockSOR::Tvmult(help, src); Vector cell_src(this->blocksize); Vector cell_dst(this->blocksize); @@ -660,7 +672,7 @@ void PreconditionBlockSSOR::Tvmult (Vector & help(row+row_cell)=this->relaxation*cell_dst(row_cell); } - PreconditionBlockSOR::vmult(dst, help); + PreconditionBlockSOR::vmult(dst, help); } #endif diff --git a/deal.II/lac/include/lac/sparse_matrix_ez.h b/deal.II/lac/include/lac/sparse_matrix_ez.h index 20053ebdc4..1e9f9cd3d9 100644 --- a/deal.II/lac/include/lac/sparse_matrix_ez.h +++ b/deal.II/lac/include/lac/sparse_matrix_ez.h @@ -1292,7 +1292,7 @@ inline number SparseMatrixEZ::operator() (const unsigned int i, const unsigned int j) const { - Entry* entry = locate(i,j); + const Entry* entry = locate(i,j); if (entry) return entry->value; Assert(false, ExcInvalidEntry(i,j)); diff --git a/deal.II/lac/source/precondition_block.cc b/deal.II/lac/source/precondition_block.cc index 3c5a665871..744547422d 100644 --- a/deal.II/lac/source/precondition_block.cc +++ b/deal.II/lac/source/precondition_block.cc @@ -16,123 +16,123 @@ // explicit instantiations for "float" PreconditionBlock -template class PreconditionBlock; +template class PreconditionBlock, float>; -// the instantiation for class PreconditionBlock is skipped +// the instantiation for class PreconditionBlock, double> is skipped // because it does not make sense to have inverse block matrices with // higher precision than the matrix itself // explicit instantiations for "double" PreconditionBlock -template class PreconditionBlock; +template class PreconditionBlock, float>; -template class PreconditionBlock; +template class PreconditionBlock, double>; /*--------------------- PreconditionBlockJacobi -----------------------*/ // explicit instantiations for "float" PreconditionBlock -template class PreconditionBlockJacobi; +template class PreconditionBlockJacobi, float>; -template void PreconditionBlockJacobi::vmult +template void PreconditionBlockJacobi, float>::vmult (Vector &, const Vector &) const; -template void PreconditionBlockJacobi::vmult +template void PreconditionBlockJacobi, float>::vmult (Vector &, const Vector &) const; -template void PreconditionBlockJacobi::Tvmult +template void PreconditionBlockJacobi, float>::Tvmult (Vector &, const Vector &) const; -template void PreconditionBlockJacobi::Tvmult +template void PreconditionBlockJacobi, float>::Tvmult (Vector &, const Vector &) const; -template void PreconditionBlockJacobi::vmult_add +template void PreconditionBlockJacobi, float>::vmult_add (Vector &, const Vector &) const; -template void PreconditionBlockJacobi::vmult_add +template void PreconditionBlockJacobi, float>::vmult_add (Vector &, const Vector &) const; -template void PreconditionBlockJacobi::Tvmult_add +template void PreconditionBlockJacobi, float>::Tvmult_add (Vector &, const Vector &) const; -template void PreconditionBlockJacobi::Tvmult_add +template void PreconditionBlockJacobi, float>::Tvmult_add (Vector &, const Vector &) const; -template class PreconditionBlockJacobi; +template class PreconditionBlockJacobi, float>; -template void PreconditionBlockJacobi::vmult +template void PreconditionBlockJacobi, float>::vmult (Vector &, const Vector &) const; -template void PreconditionBlockJacobi::vmult +template void PreconditionBlockJacobi, float>::vmult (Vector &, const Vector &) const; -template void PreconditionBlockJacobi::Tvmult +template void PreconditionBlockJacobi, float>::Tvmult (Vector &, const Vector &) const; -template void PreconditionBlockJacobi::Tvmult +template void PreconditionBlockJacobi, float>::Tvmult (Vector &, const Vector &) const; -template void PreconditionBlockJacobi::vmult_add +template void PreconditionBlockJacobi, float>::vmult_add (Vector &, const Vector &) const; -template void PreconditionBlockJacobi::vmult_add +template void PreconditionBlockJacobi, float>::vmult_add (Vector &, const Vector &) const; -template void PreconditionBlockJacobi::Tvmult_add +template void PreconditionBlockJacobi, float>::Tvmult_add (Vector &, const Vector &) const; -template void PreconditionBlockJacobi::Tvmult_add +template void PreconditionBlockJacobi, float>::Tvmult_add (Vector &, const Vector &) const; -template class PreconditionBlockJacobi; +template class PreconditionBlockJacobi, double>; -template void PreconditionBlockJacobi::vmult +template void PreconditionBlockJacobi, double>::vmult (Vector &, const Vector &) const; -template void PreconditionBlockJacobi::vmult +template void PreconditionBlockJacobi, double>::vmult (Vector &, const Vector &) const; -template void PreconditionBlockJacobi::Tvmult +template void PreconditionBlockJacobi, double>::Tvmult (Vector &, const Vector &) const; -template void PreconditionBlockJacobi::Tvmult +template void PreconditionBlockJacobi, double>::Tvmult (Vector &, const Vector &) const; -template void PreconditionBlockJacobi::vmult_add +template void PreconditionBlockJacobi, double>::vmult_add (Vector &, const Vector &) const; -template void PreconditionBlockJacobi::vmult_add +template void PreconditionBlockJacobi, double>::vmult_add (Vector &, const Vector &) const; -template void PreconditionBlockJacobi::Tvmult_add +template void PreconditionBlockJacobi, double>::Tvmult_add (Vector &, const Vector &) const; -template void PreconditionBlockJacobi::Tvmult_add +template void PreconditionBlockJacobi, double>::Tvmult_add (Vector &, const Vector &) const; /*--------------------- PreconditionBlockGaussSeidel -----------------------*/ // explicit instantiations for "float" PreconditionBlock -template class PreconditionBlockSOR; +template class PreconditionBlockSOR, float>; -template void PreconditionBlockSOR::vmult ( +template void PreconditionBlockSOR, float>::vmult ( Vector &, const Vector &) const; -template void PreconditionBlockSOR::vmult ( +template void PreconditionBlockSOR, float>::vmult ( Vector &, const Vector &) const; -template void PreconditionBlockSOR::Tvmult ( +template void PreconditionBlockSOR, float>::Tvmult ( Vector &, const Vector &) const; -template void PreconditionBlockSOR::Tvmult ( +template void PreconditionBlockSOR, float>::Tvmult ( Vector &, const Vector &) const; -// the instantiation for class PreconditionBlockSOR is skipped +// the instantiation for class PreconditionBlockSOR, double> is skipped // because it does not make sense to have inverse block matrices with // higher precision than the matrix itself // explicit instantiations for "double" PreconditionBlockSOR -template class PreconditionBlockSOR; +template class PreconditionBlockSOR, float>; -template void PreconditionBlockSOR::vmult ( +template void PreconditionBlockSOR, float>::vmult ( Vector &, const Vector &) const; -template void PreconditionBlockSOR::vmult ( +template void PreconditionBlockSOR, float>::vmult ( Vector &, const Vector &) const; -template void PreconditionBlockSOR::Tvmult ( +template void PreconditionBlockSOR, float>::Tvmult ( Vector &, const Vector &) const; -template void PreconditionBlockSOR::Tvmult ( +template void PreconditionBlockSOR, float>::Tvmult ( Vector &, const Vector &) const; -template class PreconditionBlockSOR; +template class PreconditionBlockSOR, double>; -template void PreconditionBlockSOR::vmult ( +template void PreconditionBlockSOR, double>::vmult ( Vector &, const Vector &) const; -template void PreconditionBlockSOR::vmult ( +template void PreconditionBlockSOR, double>::vmult ( Vector &, const Vector &) const; -template void PreconditionBlockSOR::Tvmult ( +template void PreconditionBlockSOR, double>::Tvmult ( Vector &, const Vector &) const; -template void PreconditionBlockSOR::Tvmult ( +template void PreconditionBlockSOR, double>::Tvmult ( Vector &, const Vector &) const; @@ -140,43 +140,43 @@ template void PreconditionBlockSOR::Tvmult ( // explicit instantiations for "float" PreconditionBlock -template class PreconditionBlockSSOR; +template class PreconditionBlockSSOR, float>; -template void PreconditionBlockSSOR::vmult ( +template void PreconditionBlockSSOR, float>::vmult ( Vector &, const Vector &) const; -template void PreconditionBlockSSOR::vmult ( +template void PreconditionBlockSSOR, float>::vmult ( Vector &, const Vector &) const; -template void PreconditionBlockSSOR::Tvmult ( +template void PreconditionBlockSSOR, float>::Tvmult ( Vector &, const Vector &) const; -template void PreconditionBlockSSOR::Tvmult ( +template void PreconditionBlockSSOR, float>::Tvmult ( Vector &, const Vector &) const; -// the instantiation for class PreconditionBlockSSOR is skipped +// the instantiation for class PreconditionBlockSSOR, double> is skipped // because it does not make sense to have inverse block matrices with // higher precision than the matrix itself // explicit instantiations for "double" PreconditionBlockSSOR -template class PreconditionBlockSSOR; +template class PreconditionBlockSSOR, float>; -template void PreconditionBlockSSOR::vmult ( +template void PreconditionBlockSSOR, float>::vmult ( Vector &, const Vector &) const; -template void PreconditionBlockSSOR::vmult ( +template void PreconditionBlockSSOR, float>::vmult ( Vector &, const Vector &) const; -template void PreconditionBlockSSOR::Tvmult ( +template void PreconditionBlockSSOR, float>::Tvmult ( Vector &, const Vector &) const; -template void PreconditionBlockSSOR::Tvmult ( +template void PreconditionBlockSSOR, float>::Tvmult ( Vector &, const Vector &) const; -template class PreconditionBlockSSOR; +template class PreconditionBlockSSOR, double>; -template void PreconditionBlockSSOR::vmult ( +template void PreconditionBlockSSOR, double>::vmult ( Vector &, const Vector &) const; -template void PreconditionBlockSSOR::vmult ( +template void PreconditionBlockSSOR, double>::vmult ( Vector &, const Vector &) const; -template void PreconditionBlockSSOR::Tvmult ( +template void PreconditionBlockSSOR, double>::Tvmult ( Vector &, const Vector &) const; -template void PreconditionBlockSSOR::Tvmult ( +template void PreconditionBlockSSOR, double>::Tvmult ( Vector &, const Vector &) const; diff --git a/tests/lac/sparse_matrices.cc b/tests/lac/sparse_matrices.cc index 241aadbd41..ccfab17f93 100644 --- a/tests/lac/sparse_matrices.cc +++ b/tests/lac/sparse_matrices.cc @@ -21,8 +21,10 @@ #include #include #include +#include #include +#include #include #include @@ -54,6 +56,15 @@ check_vmult_quadratic(std::vector& residuals, sor.initialize(A, 1.2); PreconditionSSOR ssor; ssor.initialize(A, 1.2); + PreconditionBlockJacobi block_jacobi; + block_jacobi.initialize(A, (unsigned int) (sqrt(A.n()+.3)), 1.2); + block_jacobi.invert_diagblocks(); + PreconditionBlockSSOR block_ssor; + block_ssor.initialize(A, (unsigned int) (sqrt(A.n()+.3)), 1.2); + block_ssor.invert_diagblocks(); + PreconditionBlockSOR block_sor; + block_sor.initialize(A, (unsigned int) (sqrt(A.n()+.3)), 1.2); + block_sor.invert_diagblocks(); u = 0.; f = 1.; @@ -62,6 +73,10 @@ check_vmult_quadratic(std::vector& residuals, PREC_CHECK(prich, solve, jacobi); PREC_CHECK(prich, solve, ssor); PREC_CHECK(prich, solve, sor); + u = 0.; + PREC_CHECK(prich, solve, block_jacobi); + PREC_CHECK(prich, solve, block_ssor); + PREC_CHECK(prich, solve, block_sor); u = 0.; deallog << "Transpose" << std::endl; @@ -69,6 +84,10 @@ check_vmult_quadratic(std::vector& residuals, PREC_CHECK(prich, Tsolve, jacobi); PREC_CHECK(prich, Tsolve, ssor); PREC_CHECK(prich, Tsolve, sor); + u = 0.; + PREC_CHECK(prich, Tsolve, block_jacobi); + PREC_CHECK(prich, Tsolve, block_ssor); + PREC_CHECK(prich, Tsolve, block_sor); deallog.pop(); } @@ -110,7 +129,7 @@ int main() deallog.depth_console(1000); deallog.log_execution_time(true); deallog.log_time_differences(true); - const unsigned int size = 500; + const unsigned int size = 50; const unsigned int row_length = 9; #endif @@ -156,7 +175,7 @@ int main() testproblem.nine_point(E); check_vmult_quadratic(E_res, E, "9-SparseMatrixEZ"); - for (unsigned int i=0;i 1.e-14) deallog << "SparseMatrix and SparseMatrixEZ differ!!!" << std::endl;