From 79867d06c34994d825689e9ee9fa26c3337a9d2e Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Mon, 9 Mar 1998 16:32:30 +0000 Subject: [PATCH] Kommentare von Klaus git-svn-id: https://svn.dealii.org/trunk@44 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/dfmatrix.h | 406 +++++++++++++++++++---------- 1 file changed, 275 insertions(+), 131 deletions(-) diff --git a/deal.II/lac/include/lac/dfmatrix.h b/deal.II/lac/include/lac/dfmatrix.h index 0f355658d8..325ce5082b 100644 --- a/deal.II/lac/include/lac/dfmatrix.h +++ b/deal.II/lac/include/lac/dfmatrix.h @@ -13,142 +13,286 @@ #include #endif -/// Double precision full matrix +/** + * Double precision full Matrix. + * Memory for Components is supplied explicitly

+ * ( ! Amount of memory needs not to comply with actual dimension due to reinitializations ! )

+ * - all necessary methods for matrices are supplied

+ * - operators available are `=` and `( )`

+ * 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 + * + */ class dFMatrix { - double* val; - int dim_range, dim_image, val_size; - void init(int n, int m); -// dFMatrix(const dFMatrix&); - - double& el(int i, int j) { return val[i*dim_range+j]; } - double el(int i, int j) const { return val[i*dim_range+j]; } - double el(int i) const { return val[i]; } - public: - /// Number of rows - int m() const { return dim_image; } - /// Number of columns - int n() const { return dim_range; } - - - /// copy constructor. Be very careful with this constructor, since - // it may take a hige amount of computing time for large matrices!! - dFMatrix(const dFMatrix&); - /// Constructor for quadratic n x n matrices - dFMatrix(int n = 1) { init(n,n); } - /// Constructor for rectangular matrices with m rows and n columns. - dFMatrix(int m,int n) { init(m,n); } - /// Destructor. - ~dFMatrix(); - - /// Reinitialization of rectangular matrix. - void reinit(int m, int n); - /// Reinitialization of quadratic matrix - void reinit(int n) { reinit(n,n); } - /// Reinitialization to the same dimensions of another matrix. - void reinit(const dFMatrix& A) { reinit(A.m(), A.n()); } - - /// Read access to matrix elements. - double operator() (int i, int j) const - { - //THROW2((i<0) || (i>=dim_image), IntError(IntError::Range,i)); - //THROW2((j<0) || (j>=dim_range), IntError(IntError::Range,j)); - return el(i,j); - } - - /// Read-Write access to matrix elements - double& operator() (int i, int j) - { - //THROW2((i<0) || (i>=dim_image), IntError(IntError::Range,i)); - //THROW2((j<0) || (j>=dim_range), IntError(IntError::Range,j)); - return el(i,j); - } - - /// Assignment operator - dFMatrix& operator = (const dFMatrix&); - /** Copying matrix entries. - * Copy the elements of matrix src into this beginning at - * element (i,j) - */ - - void fill(const dFMatrix& src, int i = 0, int j = 0); - void add(double s, const dFMatrix& src); - void add_diag(double s, const dFMatrix& src); - void Tadd(double s, const dFMatrix& src); - - /// - /* Add different rows of a matrix - * a(i,.) += s * a(j,.) - */ - void add_row(int i, double s, int j); - /// - /** - * Add different rows of a matrix - * a(i,.) += s * a(j,.) + t * a(k,.) - */ - void add_row(int i, double s, int j, double t, int k); - /// - /** - * Add different columns of a matrix - * a(.,i) += s * a(.,j) - */ - void add_col(int i, double s, int j); - /// - /** - * Add different columns of a matrix - * a(.,i) += s * a(.,j) + t * a(.,k) - */ - void add_col(int i, double s, int j, double t, int k); - /// Exchange contents of rows i and j - void swap_row(int i, int j); - /// Exchange contents of columns i and j - void swap_col(int i, int j); - /// Adding a scalar value on the diagonal - void diagadd(const double& src); - - /// Matrix-matrix-multiplication dst = this * src - void mmult(dFMatrix& dst, const dFMatrix& src) const; - /// - void Tmmult(dFMatrix& dst, const dFMatrix& src) const; - /// - /* - * Application of a matrix to a vector. - * - * [flag] adding determines if the result is copied to dst or added. - * - * dst (+)= this * src - */ - void vmult(dVector& dst, const dVector& src,const int adding = 0) const; - /// - void gsmult(dVector& dst, const dVector& src,const iVector& gl) const; - /// Application of the transpose matrix to a vector. - void Tvmult(dVector& dst,const dVector& src,const int adding=0) const; - /// - double residual(dVector& dst, const dVector& src, const dVector& right) const; - /// Inversion of lower triangle - void forward(dVector& dst, const dVector& src) const; - /// Inversion of upper triangle - void backward(dVector& dst, const dVector& src) const; - /// Replace this by its inverse matrix calculated with Gauß-Jordan algorithm. - void gauss_jordan(); - - /// - /* - * QR - factorization of a matrix. - * The orthogonal transformation Q is applied to the vector y and the - * matrix. After execution of householder, the upper triangle contains - * the resulting matrix R, the lower the incomplete factorization matrices. - */ - void householder(dVector& y); - - /// Least - Squares - Approximation by QR-factorization. - double least_squares(dVector& dst, dVector& src); - - /// Output of the matrix in user-defined format. - void print(FILE* fp, const char* format = 0) const; + /// Component-array. +double* val; + /// Dimension. Actual number of Columns +int dim_range; + /// Dimension. Actual number of Rows +int dim_image; + /// Dimension. Determines amount of reserved memory +int val_size; + /** + * Initialization . initialize memory for Matrix

+ * ( m rows , n columns ) + */ + void init(int m, int n); /** + * Access Elements. returns A(i,j) + */ + double& el(int i, int j) { return val[i*dim_range+j]; } + + /** + * Access Elements. returns A(i,j) + */ + double el(int i, int j) const { return val[i*dim_range+j]; } + + + public: + /// copy constructor. Be very careful with this constructor, since + // it may take a hige amount of computing time for large matrices!! + + /**@name 1: Basic Object-handling */ +//@{ + /** + * Constructor. Dimension = (n,n)

+ * -> quadratic matrix (n rows , n columns) + */ + dFMatrix(int n = 1) { init(n,n); } + + /** + * Constructor. Dimension = (m,n)

+ * -> rectangular matrix (m rows , n columns) + */ + dFMatrix(int m,int n) { init(m,n); } + + /** + * Copy constructor. Be very careful with this constructor, since + * it may take a high amount of computing time for large matrices!! + */ + dFMatrix(const dFMatrix&); + + /** + * Destructor. Clears memory + */ + ~dFMatrix(); + + /** + * A = B . Copy all elements + */ + dFMatrix& operator = (const dFMatrix& B) + + /** + * U(0-m,0-n) = s . Fill all elements + */ +void fill(const dFMatrix& src, int i = 0, int j = 0); + + /** + * Change Dimension. + * Set dimension to (m,n)

+ * ( reinit rectangular matrix ) + */ + void reinit(int m, int n); + + /** + * Change Dimension. + * Set dimension to (n,n)

+ * ( reinit quadratic matrix ) + */ + void reinit(int n) { reinit(n,n); } + + /** + * Adjust Dimension. + * Set dimension to ( m(B),n(B) )

+ * ( adjust to dimensions of another matrix B ) + */ + void reinit(const dFMatrix& B) { reinit(B.m(), B.n()); } + + /** + * Inquire Dimension (Row) . returns Number of Rows + */ + int m() const { return dim_image; } + + /** + * Inquire Dimension (Col) . returns Number of Columns + */ + int n() const { return dim_range; } +//@} + + + /**@name 2: Data-Access + */ +//@{ + /** + * Access Elements. returns element at relative 'address' i

+ * ( -> access to A(i/n , i mod n) ) + */ + double el(int i) const { return val[i]; } + + /** + * Access Elements. returns A(i,j) + */ + double operator() (int i, int j) const + { + //THROW2((i<0) || (i>=dim_image), IntError(IntError::Range,i)); + //THROW2((j<0) || (j>=dim_range), IntError(IntError::Range,i)); + return el(i,j); + } + + /** + * Access Elements. returns A(i,j) + */ + double& operator() (int i, int j) + { + //THROW2((i<0) || (i>=dim_image), IntError(IntError::Range,i)); + //THROW2((j<0) || (j>=dim_range), IntError(IntError::Range,i)); + return el(i,j); + } +//@} + + + /**@name 3: Basic applications on matrices + */ +//@{ + /** + * A+=B . Simple addition + */ + void add(double s, const dFMatrix& B); + + /** + * A+=Transp(B) . Simple addition of the transpose of B to this + */ + void Tadd(double s, const dFMatrix& B); + + /** + * C=A*B . Matrix-matrix-multiplication + */ + void mmult(dFMatrix& C, const dFMatrix& B) const; + + /** + * C=Transp(A)*B . Matrix-matrix-multiplication using transpose of this + */ + void Tmmult(dFMatrix& C, const dFMatrix& B) const; + + /** + * w (+)= A*v . Matrix-vector-multiplication ;

+ * ( application of this to a vector v ) + * flag adding=true : w=+A*v + */ + void vmult(dVector& w, const dVector& v,const int adding = 0) const; + + /** + * w (+)= Transp(A)*v . Matrix-vector-multiplication ;

+ * ( application of transpose of this to a vector v ) + * flag adding=true : w=+A*v + */ + void Tvmult(dVector& w,const dVector& v,const int adding=0) const; + + /** + * A=Inverse(A). Inversion of this by Gauss-Jordan-algorithm + */ + void gauss_jordan(); + //@} + + + /**@name 4: Basic applications on Rows or Columns + */ + //@{ + /** + * A(i,1-n)+=s*A(j,1-n). + * Simple addition of rows of this + */ + void add_row(int i, double s, int j); + + /** + * A(i,1-n)+=s*A(j,1-n)+t*A(k,1-n) . Multiple addition of rows of this + */ + void add_row(int i, double s, int j, double t, int k); + + /** + * A(1-n,i)+=s*A(1-n,j) . Simple addition of columns of this + */ + void add_col(int i, double s, int j); + + /** + * A(1-n,i)+=s*A(1-n,j)+t*A(1-n,k) . Multiple addition of columns of this + */ + void add_col(int i, double s, int j, double t, int k); + + /** + * Swap A(i,1-n) <-> A(j,1-n) . Swap rows i and j of this + */ + void swap_row(int i, int j); + + /** + * Swap A(1-n,i) <-> A(1-n,j) . Swap columns i and j of this + */ + void swap_col(int i, int j); +//@} + + + /**@name 5: Mixed stuff . Including more applications on matrices + */ + //@{ + /** + * w=b-A*v . Residual calculation , returns |w| + */ + double residual(dVector& w, const dVector& v, const dVector& b) const; + + /** + * Inversion of lower triangle . + */ + void forward(dVector& dst, const dVector& src) const; + + /** + * Inversion of upper triangle . + */ + void backward(dVector& dst, const dVector& src) const; + + /** + * QR - factorization of a matrix. + * The orthogonal transformation Q is 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. + */ + void householder(dVector& y); + + /** + * Least - Squares - Approximation by QR-factorization. + */ + double least_squares(dVector& dst, dVector& src); + + /** + * A(i,i)+=B(i,1-n). Addition of complete rows of B to diagonal-elements of this ;

+ * ( i = 1 ... m ) + */ + void add_diag(double s, const dFMatrix& B); + + /** + * A(i,i)+=s i=1-m. + * Add constant to diagonal elements of this + */ + void diagadd(const double& src); + + /** + * w+=part(A)*v . Conditional partial Matrix-vector-multiplication

+ * ( used elements of v determined by x ) + */ + void gsmult(dVector& w, const dVector& v,const iVector& x) const; + + + /** + * Output of the matrix in user-defined format. + */ + void print(FILE* fp, const char* format = 0) const; +//@} + + /** * Comparison operator. Be careful with * this thing, it may eat up huge amounts * of computing time! -- 2.39.5