innermost loop (about the better use of the index j).
Change in-class defined SolverSelector functions out-of-class to
- conform to general style guide. (Ralf?)
\ No newline at end of file
+ conform to general style guide. (Ralf?)
+
+Eliminate #vec_len# from the SparseMatrixStruct class. This should
+ be rather simple since as far as I can see it is only used once
+ (at all other places, it is only written) and also its value
+ should be equal to rowstart[row].
+
+Why not replace Vector::scale by Vector::operator *= ?
+
+Use the commented-out version in PreconditionBlock::invert_diagblocks
+ using the try-catch clauses
\ No newline at end of file
* Rectangular/quadratic full matrix.
*
* Implementation of a classical rectangular scheme of numbers. The
- * data type of the entries is provided in the template argument #number#.
- * The interface is quite fat and in fact has grown every time a new
- * feature was needed. So, a lot of functions are provided.
+ * data type of the entries is provided in the template argument
+ * #number#. The interface is quite fat and in fact has grown every
+ * time a new feature was needed. So, a lot of functions are provided.
*
* Since the instantiation of this template is quite an effort,
- * standard versions are precompiled into the library. These include all
- * combinations of 'float' and 'double' for matrices and vectors. If you need more
- * data types, the implementation of non-inline functions is in
- * "fullmatrix.templates.h". Driver files are in the source tree.
+ * standard versions are precompiled into the library. These include
+ * all combinations of #float# and #double# for matrices and
+ * vectors. If you need more data types, the implementation of
+ * non-inline functions is in #fullmatrix.templates.h#. Driver files
+ * are in the source tree.
*
- * Internal calculations are usually done with the accuracy of the vector argument to
- * functions. If there is no argument with a number type, the matrix number type is used.
+ * Internal calculations are usually done with the accuracy of the
+ * vector argument to functions. If there is no argument with a number
+ * type, the matrix number type is used.
*
* <TABLE BORDER=1>
* <TR><TH ALIGN=CENTER><B>this</B><TH ALIGN=CENTER><B>other
* <CAPTION>Instantiations provided in the library</CAPTION>
* </TABLE>
*
- * CONVENTIONS for used 'equations' : <p>
- * - THIS matrix is always named 'A' <p>
- * - matrices are always uppercase , vectors and scalars are lowercase <p>
- * - Transp(A) used for transpose of matrix A
+ * In the documentation of member functions, the following conventions are adopted:
+ * \begin{itemize}
+ * \item THIS matrix is always named #A#.
+ * \item Matrices are always uppercase, vectors and scalars are lowercase.
+ * \item #Transp(A)# denotes the transpose of matrix A.
+ * \end{itemize}
*
* @author Guido Kanschat, Franz-Theo Suttmeier, Wolfgang Bangerth
*/
template<typename number>
class FullMatrix : public Subscriptor
{
- private:
- /**
- * Component-array.
- */
- number* val;
- /**
- * Dimension of range. Actual number of Columns
- */
- unsigned int dim_range;
- /**
- * Dimension of image. Actual number of Rows
- */
- unsigned int dim_image;
- /**
- * Dimension. Determines amount of reserved memory
- */
- unsigned int val_size;
-
- /**
- * Initialization. Initialize
- * memory for a #FullMatrix#
- * of #m# rows and #n#
- * columns to zero.
- */
- void init (const unsigned int m, const unsigned int n);
-
- /**
- * Return a read-write reference to the
- * element #(i,j)#.
- *
- * This function does no bounds
- * checking and is only to be used
- * internally and in functions
- * already checked.
- */
- number& el (const unsigned int i, const unsigned int j);
-
- /**
- * Return the value of the element #(i,j)#.
- *
- * This function does no bounds checking and is only to be used
- * internally and in functions
- * already checked.
- */
- number el (const unsigned int i, const unsigned int j) const;
-
-
public:
/**
* Constructor. Initialize the matrix as
* a square matrix with dimension #n#.
+ *
+ * In order to avoid the implicit
+ * conversion of integers and other types
+ * to a matrix, this constructor is
+ * declared #explicit#.
*/
explicit FullMatrix (const unsigned int n = 1);
* Fill rectangular block.
*
* The matrix #src# is copied
- into the target. The optional
- values #i# and #j# determine the
- upper left corner of the image
- of #src#.
+ * into the target. The optional
+ * values #i# and #j# determine the
+ * upper left corner of the image
+ * of #src#.
*
* This function requires that
- #i+src.m()<=m()# and
- #j+src.n()<=n()#, that is, the
- image fits into the space of #this#.
+ * #i+src.m()<=m()# and
+ * #j+src.n()<=n()#, that is, the
+ * image fits into the space of #this#.
*/
template<typename number2>
- void fill (const FullMatrix<number2>& src,
- const unsigned int i=0, const unsigned int j=0);
+ void fill (const FullMatrix<number2> &src,
+ const unsigned int i=0,
+ const unsigned int j=0);
/**
- * Change Dimensions.
- * Set dimension to (m,n) <p>
- * ( reinit rectangular matrix )
+ * Set dimension to $m\times n$ and
+ * allocate memory if necessary. Forget
+ * the previous content of the matrix.
*/
- void reinit (const unsigned int m, const unsigned int n);
+ void reinit (const unsigned int m,
+ const unsigned int n);
/**
- * Change Dimensions.
- * Set dimension to (n,n) <p>
- * ( reinit quadratic matrix )
+ * Set dimension to $n\times n$ and
+ * allocate memory if necessary. Forget
+ * the previous content of the matrix.
*/
void reinit (const unsigned int n);
/**
- * Adjust Dimension.
- * Set dimension to ( m(B),n(B) ) <p>
- * ( adjust to dimensions of another matrix B )
+ * Set dimension to $m(B)\times n(B)$ and
+ * allocate memory if necessary. Forget
+ * the previous content of the matrix.
*/
template<typename number2>
void reinit (const FullMatrix<number2> &B);
* Return whether the matrix contains only
* elements with value zero. This function
* is mainly for internal consistency
- * check and should seldomly be used when
+ * checks and should seldomly be used when
* not in debug mode since it uses quite
* some time.
*/
bool all_zero () const;
- /*
- * Access Elements. returns element at relative 'address' i <p>
- * ( -> access to A(i/n , i mod n) )
- */
-// number el (const unsigned int i) const;
-
/**
* Return the value of the element #(i,j)#.
- * Does the same as the private #el(i,j)# function
- * but does bounds checking in
+ * Does the same as the private #el(i,j)#
+ * function but does bounds checking in
* debug mode.
*/
- number operator() (const unsigned int i, const unsigned int j) const;
+ number operator() (const unsigned int i,
+ const unsigned int j) const;
/**
* Return a read-write reference to
* the element #(i,j)#.
- * Does the same as the private #el(i,j)# function
- * but does bounds checking in
+ * Does the same as the private #el(i,j)#
+ * function but does bounds checking in
* debug mode.
*/
- number& operator() (const unsigned int i, const unsigned int j);
+ number& operator() (const unsigned int i,
+ const unsigned int j);
/**
* Set all entries in the matrix to
- * zero.
+ * zero. Do not resize the matrix.
*/
void clear ();
/**
- * Weighted addition. The matrix
- #s*B# is added to #this#.
+ * Weighted addition. The matrix
+ * #s*B# is added to #this#.
*
* $A += sB$
*/
template<typename number2>
- void add (const number s, const FullMatrix<number2>& B);
+ void add (const number s,
+ const FullMatrix<number2> &B);
/**
* Weighted addition of the
- transpose of #B# to #this#.
+ * transpose of #B# to #this#.
*
* $A += s B^T$
*/
template<typename number2>
- void Tadd (const number s, const FullMatrix<number2>& B);
+ void Tadd (const number s,
+ const FullMatrix<number2> &B);
/**
* Matrix-matrix-multiplication.
* $C=A*B$.
*/
-
template<typename number2>
- void mmult (FullMatrix<number2>& C, const FullMatrix<number2>& B) const;
+ void mmult (FullMatrix<number2> &C,
+ const FullMatrix<number2> &B) const;
/**
* Matrix-matrix-multiplication using
* $C=A^T*B$.
*/
template<typename number2>
- void Tmmult (FullMatrix<number2>& C, const FullMatrix<number2>& B) const;
+ void Tmmult (FullMatrix<number2> &C,
+ const FullMatrix<number2> &B) const;
/**
* Matrix-vector-multiplication.
*
* The optional parameter
* #adding# determines, whether the
- * result is stored in #w# or addet
+ * result is stored in #w# or added
* to #w#.
*
* if (adding)
- * w += A*v
+ * $w += A*v$
*
* if (!adding)
- * w = A*v
+ * $w = A*v$
*/
template<typename number2>
- void vmult (Vector<number2>& w, const Vector<number2>& v, const bool adding=false) const;
+ void vmult (Vector<number2> &w,
+ const Vector<number2> &v,
+ const bool adding=false) const;
/**
- * Transpose matrix-vector-multiplication. See #vmult# above.
+ * Transpose matrix-vector-multiplication.
+ * See #vmult# above.
*/
template<typename number2>
- void Tvmult (Vector<number2>& w, const Vector<number2>& v, const bool adding=false) const;
+ void Tvmult (Vector<number2> &w,
+ const Vector<number2> &v,
+ const bool adding=false) const;
/**
* Return the norm of the vector #v# with
* the finite element context.
*/
template<typename number2>
- double matrix_scalar_product (const Vector<number2> &u, const Vector<number2> &v) const;
+ double matrix_scalar_product (const Vector<number2> &u,
+ const Vector<number2> &v) const;
/**
- * Return the l1-norm of the matrix, i.e.
+ * Return the $l_1$-norm of the matrix, i.e.
* $|M|_1=max_{all columns j}\sum_{all
* rows i} |M_ij|$,
- * (max. sum of columns).
- * This is the
+ * (max. sum of columns). This is the
* natural matrix norm that is compatible
- * to the l1-norm for vectors, i.e.
+ * to the $l_1$-norm for vectors, i.e.
* $|Mv|_1\leq |M|_1 |v|_1$.
* (cf. Rannacher Numerik0)
*/
number l1_norm () const;
/**
- * Return the linfty-norm of the
+ * Return the $l_\infty$-norm of the
* matrix, i.e.
- * $|M|_infty=max_{all rows i}\sum_{all
- * columns j} |M_ij|$,
+ * $|M|_\infty=\max_{all rows i}\sum_{all
+ * columns j} |M_{ij}|$,
* (max. sum of rows).
* This is the
* natural matrix norm that is compatible
- * to the linfty-norm of vectors, i.e.
- * $|Mv|_infty \leq |M|_infty |v|_infty$.
+ * to the $l_\infty$-norm of vectors, i.e.
+ * $|Mv|_\infty \leq |M|_\infty |v|_\infty$.
* (cf. Rannacher Numerik0)
*/
number linfty_norm () const;
+ /**
+ * Compute the quadratic matrix norm.
+ * Return value is the root of the square
+ * sum of all matrix entries.
+ */
+ number norm2 () const;
+
/**
* A=Inverse(A). Inversion of this by
- * Gauss-Jordan-algorithm
+ * Gauss-Jordan-algorithm. Note that this
+ * is a rather expensive operation, so
+ * you may not want to use it for
+ * larger matrices if not necessary.
*/
void gauss_jordan ();
*/
double determinant () const;
- /**
- * Compute the quadratic matrix norm.
- * Return value is the root of the square
- * sum of all matrix entries.
- */
- double norm2 () const;
-
/**
* Assign the inverse of the given
* matrix to #*this#. This function is
* only implemented (hardcoded) for
* square matrices of dimension one,
- * two and three.
+ * two, three and four, since the
+ * amount of code needed grows quickly.
+ * The implementation does not use
+ * an elimination method like the
+ * Gauss-Jordan one, but rather sets
+ * the element directly; their values
+ * are precomputed symbolically using
+ * Maple. This way, we can avoid the
+ * overhead of loops and local variables
+ * but the number of lines of code
+ * grows rapidly.
+ *
+ * For all other sizes than the ones given
+ * above, an exception of type
+ * #ExcNotImplemented(dim_range)# is
+ * thrown, which you can catch and use
+ * some other method to invert the matrix,
+ * e.g. the #gauss_jordan# function.
*/
void invert (const FullMatrix<number> &M);
/**
- * A(i,1-n)+=s*A(j,1-n).
+ * $A(i,1-n)+=s*A(j,1-n)$.
* Simple addition of rows of this
*/
- void add_row (const unsigned int i, const number s, const unsigned int j);
+ void add_row (const unsigned int i,
+ const number s,
+ const unsigned int j);
/**
- * A(i,1-n)+=s*A(j,1-n)+t*A(k,1-n).
- * Multiple addition of rows of this
+ * $A(i,1-n)+=s*A(j,1-n)+t*A(k,1-n)$.
+ * Multiple addition of rows of this.
*/
void add_row (const unsigned int i,
const number s, const unsigned int j,
const number t, const unsigned int k);
/**
- * A(1-n,i)+=s*A(1-n,j).
- * Simple addition of columns of this
+ * $A(1-n,i)+=s*A(1-n,j)$.
+ * Simple addition of columns of this.
*/
- void add_col (const unsigned int i, const number s, const unsigned int j);
+ void add_col (const unsigned int i,
+ const number s,
+ const unsigned int j);
/**
- * A(1-n,i)+=s*A(1-n,j)+t*A(1-n,k).
- * Multiple addition of columns of this
+ * $A(1-n,i)+=s*A(1-n,j)+t*A(1-n,k)$.
+ * Multiple addition of columns of this.
*/
void add_col (const unsigned int i,
const number s, const unsigned int j,
void swap_col (const unsigned int i, const unsigned int j);
/**
- * w=b-A*v.
- * Residual calculation , returns |w|
+ * $w=b-A*v$.
+ * Residual calculation , returns
+ * the $l_2$-norm $|w|$
*/
template<typename number2, typename number3>
- double residual (Vector<number2>& w, const Vector<number2>& v, const Vector<number3>& b) const;
+ double residual (Vector<number2> & w,
+ const Vector<number2>& v,
+ const Vector<number3>& b) const;
/**
* Forward elimination of lower triangle.
* is considered
*/
template<typename number2>
- void forward (Vector<number2>& dst, const Vector<number2>& src) const;
+ void forward (Vector<number2> &dst,
+ const Vector<number2> &src) const;
/**
* Backward elimination of upper triangle.
* @see forward
*/
template<typename number2>
- void backward (Vector<number2>& dst, const Vector<number2>& src) const;
+ void backward (Vector<number2> &dst,
+ const Vector<number2> &src) const;
/**
- * QR - factorization of a matrix.
+ * QR-factorization of a matrix.
* The orthogonal transformation Q is
- * applied to the vector y and this matrix. <p>
+ * applied to the vector y and this matrix.
+ *
* After execution of householder, the upper
- * triangle contains the resulting matrix R, <p>
- * the lower the incomplete factorization matrices.
+ * triangle contains the resulting matrix R,
+ * the lower the incomplete factorization
+ * matrices.
*/
template<typename number2>
- void householder (Vector<number2>& y);
+ void householder (Vector<number2> &y);
/**
- * Least - Squares - Approximation by QR-factorization.
+ * Least-Squares-Approximation by
+ * QR-factorization.
*/
template<typename number2>
- double least_squares (Vector<number2>& dst, Vector<number2>& src);
+ double least_squares (Vector<number2> &dst,
+ Vector<number2> &src);
/**
* A(i,i)+=B(i,1-n). Addition of complete
* ( i = 1 ... m )
*/
template<typename number2>
- void add_diag (const number s, const FullMatrix<number2>& B);
+ void add_diag (const number s,
+ const FullMatrix<number2> &B);
/**
* A(i,i)+=s i=1-m.
* Exception
*/
DeclException0 (ExcIO);
+
+ private:
+ /**
+ * Component-array.
+ */
+ number *val;
+
+ /**
+ * Dimension of range.
+ * Actual number of Columns
+ */
+ unsigned int dim_range;
+
+ /**
+ * Dimension of image. Actual number of Rows
+ */
+ unsigned int dim_image;
+
+ /**
+ * Dimension of the array
+ * holding the values of the
+ * matrix elements. Determines
+ * amount of reserved memory.
+ *
+ * Actually, the allocated array may
+ * not have a size equal to the number
+ * of elements of this matrix, since
+ * reallocation only happens when the
+ * size of the matrix is increased.
+ * Therefore, if the matrix size was
+ * decreased somewhen in the past,
+ * #val_size# will be larger than
+ * #dim_range * dim_image#.
+ */
+ unsigned int val_size;
+
+ /**
+ * Initialization. Initialize
+ * memory for a #FullMatrix#
+ * of #m# rows and #n#
+ * columns to zero.
+ */
+ void init (const unsigned int m, const unsigned int n);
+
+ /**
+ * Return a read-write reference to the
+ * element #(i,j)#.
+ *
+ * This function does no bounds
+ * checking and is only to be used
+ * internally and in functions
+ * already checked.
+ */
+ number & el (const unsigned int i, const unsigned int j);
+
+ /**
+ * Return the value of the element #(i,j)#.
+ *
+ * This function does no bounds checking
+ * and is only to be used
+ * internally and in functions
+ * already checked.
+ */
+ number el (const unsigned int i, const unsigned int j) const;
};
template <typename number>
void
-FullMatrix<number>::reinit (const unsigned int mm, const unsigned int nn)
+FullMatrix<number>::reinit (const unsigned int mm,
+ const unsigned int nn)
{
if (val_size<nn*mm)
{
Assert(src.size() == n(), ExcDimensionMismatch(src.size(), n()));
Assert(right.size() == m(), ExcDimensionMismatch(right.size(), m()));
- unsigned int i,j;
number2 s, res = 0.;
const unsigned int size_m = m(),
size_n = n();
- for (i=0; i<size_n; ++i)
+ for (unsigned int i=0; i<size_n; ++i)
{
s = right(i);
- for (j=0; j<size_m; ++j)
+ for (unsigned int j=0; j<size_m; ++j)
s -= src(j) * el(i,j);
dst(i) = s;
res += s*s;
template <typename number>
-double
+number
FullMatrix<number>::norm2 () const
{
number s = 0.;
{
friend class dFMatrix;
-protected:
+ protected:
/// Dimension. Actual number of components
unsigned int dim;
* Base class for #PreconditionBlockJacobi#, #PreconditionBlockSOR#, ...
* This class assumes the #SparseMatrix<number># consisting of invertible blocks
* of #blocksize# on the diagonal and provides the inversion of the diagonal blocks
- * of the matrix. NOT only diagonal block matrices are allowed but all
- * matrices of arbitrary structure with the minimal property of having got
+ * of the matrix. NOT only block diagonal matrices are allowed but all
+ * matrices of arbitrary structure with the minimal property of having
* invertible blocks on the diagonal!
*
* This block matrix structure is given e.g. for the DG method
* for the transport equation. For a downstream numbering the matrices
- * even have got a lower block diagonal matrix structure, i.e. the matrices
+ * even have got a block lower left matrix structure, i.e. the matrices
* are empty above the diagonal blocks.
*
* For all matrices that are empty above and below the diagonal
* blocks (i.e. for all block diagonal matrices) the #BlockJacobi# preconditioner
* is a direct solver. For all matrices that are empty only above the diagonal blocks
- * (e.g. the matrices one get by the DG method with downstream numbering) the
+ * (e.g. the matrices one gets by the DG method with downstream numbering) the
* #BlockSOR# is a direct solver.
*
* This first implementation of the #PreconditionBlock# assumes the
- * matrix having blocks each of the same block size. Varying
+ * matrix has blocks each of the same block size. Varying
* block sizes within the matrix must still be implemented if needed.
*
* The first template parameter denotes the type of number representation in
* the sparse matrix, the second denotes the type of number representation in
- * which the inverse diagonal block
- * matrices are stored by #invert_diagblocks()#.
+ * which the inverted diagonal block matrices are stored within this class
+ * by #invert_diagblocks()#. If you don't want to use the block inversion as
+ * an exact solver, but rather as a preconditioner, you may probably want to
+ * store the inverted blocks with less accuracy than the original matrix;
+ * for example, #number==double, inverse_type=float# might be a viable choice.
*/
template<typename number, typename inverse_type>
class PreconditionBlock: public Subscriptor
/**
* Takes the matrix that should be used
- * for the preconditioning.
+ * for the preconditioning. A reference
+ * to it is stored within this class,
+ * but ownership of the matrix remains
+ * with the caller of this function.
*/
void use_matrix(const SparseMatrix<number> &M);
/**
* Stores the inverse of
- * the diagonal blocks matrices
+ * the diagonal blocks
* in #inverse#. This costs some
* additional memory - for DG
* methods about 1/3 (for double inverses)
*
* It is not allowed to call this function
* twice (will produce an error) before
- * a call of #reinit(..)#
+ * a call of #clear(..)#
* because at the second time there already
* exist the inverse matrices.
+ *
+ * After this function is called, the
+ * lock on the matrix given through the
+ * #use_matrix# function is released,
+ * i.e. you may overwrite of delete it.
+ * You may want to do this in case
+ * you use this matrix to precondition
+ * another matrix.
*/
void invert_diagblocks();
/**
- * Gives back the size of the blocks.
+ * Return the size of the blocks.
*/
unsigned int block_size () const;
* Constructor. Takes the damping
* Parameter as
*/
- PreconditionBlockSOR(number omega=1.);
+ PreconditionBlockSOR (const number omega=1.);
/**
* Destructor.
void operator() (Vector<number2>&, const Vector<number2>&) const;
private:
- number omega;
+ /**
+ * Damping parameter.
+ */
+ const number omega;
};
/* end of #ifndef __precondition_block_H */
#endif
/*------------------------ precondition_block.h ---------------------------*/
+
void PreconditionBlock<number, inverse_type>::use_matrix(
const SparseMatrix<number> &M)
{
- A=0;
- A=&M;
+ A = &M;
}
Assert (blocksize!=0, ExcBlockSizeNotSet());
Assert (M.m()%blocksize==0, ExcWrongBlockSize(blocksize, M.m()));
- unsigned int n_cells=M.m()/blocksize;
+ const unsigned int n_cells = M.m()/blocksize;
// cell_row, cell_column are the
// numbering of the blocks (cells).
// row, column are the global numbering
// of the unkowns.
- inverse = vector<FullMatrix<inverse_type> > (
- n_cells, FullMatrix<inverse_type>(blocksize));
+ // set the #inverse# array to the right
+ // size. we could do it like this:
+ // inverse = vector<>(n_cells,FullMatrix<>())
+ // but this would involve copying many
+ // FullMatrix objects.
+ //
+ // the following is a neat trick which
+ // avoids copying
+ if (true)
+ {
+ vector<FullMatrix<inverse_type> > tmp(n_cells,
+ FullMatrix<inverse_type>(blocksize));
+ inverse.swap (tmp);
+ };
+
FullMatrix<inverse_type> M_cell(blocksize);
for (unsigned int cell=0, row=0; cell<n_cells; ++cell)
template <typename number, typename inverse_type>
template <typename number2>
-void PreconditionBlockSOR<number, inverse_type>::operator() (Vector<number2> &dst, const Vector<number2> &src) const
+void PreconditionBlockSOR<number, inverse_type>::operator() (Vector<number2> &dst,
+ const Vector<number2> &src) const
{
Assert(A!=0, ExcNoMatrixGivenToUse());
const SparseMatrix<number> &M=*A;
Assert (M.m() == M.n(), ExcMatrixNotSquare());
Assert (blocksize!=0, ExcBlockSizeNotSet());
Assert (M.m()%blocksize==0, ExcWrongBlockSize(blocksize, M.m()));
- unsigned int n_cells=M.m()/blocksize;
+ const unsigned int n_cells=M.m()/blocksize;
Assert (inverse.size()==0 || inverse.size()==n_cells,
ExcWrongNumberOfInverses(inverse.size(), n_cells));
/*
-
- @author Original version by Roland Becker, Guido Kanschat, Franz-Theo Suttmeier; lots of enhancements, reorganisation and documentation by Wolfgang Bangerth
- */
+ * Structure representing the sparsity pattern of a sparse matrix.
+ *
+ * @author Original version by Roland Becker, Guido Kanschat, Franz-Theo Suttmeier; lots of enhancements, reorganisation and documentation by Wolfgang Bangerth
+ */
class SparseMatrixStruct : public Subscriptor
{
public:
* readonly.
*
* Though the return value is declared
- * #const#, you shoudl be aware that it
+ * #const#, you should be aware that it
* may change if you call any nonconstant
* function of objects which operate on
* it.
DeclException0 (ExcNotSquare);
private:
+ /**
+ * Maximum number of rows that can
+ * be stored in the #row_start# array.
+ * Since reallocation of that array
+ * only happens if the present one is
+ * too small, but never when the size
+ * of this matrix structure shrinks,
+ * #max_dim# might be larger than
+ * #rows# and in this case #row_start#
+ * has more elements than are used.
+ */
unsigned int max_dim;
- unsigned int rows, cols;
- unsigned int vec_len, max_vec_len;
+
+ /**
+ * Number of rows that this sparsity
+ * structure shall represent.
+ */
+ unsigned int rows;
+
+ /**
+ * Number of columns that this sparsity
+ * structure shall represent.
+ */
+ unsigned int cols;
+
+ /**
+ * Size of the used part of the
+ * #colnums# array. Might be lower than
+ * #max_vec_len# if the size was reduced
+ * somewhen in the past.
+ */
+ unsigned int vec_len;
+
+ /**
+ * Size of the actually allocated array
+ * #colnums#. Here, the same applies as
+ * for the #rowstart# array, i.e. it
+ * may be larger than the actually used
+ * part of the array.
+ */
+ unsigned int max_vec_len;
+
+ /**
+ * Maximum number of elements per
+ * row. This is set to the value given
+ * to the #reinit# function (or to the
+ * constructor). Its value is more
+ * or less meaningsless after #compress()#
+ * has been called.
+ */
unsigned int max_row_len;
- unsigned int* rowstart;
- int* colnums;
+
+ /**
+ * Array which hold for each row which
+ * is the first element in #colnums#
+ * belonging to that row. Note that
+ * the size of the array is one larger
+ * than the number of rows, because
+ * the last element is used for
+ * #row=rows#, i.e. the row past the
+ * last used one. The value of
+ * #rowstart[rows]# equals the index
+ * of the element past the end in
+ * #colnums#; this way, we are able to
+ * write loops like
+ * #for (i=rowstart[k]; i<rowstart[k+1]; ++i)#
+ * also for the last row.
+ *
+ * Note that the actual size of the
+ * allocated memory may be larger than
+ * the region that is used. The actual
+ * number of elements that was allocated
+ * is stored in #max_dim#.
+ */
+ unsigned int *rowstart;
+
+ /**
+ * Array of column numbers. In this array,
+ * we store for each non-zero element its
+ * column number. The column numbers for
+ * the elements in row #r# are stored
+ * within the index range
+ * #rowstart[r]...rowstart[r+1]#. Therefore
+ * to find out whether a given element
+ * #(r,c)# exists, we have to check
+ * whether the column number #c# exists in
+ * the abovementioned range within this
+ * array. If it exists, say at position
+ * #p# within this array, the value of
+ * the respective element in the sparse
+ * matrix will also be at position #p#
+ * of the values array of that class.
+ *
+ * At the beginning, all elements of
+ * this array are set to #-1# indicating
+ * invalid (unused) column numbers
+ * (however, note that if this object
+ * refers to a square matrix, the diagonal
+ * elements are preset, see below). Now, if
+ * nonzero elements are added, one column
+ * number in the row's respective range
+ * after the other is set to the column
+ * number of the added element. When
+ * compress is called, unused elements
+ * (indicated by column numbers #-1#)
+ * are eliminated by copying the column
+ * number of subsequent rows and the
+ * column numbers within each row (with
+ * the exception of the diagonal element)
+ * are sorted, such that finding whether
+ * an element exists and determining its
+ * position can be done by a binary search.
+ *
+ * If this object represents a square
+ * matrix, the first element in each
+ * row always denotes the diagonal
+ * element, i.e.
+ * #colnums[rowstart[r]]==r#.
+ */
+ int *colnums;
/**
* Store whether the #compress# function
/*
-CLASS
- SparseMatrix
-
- @author Original version by Roland Becker, Guido Kanschat, Franz-Theo Suttmeier; lots of enhancements, reorganisation and documentation by Wolfgang Bangerth 1998
- */
+ * Sparse matrix.
+ *
+ * @author Original version by Roland Becker, Guido Kanschat, Franz-Theo Suttmeier; lots of enhancements, reorganisation and documentation by Wolfgang Bangerth 1998
+ */
template <typename number>
class SparseMatrix : public Subscriptor
{
public:
/**
- * Type of matrix entries.
+ * Type of matrix entries. In analogy to
+ * the STL container classes.
*/
typedef number value_type;
* #dst#.
*/
template <typename somenumber>
- somenumber residual (Vector<somenumber>& dst, const Vector<somenumber>& x,
- const Vector<somenumber>& b) const;
+ somenumber residual (Vector<somenumber> &dst,
+ const Vector<somenumber> &x,
+ const Vector<somenumber> &b) const;
//
template <typename somenumber>
- void precondition_Jacobi (Vector<somenumber>& dst, const Vector<somenumber>& src,
- const number om = 1.) const;
+ void precondition_Jacobi (Vector<somenumber> &dst,
+ const Vector<somenumber> &src,
+ const number om = 1.) const;
//
template <typename somenumber>
- void precondition_SSOR (Vector<somenumber>& dst, const Vector<somenumber>& src,
- const number om = 1.) const;
+ void precondition_SSOR (Vector<somenumber> &dst,
+ const Vector<somenumber> &src,
+ const number om = 1.) const;
//
template <typename somenumber>
- void precondition_SOR (Vector<somenumber>& dst, const Vector<somenumber>& src,
- const number om = 1.) const;
- //
+ void precondition_SOR (Vector<somenumber> &dst,
+ const Vector<somenumber> &src,
+ const number om = 1.) const;
+
+ /**
+ * Perform an SSOR step in-place, i.e.
+ * without copying to a second vector.
+ * #omega# is the damping parameter.
+ */
template <typename somenumber>
- void SSOR (Vector<somenumber>& dst, const number om = 1.) const;
- //
+ void SSOR (Vector<somenumber> &dst,
+ const number omega = 1.) const;
+
+ /**
+ * Perform an SOR step in-place, i.e.
+ * without copying to a second vector.
+ * #omega# is the damping parameter.
+ */
template <typename somenumber>
- void SOR (Vector<somenumber>& dst, const number om = 1.) const;
+ void SOR (Vector<somenumber> &dst,
+ const number om = 1.) const;
/**
* Return a (constant) reference to the
* matrix.
*
* Though the return value is declared
- * #const#, you shoudl be aware that it
+ * #const#, you should be aware that it
* may change if you call any nonconstant
* function of objects which operate on
* it.
DeclException0 (ExcInvalidConstructorCall);
private:
+ /**
+ * Pointer to the sparsity pattern used
+ * for this matrix. In order to guarantee
+ * that it is not deleted while still in
+ * use, we subscribe to it using the
+ * #SmartPointer# class.
+ */
SmartPointer<const SparseMatrixStruct> cols;
- number* val;
+
+ /**
+ * Array of values for all the nonzero
+ * entries. The position within the matrix,
+ * i.e. the row and column number for a
+ * given entry can only be deduced using
+ * the sparsity pattern. The same holds
+ * for the more common operation of
+ * finding an entry by its coordinates.
+ */
+ number *val;
+
+ /**
+ * Allocated size of #val#. This can
+ * be larger than the actually used part
+ * if the size of the matrix was
+ * reduced somewhen in the past by
+ * associating a sparsity pattern
+ * with a smaller size to this object
+ * somewhen, using the #reinit#
+ * function.
+ */
unsigned int max_len;
*/
template <typename Number>
class Vector {
- protected:
-
- /**
- * Dimension. Actual number of components
- * contained in the vector.
- * Get this number by calling #size()#.
- */
- unsigned int dim;
-
- /**
- * Amount of memory actually reserved for
- * this vector. This number may be greater
- * than #dim# if a #reinit# was called with
- * less memory requirements than the vector
- * needed last time. At present #reinit#
- * does not free memory when the number of
- * needed elements is reduced.
- */
- unsigned int maxdim;
-
- /**
- * Pointer to the array of components.
- */
- Number *val;
-
public:
-
- /**
- * Declare iterator types just like those
- * for the C++ standard library:
- * Data type stored by this container.
- */
- typedef Number value_type;
-
/**
* Declare standard types used in all
- * containers.
+ * containers. These types parallel
+ * those in the #C++# standard libraries
+ * #vector<...># class.
*/
+ typedef Number value_type;
typedef value_type* pointer;
typedef const value_type* const_pointer;
typedef value_type* iterator;
Vector ();
/**
- * Copy-Constructor. Dimension set to that of V , <p>
- * all components are copied from V
+ * Copy-Constructor. Dimension set to
+ * that of V, all components are copied
+ * from V
*/
Vector (const Vector<Number>& V);
// note: I disabled this function for the time being, since egcs1.1.2
// does not respect the "explicit" keyword for template constructors.
// this leads to unwanted conversions and in some places to automatically
-// generated temporaries, where this is not a good idea
+// generated temporaries, where this is not a good idea. [WB]
// /**
// * Copy constructor taking a vector of
// * another data type. This will fail if
Vector (const unsigned int n);
/**
- * Destructor. Clears memory
+ * Destructor. Clears memory
*/
~Vector ();
void clear ();
/**
- * U(0-N) = s . Fill all components
+ * $U(0-N) = s$: fill all components.
*/
- Vector<Number>& operator= (const Number s);
+ Vector<Number> & operator= (const Number s);
/**
- * U = V . Copy all components
+ * $U = V$: copy all components.
*/
- Vector<Number>& operator= (const Vector<Number>& V);
+ Vector<Number> & operator= (const Vector<Number>& V);
/**
- * U = V for different types.
+ * $U = V$ for different types.
*/
template<typename Number2>
- Vector<Number>& operator= (const Vector<Number2>& V);
+ Vector<Number> & operator= (const Vector<Number2>& V);
/**
- * U = U * V . Scalar Produkt
+ * $U = U * V$: scalar product.
*/
Number operator* (const Vector<Number>& V) const;
/**
- * Return square of the l2-norm.
+ * Return square of the $l_2$-norm.
*/
Number norm_sqr () const;
Number mean_value () const;
/**
- * Return the l1-norm of the vector, i.e.
+ * Return the $l_1$-norm of the vector, i.e.
* the sum of the absolute values.
*/
Number l1_norm () const;
/**
- * Return the l2-norm of the vector, i.e.
+ * Return the $l_2$-norm of the vector, i.e.
* the square root of the sum of the
* squares of the elements.
*/
/**
* Return the maximum absolute value of the
- * elements of this vector.
+ * elements of this vector, which is the
+ * $l_\infty$-norm of a vector.
*/
Number linfty_norm () const;
* On #fast==false#, the vector is filled by
* zeros.
*/
- void reinit (const unsigned int N, const bool fast=false);
+ void reinit (const unsigned int N,
+ const bool fast=false);
/**
* Change the dimension to that of the
* this function is the same as calling
* #reinit (V.size(), fast)#.
*/
- void reinit (const Vector<Number>& V, const bool fast=false);
+ void reinit (const Vector<Number> &V,
+ const bool fast=false);
/**
* Return dimension of the vector. This
*/
//@{
/**
- * Access Components. returns U(i) ,
- * INLINE
+ * Access components, returns U(i).
*/
Number operator() (const unsigned int i) const;
/**
- * Access Components. returns U(i) ,
- * INLINE
+ * Access components, returns U(i)
+ * as a writeable reference.
*/
Number& operator() (const unsigned int i);
//@}
void scale (const Number factor);
/**
- * U=a*V. Replacing
+ * U=a*V. Replacing.
*/
void equ (const Number a, const Vector<Number>& V);
* is undefined. No attempt is made to
* catch such situations.
*/
- void ratio (const Vector<Number> &a, const Vector<Number> &b);
+ void ratio (const Vector<Number> &a,
+ const Vector<Number> &b);
//@}
* Exception
*/
DeclException0 (ExcIO);
+
+ protected:
+
+ /**
+ * Dimension. Actual number of components
+ * contained in the vector.
+ * Get this number by calling #size()#.
+ */
+ unsigned int dim;
+
+ /**
+ * Amount of memory actually reserved for
+ * this vector. This number may be greater
+ * than #dim# if a #reinit# was called with
+ * less memory requirements than the vector
+ * needed last time. At present #reinit#
+ * does not free memory when the number of
+ * needed elements is reduced.
+ */
+ unsigned int maxdim;
+
+ /**
+ * Pointer to the array of components.
+ */
+ Number *val;
};
#include <base/subscriptor.h>
+
+
/**
* Memory management for vectors. This class is used by all
* iterative methods to allocate space for auxilliary
../lib/o/%.go :
@echo ============================ Compiling with debugging information: $<
- $(CXX) $(CXXFLAGS.g) -c $< -o $@
+ @$(CXX) $(CXXFLAGS.g) -c $< -o $@
../lib/o/%.o :
@echo ============================ Compiling with optimization: $<
@$(CXX) $(CXXFLAGS) -c $< -o $@
../lib/liblac.a: $(forward-declarations) $(o-files)
+ @echo ======================================== Updating library: $@
@ar ruv $@ $(o-files)
../lib/liblac.g.a: $(forward-declarations) $(go-files)
+ @echo ======================================== Updating library: $@
@ar ruv $@ $(go-files)
clean:
row_length = 0;
// reserve temporary storage to
- // store the entries of one wor
+ // store the entries of one row
int *tmp_entries = new int[max_row_len];
// Traverse all rows
unsigned int
SparseMatrixStruct::max_entries_per_row () const
{
+ // if compress() has not yet been
+ // called, we can get the maximum
+ // number of elements per row using
+ // the stored value
+ if (!compressed)
+ return max_row_len;
+
+ // if compress() was called, we
+ // use a better algorithm which
+ // gives us a sharp bound
unsigned int m = 0;
for (unsigned int i=1; i<rows; ++i)
m = max (m, rowstart[i]-rowstart[i-1]);