From: wolf Date: Thu, 6 May 1999 12:09:45 +0000 (+0000) Subject: Doc updates and small formatting things. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=2b0f14a8c9bdfa337a4472fc598e43d63d13b080;p=dealii-svn.git Doc updates and small formatting things. git-svn-id: https://svn.dealii.org/trunk@1287 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/Todo b/deal.II/lac/Todo index 21b97904c7..a706b31c48 100644 --- a/deal.II/lac/Todo +++ b/deal.II/lac/Todo @@ -29,4 +29,14 @@ Make ILU decomposition faster by following the comment in the 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 diff --git a/deal.II/lac/include/lac/full_matrix.h b/deal.II/lac/include/lac/full_matrix.h index 732efc0b34..523f6c6708 100644 --- a/deal.II/lac/include/lac/full_matrix.h +++ b/deal.II/lac/include/lac/full_matrix.h @@ -26,18 +26,20 @@ class iVector; * 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. * * * *
thisother @@ -52,67 +54,27 @@ class iVector; *
Instantiations provided in the library
* - * CONVENTIONS for used 'equations' :

- * - THIS matrix is always named 'A'

- * - matrices are always uppercase , vectors and scalars are lowercase

- * - 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 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); @@ -172,38 +134,40 @@ class FullMatrix : public Subscriptor * 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 - void fill (const FullMatrix& src, - const unsigned int i=0, const unsigned int j=0); + void fill (const FullMatrix &src, + const unsigned int i=0, + const unsigned int j=0); /** - * Change Dimensions. - * Set dimension to (m,n)

- * ( 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)

- * ( 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) )

- * ( 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 void reinit (const FullMatrix &B); @@ -226,66 +190,64 @@ class FullMatrix : public Subscriptor * 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

- * ( -> 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 - void add (const number s, const FullMatrix& B); + void add (const number s, + const FullMatrix &B); /** * Weighted addition of the - transpose of #B# to #this#. + * transpose of #B# to #this#. * * $A += s B^T$ */ template - void Tadd (const number s, const FullMatrix& B); + void Tadd (const number s, + const FullMatrix &B); /** * Matrix-matrix-multiplication. * $C=A*B$. */ - template - void mmult (FullMatrix& C, const FullMatrix& B) const; + void mmult (FullMatrix &C, + const FullMatrix &B) const; /** * Matrix-matrix-multiplication using @@ -293,30 +255,36 @@ class FullMatrix : public Subscriptor * $C=A^T*B$. */ template - void Tmmult (FullMatrix& C, const FullMatrix& B) const; + void Tmmult (FullMatrix &C, + const FullMatrix &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 - void vmult (Vector& w, const Vector& v, const bool adding=false) const; + void vmult (Vector &w, + const Vector &v, + const bool adding=false) const; /** - * Transpose matrix-vector-multiplication. See #vmult# above. + * Transpose matrix-vector-multiplication. + * See #vmult# above. */ template - void Tvmult (Vector& w, const Vector& v, const bool adding=false) const; + void Tvmult (Vector &w, + const Vector &v, + const bool adding=false) const; /** * Return the norm of the vector #v# with @@ -350,38 +318,48 @@ class FullMatrix : public Subscriptor * the finite element context. */ template - double matrix_scalar_product (const Vector &u, const Vector &v) const; + double matrix_scalar_product (const Vector &u, + const Vector &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 (); @@ -395,45 +373,59 @@ class FullMatrix : public Subscriptor */ 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 &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, @@ -452,11 +444,14 @@ class FullMatrix : public Subscriptor 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 - double residual (Vector& w, const Vector& v, const Vector& b) const; + double residual (Vector & w, + const Vector& v, + const Vector& b) const; /** * Forward elimination of lower triangle. @@ -470,31 +465,37 @@ class FullMatrix : public Subscriptor * is considered */ template - void forward (Vector& dst, const Vector& src) const; + void forward (Vector &dst, + const Vector &src) const; /** * Backward elimination of upper triangle. * @see forward */ template - void backward (Vector& dst, const Vector& src) const; + void backward (Vector &dst, + const Vector &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.

+ * applied to the vector y and this matrix. + * * After execution of householder, the upper - * triangle contains the resulting matrix R,

- * the lower the incomplete factorization matrices. + * triangle contains the resulting matrix R, + * the lower the incomplete factorization + * matrices. */ template - void householder (Vector& y); + void householder (Vector &y); /** - * Least - Squares - Approximation by QR-factorization. + * Least-Squares-Approximation by + * QR-factorization. */ template - double least_squares (Vector& dst, Vector& src); + double least_squares (Vector &dst, + Vector &src); /** * A(i,i)+=B(i,1-n). Addition of complete @@ -502,7 +503,8 @@ class FullMatrix : public Subscriptor * ( i = 1 ... m ) */ template - void add_diag (const number s, const FullMatrix& B); + void add_diag (const number s, + const FullMatrix &B); /** * A(i,i)+=s i=1-m. @@ -592,6 +594,70 @@ class FullMatrix : public Subscriptor * 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; }; diff --git a/deal.II/lac/include/lac/full_matrix.templates.h b/deal.II/lac/include/lac/full_matrix.templates.h index f6193a393d..998cf0ad77 100644 --- a/deal.II/lac/include/lac/full_matrix.templates.h +++ b/deal.II/lac/include/lac/full_matrix.templates.h @@ -71,7 +71,8 @@ FullMatrix::all_zero () const template void -FullMatrix::reinit (const unsigned int mm, const unsigned int nn) +FullMatrix::reinit (const unsigned int mm, + const unsigned int nn) { if (val_size::residual (Vector& dst, 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::determinant () const template -double +number FullMatrix::norm2 () const { number s = 0.; diff --git a/deal.II/lac/include/lac/ivector.h b/deal.II/lac/include/lac/ivector.h index cc2d22ee48..d263cc2b00 100644 --- a/deal.II/lac/include/lac/ivector.h +++ b/deal.II/lac/include/lac/ivector.h @@ -26,7 +26,7 @@ class iVector { friend class dFMatrix; -protected: + protected: /// Dimension. Actual number of components unsigned int dim; diff --git a/deal.II/lac/include/lac/precondition_block.h b/deal.II/lac/include/lac/precondition_block.h index 26e9e0f8ad..8eb6353cc0 100644 --- a/deal.II/lac/include/lac/precondition_block.h +++ b/deal.II/lac/include/lac/precondition_block.h @@ -16,29 +16,32 @@ * Base class for #PreconditionBlockJacobi#, #PreconditionBlockSOR#, ... * This class assumes the #SparseMatrix# 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 class PreconditionBlock: public Subscriptor @@ -66,7 +69,10 @@ 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 &M); @@ -85,7 +91,7 @@ class PreconditionBlock: public Subscriptor /** * 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) @@ -95,14 +101,22 @@ class PreconditionBlock: public Subscriptor * * 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; @@ -191,7 +205,7 @@ class PreconditionBlockSOR : public PreconditionBlock * Constructor. Takes the damping * Parameter as */ - PreconditionBlockSOR(number omega=1.); + PreconditionBlockSOR (const number omega=1.); /** * Destructor. @@ -220,7 +234,10 @@ class PreconditionBlockSOR : public PreconditionBlock void operator() (Vector&, const Vector&) const; private: - number omega; + /** + * Damping parameter. + */ + const number omega; }; @@ -229,3 +246,4 @@ class PreconditionBlockSOR : public PreconditionBlock /* end of #ifndef __precondition_block_H */ #endif /*------------------------ precondition_block.h ---------------------------*/ + diff --git a/deal.II/lac/include/lac/precondition_block.templates.h b/deal.II/lac/include/lac/precondition_block.templates.h index a15886ff34..bb747393ef 100644 --- a/deal.II/lac/include/lac/precondition_block.templates.h +++ b/deal.II/lac/include/lac/precondition_block.templates.h @@ -34,8 +34,7 @@ template void PreconditionBlock::use_matrix( const SparseMatrix &M) { - A=0; - A=&M; + A = &M; } @@ -62,7 +61,7 @@ void PreconditionBlock::invert_diagblocks() 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). @@ -72,8 +71,21 @@ void PreconditionBlock::invert_diagblocks() // row, column are the global numbering // of the unkowns. - inverse = vector > ( - n_cells, FullMatrix(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 > tmp(n_cells, + FullMatrix(blocksize)); + inverse.swap (tmp); + }; + FullMatrix M_cell(blocksize); for (unsigned int cell=0, row=0; cell::~PreconditionBlockSOR(){} template template -void PreconditionBlockSOR::operator() (Vector &dst, const Vector &src) const +void PreconditionBlockSOR::operator() (Vector &dst, + const Vector &src) const { Assert(A!=0, ExcNoMatrixGivenToUse()); const SparseMatrix &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)); diff --git a/deal.II/lac/include/lac/sparse_matrix.h b/deal.II/lac/include/lac/sparse_matrix.h index 06ff1a3bb7..88f0b6e956 100644 --- a/deal.II/lac/include/lac/sparse_matrix.h +++ b/deal.II/lac/include/lac/sparse_matrix.h @@ -26,9 +26,10 @@ class ostream; /* - - @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: @@ -337,7 +338,7 @@ class SparseMatrixStruct : public Subscriptor * 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. @@ -407,12 +408,126 @@ class SparseMatrixStruct : public Subscriptor 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 class SparseMatrix : public Subscriptor { public: /** - * Type of matrix entries. + * Type of matrix entries. In analogy to + * the STL container classes. */ typedef number value_type; @@ -765,26 +880,42 @@ class SparseMatrix : public Subscriptor * #dst#. */ template - somenumber residual (Vector& dst, const Vector& x, - const Vector& b) const; + somenumber residual (Vector &dst, + const Vector &x, + const Vector &b) const; // template - void precondition_Jacobi (Vector& dst, const Vector& src, - const number om = 1.) const; + void precondition_Jacobi (Vector &dst, + const Vector &src, + const number om = 1.) const; // template - void precondition_SSOR (Vector& dst, const Vector& src, - const number om = 1.) const; + void precondition_SSOR (Vector &dst, + const Vector &src, + const number om = 1.) const; // template - void precondition_SOR (Vector& dst, const Vector& src, - const number om = 1.) const; - // + void precondition_SOR (Vector &dst, + const Vector &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 - void SSOR (Vector& dst, const number om = 1.) const; - // + void SSOR (Vector &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 - void SOR (Vector& dst, const number om = 1.) const; + void SOR (Vector &dst, + const number om = 1.) const; /** * Return a (constant) reference to the @@ -792,7 +923,7 @@ class SparseMatrix : public Subscriptor * 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. @@ -881,8 +1012,36 @@ class SparseMatrix : public Subscriptor 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 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; diff --git a/deal.II/lac/include/lac/vector.h b/deal.II/lac/include/lac/vector.h index 5e82d24093..4b51edef3f 100644 --- a/deal.II/lac/include/lac/vector.h +++ b/deal.II/lac/include/lac/vector.h @@ -28,44 +28,14 @@ */ template 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; @@ -85,8 +55,9 @@ class Vector { Vector (); /** - * Copy-Constructor. Dimension set to that of V ,

- * all components are copied from V + * Copy-Constructor. Dimension set to + * that of V, all components are copied + * from V */ Vector (const Vector& V); @@ -94,7 +65,7 @@ class Vector { // 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 @@ -115,7 +86,7 @@ class Vector { Vector (const unsigned int n); /** - * Destructor. Clears memory + * Destructor. Clears memory */ ~Vector (); @@ -129,28 +100,28 @@ class Vector { void clear (); /** - * U(0-N) = s . Fill all components + * $U(0-N) = s$: fill all components. */ - Vector& operator= (const Number s); + Vector & operator= (const Number s); /** - * U = V . Copy all components + * $U = V$: copy all components. */ - Vector& operator= (const Vector& V); + Vector & operator= (const Vector& V); /** - * U = V for different types. + * $U = V$ for different types. */ template - Vector& operator= (const Vector& V); + Vector & operator= (const Vector& V); /** - * U = U * V . Scalar Produkt + * $U = U * V$: scalar product. */ Number operator* (const Vector& V) const; /** - * Return square of the l2-norm. + * Return square of the $l_2$-norm. */ Number norm_sqr () const; @@ -161,13 +132,13 @@ class Vector { 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. */ @@ -175,7 +146,8 @@ class Vector { /** * 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; @@ -197,7 +169,8 @@ class Vector { * 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 @@ -208,7 +181,8 @@ class Vector { * this function is the same as calling * #reinit (V.size(), fast)#. */ - void reinit (const Vector& V, const bool fast=false); + void reinit (const Vector &V, + const bool fast=false); /** * Return dimension of the vector. This @@ -263,14 +237,13 @@ class Vector { */ //@{ /** - * 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); //@} @@ -356,7 +329,7 @@ class Vector { void scale (const Number factor); /** - * U=a*V. Replacing + * U=a*V. Replacing. */ void equ (const Number a, const Vector& V); @@ -382,7 +355,8 @@ class Vector { * is undefined. No attempt is made to * catch such situations. */ - void ratio (const Vector &a, const Vector &b); + void ratio (const Vector &a, + const Vector &b); //@} @@ -465,6 +439,31 @@ class Vector { * 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; }; diff --git a/deal.II/lac/include/lac/vector_memory.h b/deal.II/lac/include/lac/vector_memory.h index 7ccb502815..f6e989c530 100644 --- a/deal.II/lac/include/lac/vector_memory.h +++ b/deal.II/lac/include/lac/vector_memory.h @@ -7,6 +7,8 @@ #include + + /** * Memory management for vectors. This class is used by all * iterative methods to allocate space for auxilliary diff --git a/deal.II/lac/source/Makefile b/deal.II/lac/source/Makefile index d59ead60ed..72204a29f6 100644 --- a/deal.II/lac/source/Makefile +++ b/deal.II/lac/source/Makefile @@ -14,7 +14,7 @@ include $D/deal.II/Make.global_options ../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 $@ @@ -24,9 +24,11 @@ lib: ../lib/liblac.g.a ../lib/liblac.a ../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: diff --git a/deal.II/lac/source/sparse_matrix.cc b/deal.II/lac/source/sparse_matrix.cc index 1bfd73193a..1027c32c9d 100644 --- a/deal.II/lac/source/sparse_matrix.cc +++ b/deal.II/lac/source/sparse_matrix.cc @@ -226,7 +226,7 @@ SparseMatrixStruct::compress () 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 @@ -324,6 +324,16 @@ SparseMatrixStruct::empty () const { 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