<h3>lac</h3>
<ol>
+ <li> <p> Changed: The classes <code
+ class="class">PreconditionBlock</code>, <code
+ class="class">PreconditionBlockJacobi</code>, <code
+ class="class">PreconditionBlockSOR</code>, and <code
+ class="class">PreconditionBlockSSOR</code> have changed their
+ template signature. The first template argument is now the matrix
+ type, not just a number type.
+ <br>
+ (GK 2002/09/18)
+ </p>
+
<li> <p> New: class <code class="class">BlockVector</code> has a
function <code class="member">collect_sizes()</code>, very much as
<code class="class">BlockSparsityPattern</code>. This allows
*
* @author Ralf Hartmann, Guido Kanschat, 1999, 2000
*/
-template<typename number, typename inverse_type = number>
+template<class MATRIX, typename inverse_type = typename MATRIX::value_type>
class PreconditionBlock : public virtual Subscriptor
{
+ private:
+ /**
+ * Define number type of matrix.
+ */
+ typedef typename MATRIX::value_type number;
+
public:
/**
* Constructor.
* parameter for derived classes
* may be provided.
*/
- void initialize (const SparseMatrix<number>& A,
+ void initialize (const MATRIX& A,
const unsigned int block_size,
const double relaxation = 1.);
* @p{vmult} function of the
* derived classes.
*/
- SmartPointer<const SparseMatrix<number> > A;
+ SmartPointer<const MATRIX> A;
/**
* Relaxation parameter to be
* used by derived classes.
* See @ref{PreconditionBlock} for requirements on the matrix.
* @author Ralf Hartmann, Guido Kanschat, 1999, 2000
*/
-template<typename number, typename inverse_type = number>
+template<class MATRIX, typename inverse_type = typename MATRIX::value_type>
class PreconditionBlockJacobi : public virtual Subscriptor,
- private PreconditionBlock<number,inverse_type>
+ private PreconditionBlock<MATRIX, inverse_type>
{
+ private:
+ /**
+ * Define number type of matrix.
+ */
+ typedef typename MATRIX::value_type number;
+
public:
/**
* Make initialization function
* publicly available.
*/
- PreconditionBlock<number,inverse_type>::initialize;
+ PreconditionBlock<MATRIX, inverse_type>::initialize;
/**
* Make function of base class public again.
*/
- PreconditionBlock<number,inverse_type>::clear;
+ PreconditionBlock<MATRIX, inverse_type>::clear;
/**
* Make function of base class public again.
*/
- PreconditionBlock<number,inverse_type>::set_same_diagonal;
+ PreconditionBlock<MATRIX, inverse_type>::set_same_diagonal;
/**
* Make function public.
*/
- PreconditionBlock<number,inverse_type>::invert_diagblocks;
+ PreconditionBlock<MATRIX, inverse_type>::invert_diagblocks;
/**
* Execute block Jacobi
* See @ref{PreconditionBlock} for requirements on the matrix.
* @author Ralf Hartmann, Guido Kanschat, 1999, 2000
*/
-template<typename number, typename inverse_type = number>
+template<class MATRIX, typename inverse_type = typename MATRIX::value_type>
class PreconditionBlockSOR : public virtual Subscriptor,
- protected PreconditionBlock<number,inverse_type>
+ protected PreconditionBlock<MATRIX, inverse_type>
{
+ private:
+ /**
+ * Define number type of matrix.
+ */
+ typedef typename MATRIX::value_type number;
+
public:
/**
* Make initialization function
* publicly available.
*/
- PreconditionBlock<number,inverse_type>::initialize;
+ PreconditionBlock<MATRIX, inverse_type>::initialize;
/**
* Make function of base class public again.
*/
- PreconditionBlock<number,inverse_type>::clear;
+ PreconditionBlock<MATRIX, inverse_type>::clear;
/**
* Make function of base class public again.
*/
- PreconditionBlock<number,inverse_type>::set_same_diagonal;
+ PreconditionBlock<MATRIX, inverse_type>::set_same_diagonal;
/**
* Make function of base class public again.
*/
- PreconditionBlock<number,inverse_type>::invert_diagblocks;
+ PreconditionBlock<MATRIX, inverse_type>::invert_diagblocks;
/**
* Execute block SOR
* See @ref{PreconditionBlock} for requirements on the matrix.
* @author Ralf Hartmann, Guido Kanschat, 1999, 2000
*/
-template<typename number, typename inverse_type = number>
+template<class MATRIX, typename inverse_type = typename MATRIX::value_type>
class PreconditionBlockSSOR : public virtual Subscriptor,
- private PreconditionBlockSOR<number,inverse_type>
+ private PreconditionBlockSOR<MATRIX, inverse_type>
{
+ private:
+ /**
+ * Define number type of matrix.
+ */
+ typedef typename MATRIX::value_type number;
+
public:
/**
* Constructor.
* Make initialization function
* publicly available.
*/
- PreconditionBlockSOR<number,inverse_type>::initialize;
+ PreconditionBlockSOR<MATRIX,inverse_type>::initialize;
/**
* Make function of base class public again.
*/
- PreconditionBlockSOR<number,inverse_type>::clear;
+ PreconditionBlockSOR<MATRIX,inverse_type>::clear;
/**
* Make function of base class public again.
*/
- PreconditionBlockSOR<number,inverse_type>::set_same_diagonal;
+ PreconditionBlockSOR<MATRIX,inverse_type>::set_same_diagonal;
/**
* Make function of base class public again.
*/
- PreconditionBlockSOR<number,inverse_type>::invert_diagblocks;
+ PreconditionBlockSOR<MATRIX,inverse_type>::invert_diagblocks;
/**
* Execute block SSOR
#include <lac/sparse_matrix.h>
-template <typename number, typename inverse_type>
-PreconditionBlock<number,inverse_type>::PreconditionBlock ():
+template <class MATRIX, typename inverse_type>
+PreconditionBlock<MATRIX,inverse_type>::PreconditionBlock ():
blocksize(0),
A(0),
store_diagonals(false),
{};
-template <typename number, typename inverse_type>
-PreconditionBlock<number,inverse_type>::~PreconditionBlock ()
+template <class MATRIX, typename inverse_type>
+PreconditionBlock<MATRIX,inverse_type>::~PreconditionBlock ()
{
if (var_inverse.size()!=0)
var_inverse.erase(var_inverse.begin(), var_inverse.end());
}
-template <typename number, typename inverse_type>
-void PreconditionBlock<number,inverse_type>::clear ()
+template <class MATRIX, typename inverse_type>
+void PreconditionBlock<MATRIX,inverse_type>::clear ()
{
if (var_inverse.size()!=0)
var_inverse.erase(var_inverse.begin(), var_inverse.end());
}
-template <typename number, typename inverse_type>
-void PreconditionBlock<number,inverse_type>::initialize (const SparseMatrix<number> &M,
+template <class MATRIX, typename inverse_type>
+void PreconditionBlock<MATRIX,inverse_type>::initialize (const MATRIX &M,
const unsigned int bsize,
const double relax)
{
}
-template <typename number, typename inverse_type>
+template <class MATRIX, typename inverse_type>
const FullMatrix<inverse_type>&
-PreconditionBlock<number,inverse_type>::inverse(unsigned int i) const
+PreconditionBlock<MATRIX,inverse_type>::inverse(unsigned int i) const
{
if (same_diagonal)
return var_inverse[0];
}
-template <typename number, typename inverse_type>
+template <class MATRIX, typename inverse_type>
const FullMatrix<inverse_type>&
-PreconditionBlock<number,inverse_type>::diagonal(unsigned int i) const
+PreconditionBlock<MATRIX,inverse_type>::diagonal(unsigned int i) const
{
Assert(store_diagonals, ExcDiagonalsNotStored());
}
-template <typename number, typename inverse_type>
-unsigned int PreconditionBlock<number,inverse_type>::block_size() const
+template <class MATRIX, typename inverse_type>
+unsigned int PreconditionBlock<MATRIX,inverse_type>::block_size() const
{
return blocksize;
}
-template <typename number, typename inverse_type>
+template <class MATRIX, typename inverse_type>
void
-PreconditionBlock<number,inverse_type>::set_same_diagonal()
+PreconditionBlock<MATRIX,inverse_type>::set_same_diagonal()
{
Assert(var_inverse.size()==0, ExcInverseMatricesAlreadyExist());
same_diagonal = true;
}
-template <typename number, typename inverse_type>
+template <class MATRIX, typename inverse_type>
bool
-PreconditionBlock<number,inverse_type>::inverses_ready() const
+PreconditionBlock<MATRIX,inverse_type>::inverses_ready() const
{
return (var_inverse.size() != 0);
}
-template <typename number, typename inverse_type>
-void PreconditionBlock<number,inverse_type>::invert_diagblocks()
+template <class MATRIX, typename inverse_type>
+void PreconditionBlock<MATRIX,inverse_type>::invert_diagblocks()
{
Assert (A!=0, ExcNotInitialized());
Assert (blocksize!=0, ExcNotInitialized());
- const SparseMatrix<number> &M=*A;
+ const MATRIX &M=*A;
Assert (var_inverse.size()==0, ExcInverseMatricesAlreadyExist());
const unsigned int n_cells = M.m()/blocksize;
for (unsigned int row_cell=0; row_cell<blocksize; ++row_cell)
{
- SparseMatrix<number>::const_iterator entry = M.begin(row_cell);
- const SparseMatrix<number>::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;
}
}
for (unsigned int row_cell=0; row_cell<blocksize; ++row_cell)
{
const unsigned int row = row_cell + cell_start;
- SparseMatrix<number>::const_iterator entry = M.begin(row);
- const SparseMatrix<number>::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()<cell_start)
- ++entry;
- while(entry != row_end)
+ for (;entry != row_end; ++entry)
{
+ if (entry->column()<cell_start)
+ continue;
+
const unsigned int column_cell = entry->column()-cell_start;
if (column_cell >= blocksize)
- break;
+ continue;
M_cell(row_cell, column_cell) = entry->value();
- ++entry;
}
}
-template <typename number, typename inverse_type>
+template <class MATRIX, typename inverse_type>
unsigned int
-PreconditionBlock<number,inverse_type>::memory_consumption () const
+PreconditionBlock<MATRIX,inverse_type>::memory_consumption () const
{
unsigned int mem = sizeof(*this);
for (unsigned int i=0; i<var_inverse.size(); ++i)
/*--------------------- PreconditionBlockJacobi -----------------------*/
-template <typename number, typename inverse_type>
+template <class MATRIX, typename inverse_type>
template <typename number2>
-void PreconditionBlockJacobi<number,inverse_type>
+void PreconditionBlockJacobi<MATRIX,inverse_type>
::do_vmult (Vector<number2> &dst,
const Vector<number2> &src,
bool adding) const
// parse error for gcc for the
// exceptions that do not take
// args...
- typedef PreconditionBlock<number,inverse_type> BaseClass;
+ typedef PreconditionBlock<MATRIX,inverse_type> BaseClass;
Assert(this->A!=0, ExcNotInitialized());
- const SparseMatrix<number> &M=*this->A;
+ const MATRIX &M=*this->A;
const unsigned int n_cells=M.m()/this->blocksize;
Vector<number2> b_cell(this->blocksize), x_cell(this->blocksize);
}
-template <typename number, typename inverse_type>
+template <class MATRIX, typename inverse_type>
template <typename number2>
-void PreconditionBlockJacobi<number,inverse_type>
+void PreconditionBlockJacobi<MATRIX,inverse_type>
::vmult (Vector<number2> &dst,
const Vector<number2> &src) const
{
}
-template <typename number, typename inverse_type>
+template <class MATRIX, typename inverse_type>
template <typename number2>
-void PreconditionBlockJacobi<number,inverse_type>
+void PreconditionBlockJacobi<MATRIX,inverse_type>
::Tvmult (Vector<number2> &dst,
const Vector<number2> &src) const
{
}
-template <typename number, typename inverse_type>
+template <class MATRIX, typename inverse_type>
template <typename number2>
-void PreconditionBlockJacobi<number,inverse_type>
+void PreconditionBlockJacobi<MATRIX,inverse_type>
::vmult_add (Vector<number2> &dst,
const Vector<number2> &src) const
{
}
-template <typename number, typename inverse_type>
+template <class MATRIX, typename inverse_type>
template <typename number2>
-void PreconditionBlockJacobi<number,inverse_type>
+void PreconditionBlockJacobi<MATRIX,inverse_type>
::Tvmult_add (Vector<number2> &dst,
const Vector<number2> &src) const
{
/*--------------------- PreconditionBlockSOR -----------------------*/
-template <typename number, typename inverse_type>
+template <class MATRIX, typename inverse_type>
template <typename number2>
-void PreconditionBlockSOR<number,inverse_type>::vmult (Vector<number2> &dst,
+void PreconditionBlockSOR<MATRIX,inverse_type>::vmult (Vector<number2> &dst,
const Vector<number2> &src) const
{
// introduce the following typedef
// parse error for gcc for the
// exceptions that do not take
// args...
- typedef PreconditionBlock<number,inverse_type> BaseClass;
+ typedef PreconditionBlock<MATRIX,inverse_type> BaseClass;
Assert (this->A!=0, ExcNotInitialized());
- const SparseMatrix<number> &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<number2> b_cell(this->blocksize), x_cell(this->blocksize);
// cell_row, cell_column are the
// 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())
FullMatrix<number> M_cell(this->blocksize);
for (unsigned int cell=0; cell<n_cells; ++cell)
{
- for (row=cell*this->blocksize, row_cell=0;
- row_cell<this->blocksize;
+ 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]; j<rowstart[row+1]; ++j)
- if ((column=columns[j]) < begin_diag_block)
- b_cell_row -= M.global_entry(j) * dst(column);
+ for (; entry != row_end; ++entry)
+ {
+ const unsigned int column = entry->column();
+ 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_cell<this->blocksize;
- ++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_cell<this->blocksize;
++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; cell<n_cells; ++cell)
{
- for (row=cell*this->blocksize, row_cell=0;
+ 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]; j<rowstart[row+1]; ++j)
- if ((column=columns[j]) < begin_diag_block)
- {
- b_cell_row -= M.global_entry(j) * dst(column);
- }
+ for (; entry != row_end; ++entry)
+ {
+ const unsigned int column = entry->column();
+ 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);
++row_cell, ++row)
dst(row)=this->relaxation*x_cell(row_cell);
- begin_diag_block+=this->blocksize;
+ block_start+=this->blocksize;
}
}
-template <typename number, typename inverse_type>
+template <class MATRIX, typename inverse_type>
template <typename number2>
-void PreconditionBlockSOR<number,inverse_type>::Tvmult (Vector<number2> &dst,
+void PreconditionBlockSOR<MATRIX,inverse_type>::Tvmult (Vector<number2> &dst,
const Vector<number2> &src) const
{
// introduce the following typedef
// parse error for gcc for the
// exceptions that do not take
// args...
- typedef PreconditionBlock<number,inverse_type> BaseClass;
+ typedef PreconditionBlock<MATRIX,inverse_type> BaseClass;
Assert (this->A!=0, ExcNotInitialized());
- const SparseMatrix<number> &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<number2> b_cell(this->blocksize), x_cell(this->blocksize);
// cell_row, cell_column are the
// 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())
{
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_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]; j<rowstart[row+1]; ++j)
+ for (; entry != row_end; ++entry)
{
- column=columns[j];
- if (column >= 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_cell<this->blocksize; ++column_cell, ++column)
- M_cell(row_cell,column_cell)=M(row,column);
}
M_cell.householder(b_cell);
M_cell.backward(x_cell,b_cell);
++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
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]; j<rowstart[row+1]; ++j)
+ for (; entry != row_end; ++entry)
{
- column=columns[j];
- if (column >= 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);
++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 <typename number, typename inverse_type>
-PreconditionBlockSSOR<number,inverse_type>::PreconditionBlockSSOR ()
+template <class MATRIX, typename inverse_type>
+PreconditionBlockSSOR<MATRIX,inverse_type>::PreconditionBlockSSOR ()
{
this->store_diagonals = 1;
}
-template <typename number, typename inverse_type>
+template <class MATRIX, typename inverse_type>
template <typename number2>
-void PreconditionBlockSSOR<number,inverse_type>::vmult (Vector<number2> &dst,
+void PreconditionBlockSSOR<MATRIX,inverse_type>::vmult (Vector<number2> &dst,
const Vector<number2> &src) const
{
Vector<number2> help;
help.reinit(dst);
- PreconditionBlockSOR<number,inverse_type>::vmult(help, src);
+ PreconditionBlockSOR<MATRIX,inverse_type>::vmult(help, src);
Vector<inverse_type> cell_src(this->blocksize);
Vector<inverse_type> cell_dst(this->blocksize);
help(row+row_cell)=this->relaxation*cell_dst(row_cell);
}
- PreconditionBlockSOR<number,inverse_type>::Tvmult(dst, help);
+ PreconditionBlockSOR<MATRIX,inverse_type>::Tvmult(dst, help);
}
-template <typename number, typename inverse_type>
+template <class MATRIX, typename inverse_type>
template <typename number2>
-void PreconditionBlockSSOR<number,inverse_type>::Tvmult (Vector<number2> &dst,
+void PreconditionBlockSSOR<MATRIX,inverse_type>::Tvmult (Vector<number2> &dst,
const Vector<number2> &src) const
{
Vector<number2> help;
help.reinit(dst);
- PreconditionBlockSOR<number,inverse_type>::Tvmult(help, src);
+ PreconditionBlockSOR<MATRIX,inverse_type>::Tvmult(help, src);
Vector<inverse_type> cell_src(this->blocksize);
Vector<inverse_type> cell_dst(this->blocksize);
help(row+row_cell)=this->relaxation*cell_dst(row_cell);
}
- PreconditionBlockSOR<number,inverse_type>::vmult(dst, help);
+ PreconditionBlockSOR<MATRIX,inverse_type>::vmult(dst, help);
}
#endif
number SparseMatrixEZ<number>::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));
// explicit instantiations for "float" PreconditionBlock
-template class PreconditionBlock<float, float>;
+template class PreconditionBlock<SparseMatrix<float>, float>;
-// the instantiation for class PreconditionBlock<float, double> is skipped
+// the instantiation for class PreconditionBlock<SparseMatrix<float>, 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<double, float>;
+template class PreconditionBlock<SparseMatrix<double>, float>;
-template class PreconditionBlock<double, double>;
+template class PreconditionBlock<SparseMatrix<double>, double>;
/*--------------------- PreconditionBlockJacobi -----------------------*/
// explicit instantiations for "float" PreconditionBlock
-template class PreconditionBlockJacobi<float, float>;
+template class PreconditionBlockJacobi<SparseMatrix<float>, float>;
-template void PreconditionBlockJacobi<float, float>::vmult<float>
+template void PreconditionBlockJacobi<SparseMatrix<float>, float>::vmult<float>
(Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockJacobi<float, float>::vmult<double>
+template void PreconditionBlockJacobi<SparseMatrix<float>, float>::vmult<double>
(Vector<double> &, const Vector<double> &) const;
-template void PreconditionBlockJacobi<float, float>::Tvmult<float>
+template void PreconditionBlockJacobi<SparseMatrix<float>, float>::Tvmult<float>
(Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockJacobi<float, float>::Tvmult<double>
+template void PreconditionBlockJacobi<SparseMatrix<float>, float>::Tvmult<double>
(Vector<double> &, const Vector<double> &) const;
-template void PreconditionBlockJacobi<float, float>::vmult_add<float>
+template void PreconditionBlockJacobi<SparseMatrix<float>, float>::vmult_add<float>
(Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockJacobi<float, float>::vmult_add<double>
+template void PreconditionBlockJacobi<SparseMatrix<float>, float>::vmult_add<double>
(Vector<double> &, const Vector<double> &) const;
-template void PreconditionBlockJacobi<float, float>::Tvmult_add<float>
+template void PreconditionBlockJacobi<SparseMatrix<float>, float>::Tvmult_add<float>
(Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockJacobi<float, float>::Tvmult_add<double>
+template void PreconditionBlockJacobi<SparseMatrix<float>, float>::Tvmult_add<double>
(Vector<double> &, const Vector<double> &) const;
-template class PreconditionBlockJacobi<double, float>;
+template class PreconditionBlockJacobi<SparseMatrix<double>, float>;
-template void PreconditionBlockJacobi<double, float>::vmult<float>
+template void PreconditionBlockJacobi<SparseMatrix<double>, float>::vmult<float>
(Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockJacobi<double, float>::vmult<double>
+template void PreconditionBlockJacobi<SparseMatrix<double>, float>::vmult<double>
(Vector<double> &, const Vector<double> &) const;
-template void PreconditionBlockJacobi<double, float>::Tvmult<float>
+template void PreconditionBlockJacobi<SparseMatrix<double>, float>::Tvmult<float>
(Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockJacobi<double, float>::Tvmult<double>
+template void PreconditionBlockJacobi<SparseMatrix<double>, float>::Tvmult<double>
(Vector<double> &, const Vector<double> &) const;
-template void PreconditionBlockJacobi<double, float>::vmult_add<float>
+template void PreconditionBlockJacobi<SparseMatrix<double>, float>::vmult_add<float>
(Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockJacobi<double, float>::vmult_add<double>
+template void PreconditionBlockJacobi<SparseMatrix<double>, float>::vmult_add<double>
(Vector<double> &, const Vector<double> &) const;
-template void PreconditionBlockJacobi<double, float>::Tvmult_add<float>
+template void PreconditionBlockJacobi<SparseMatrix<double>, float>::Tvmult_add<float>
(Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockJacobi<double, float>::Tvmult_add<double>
+template void PreconditionBlockJacobi<SparseMatrix<double>, float>::Tvmult_add<double>
(Vector<double> &, const Vector<double> &) const;
-template class PreconditionBlockJacobi<double, double>;
+template class PreconditionBlockJacobi<SparseMatrix<double>, double>;
-template void PreconditionBlockJacobi<double, double>::vmult<float>
+template void PreconditionBlockJacobi<SparseMatrix<double>, double>::vmult<float>
(Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockJacobi<double, double>::vmult<double>
+template void PreconditionBlockJacobi<SparseMatrix<double>, double>::vmult<double>
(Vector<double> &, const Vector<double> &) const;
-template void PreconditionBlockJacobi<double, double>::Tvmult<float>
+template void PreconditionBlockJacobi<SparseMatrix<double>, double>::Tvmult<float>
(Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockJacobi<double, double>::Tvmult<double>
+template void PreconditionBlockJacobi<SparseMatrix<double>, double>::Tvmult<double>
(Vector<double> &, const Vector<double> &) const;
-template void PreconditionBlockJacobi<double, double>::vmult_add<float>
+template void PreconditionBlockJacobi<SparseMatrix<double>, double>::vmult_add<float>
(Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockJacobi<double, double>::vmult_add<double>
+template void PreconditionBlockJacobi<SparseMatrix<double>, double>::vmult_add<double>
(Vector<double> &, const Vector<double> &) const;
-template void PreconditionBlockJacobi<double, double>::Tvmult_add<float>
+template void PreconditionBlockJacobi<SparseMatrix<double>, double>::Tvmult_add<float>
(Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockJacobi<double, double>::Tvmult_add<double>
+template void PreconditionBlockJacobi<SparseMatrix<double>, double>::Tvmult_add<double>
(Vector<double> &, const Vector<double> &) const;
/*--------------------- PreconditionBlockGaussSeidel -----------------------*/
// explicit instantiations for "float" PreconditionBlock
-template class PreconditionBlockSOR<float, float>;
+template class PreconditionBlockSOR<SparseMatrix<float>, float>;
-template void PreconditionBlockSOR<float, float>::vmult<float> (
+template void PreconditionBlockSOR<SparseMatrix<float>, float>::vmult<float> (
Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockSOR<float, float>::vmult<double> (
+template void PreconditionBlockSOR<SparseMatrix<float>, float>::vmult<double> (
Vector<double> &, const Vector<double> &) const;
-template void PreconditionBlockSOR<float, float>::Tvmult<float> (
+template void PreconditionBlockSOR<SparseMatrix<float>, float>::Tvmult<float> (
Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockSOR<float, float>::Tvmult<double> (
+template void PreconditionBlockSOR<SparseMatrix<float>, float>::Tvmult<double> (
Vector<double> &, const Vector<double> &) const;
-// the instantiation for class PreconditionBlockSOR<float, double> is skipped
+// the instantiation for class PreconditionBlockSOR<SparseMatrix<float>, 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<double, float>;
+template class PreconditionBlockSOR<SparseMatrix<double>, float>;
-template void PreconditionBlockSOR<double, float>::vmult<float> (
+template void PreconditionBlockSOR<SparseMatrix<double>, float>::vmult<float> (
Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockSOR<double, float>::vmult<double> (
+template void PreconditionBlockSOR<SparseMatrix<double>, float>::vmult<double> (
Vector<double> &, const Vector<double> &) const;
-template void PreconditionBlockSOR<double, float>::Tvmult<float> (
+template void PreconditionBlockSOR<SparseMatrix<double>, float>::Tvmult<float> (
Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockSOR<double, float>::Tvmult<double> (
+template void PreconditionBlockSOR<SparseMatrix<double>, float>::Tvmult<double> (
Vector<double> &, const Vector<double> &) const;
-template class PreconditionBlockSOR<double, double>;
+template class PreconditionBlockSOR<SparseMatrix<double>, double>;
-template void PreconditionBlockSOR<double, double>::vmult<float> (
+template void PreconditionBlockSOR<SparseMatrix<double>, double>::vmult<float> (
Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockSOR<double, double>::vmult<double> (
+template void PreconditionBlockSOR<SparseMatrix<double>, double>::vmult<double> (
Vector<double> &, const Vector<double> &) const;
-template void PreconditionBlockSOR<double, double>::Tvmult<float> (
+template void PreconditionBlockSOR<SparseMatrix<double>, double>::Tvmult<float> (
Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockSOR<double, double>::Tvmult<double> (
+template void PreconditionBlockSOR<SparseMatrix<double>, double>::Tvmult<double> (
Vector<double> &, const Vector<double> &) const;
// explicit instantiations for "float" PreconditionBlock
-template class PreconditionBlockSSOR<float, float>;
+template class PreconditionBlockSSOR<SparseMatrix<float>, float>;
-template void PreconditionBlockSSOR<float, float>::vmult<float> (
+template void PreconditionBlockSSOR<SparseMatrix<float>, float>::vmult<float> (
Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockSSOR<float, float>::vmult<double> (
+template void PreconditionBlockSSOR<SparseMatrix<float>, float>::vmult<double> (
Vector<double> &, const Vector<double> &) const;
-template void PreconditionBlockSSOR<float, float>::Tvmult<float> (
+template void PreconditionBlockSSOR<SparseMatrix<float>, float>::Tvmult<float> (
Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockSSOR<float, float>::Tvmult<double> (
+template void PreconditionBlockSSOR<SparseMatrix<float>, float>::Tvmult<double> (
Vector<double> &, const Vector<double> &) const;
-// the instantiation for class PreconditionBlockSSOR<float, double> is skipped
+// the instantiation for class PreconditionBlockSSOR<SparseMatrix<float>, 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<double, float>;
+template class PreconditionBlockSSOR<SparseMatrix<double>, float>;
-template void PreconditionBlockSSOR<double, float>::vmult<float> (
+template void PreconditionBlockSSOR<SparseMatrix<double>, float>::vmult<float> (
Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockSSOR<double, float>::vmult<double> (
+template void PreconditionBlockSSOR<SparseMatrix<double>, float>::vmult<double> (
Vector<double> &, const Vector<double> &) const;
-template void PreconditionBlockSSOR<double, float>::Tvmult<float> (
+template void PreconditionBlockSSOR<SparseMatrix<double>, float>::Tvmult<float> (
Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockSSOR<double, float>::Tvmult<double> (
+template void PreconditionBlockSSOR<SparseMatrix<double>, float>::Tvmult<double> (
Vector<double> &, const Vector<double> &) const;
-template class PreconditionBlockSSOR<double, double>;
+template class PreconditionBlockSSOR<SparseMatrix<double>, double>;
-template void PreconditionBlockSSOR<double, double>::vmult<float> (
+template void PreconditionBlockSSOR<SparseMatrix<double>, double>::vmult<float> (
Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockSSOR<double, double>::vmult<double> (
+template void PreconditionBlockSSOR<SparseMatrix<double>, double>::vmult<double> (
Vector<double> &, const Vector<double> &) const;
-template void PreconditionBlockSSOR<double, double>::Tvmult<float> (
+template void PreconditionBlockSSOR<SparseMatrix<double>, double>::Tvmult<float> (
Vector<float> &, const Vector<float> &) const;
-template void PreconditionBlockSSOR<double, double>::Tvmult<double> (
+template void PreconditionBlockSSOR<SparseMatrix<double>, double>::Tvmult<double> (
Vector<double> &, const Vector<double> &) const;
#include <lac/vector.h>
#include <lac/solver_richardson.h>
#include <lac/precondition.h>
+#include <lac/precondition_block.h>
#include <lac/sparse_matrix_ez.templates.h>
+#include <lac/precondition_block.templates.h>
#include <fstream>
#include <cstdio>
sor.initialize(A, 1.2);
PreconditionSSOR<MATRIX> ssor;
ssor.initialize(A, 1.2);
+ PreconditionBlockJacobi<MATRIX, float> block_jacobi;
+ block_jacobi.initialize(A, (unsigned int) (sqrt(A.n()+.3)), 1.2);
+ block_jacobi.invert_diagblocks();
+ PreconditionBlockSSOR<MATRIX, float> block_ssor;
+ block_ssor.initialize(A, (unsigned int) (sqrt(A.n()+.3)), 1.2);
+ block_ssor.invert_diagblocks();
+ PreconditionBlockSOR<MATRIX, float> block_sor;
+ block_sor.initialize(A, (unsigned int) (sqrt(A.n()+.3)), 1.2);
+ block_sor.invert_diagblocks();
u = 0.;
f = 1.;
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;
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();
}
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
testproblem.nine_point(E);
check_vmult_quadratic(E_res, E, "9-SparseMatrixEZ<double>");
- for (unsigned int i=0;i<A_res.size();++i)
+ for (unsigned int i=0;i<E_res.size();++i)
if (std::fabs(A_res[i] - E_res[i]) > 1.e-14)
deallog << "SparseMatrix and SparseMatrixEZ differ!!!"
<< std::endl;