*/
template <typename number2>
void Tvmult (Vector<number2>&, const Vector<number2>&) const;
+ /**
+ * Execute block Jacobi
+ * preconditioning, adding to @p{dst}.
+ *
+ * This function will
+ * automatically use the inverse
+ * matrices if they exist, if not
+ * then BlockJacobi will need
+ * much time inverting the
+ * diagonal block matrices in
+ * each preconditioning step.
+ */
+ template <typename number2>
+ void vmult_add (Vector<number2>&, const Vector<number2>&) const;
+
+ /**
+ * Same as @p{vmult_add}, since Jacobi is symmetric.
+ */
+ template <typename number2>
+ void Tvmult_add (Vector<number2>&, const Vector<number2>&) const;
+
+ private:
+ /**
+ * Actual implementation of the
+ * preconditioner.
+ */
+ template <typename number2>
+ void do_vmult (Vector<number2>&,
+ const Vector<number2>&,
+ bool adding) const;
};
template <typename number, typename inverse_type>
template <typename number2>
void PreconditionBlockJacobi<number,inverse_type>
-::vmult (Vector<number2> &dst,
- const Vector<number2> &src) const
+::do_vmult (Vector<number2> &dst,
+ const Vector<number2> &src,
+ bool adding) const
{
// introduce the following typedef
// since in the use of exceptions,
M_cell.backward(x_cell,b_cell);
// distribute x_cell to dst
for (row=cell*blocksize, row_cell=0; row_cell<blocksize; ++row_cell, ++row)
- dst(row)=x_cell(row_cell);
+ if (adding)
+ dst(row)+=x_cell(row_cell);
+ else
+ dst(row)=x_cell(row_cell);
begin_diag_block+=blocksize;
}
inverse(cell).vmult(x_cell, b_cell);
// distribute x_cell to dst
for (row=cell*blocksize, row_cell=0; row_cell<blocksize; ++row_cell, ++row)
- dst(row)=x_cell(row_cell);
+ if (adding)
+ dst(row)+=x_cell(row_cell);
+ else
+ dst(row)=x_cell(row_cell);
begin_diag_block+=blocksize;
}
template <typename number, typename inverse_type>
template <typename number2>
void PreconditionBlockJacobi<number,inverse_type>
-::Tvmult (Vector<number2> &dst,
+::vmult (Vector<number2> &dst,
const Vector<number2> &src) const
{
- vmult(dst, src);
+ do_vmult(dst, src, false);
+}
+
+
+template <typename number, typename inverse_type>
+template <typename number2>
+void PreconditionBlockJacobi<number,inverse_type>
+::Tvmult (Vector<number2> &dst,
+ const Vector<number2> &src) const
+{
+ do_vmult(dst, src, false);
+}
+
+
+template <typename number, typename inverse_type>
+template <typename number2>
+void PreconditionBlockJacobi<number,inverse_type>
+::vmult_add (Vector<number2> &dst,
+ const Vector<number2> &src) const
+{
+ do_vmult(dst, src, true);
+}
+
+
+template <typename number, typename inverse_type>
+template <typename number2>
+void PreconditionBlockJacobi<number,inverse_type>
+::Tvmult_add (Vector<number2> &dst,
+ const Vector<number2> &src) const
+{
+ do_vmult(dst, src, true);
}
(Vector<float> &, const Vector<float> &) const;
template void PreconditionBlockJacobi<float, float>::Tvmult
(Vector<double> &, const Vector<double> &) const;
+template void PreconditionBlockJacobi<float, float>::vmult_add
+(Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockJacobi<float, float>::vmult_add
+(Vector<double> &, const Vector<double> &) const;
+template void PreconditionBlockJacobi<float, float>::Tvmult_add
+(Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockJacobi<float, float>::Tvmult_add
+(Vector<double> &, const Vector<double> &) const;
template class PreconditionBlockJacobi<double, float>;
(Vector<float> &, const Vector<float> &) const;
template void PreconditionBlockJacobi<double, float>::Tvmult
(Vector<double> &, const Vector<double> &) const;
+template void PreconditionBlockJacobi<double, float>::vmult_add
+(Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockJacobi<double, float>::vmult_add
+(Vector<double> &, const Vector<double> &) const;
+template void PreconditionBlockJacobi<double, float>::Tvmult_add
+(Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockJacobi<double, float>::Tvmult_add
+(Vector<double> &, const Vector<double> &) const;
template class PreconditionBlockJacobi<double, double>;
(Vector<float> &, const Vector<float> &) const;
template void PreconditionBlockJacobi<double, double>::Tvmult
(Vector<double> &, const Vector<double> &) const;
+template void PreconditionBlockJacobi<double, double>::vmult_add
+(Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockJacobi<double, double>::vmult_add
+(Vector<double> &, const Vector<double> &) const;
+template void PreconditionBlockJacobi<double, double>::Tvmult_add
+(Vector<float> &, const Vector<float> &) const;
+template void PreconditionBlockJacobi<double, double>::Tvmult_add
+(Vector<double> &, const Vector<double> &) const;
/*--------------------- PreconditionBlockGaussSeidel -----------------------*/