double residual (Vector<number2>& w, const Vector<number2>& v, const Vector<number3>& b) const;
/**
- * Inversion of lower triangle .
+ * Forward elimination of lower triangle.
+ * Inverts the lower triangle of a
+ * quadratic matrix.
+ *
+ * If the matrix has more columns than rows,
+ * this function only operates on the left
+ * square submatrix. If there are more rows,
+ * the upper square part of the matrix
+ * is considered
*/
template<typename number2>
void forward (Vector<number2>& dst, const Vector<number2>& src) const;
/**
- * Inversion of upper triangle .
+ * Backward elimination of upper triangle.
+ * @see forward
*/
template<typename number2>
void backward (Vector<number2>& dst, const Vector<number2>& src) const;
template <typename number>
template <typename number2>
-void FullMatrix<number>::forward (Vector<number2>& dst, const Vector<number2>& src) const
+void FullMatrix<number>::forward (Vector<number2>& dst,
+ const Vector<number2>& src) const
{
- Assert(n() == m(), ExcNotQuadratic());
- Assert(dst.size() == n(), ExcDimensionMismatch(dst.size(), n()));
+ Assert(dst.size() == m(), ExcDimensionMismatch(dst.size(), m()));
Assert(src.size() == n(), ExcDimensionMismatch(src.size(), n()));
unsigned int i,j;
- unsigned int nu = (m()<n() ? m() : n());
+ unsigned int nu = ( (m()<n()) ? m() : n());
double s;
for (i=0; i<nu; ++i)
{
template <typename number>
template <typename number2>
-void FullMatrix<number>::backward (Vector<number2>& dst, const Vector<number2>& src) const
+void FullMatrix<number>::backward (Vector<number2>& dst,
+ const Vector<number2>& src) const
{
unsigned int j;
unsigned int nu = (m()<n() ? m() : n());
/**
* Initialize a rectangular matrix with
- * #m# rows and #n# columns,
- * with at most #max_per_row#
+ * #m# rows and #n# columns.
+ * The matrix may contain at most #max_per_row#
* nonzero entries per row.
*/
SparseMatrixStruct (const unsigned int m,
SparseMatrixStruct (const unsigned int n,
const unsigned int max_per_row);
+ /**
+ * Make a copy with extra off-diagonals.
+ *
+ * This constructs objects intended for
+ * the application of the ILU(n)-method.
+ * Therefore, additional to the original
+ * entry structure, space for #extra_cols#
+ * side-diagonals is provided.
+ */
+ SparseMatrixStruct(const SparseMatrixStruct& original,
+ unsigned int extra_cols);
+
/**
* Destructor.
*/
* #this#.
*/
template <typename somenumber>
- SparseMatrix<number> & copy_from (const SparseMatrix<somenumber> &);
+ SparseMatrix<number> & copy_from (const SparseMatrix<somenumber> &source);
+
+ /**
+ * Generate ILU.
+ * The matrix will entries contain the
+ * incomplete LU-factorization of
+ * the matrix #source#. Having a matrix
+ * #source# and a matrix structure object
+ * #struct#, the code for generating
+ * an ILU factorization reads
+ * \begin{verbatim}
+ * SparseMatrix<float> ilu(struct);
+ * ilu.ILU(source);
+ * \end{verbatim}
+ *
+ * If additional side diagonals are
+ * needed (ILU(n)-algorithm), you have to
+ * construct a second matrix structure:
+ * \begin{verbatim}
+ * SparseMatrixStruct ilustruct(struct,n);
+ * SparseMatrix<float> ilu(ilustruct);
+ * ilu.ILU(source);
+ * \end{verbatim}
+ *
+ * After generating the ILU-decomposition,
+ * it can be applied to a vector by
+ * #backward_forward#.
+ */
+ template <typename somenumber>
+ SparseMatrix<number> & ILU (const SparseMatrix<somenumber> &source);
/**
* Add #matrix# scaled by #factor# to this
* data type of this matrix.
*/
template <typename somenumber>
- void add_scaled (const number factor, const SparseMatrix<somenumber> &matrix);
+ void add_scaled (const number factor,
+ const SparseMatrix<somenumber> &matrix);
/**
* Return the value of the entry (i,j).
template <typename somenumber>
void Tvmult (Vector<somenumber>& dst, const Vector<somenumber>& src) const;
-
+ /**
+ * Do backward-forward solution of a
+ * previously generated ILU.
+ */
+ template <typename somenumber>
+ void backward_forward (Vector<somenumber>& v);
+
/**
* Return the norm of the vector #v# with
* respect to the norm induced by this
/**
* Exception
*/
+ DeclException0 (ExcNoILU);
+ /**
+ * Exception
+ */
DeclException0 (ExcInvalidConstructorCall);
private:
const SparseMatrixStruct * cols;
number* val;
unsigned int max_len;
+ bool is_ilu;
+
// make all other sparse matrices
// friends
SparseMatrix<number>::SparseMatrix () :
cols(0),
val(0),
- max_len(0)
+ max_len(0),
+ is_ilu(false)
{};
-
template <typename number>
SparseMatrix<number>::SparseMatrix (const SparseMatrix &m) :
cols(0),
template <typename number>
SparseMatrix<number>::SparseMatrix (const SparseMatrixStruct &c)
- : cols(&c), val(0), max_len(0)
+ : cols(&c),
+ val(0),
+ max_len(0),
+ is_ilu(false)
{
reinit();
};
if (val)
fill_n (&val[0], cols->vec_len, 0);
+
+ is_ilu = false;
}
template <typename number>
void
-SparseMatrix<number>::reinit (const SparseMatrixStruct &sparsity) {
+SparseMatrix<number>::reinit (const SparseMatrixStruct &sparsity)
+{
cols = &sparsity;
reinit ();
};
template <typename number>
void
-SparseMatrix<number>::clear () {
+SparseMatrix<number>::clear ()
+{
cols = 0;
if (val) delete[] val;
val = 0;
max_len = 0;
+ is_ilu = false;
};
template <typename number>
unsigned int
-SparseMatrix<number>::n_nonzero_elements () const {
+SparseMatrix<number>::n_nonzero_elements () const
+{
Assert (cols != 0, ExcMatrixNotInitialized());
return cols->n_nonzero_elements ();
};
template <typename number>
template <typename somenumber>
SparseMatrix<number> &
-SparseMatrix<number>::copy_from (const SparseMatrix<somenumber> &matrix) {
+SparseMatrix<number>::copy_from (const SparseMatrix<somenumber> &matrix)
+{
Assert (cols != 0, ExcMatrixNotInitialized());
Assert (val != 0, ExcMatrixNotInitialized());
Assert (cols == matrix.cols, ExcDifferentSparsityPatterns());
while (val_ptr != end_ptr)
*val_ptr++ = *matrix_ptr++;
+ is_ilu = false;
return *this;
};
template <typename somenumber>
void
SparseMatrix<number>::add_scaled (const number factor,
- const SparseMatrix<somenumber> &matrix) {
+ const SparseMatrix<somenumber> &matrix)
+{
Assert (cols != 0, ExcMatrixNotInitialized());
Assert (val != 0, ExcMatrixNotInitialized());
Assert (cols == matrix.cols, ExcDifferentSparsityPatterns());
while (val_ptr != end_ptr)
*val_ptr++ += factor * *matrix_ptr++;
+ is_ilu = false;
};
template <typename number>
template <typename somenumber>
void
-SparseMatrix<number>::precondition_Jacobi (Vector<somenumber>& dst, const Vector<somenumber>& src,
+SparseMatrix<number>::precondition_Jacobi (Vector<somenumber>& dst,
+ const Vector<somenumber>& src,
const number om) const
{
Assert (cols != 0, ExcMatrixNotInitialized());
template <typename number>
template <typename somenumber>
void
-SparseMatrix<number>::precondition_SSOR (Vector<somenumber>& dst, const Vector<somenumber>& src,
+SparseMatrix<number>::precondition_SSOR (Vector<somenumber>& dst,
+ const Vector<somenumber>& src,
const number om) const
{
// to understand how this function works