From: wolf Date: Fri, 26 Feb 1999 16:03:55 +0000 (+0000) Subject: Remove the old-style matrices and vectors which were not templated. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=0b41be34f5300fd4b6358cfef874a2b1e467b41d;p=dealii-svn.git Remove the old-style matrices and vectors which were not templated. git-svn-id: https://svn.dealii.org/trunk@919 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/dblocksmatrix.h b/deal.II/lac/include/lac/dblocksmatrix.h deleted file mode 100644 index 5501acc8d1..0000000000 --- a/deal.II/lac/include/lac/dblocksmatrix.h +++ /dev/null @@ -1,154 +0,0 @@ -/*---------------------------- dblocksmatrix.h ---------------------------*/ -/* $Id$ */ -#ifndef __dblocksmatrix_H -#define __dblocksmatrix_H -/*---------------------------- dblocksmatrix.h ---------------------------*/ - - -#include -#include -#include - -/** - * Double precision block sparse matrix. - * The block matrix assumes the matrix consisting of blocks on - * the diagonal. These diagonal blocks and the elements below the - * diagonal blocks are used in the #precondition_BlockSOR#. - * - * This block matrix structure is given e.g. for the DG method - * for the transport equation and a downstream numbering. - * If (as for this DG method) the matrix is empty above the - * diagonal blocks BlockSOR is a direct solver. - * - * This first implementation of the BlockMatrix assumes the - * matrix having blocks each of the same block size. Varying - * block sizes within the matrix must still be implemented if needed. - * @author Ralf Hartmann, 1999 - */ -class dBlockSMatrix: public dSMatrix -{ - public: - /** - * Constructor - */ - dBlockSMatrix(); - - /** - * Destructor - */ - virtual ~dBlockSMatrix(); - - /** - * Call #dSMatrix::reinit()# and - * delete the inverse matrices if existent. - */ - - virtual void reinit(); - - /** - * Call #dSMatrix::reinit - * (const dSMatrixStruct &sparsity)# and - * delete the inverse matrices if existent. - */ - virtual void reinit (const dSMatrixStruct &sparsity); - - /** - * Call #dSMatrix::clear# and - * delete the inverse matrices if existent. - */ - virtual void clear (); - - /** - * Stores the inverse matrices of - * the diagonal blocks matrices - * in #inverse#. This costs some - * additional memory (for DG - * methods about 1/3 of that used for - * the matrix) but it - * makes the preconditioning much faster. - */ - void invert_diagblocks(); - - /** - * Block SOR. Make sure that the right block size - * of the matrix is set by #set_block_size# - * before calling this function. - * - * BlockSOR will automatically use the - * inverse matrices if they exist, if not - * then BlockSOR will waste much time - * inverting the diagonal block - * matrices in each preconditioning step. - * - * For matrices which are - * empty above the diagonal blocks - * BlockSOR is a direct solver. - */ - void precondition_BlockSOR (dVector &dst, const dVector &src) const; - - /** - * Set the right block size before calling - * #precondition_BlockSOR#. - * If block_size==1 BlockSOR is the same as SOR. - */ - void set_block_size (const unsigned int bsize); - - /** - * Gives back the size of the blocks. - */ - unsigned int block_size() const; - - /** - * Exception - */ - DeclException2 (ExcWrongBlockSize, - int, int, - << "The blocksize " << arg1 - << " and the size of the matrix " << arg2 - << " do not match."); - - DeclException2 (ExcWrongInverses, - int, int, - << "There are " << arg1 - << " inverse matrices but " << arg2 - << " cells."); - - /** - * Exception - */ - DeclException0 (ExcInverseMatricesDoNotExist); - - /** - * Exception - */ - DeclException0 (ExcInverseMatricesAlreadyExist); - - /** - * Exception - */ - DeclException0 (ExcBlockSizeNotSet); - - /** - * Exception - */ - DeclException0 (ExcInternalError); - - private: - /** - * size of the blocks. - */ - unsigned int blocksize; - - /** - * stores the inverse matrices of - * the diagonal blocks matrices - */ - vector inverse; -}; - - - -/*---------------------------- dblocksmatrix.h ---------------------------*/ -/* end of #ifndef __dblocksmatrix_H */ -#endif -/*---------------------------- dblocksmatrix.h ---------------------------*/ diff --git a/deal.II/lac/include/lac/dfmatrix.h b/deal.II/lac/include/lac/dfmatrix.h deleted file mode 100644 index a5fda7a698..0000000000 --- a/deal.II/lac/include/lac/dfmatrix.h +++ /dev/null @@ -1,568 +0,0 @@ -/*---------------------------- dfmatrix.h ---------------------------*/ -/* $Id$ */ -#ifndef __dfmatrix_H -#define __dfmatrix_H -/*---------------------------- dfmatrix.h ---------------------------*/ - -// This file is part of the DEAL Library -// DEAL is Copyright(1995) by -// Roland Becker, Guido Kanschat, Franz-Theo Suttmeier -// Revised by Wolfgang Bangerth - - -#include - - -// forward declarations -class dVector; -class iVector; - - - -/** - * 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 -{ - /** - * Component-array. - */ - double* val; - /** - * Dimension. Actual number of Columns - */ - unsigned int dim_range; - /** - * Dimension. Actual number of Rows - */ - unsigned int dim_image; - /** - * Dimension. Determines amount of reserved memory - */ - unsigned int val_size; - - /** - * Initialization . initialize memory for Matrix

- * ( m rows , n columns ) - */ - 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. - */ - double& el (const unsigned int i, const unsigned int j); - - /** - * Return the value of the element #(i,j)#. - * - * This function does no bounds checking. - */ - double el (const unsigned int i, const unsigned int j) const; - - - public: - /**@name 1: Basic Object-handling */ - //@{ - /** - * Constructor. Initialize the matrix as - * a square matrix with dimension #n#. - */ - explicit dFMatrix (const unsigned int n = 1); - - /** - * Constructor. Initialize the matrix as - * a rectangular #m# times #n# matrix. - */ - dFMatrix (const unsigned int m, const unsigned int n); - - /** - * Copy constructor. Be very careful with - * this constructor, since it may take a - * huge amount of computing time for large - * matrices!! - */ - dFMatrix (const dFMatrix&); - - /** - * Destructor. Release all memory. - */ - ~dFMatrix(); - - /** - * Comparison operator. Be careful with - * this thing, it may eat up huge amounts - * of computing time! It is most commonly - * used for internal consistency checks - * of programs. - */ - bool operator == (const dFMatrix &) const; - - /** - * A = B . Copy all elements - */ - dFMatrix& operator = (const dFMatrix& B); - - - /** - * U(0-m,0-n) = s . Fill all elements - */ - void fill (const dFMatrix& src, - const unsigned int i=0, const unsigned int j=0); - - /** - * Change Dimension. - * Set dimension to (m,n)

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

- * ( reinit quadratic matrix ) - */ - void reinit (const unsigned int n); - - /** - * Adjust Dimension. - * Set dimension to ( m(B),n(B) )

- * ( adjust to dimensions of another matrix B ) - */ - void reinit (const dFMatrix &B); - - /** - * Return number of rows of this matrix. - * To remember: this matrix is an - * $m \times n$-matrix. - */ - unsigned int m () const; - - /** - * Return number of columns of this matrix. - * To remember: this matrix is an - * $m \times n$-matrix. - */ - unsigned int n () const; - - /** - * Return whether the matrix contains only - * elements with value zero. This function - * is mainly for internal consistency - * check and should seldomly be used when - * not in debug mode since it uses quite - * some time. - */ - bool all_zero () const; - - //@} - - - /**@name 2: Data-Access - */ - //@{ - /** - * Access Elements. returns element at relative 'address' i

- * ( -> access to A(i/n , i mod n) ) - */ - double el (const unsigned int i) const; - - /** - * Return the value of the element #(i,j)#. - * Does the same as the #el(i,j)# function - * but does bounds checking. - */ - double 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 #el(i,j)# function - * but does bounds checking. - */ - double& operator() (const unsigned int i, const unsigned int j); - - /** - * Set all entries in the matrix to - * zero. - */ - void clear (); - //@} - - - /**@name 3: Basic applications on matrices - */ - //@{ - /** - * A+=B . Simple addition - */ - void add (const double s, const dFMatrix& B); - - /** - * A+=Transp(B). - * Simple addition of the transpose of B to this - */ - void Tadd (const 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 bool adding=false) 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 bool adding=false) const; - - /** - * Return the norm of the vector #v# with - * respect to the norm induced by this - * matrix, i.e. $\left$. This - * is useful, e.g. in the finite element - * context, where the $L_2$ norm of a - * function equals the matrix norm with - * respect to the mass matrix of the vector - * representing the nodal values of the - * finite element function. - * - * Note the order in which the matrix - * appears. For non-symmetric matrices - * there is a difference whether the - * matrix operates on the first - * or on the second operand of the - * scalar product. - * - * Obviously, the matrix needs to be square - * for this operation. - */ - double matrix_norm (const dVector &v) const; - - /** - * Build the matrix scalar product - * #u^T M v#. This function is mostly - * useful when building the cellwise - * scalar product of two functions in - * the finite element context. - */ - double matrix_scalar_product (const dVector &u, const dVector &v) const; - - /** - * A=Inverse(A). Inversion of (*this) by - * Gauss-Jordan-algorithm - */ - void gauss_jordan (); - - /** - * Computes the determinant of a matrix. - * This is only implemented for one two and - * three dimensions, since for higher - * dimensions the numerical work explodes. - * Obviously, the matrix needs to be square - * for this function. - */ - 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. - */ - void invert (const dFMatrix &M); - //@} - - - /**@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 (const unsigned int i, const double 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 - */ - void add_row (const unsigned int i, - const double s, const unsigned int j, - const double t, const unsigned int k); - - /** - * A(1-n,i)+=s*A(1-n,j). - * Simple addition of columns of this - */ - void add_col (const unsigned int i, const double 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 - */ - void add_col (const unsigned int i, - const double s, const unsigned int j, - const double t, const unsigned int k); - - /** - * Swap A(i,1-n) <-> A(j,1-n). - * Swap rows i and j of this - */ - void swap_row (const unsigned int i, const unsigned int j); - - /** - * Swap A(1-n,i) <-> A(1-n,j). - * Swap columns i and j of this - */ - void swap_col (const unsigned int i, const unsigned 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. - * - * #householder(src); backward(dst, src);# gives - * the solution #dst# of the linear system - * #(*this)dst=src#. - * - * Note that #src# and #(*this)# (i.e. the - * matrix itself) is changed in - * the process of the #householder(src)# function!! - */ - void householder (dVector& src); - - /** - * Least - Squares - Approximation by QR-factorization. - * - * Note that #src# and #(*this)# (i.e. the - * matrix itself) is changed in - * the process of this function!! - */ - 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 (const 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 (ostream& s, int width=5, int precision=2) const; - - /** - * Print the matrix in the usual format, - * i.e. as a matrix and not as a list of - * nonzero elements. For better - * readability, zero elements - * are displayed as empty space. - * - * Each entry is printed in scientific - * format, with one pre-comma digit and - * the number of digits given by - * #precision# after the comma, with one - * space following. - * The precision defaults to four, which - * suffices for most cases. The precision - * and output format are {\it not} - * properly reset to the old values - * when the function exits. - * - * You should be aware that this function - * may produce {\bf large} amounts of - * output if applied to a large matrix! - * Be careful with it. - */ - void print_formatted (ostream &out, - const unsigned int presicion=3) const; - //@} - - /** - * Exception - */ - DeclException2 (ExcInvalidIndex, - int, int, - << "The given index " << arg1 - << " should be less than " << arg2 << "."); - /** - * Exception - */ - DeclException2 (ExcDimensionMismatch, - int, int, - << "The two dimensions " << arg1 << " and " << arg2 - << " do not match here."); - /** - * Exception - */ - DeclException0 (ExcNotQuadratic); - /** - * Exception - */ - DeclException0 (ExcNotRegular); - /** - * Exception - */ - DeclException0 (ExcInternalError); - /** - * Exception - */ - DeclException3 (ExcInvalidDestination, - int, int, int, - << "Target region not in matrix: size in this direction=" - << arg1 << ", size of new matrix=" << arg2 - << ", offset=" << arg3); - /** - * Exception - */ - DeclException1 (ExcNotImplemented, - int, - << "This function is not implemented for the given" - << " matrix dimension " << arg1); - /** - * Exception - */ - DeclException0 (ExcIO); -}; - - - - - -/*-------------------------Inline functions -------------------------------*/ - - -inline -double & dFMatrix::el (const unsigned int i, const unsigned int j) { - return val[i*dim_range+j]; -}; - - - -inline -double dFMatrix::el (const unsigned int i, const unsigned int j) const { - return val[i*dim_range+j]; -}; - - - -inline -unsigned int dFMatrix::m() const { - return dim_image; -}; - - - -inline -unsigned int dFMatrix::n() const { - return dim_range; -}; - - - -inline -double dFMatrix::el (const unsigned int i) const { - return val[i]; -}; - - - -inline -double dFMatrix::operator() (const unsigned int i, const unsigned int j) const { - Assert (i - - -//forward declarations -class dVector; -class iVector; -class ostream; - - - -/* -CLASS - dSMatrixStruct - - @author Original version by Roland Becker, Guido Kanschat, Franz-Theo Suttmeier; lots of enhancements, reorganisation and documentation by Wolfgang Bangerth - */ -class dSMatrixStruct -{ - private: - /** - * Copy constructor, made private in order to - * prevent copying such an object which does - * not make much sense because you can use - * a structure like this for more than one - * matrix. - * - * Because it is not needed, this function - * is not implemented. - */ - dSMatrixStruct (const dSMatrixStruct &); - - public: - /** - * Initialize the matrix empty, i.e. with - * no memory allocated. This is useful if - * you want such objects as member - * variables in other classes. You can make - * the structure usable by calling the - * #reinit# function. - */ - dSMatrixStruct (); - - /** - * Initialize a rectangular matrix with - * #m# rows and #n# columns, - * with at most #max_per_row# - * nonzero entries per row. - */ - dSMatrixStruct (const unsigned int m, - const unsigned int n, - const unsigned int max_per_row); - - /** - * Initialize a square matrix of dimension - * #n# with at most #max_per_row# - * nonzero entries per row. - */ - dSMatrixStruct (const unsigned int n, - const unsigned int max_per_row); - - /** - * Destructor. - */ - ~dSMatrixStruct (); - - /** - * Reallocate memory and set up data - * structures for a new matrix with - * #m# rows and #n# columns, - * with at most #max_per_row# - * nonzero entries per row. - * - * If #m*n==0# all memory is freed, - * resulting in a total reinitialization - * of the object. If it is nonzero, new - * memory is only allocated if the new - * size extends the old one. This is done - * to save time and to avoid fragmentation - * of the heap. - */ - void reinit (const unsigned int m, - const unsigned int n, - const unsigned int max_per_row); - - /** - * This function compresses the sparsity - * structure that this object represents. - * It does so by eliminating unused - * entries and sorting the remaining - * ones to allow faster access by usage - * of binary search algorithms. A special - * sorting scheme is used for the diagonal - * entry of square matrices, which is - * always the first entry of each row. - * - * #dSMatrix# objects require the - * #dSMatrixStruct# objects they are - * initialized with to be compressed, to - * reduce memory requirements. - */ - void compress (); - - /** - * Return whether the object is empty. It - * is empty if no memory is allocated, - * which is the same as that both - * dimensions are zero. - */ - bool empty () const; - - - /** - * Return the index of the matrix - * element with row number #i# and - * column number #j#. If the matrix - * element is not a nonzero one, - * return -1. - * - * This function is usually called - * by the #operator()# of the - * #dSMatrix#. It shall only be - * called for compressed sparsity - * patterns, since in this case - * searching whether the entry - * exists can be done quite fast - * with a binary sort algorithm - * because the column numbers are - * sorted. - */ - int operator() (const unsigned int i, const unsigned int j) const; - - /** - * Add a nonzero entry to the matrix. - * This function may only be called - * for non-compressed sparsity patterns. - * - * If the entry already exists, nothing - * bad happens. - */ - void add (const unsigned int i, const unsigned int j); - - /** - * This matrix adds a whole connectivity - * list to the sparsity structure - * respresented by this object. It assumes - * the #rowcols# array to be a list of - * indices which are all linked together, - * i.e. all entries - * #(rowcols[i], rowcols[j])# for all - * #i,j=0...n# for this sparsity pattern - * are created. #n# is assumed to be the - * number of elements pointed to by - * #rowcols#. - */ - void add_matrix (const unsigned int n, const int* rowcols); - - ////////// - void add_matrix (const unsigned int m, const unsigned int n, - const int* rows, const int* cols); - ////////// - void add_matrix (const iVector& rowcols); - ////////// - void add_matrix (const iVector& rows, const iVector& cols); - - /** - * Print the sparsity of the matrix - * in a format that #gnuplot# understands - * and which can be used to plot the - * sparsity pattern in a graphical - * way. The format consists of pairs - * #i j# of nonzero elements, each - * representing one entry of this - * matrix, one per line of the output - * file. Indices are counted from - * zero on, as usual. Since sparsity - * patterns are printed in the same - * way as matrices are displayed, we - * print the negative of the column - * index, which means that the - * #(0,0)# element is in the top left - * rather than in the bottom left - * corner. - * - * Print the sparsity pattern in - * gnuplot by setting the data style - * to dots or points and use the - * #plot# command. - */ - void print_gnuplot (ostream &out) const; - - /** - * Return number of rows of this - * matrix, which equals the dimension - * of the image space. - */ - unsigned int n_rows () const; - - /** - * Return number of columns of this - * matrix, which equals the dimension - * of the range space. - */ - unsigned int n_cols () const; - - /** - * Compute the bandwidth of the matrix - * represented by this structure. The - * bandwidth is the maximum of - * $|i-j|$ for which the index pair - * $(i,j)$ represents a nonzero entry - * of the matrix. - */ - unsigned int bandwidth () const; - - /** - * Return the number of nonzero elements of - * this matrix. Actually, it returns the - * number of entries in the sparsity - * pattern; if any of the entries should - * happen to be zero, it is counted - * anyway. - * - * This function may only be called if the - * matrix struct is compressed. It does not - * make too much sense otherwise anyway. - */ - unsigned int n_nonzero_elements () const; - - /** - * Return whether the structure is - * compressed or not. - */ - bool is_compressed () const; - - /** - * This is kind of an expert mode: get - * access to the rowstart array, but - * readonly. - * - * Though the return value is declared - * #const#, you should be aware that it - * may change if you call any nonconstant - * function of objects which operate on - * it. - * - * You should use this interface very - * carefully and only if you are absolutely - * sure to know what you do. You should - * also note that the structure of these - * arrays may change over time. - * If you change the layout yourself, you - * should also rename this function to - * avoid programs relying on outdated - * information! - */ - const unsigned int * get_rowstart_indices () const; - - /** - * This is kind of an expert mode: get - * access to the colnums array, but - * readonly. - * - * Though the return value is declared - * #const#, you shoudl be aware that it - * may change if you call any nonconstant - * function of objects which operate on - * it. - * - * You should use this interface very - * carefully and only if you are absolutely - * sure to know what you do. You should - * also note that the structure of these - * arrays may change over time. - * If you change the layout yourself, you - * should also rename this function to - * avoid programs relying on outdated - * information! - */ - const int * get_column_numbers () const; - - - /** - * Exception - */ - DeclException1 (ExcInvalidNumber, - int, - << "The provided number is invalid here: " << arg1); - /** - * Exception - */ - DeclException2 (ExcInvalidIndex, - int, int, - << "The given index " << arg1 - << " should be less than " << arg2 << "."); - /** - * Exception - */ - DeclException2 (ExcNotEnoughSpace, - int, int, - << "Upon entering a new entry to row " << arg1 - << ": there was no free entry any more. " << endl - << "(Maximum number of entries for this row: " - << arg2 << "; maybe the matrix is already compressed?)"); - /** - * Exception - */ - DeclException0 (ExcNotCompressed); - /** - * Exception - */ - DeclException0 (ExcMatrixIsCompressed); - /** - * Exception - */ - DeclException0 (ExcEmptyObject); - /** - * Exception - */ - DeclException0 (ExcInternalError); - /** - * Exception - */ - DeclException0 (ExcIO); - - private: - unsigned int max_dim; - unsigned int rows, cols; - unsigned int vec_len, max_vec_len; - unsigned int max_row_len; - unsigned int* rowstart; - int* colnums; - - /** - * Store whether the #compress# function - * was called for this object. - */ - bool compressed; - - friend class dSMatrix; -}; - - - - -/* -CLASS - dSMatrix - - @author Original version by Roland Becker, Guido Kanschat, Franz-Theo Suttmeier; lots of enhancements, reorganisation and documentation by Wolfgang Bangerth 1998 - */ -class dSMatrix -{ - public: - - /** - * Constructor; initializes the matrix to - * be empty, without any structure, i.e. - * the matrix is not usable at all. This - * constructor is therefore only useful - * for matrices which are members of a - * class. All other matrices should be - * created at a point in the data flow - * where all necessary information is - * available. - * - * You have to initialize - * the matrix before usage with - * #reinit(dSMatrixStruct)#. - */ - dSMatrix (); - - /** - * Constructor. Takes the given matrix - * sparisty structure to represent the - * sparsity pattern of this matrix. You - * can change the sparsity pattern later - * on by calling the #reinit# function. - * - * You have to make sure that the lifetime - * of the sparsity structure is at least - * as long as that of this matrix or as - * long as #reinit# is not called with a - * new sparsity structure. - */ - dSMatrix (const dSMatrixStruct &sparsity); - - /** - * Destructor. Free all memory, but do not - * release the memory of the sparsity - * structure. - */ - virtual ~dSMatrix (); - - - /** - * Reinitialize the object but keep to - * the sparsity pattern previously used. - * This may be necessary if you #reinit#'d - * the sparsity structure and want to - * update the size of the matrix. - * - * Note that memory is only reallocated if - * the new size exceeds the old size. If - * that is not the case, the allocated - * memory is not reduced. However, if the - * sparsity structure is empty (i.e. the - * dimensions are zero), then all memory - * is freed. - */ - virtual void reinit (); - - /** - * Reinitialize the sparse matrix with the - * given sparsity pattern. The latter tells - * the matrix how many nonzero elements - * there need to be reserved. - * - * Regarding memory allocation, the same - * applies as said above. - * - * You have to make sure that the lifetime - * of the sparsity structure is at least - * as long as that of this matrix or as - * long as #reinit# is not called with a - * new sparsity structure. - */ - virtual void reinit (const dSMatrixStruct &sparsity); - - /** - * Release all memory and return to a state - * just like after having called the - * default constructor. It also forgets the - * sparsity pattern it was previously tied - * to. - */ - virtual void clear (); - - /** - * Return the dimension of the image space. - * To remember: the matrix is of dimension - * $m \times n$. - */ - unsigned int m () const; - - /** - * Return the dimension of the range space. - * To remember: the matrix is of dimension - * $m \times n$. - */ - unsigned int n () const; - - /** - * Return the number of nonzero elements of - * this matrix. Actually, it returns the - * number of entries in the sparsity - * pattern; if any of the entries should - * happen to be zero, it is counted - * anyway. - */ - unsigned int n_nonzero_elements () const; - - /** - * Set the element #(i,j)# to #value#. - * Throws an error if the entry does - * not exist. Still, it is allowed to store - * zero values in non-existent fields. - */ - void set (const unsigned int i, const unsigned int j, - const double value); - - /** - * Add #value# to the element #(i,j)#. - * Throws an error if the entry does - * not exist. Still, it is allowed to store - * zero values in non-existent fields. - */ - void add (const unsigned int i, const unsigned int j, - const double value); - - /** - * Copy the given matrix to this one. - * The operation throws an error if the - * sparsity patterns of the two involved - * matrices do not point to the same - * object, since in this case the copy - * operation is cheaper. Since this - * operation is notheless not for free, - * we do not make it available through - * #operator =#, since this may lead - * to unwanted usage, e.g. in copy - * arguments to functions, which should - * really be arguments by reference. - * - * The function returns a reference to - * #this#. - */ - dSMatrix & copy_from (const dSMatrix &); - - /** - * Add #matrix# scaled by #factor# to this - * matrix. The function throws an error - * if the sparsity patterns of the two - * involved matrices do not point to the - * same object, since in this case the - * operation is cheaper. - */ - void add_scaled (const double factor, const dSMatrix &matrix); - - /** - * Return the value of the entry (i,j). - * This may be an expensive operation - * and you should always take care - * where to call this function. - * In order to avoid abuse, this function - * throws an exception if the wanted - * element does not exist in the matrix. - */ - double operator () (const unsigned int i, const unsigned int j) const; - - /** - * Return the main diagonal element in - * the #i#th row. This function throws an - * error if the matrix is not square. - * - * This function is considerably faster - * than the #operator()#, since for - * square matrices, the diagonal entry is - * always the first to be stored in each - * row and access therefore does not - * involve searching for the right column - * number. - */ - double diag_element (const unsigned int i) const; - - /** - * This is kind of an expert mode: get - * access to the #i#th element of this - * matrix. The elements are stored in - * a consecutive way, refer to the - * #dSMatrixStruct# class for more details. - * - * You should use this interface very - * carefully and only if you are absolutely - * sure to know what you do. You should - * also note that the structure of these - * arrays may change over time. - * If you change the layout yourself, you - * should also rename this function to - * avoid programs relying on outdated - * information! - */ - double global_entry (const unsigned int i) const; - - /** - * Same as above, but with write access. - * You certainly know what you do? - */ - double & global_entry (const unsigned int i); - - /** - * Matrix-vector multiplication: let - * #dst = M*src# with #M# being this matrix. - */ - void vmult (dVector& dst, const dVector& src) const; - - /** - * Matrix-vector multiplication: let - * #dst = M^T*src# with #M# being this - * matrix. This function does the same as - * #vmult# but takes the transposed matrix. - */ - void Tvmult (dVector& dst, const dVector& src) const; - - - /** - * Return the norm of the vector #v# with - * respect to the norm induced by this - * matrix, i.e. $\left$. This - * is useful, e.g. in the finite element - * context, where the $L_2$ norm of a - * function equals the matrix norm with - * respect to the mass matrix of the vector - * representing the nodal values of the - * finite element function. - * - * Note the order in which the matrix - * appears. For non-symmetric matrices - * there is a difference whether the - * matrix operates on the first - * or on the second operand of the - * scalar product. - * - * Obviously, the matrix needs to be square - * for this operation. - */ - double matrix_norm (const dVector &v) const; - - // - double residual (dVector& dst, const dVector& x, - const dVector& b) const; - // - void precondition_Jacobi (dVector& dst, const dVector& src, - const double om = 1.) const; - // - void precondition_SSOR (dVector& dst, const dVector& src, - const double om = 1.) const; - // - void precondition_SOR (dVector& dst, const dVector& src, - const double om = 1.) const; - // - void SSOR (dVector& dst, const double om = 1.) const; - // - void SOR (dVector& dst, const double om = 1.) const; - // - void precondition (dVector& dst, const dVector& src) const; - - /** - * Return a (constant) reference to the - * underlying sparsity pattern of this - * matrix. - * - * Though the return value is declared - * #const#, you shoudl be aware that it - * may change if you call any nonconstant - * function of objects which operate on - * it. - */ - const dSMatrixStruct & get_sparsity_pattern () const; - - /** - * Print the matrix to the given stream, - * using the format - * #(line,col) value#, i.e. one - * nonzero entry of the matrix per line. - */ - void print (ostream &out) const; - - /** - * Print the matrix in the usual format, - * i.e. as a matrix and not as a list of - * nonzero elements. For better - * readability, elements not in the matrix - * are displayed as empty space, while - * matrix elements which are explicitely - * set to zero are displayed as such. - * - * Each entry is printed in scientific - * format, with one pre-comma digit and - * the number of digits given by - * #precision# after the comma, with one - * space following. - * The precision defaults to four, which - * suffices for most cases. The precision - * and output format are {\it not} - * properly reset to the old values - * when the function exits. - * - * You should be aware that this function - * may produce {\bf large} amounts of - * output if applied to a large matrix! - * Be careful with it. - */ - void print_formatted (ostream &out, - const unsigned int presicion=3) const; - - /** - * Exception - */ - DeclException0 (ExcNotCompressed); - /** - * Exception - */ - DeclException0 (ExcMatrixNotInitialized); - /** - * Exception - */ - DeclException2 (ExcDimensionsDontMatch, - int, int, - << "The dimensions " << arg1 << " and " << arg2 - << " do not match properly."); - /** - * Exception - */ - DeclException2 (ExcInvalidIndex, - int, int, - << "The entry with index <" << arg1 << ',' << arg2 - << "> does not exist."); - /** - * Exception - */ - DeclException1 (ExcInvalidIndex1, - int, - << "The index " << arg1 << " is not in the allowed range."); - /** - * Exception - */ - DeclException0 (ExcMatrixNotSquare); - /** - * Exception - */ - DeclException0 (ExcDifferentSparsityPatterns); - /** - * Exception - */ - DeclException0 (ExcIO); - - private: - const dSMatrixStruct * cols; - double* val; - unsigned int max_len; -}; - - - - - -/*---------------------- Inline functions -----------------------------------*/ - -inline -unsigned int dSMatrixStruct::n_rows () const { - return rows; -}; - - - -inline -unsigned int dSMatrixStruct::n_cols () const { - return cols; -}; - - - -inline -bool dSMatrixStruct::is_compressed () const { - return compressed; -}; - - - -inline -const unsigned int * dSMatrixStruct::get_rowstart_indices () const { - return rowstart; -}; - - - -inline -const int * dSMatrixStruct::get_column_numbers () const { - return colnums; -}; - - - -inline -unsigned int dSMatrix::m () const -{ - Assert (cols != 0, ExcMatrixNotInitialized()); - return cols->rows; -}; - - - -inline -unsigned int dSMatrix::n () const -{ - Assert (cols != 0, ExcMatrixNotInitialized()); - return cols->cols; -}; - - - -inline -void dSMatrix::set (const unsigned int i, const unsigned int j, - const double value) { - Assert (cols != 0, ExcMatrixNotInitialized()); - Assert ((cols->operator()(i,j) != -1) || (value == 0.), - ExcInvalidIndex(i,j)); - - const int index = cols->operator()(i,j); - - if (index >= 0) val[index] = value; -}; - - - -inline -void dSMatrix::add (const unsigned int i, const unsigned int j, - const double value) { - Assert (cols != 0, ExcMatrixNotInitialized()); - Assert ((cols->operator()(i,j) != -1) || (value == 0.), - ExcInvalidIndex(i,j)); - - const int index = cols->operator()(i,j); - - if (index >= 0) val[index] += value; -}; - - - - - -inline -double dSMatrix::operator () (const unsigned int i, const unsigned int j) const { - Assert (cols != 0, ExcMatrixNotInitialized()); - Assert (cols->operator()(i,j) != -1, - ExcInvalidIndex(i,j)); - return val[cols->operator()(i,j)]; -}; - - - -inline -double dSMatrix::diag_element (const unsigned int i) const { - Assert (cols != 0, ExcMatrixNotInitialized()); - Assert (m() == n(), ExcMatrixNotSquare()); - Assert (irowstart[i]]; -}; - - - -inline -double dSMatrix::global_entry (const unsigned int j) const { - Assert (cols != 0, ExcMatrixNotInitialized()); - return val[j]; -}; - - - -inline -double & dSMatrix::global_entry (const unsigned int j) { - Assert (cols != 0, ExcMatrixNotInitialized()); - return val[j]; -}; - - - -/*---------------------------- dsmatrix.h ---------------------------*/ -/* end of #ifndef __dsmatrix_H */ -#endif -/*---------------------------- dsmatrix.h ---------------------------*/ - - diff --git a/deal.II/lac/include/lac/dvector.h b/deal.II/lac/include/lac/dvector.h deleted file mode 100644 index 39f61eaf6e..0000000000 --- a/deal.II/lac/include/lac/dvector.h +++ /dev/null @@ -1,483 +0,0 @@ -/*---------------------------- dvector.h ---------------------------*/ -/* $Id$ */ -#ifndef __dvector_H -#define __dvector_H -/*---------------------------- dvector.h ---------------------------*/ - -// This file is part of the DEAL Library -// DEAL is Copyright(1995) by -// Roland Becker, Guido Kanschat, Franz-Theo Suttmeier -// Revised by Wolfgang Bangerth - -#include -#include - - - -/** - * Double precision Vector. - * Memory for Components is supplied explicitly

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

- * - all necessary methods for Vectors are supplied

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

- * CONVENTIONS for used `equations` :

- * - THIS vector is always named `U`

- * - vectors are always uppercase , scalars are lowercase - * - * @author Roland Becker, Guido Kanschat, Franz-Theo Suttmeier, revised and extended by Wolfgang Bangerth, documented by Klaus Mampel and Wolfgang Bangerth - */ -class dVector { - friend class dFMatrix; - - 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. - */ - double *val; - - public: - - /** - * Declare iterator types just like those - * for the C++ standard library: - * - * Data type stored by this container. - */ - typedef double value_type; - - /** - * Declare standard types used in all - * containers. - */ - typedef value_type* pointer; - typedef const value_type* const_pointer; - typedef value_type* iterator; - typedef const value_type* const_iterator; - typedef value_type& reference; - typedef const value_type& const_reference; - typedef size_t size_type; - - - /** - * @name 1: Basic Object-handling - */ - //@{ - /** - * Dummy-Constructor. Dimension=0 - */ - dVector (); - - /** - * Copy-Constructor. Dimension set to that of V ,

- * all components are copied from V - */ - dVector (const dVector& V); - - /** - * Constructor. Set dimension to #n# and - * initialize all elements with zero. - */ - dVector (const unsigned int n); - - /** - * Destructor. Clears memory - */ - ~dVector (); - - /** - * Set all entries to zero. Equivalent to - * #v = 0#, but more obvious and faster. - * Note that this function does not change - * the size of the vector, unlike the - * STL's #vector<>::clear# function. - */ - void clear (); - - /** - * U(0-N) = s . Fill all components - */ - dVector& operator= (const double s); - - /** - * U = V . Copy all components - */ - dVector& operator= (const dVector& V); - - /** - * U = U * V . Scalar Produkt - */ - double operator* (const dVector& V) const; - - /** - * Return square of the l2-norm. - */ - double norm_sqr () const; - - /** - * Return the mean value of the elements of - * this vector. - */ - double mean_value () const; - - /** - * Return the l1-norm of the vector, i.e. - * the sum of the absolute values. - */ - double l1_norm () const; - - /** - * Return the l2-norm of the vector, i.e. - * the square root of the sum of the - * squares of the elements. - */ - double l2_norm () const; - - /** - * Return the maximum absolute value of the - * elements of this vector. - */ - double linfty_norm () const; - - - /** - * Change the dimension of the vector to - * #N#. The reserved memory for this vector - * remains unchanged if possible, to make - * things faster, but this may waste some - * memory, so take this in the back of your - * head. - * However, if #N==0# all memory is freed, - * i.e. if you want to resize the vector - * and release the memory not needed, you - * have to first call #reinit(0)# and then - * #reinit(N)#. This cited behaviour is - * analogous to that of the STL containers. - * - * On #fast==false#, the vector is filled by - * zeros. - */ - void reinit (const unsigned int N, const bool fast=false); - - /** - * Change the dimension to that of the - * vector #V#. The same applies as for - * the other #reinit# function. - * - * The elements of #V# are not copied, i.e. - * this function is the same as calling - * #reinit (V.size(), fast)#. - */ - void reinit (const dVector& V, const bool fast=false); - - /** - * Return dimension of the vector. This - * function was formerly called #n()#, but - * was renamed to get the #dVector# class - * closer to the C++ standard library's - * #vector# container. - */ - unsigned int size () const; - - /** - * Return whether the vector contains only - * elements with value zero. This function - * is mainly for internal consistency - * check and should seldomly be used when - * not in debug mode since it uses quite - * some time. - */ - bool all_zero () const; - - /** - * Make the #dVector# class a bit like the - * #vector<># class of the C++ standard - * library by returning iterators to - * the start and end of the elements of this - * vector. - */ - iterator begin (); - - /** - * Return constant iterator to the start of - * the vectors. - */ - const_iterator begin () const; - - /** - * Return an iterator pointing to the - * element past the end of the array. - */ - iterator end (); - - /** - * Return a constant iterator pointing to - * the element past the end of the array. - */ - const_iterator end () const; - //@} - - - /** - * @name 2: Data-Access - */ - //@{ - /** - * Access Components. returns U(i) , - * INLINE - */ - double operator() (const unsigned int i) const; - - /** - * Access Components. returns U(i) , - * INLINE - */ - double& operator() (const unsigned int i); - //@} - - - /** - * @name 3: Modification of vectors - */ - //@{ - /** - * Fast equivalent to #U.add(1, V)#. - */ - dVector & operator += (const dVector &V); - - /** - * Fast equivalent to #U.add(-1, V)#. - */ - dVector & operator -= (const dVector &V); - - /** - * U(0-DIM)+=s. - * Addition of #s# to all components. Note - * that #s# is a scalar and not a vector. - */ - void add (const double s); - - /** - * U+=V. - * Simple vector addition, equal to the - * #operator +=#. - */ - void add (const dVector& V); - - /** - * U+=a*V. - * Simple addition of a scaled vector. - */ - void add (const double a, const dVector& V); - - /** - * U+=a*V+b*W. - * Multiple addition of scaled vectors. - */ - void add (const double a, const dVector& V, - const double b, const dVector& W); - - /** - * U=s*U+V. - * Scaling and simple vector addition. - */ - void sadd (const double s, const dVector& V); - - /** - * U=s*U+a*V. - * Scaling and simple addition. - */ - void sadd (const double s, const double a, const dVector& V); - - /** - * U=s*U+a*V+b*W. - * Scaling and multiple addition. - */ - void sadd (const double s, const double a, - const dVector& V, const double b, const dVector& W); - - /** - * U=s*U+a*V+b*W+c*X. - * Scaling and multiple addition. - */ - void sadd (const double s, const double a, - const dVector& V, const double b, const dVector& W, - const double c, const dVector& X); - - /** - * Scale each element of the vector by the - * given factor. This function was - * previously called #equ(double)#, which - * in my eyes is an extremely unintuitive - * naming and was thus replaced. - */ - void scale (const double factor); - - /** - * U=a*V. Replacing - */ - void equ (const double a, const dVector& V); - - /** - * U=a*V+b*W. - * Replacing by sum. - */ - void equ (const double a, const dVector& V, - const double b, const dVector& W); - - /** - * Compute the elementwise ratio of the - * two given vectors, that is let - * #this[i] = a[i]/b[i]#. This is useful - * for example if you want to compute - * the cellwise ratio of true to estimated - * error. - * - * This vector is appropriately scaled to - * hold the result. - * - * If any of the #b[i]# is zero, the result - * is undefined. No attempt is made to - * catch such situations. - */ - void ratio (const dVector &a, const dVector &b); - //@} - - - /** - * @name 5: Mixed stuff - */ - //@{ - /** - * Output of vector in user-defined format. - */ - void print (FILE* fp, const char* format = 0) const; - - /** - * Output of vector in user-defined format. - */ - void print (const char* format = 0) const; - - /** - * Print to given stream, one element per line. - */ - void print (ostream &) const; - //@} - - /** - * Exception - */ - DeclException2 (ExcDimensionsDontMatch, - int, int, - << "The dimensions " << arg1 << " and " << arg2 - << " do not match here."); - /** - * Exception - */ - DeclException2 (ExcInvalidIndex, - int, int, - << "The given index " << arg1 - << " should be less than " << arg2 << "."); - /** - * Exception - */ - DeclException1 (ExcInvalidNumber, - int, - << "The provided number is invalid here: " << arg1); - /** - * Exception - */ - DeclException0 (ExcOutOfMemory); - /** - * Exception - */ - DeclException0 (ExcEmptyVector); - /** - * Exception - */ - DeclException0 (ExcIO); -}; - - - - - - -/*----------------------- Inline functions ----------------------------------*/ - - -inline unsigned int dVector::size () const -{ - return dim; -} - - - -inline -dVector::iterator dVector::begin () { - return &val[0]; -}; - - - -inline -dVector::const_iterator dVector::begin () const { - return &val[0]; -}; - - - -inline -dVector::iterator dVector::end () { - return &val[dim]; -}; - - - -inline -dVector::const_iterator dVector::end () const { - return &val[dim]; -}; - - - -inline double dVector::operator() (const unsigned int i) const -{ - Assert (i -#include - - -dBlockSMatrix::dBlockSMatrix (): - blocksize(0) {}; - -dBlockSMatrix::~dBlockSMatrix () -{ - if (inverse.size()!=0) - inverse.erase(inverse.begin(), inverse.end()); -} - - -void dBlockSMatrix::reinit () -{ - if (inverse.size()!=0) - inverse.erase(inverse.begin(), inverse.end()); - blocksize=0; - dSMatrix::reinit (); -} - - -void dBlockSMatrix::reinit (const dSMatrixStruct &sparsity) -{ - if (inverse.size()!=0) - inverse.erase(inverse.begin(), inverse.end()); - blocksize=0; - dSMatrix::reinit (sparsity); -} - - -void dBlockSMatrix::clear () -{ - dSMatrix::clear(); - if (inverse.size()!=0) - inverse.erase(inverse.begin(), inverse.end()); - blocksize=0; -} - - -void dBlockSMatrix::set_block_size(unsigned int bsize) { - blocksize=bsize; -} - - - -unsigned int dBlockSMatrix::block_size() const { - return blocksize; -} - - - -void dBlockSMatrix::precondition_BlockSOR (dVector &dst, const dVector &src) const -{ - Assert (m() == n(), ExcMatrixNotSquare()); - Assert (blocksize!=0, ExcBlockSizeNotSet()); - Assert (m()%blocksize==0, ExcWrongBlockSize(blocksize, m())); - unsigned int n_cells=m()/blocksize; - Assert (inverse.size()==0 || inverse.size()==n_cells, - ExcWrongInverses(inverse.size(), n_cells)); - - const dSMatrixStruct &spars=get_sparsity_pattern(); - const unsigned int *rowstart = spars.get_rowstart_indices(); - const int *columns = spars.get_column_numbers(); - - dVector b_cell(blocksize), x_cell(blocksize); - - // cell_row, cell_column are the - // numbering of the blocks (cells). - // row_cell, column_cell are the local - // numbering of the unknowns in the - // blocks. - // row, column are the global numbering - // of the unkowns. - unsigned int row, column, row_cell, begin_diag_block=0; - double b_cell_row; - - if (inverse.size()==0) - { - dFMatrix M_cell(blocksize); - for (unsigned int cell=0; cell(columns[j])) - < begin_diag_block) - b_cell_row -= global_entry(j) * dst(column); - b_cell(row_cell)=b_cell_row; - for (unsigned int column_cell=0, column=cell*blocksize; - column_cell(columns[j])) < begin_diag_block) - { - b_cell_row -= global_entry(j) * dst(column); - } - b_cell(row_cell)=b_cell_row; - } - inverse[cell].vmult(x_cell, b_cell); - // distribute x_cell to dst - for (row=cell*blocksize, row_cell=0; row_cell 4 - if (blocksize<=4) - inverse[cell].invert(M_cell); - else - { - M_cell.gauss_jordan(); - inverse[cell]=M_cell; - } - } -} - - - -/*---------------------------- dblocksmatrix.cc ---------------------------*/ diff --git a/deal.II/lac/source/dfmatrix.cc b/deal.II/lac/source/dfmatrix.cc deleted file mode 100644 index 040141a95c..0000000000 --- a/deal.II/lac/source/dfmatrix.cc +++ /dev/null @@ -1,1333 +0,0 @@ -// $Id$ - -#include -#include -#include - -#include -#include -#include -#include - - - -dFMatrix::dFMatrix (const unsigned int n) { - init (n,n); -}; - - - -dFMatrix::dFMatrix (const unsigned int m, const unsigned int n) { - init (m,n); -}; - - - -dFMatrix::dFMatrix (const dFMatrix &m) -{ - init (m.dim_image, m.dim_range); - double * p = &val[0]; - const double * vp = &m.val[0]; - const double * const e = &val[dim_image*dim_range]; - - while (p!=e) - *p++ = *vp++; -}; - - - -void dFMatrix::init (const unsigned int mm, const unsigned int nn) -{ - val_size = nn*mm; - val = new double[val_size]; - dim_range = nn; - dim_image = mm; - clear (); -}; - - - -dFMatrix::~dFMatrix () { - delete[] val; -}; - - - -bool dFMatrix::all_zero () const { - const double *p = &val[0], - *e = &val[n()*m()]; - while (p!=e) - if (*p++ != 0.0) - return false; - - return true; -}; - - - -void dFMatrix::reinit (const unsigned int mm, const unsigned int nn) -{ - if (val_size=0; --i) - { - s = src(i); - for (j=i+1; j= src.n() + j, ExcInvalidDestination(n(), src.n(), j)); - Assert (m() >= src.m() + i, ExcInvalidDestination(m(), src.m(), i)); - - for (unsigned int ii=0; ii=1) && (dim_range<=3), ExcNotImplemented(dim_range)); - - switch (dim_range) - { - case 1: - return el(0,0); - case 2: - return el(0,0)*el(1,1) - el(1,0)*el(0,1); - case 3: - return (el(0,0)*el(1,1)*el(2,2) - -el(0,0)*el(1,2)*el(2,1) - -el(1,0)*el(0,1)*el(2,2) - +el(1,0)*el(0,2)*el(2,1) - +el(2,0)*el(0,1)*el(1,2) - -el(2,0)*el(0,2)*el(1,1)); - default: - return 0; - }; -}; - -double dFMatrix::norm2 () const -{ - double s = 0.; - for (unsigned int i=0;i=1) && (dim_range<=4), ExcNotImplemented(dim_range)); - Assert (dim_range == M.dim_range, - ExcDimensionMismatch(dim_range,M.dim_range)); - Assert (dim_image == M.dim_image, - ExcDimensionMismatch(dim_image,M.dim_image)); - - switch (dim_range) - { - case 1: - val[0] = 1.0/M.val[0]; - return; - case 2: - // this is Maple output, - // thus a bit unstructured - { - const double t4 = 1.0/(M.el(0,0)*M.el(1,1)-M.el(0,1)*M.el(1,0)); - el(0,0) = M.el(1,1)*t4; - el(0,1) = -M.el(0,1)*t4; - el(1,0) = -M.el(1,0)*t4; - el(1,1) = M.el(0,0)*t4; - return; - }; - - case 3: - { - const double t4 = M.el(0,0)*M.el(1,1), - t6 = M.el(0,0)*M.el(1,2), - t8 = M.el(0,1)*M.el(1,0), - t00 = M.el(0,2)*M.el(1,0), - t01 = M.el(0,1)*M.el(2,0), - t04 = M.el(0,2)*M.el(2,0), - t07 = 1.0/(t4*M.el(2,2)-t6*M.el(2,1)-t8*M.el(2,2)+ - t00*M.el(2,1)+t01*M.el(1,2)-t04*M.el(1,1)); - el(0,0) = (M.el(1,1)*M.el(2,2)-M.el(1,2)*M.el(2,1))*t07; - el(0,1) = -(M.el(0,1)*M.el(2,2)-M.el(0,2)*M.el(2,1))*t07; - el(0,2) = -(-M.el(0,1)*M.el(1,2)+M.el(0,2)*M.el(1,1))*t07; - el(1,0) = -(M.el(1,0)*M.el(2,2)-M.el(1,2)*M.el(2,0))*t07; - el(1,1) = (M.el(0,0)*M.el(2,2)-t04)*t07; - el(1,2) = -(t6-t00)*t07; - el(2,0) = -(-M.el(1,0)*M.el(2,1)+M.el(1,1)*M.el(2,0))*t07; - el(2,1) = -(M.el(0,0)*M.el(2,1)-t01)*t07; - el(2,2) = (t4-t8)*t07; - return; - }; - - case 4: - { - // with (linalg); - // a:=matrix(4,4); - // evalm(a); - // ai:=inverse(a); - // readlib(C); - // C(ai,optimized,filename=x4); - - const double t14 = M.el(0,0)*M.el(1,1); - const double t15 = M.el(2,2)*M.el(3,3); - const double t17 = M.el(2,3)*M.el(3,2); - const double t19 = M.el(0,0)*M.el(2,1); - const double t20 = M.el(1,2)*M.el(3,3); - const double t22 = M.el(1,3)*M.el(3,2); - const double t24 = M.el(0,0)*M.el(3,1); - const double t25 = M.el(1,2)*M.el(2,3); - const double t27 = M.el(1,3)*M.el(2,2); - const double t29 = M.el(1,0)*M.el(0,1); - const double t32 = M.el(1,0)*M.el(2,1); - const double t33 = M.el(0,2)*M.el(3,3); - const double t35 = M.el(0,3)*M.el(3,2); - const double t37 = M.el(1,0)*M.el(3,1); - const double t38 = M.el(0,2)*M.el(2,3); - const double t40 = M.el(0,3)*M.el(2,2); - const double t42 = t14*t15-t14*t17-t19*t20+t19*t22+ - t24*t25-t24*t27-t29*t15+t29*t17+ - t32*t33-t32*t35-t37*t38+t37*t40; - const double t43 = M.el(2,0)*M.el(0,1); - const double t46 = M.el(2,0)*M.el(1,1); - const double t49 = M.el(2,0)*M.el(3,1); - const double t50 = M.el(0,2)*M.el(1,3); - const double t52 = M.el(0,3)*M.el(1,2); - const double t54 = M.el(3,0)*M.el(0,1); - const double t57 = M.el(3,0)*M.el(1,1); - const double t60 = M.el(3,0)*M.el(2,1); - const double t63 = t43*t20-t43*t22-t46*t33+t46*t35+ - t49*t50-t49*t52-t54*t25+t54*t27+ - t57*t38-t57*t40-t60*t50+t60*t52; - const double t65 = 1/(t42+t63); - const double t71 = M.el(0,2)*M.el(2,1); - const double t73 = M.el(0,3)*M.el(2,1); - const double t75 = M.el(0,2)*M.el(3,1); - const double t77 = M.el(0,3)*M.el(3,1); - const double t81 = M.el(0,1)*M.el(1,2); - const double t83 = M.el(0,1)*M.el(1,3); - const double t85 = M.el(0,2)*M.el(1,1); - const double t87 = M.el(0,3)*M.el(1,1); - const double t101 = M.el(1,0)*M.el(2,2); - const double t103 = M.el(1,0)*M.el(2,3); - const double t105 = M.el(2,0)*M.el(1,2); - const double t107 = M.el(2,0)*M.el(1,3); - const double t109 = M.el(3,0)*M.el(1,2); - const double t111 = M.el(3,0)*M.el(1,3); - const double t115 = M.el(0,0)*M.el(2,2); - const double t117 = M.el(0,0)*M.el(2,3); - const double t119 = M.el(2,0)*M.el(0,2); - const double t121 = M.el(2,0)*M.el(0,3); - const double t123 = M.el(3,0)*M.el(0,2); - const double t125 = M.el(3,0)*M.el(0,3); - const double t129 = M.el(0,0)*M.el(1,2); - const double t131 = M.el(0,0)*M.el(1,3); - const double t133 = M.el(1,0)*M.el(0,2); - const double t135 = M.el(1,0)*M.el(0,3); - el(0,0) = (M.el(1,1)*M.el(2,2)*M.el(3,3)-M.el(1,1)*M.el(2,3)*M.el(3,2)- - M.el(2,1)*M.el(1,2)*M.el(3,3)+M.el(2,1)*M.el(1,3)*M.el(3,2)+ - M.el(3,1)*M.el(1,2)*M.el(2,3)-M.el(3,1)*M.el(1,3)*M.el(2,2))*t65; - el(0,1) = -(M.el(0,1)*M.el(2,2)*M.el(3,3)-M.el(0,1)*M.el(2,3)*M.el(3,2)- - t71*M.el(3,3)+t73*M.el(3,2)+t75*M.el(2,3)-t77*M.el(2,2))*t65; - el(0,2) = (t81*M.el(3,3)-t83*M.el(3,2)-t85*M.el(3,3)+t87*M.el(3,2)+ - t75*M.el(1,3)-t77*M.el(1,2))*t65; - el(0,3) = -(t81*M.el(2,3)-t83*M.el(2,2)-t85*M.el(2,3)+t87*M.el(2,2)+ - t71*M.el(1,3)-t73*M.el(1,2))*t65; - el(1,0) = -(t101*M.el(3,3)-t103*M.el(3,2)-t105*M.el(3,3)+t107*M.el(3,2)+ - t109*M.el(2,3)-t111*M.el(2,2))*t65; - el(1,1) = (t115*M.el(3,3)-t117*M.el(3,2)-t119*M.el(3,3)+t121*M.el(3,2)+ - t123*M.el(2,3)-t125*M.el(2,2))*t65; - el(1,2) = -(t129*M.el(3,3)-t131*M.el(3,2)-t133*M.el(3,3)+t135*M.el(3,2)+ - t123*M.el(1,3)-t125*M.el(1,2))*t65; - el(1,3) = (t129*M.el(2,3)-t131*M.el(2,2)-t133*M.el(2,3)+t135*M.el(2,2)+ - t119*M.el(1,3)-t121*M.el(1,2))*t65; - el(2,0) = (t32*M.el(3,3)-t103*M.el(3,1)-t46*M.el(3,3)+t107*M.el(3,1)+ - t57*M.el(2,3)-t111*M.el(2,1))*t65; - el(2,1) = -(t19*M.el(3,3)-t117*M.el(3,1)-t43*M.el(3,3)+t121*M.el(3,1)+ - t54*M.el(2,3)-t125*M.el(2,1))*t65; - el(2,2) = (t14*M.el(3,3)-t131*M.el(3,1)-t29*M.el(3,3)+t135*M.el(3,1)+ - t54*M.el(1,3)-t125*M.el(1,1))*t65; - el(2,3) = -(t14*M.el(2,3)-t131*M.el(2,1)-t29*M.el(2,3)+t135*M.el(2,1)+ - t43*M.el(1,3)-t121*M.el(1,1))*t65; - el(3,0) = -(t32*M.el(3,2)-t101*M.el(3,1)-t46*M.el(3,2)+t105*M.el(3,1)+ - t57*M.el(2,2)-t109*M.el(2,1))*t65; - el(3,1) = (t19*M.el(3,2)-t115*M.el(3,1)-t43*M.el(3,2)+t119*M.el(3,1)+ - t54*M.el(2,2)-t123*M.el(2,1))*t65; - el(3,2) = -(t14*M.el(3,2)-t129*M.el(3,1)-t29*M.el(3,2)+t133*M.el(3,1)+ - t54*M.el(1,2)-t123*M.el(1,1))*t65; - el(3,3) = (t14*M.el(2,2)-t129*M.el(2,1)-t29*M.el(2,2)+t133*M.el(2,1)+ - t43*M.el(1,2)-t119*M.el(1,1))*t65; - } - }; -}; - - - -void dFMatrix::print_formatted (ostream &out, const unsigned int precision) const { - out.precision (precision); - out.setf (ios::scientific, ios::floatfield); // set output format - - for (unsigned int i=0; i max) - { - max = fabs(el(i,j)); - r = i; - } - } - Assert(max>1.e-16, ExcNotRegular()); - // rowinterchange - if (r>j) - { - for (k=0; k n, src.n() = m - Assert (dim_range <= dim_image, ExcDimensionMismatch(dim_range, dim_image)); - Assert (src.size() == dim_range, ExcDimensionMismatch(src.size(), dim_range)); - - for (unsigned int j=0 ; j n, m = src.n, n = dst.n - - householder(src); - backward(dst, src); - - double sum = 0.; - for (unsigned int i=n() ; i -#include -#include - -#include -#include -#include - - - -dSMatrixStruct::dSMatrixStruct () : - max_dim(0), - max_vec_len(0), - rowstart(0), - colnums(0) -{ - reinit (0,0,0); -}; - - - -dSMatrixStruct::dSMatrixStruct (const unsigned int m, const unsigned int n, - const unsigned int max_per_row) - : max_dim(0), - max_vec_len(0), - rowstart(0), - colnums(0) -{ - reinit (m,n,max_per_row); -}; - - - -dSMatrixStruct::dSMatrixStruct (const unsigned int n, - const unsigned int max_per_row) - : max_dim(0), - max_vec_len(0), - rowstart(0), - colnums(0) -{ - reinit (n,n,max_per_row); -}; - - - -dSMatrixStruct::~dSMatrixStruct () -{ - if (rowstart != 0) delete[] rowstart; - if (colnums != 0) delete[] colnums; -} - - - - -void -dSMatrixStruct::reinit (const unsigned int m, const unsigned int n, - const unsigned int max_per_row) -{ - Assert ((max_per_row>0) || ((m==0) && (n==0)), ExcInvalidNumber(max_per_row)); - rows = m; - cols = n; - vec_len = m * max_per_row; - max_row_len = max_per_row; - - // delete empty matrices - if ((m==0) || (n==0)) - { - if (rowstart) delete[] rowstart; - if (colnums) delete[] colnums; - rowstart = 0; - colnums = 0; - max_vec_len = vec_len = max_dim = rows = cols = 0; - compressed = false; - return; - }; - - if (rows > max_dim) - { - if (rowstart) delete[] rowstart; - max_dim = rows; - rowstart = new unsigned int[max_dim+1]; - }; - - if (vec_len > max_vec_len) - { - if (colnums) delete[] colnums; - max_vec_len = vec_len; - colnums = new int[max_vec_len]; - }; - - for (unsigned int i=0; i<=rows; i++) - rowstart[i] = i * max_per_row; - fill_n (&colnums[0], vec_len, -1); - - if (rows == cols) - for (unsigned int i=0;i(line)), - ExcInternalError()); - // assert that the first entry - // does not show up in - // the remaining ones and that - // the remaining ones are unique - // among themselves (this handles - // both cases, quadratic and - // rectangular matrices) - Assert (find (&colnums[rowstart[line]+1], - &colnums[next_row_start], - colnums[rowstart[line]]) == - &colnums[next_row_start], - ExcInternalError()); - Assert (adjacent_find(&colnums[rowstart[line]+1], - &colnums[next_row_start]) == - &colnums[next_row_start], - ExcInternalError()); - }; - - vec_len = rowstart[rows] = next_row_start; - compressed = true; - - delete[] tmp_entries; -}; - - - -bool -dSMatrixStruct::empty () const { - // let's try to be on the safe side of - // life by using multiple possibilities in - // the check for emptiness... (sorry for - // this kludge -- emptying matrices and - // freeing memory was not present in the - // original implementation and I don't - // know at how many places I missed - // something in adding it, so I try to - // be cautious. wb) - if ((rowstart==0) || (rows==0) || (cols==0)) - { - Assert (rowstart==0, ExcInternalError()); - Assert (rows==0, ExcInternalError()); - Assert (cols==0, ExcInternalError()); - Assert (colnums==0, ExcInternalError()); - Assert (vec_len==0, ExcInternalError()); - Assert (max_vec_len==0, ExcInternalError()); - Assert (vec_len==0, ExcInternalError()); - - return true; - }; - return false; -}; - - - -int -dSMatrixStruct::operator () (const unsigned int i, const unsigned int j) const -{ - Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject()); - Assert (i(j) == colnums[rowstart[i]]) - return rowstart[i]; - } - else - // no first entry exists for this - // line - return -1; - - // all other entries are sorted, so - // we can use a binary seach algorithm - const int* p = lower_bound (&colnums[rowstart[i]+1], - &colnums[rowstart[i+1]], - static_cast(j)); - if ((*p == static_cast(j)) && - (p != &colnums[rowstart[i+1]])) - return (p - &colnums[0]); - else - return -1; -} - - -void -dSMatrixStruct::add (const unsigned int i, const unsigned int j) -{ - Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject()); - Assert (i=0) - out << i << " " << -colnums[j] << endl; - - AssertThrow (out, ExcIO()); -} - - - -unsigned int -dSMatrixStruct::bandwidth () const -{ - Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject()); - unsigned int b=0; - for (unsigned int i=0; i=0) - { - if (static_cast(abs(static_cast(i-colnums[j]))) > b) - b = abs(static_cast(i-colnums[j])); - } - else - // leave if at the end of - // the entries of this line - break; - return b; -}; - - - -unsigned int -dSMatrixStruct::n_nonzero_elements () const { - Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject()); - Assert (compressed, ExcNotCompressed()); - return colnums[rows]-colnums[0]; -}; - - - - - -/*-------------------------------------------------------------------------*/ - - -dSMatrix::dSMatrix () : - cols(0), - val(0), - max_len(0) {}; - - - -dSMatrix::dSMatrix (const dSMatrixStruct &c) - : cols(&c), val(0), max_len(0) -{ - reinit(); -}; - - - -dSMatrix::~dSMatrix () -{ - delete[] val; -}; - - - -void -dSMatrix::reinit () -{ - Assert (cols != 0, ExcMatrixNotInitialized()); - Assert (cols->compressed || cols->empty(), ExcNotCompressed()); - - if (cols->empty()) - { - if (val) delete[] val; - val = 0; - max_len = 0; - return; - }; - - if (max_lenvec_len) - { - if (val) delete[] val; - val = new double[cols->vec_len]; - max_len = cols->vec_len; - }; - - if (val) - fill_n (&val[0], cols->vec_len, 0); -} - - - -void -dSMatrix::reinit (const dSMatrixStruct &sparsity) { - cols = &sparsity; - reinit (); -}; - - - -void -dSMatrix::clear () { - cols = 0; - if (val) delete[] val; - val = 0; - max_len = 0; -}; - - - -unsigned int -dSMatrix::n_nonzero_elements () const { - Assert (cols != 0, ExcMatrixNotInitialized()); - return cols->n_nonzero_elements (); -}; - - - -dSMatrix & -dSMatrix::copy_from (const dSMatrix &matrix) { - Assert (cols != 0, ExcMatrixNotInitialized()); - Assert (val != 0, ExcMatrixNotInitialized()); - Assert (cols == matrix.cols, ExcDifferentSparsityPatterns()); - - double *val_ptr = &val[0]; - const double *matrix_ptr = &matrix.val[0]; - const double *const end_ptr = &val[cols->vec_len]; - - while (val_ptr != end_ptr) - *val_ptr++ = *matrix_ptr++; - - return *this; -}; - - - -void -dSMatrix::add_scaled (const double factor, const dSMatrix &matrix) { - Assert (cols != 0, ExcMatrixNotInitialized()); - Assert (val != 0, ExcMatrixNotInitialized()); - Assert (cols == matrix.cols, ExcDifferentSparsityPatterns()); - - double *val_ptr = &val[0]; - const double *matrix_ptr = &matrix.val[0]; - const double *const end_ptr = &val[cols->vec_len]; - - while (val_ptr != end_ptr) - *val_ptr++ += factor * *matrix_ptr++; -}; - - - -void -dSMatrix::vmult (dVector& dst, const dVector& src) const -{ - Assert (cols != 0, ExcMatrixNotInitialized()); - Assert (val != 0, ExcMatrixNotInitialized()); - Assert(m() == dst.size(), ExcDimensionsDontMatch(m(),dst.size())); - Assert(n() == src.size(), ExcDimensionsDontMatch(n(),src.size())); - - const unsigned int n_rows = m(); - const double *val_ptr = &val[cols->rowstart[0]]; - const int *colnum_ptr = &cols->colnums[cols->rowstart[0]]; - double *dst_ptr = &dst(0); - for (unsigned int row=0; rowrowstart[row+1]]; - while (val_ptr != val_end_of_row) - s += *val_ptr++ * src(*colnum_ptr++); - *dst_ptr++ = s; - }; -}; - - -void -dSMatrix::Tvmult (dVector& dst, const dVector& src) const -{ - Assert (val != 0, ExcMatrixNotInitialized()); - Assert (cols != 0, ExcMatrixNotInitialized()); - Assert(n() == dst.size(), ExcDimensionsDontMatch(n(),dst.size())); - Assert(m() == src.size(), ExcDimensionsDontMatch(m(),src.size())); - - dst.clear (); - - for (unsigned int i=0;irowstart[i]; jrowstart[i+1] ;j++) - { - int p = cols->colnums[j]; - dst(p) += val[j] * src(i); - } - } -} - - - -double -dSMatrix::matrix_norm (const dVector& v) const -{ - Assert (cols != 0, ExcMatrixNotInitialized()); - Assert (val != 0, ExcMatrixNotInitialized()); - Assert(m() == v.size(), ExcDimensionsDontMatch(m(),v.size())); - Assert(n() == v.size(), ExcDimensionsDontMatch(n(),v.size())); - - double sum = 0.; - const unsigned int n_rows = m(); - const double *val_ptr = &val[cols->rowstart[0]]; - const int *colnum_ptr = &cols->colnums[cols->rowstart[0]]; - for (unsigned int row=0; rowrowstart[row+1]]; - while (val_ptr != val_end_of_row) - s += *val_ptr++ * v(*colnum_ptr++); - - sum += s* v(row); - }; - - return sum; -}; - - - -double -dSMatrix::residual (dVector& dst, const dVector& u, const dVector& b) const -{ - Assert (cols != 0, ExcMatrixNotInitialized()); - Assert (val != 0, ExcMatrixNotInitialized()); - Assert(m() == dst.size(), ExcDimensionsDontMatch(m(),dst.size())); - Assert(m() == b.size(), ExcDimensionsDontMatch(m(),b.size())); - Assert(n() == u.size(), ExcDimensionsDontMatch(n(),u.size())); - - double s,norm=0.; - - for (unsigned int i=0;irowstart[i]; jrowstart[i+1] ;j++) - { - int p = cols->colnums[j]; - s -= val[j] * u(p); - } - dst(i) = s; - norm += dst(i)*dst(i); - } - return sqrt(norm); -} - -void -dSMatrix::precondition_Jacobi (dVector& dst, const dVector& src, - const double om) const -{ - Assert (cols != 0, ExcMatrixNotInitialized()); - Assert (val != 0, ExcMatrixNotInitialized()); - Assert (m() == n(), ExcMatrixNotSquare()); - - const unsigned int n = src.size(); - double *dst_ptr = dst.begin(); - const double *src_ptr = src.begin(); - const unsigned int *rowstart_ptr = &cols->rowstart[0]; - - for (unsigned int i=0; irowstart[0]; - double *dst_ptr = &dst(0); - - for (unsigned int row=0; rowcolnums[*rowstart_ptr+1], - &cols->colnums[*(rowstart_ptr+1)], - static_cast(row)) - - &cols->colnums[0]); - - for (unsigned int j=(*rowstart_ptr)+1; jcolnums[j]); - *dst_ptr /= val[*rowstart_ptr]; - }; - - rowstart_ptr = &cols->rowstart[0]; - dst_ptr = &dst(0); - for (unsigned int row=0; rowrowstart[n-1]; - dst_ptr = &dst(n-1); - for (int row=n-1; row>=0; --row, --rowstart_ptr, --dst_ptr) - { - const unsigned int first_right_of_diagonal_index - = (lower_bound (&cols->colnums[*rowstart_ptr+1], - &cols->colnums[*(rowstart_ptr+1)], - static_cast(row)) - - &cols->colnums[0]); - for (unsigned int j=first_right_of_diagonal_index; j<*(rowstart_ptr+1); ++j) - if (cols->colnums[j] > row) - *dst_ptr -= om* val[j] * dst(cols->colnums[j]); - - *dst_ptr /= val[*rowstart_ptr]; - }; -} - -void -dSMatrix::precondition_SOR (dVector& dst, const dVector& src, - const double om) const -{ - Assert (cols != 0, ExcMatrixNotInitialized()); - Assert (val != 0, ExcMatrixNotInitialized()); - Assert (m() == n(), ExcMatrixNotSquare()); - - dst = src; - SOR(dst,om); -}; - - -void dSMatrix::precondition (dVector &dst, const dVector &src) const { - Assert (m() == n(), ExcMatrixNotSquare()); - dst=src; -}; - - -void -dSMatrix::SOR (dVector& dst, const double om) const -{ - Assert (cols != 0, ExcMatrixNotInitialized()); - Assert (val != 0, ExcMatrixNotInitialized()); - Assert (m() == n(), ExcMatrixNotSquare()); - Assert (m() == dst.size(), ExcDimensionsDontMatch(m(),dst.size())); - - for (unsigned int row=0; rowrowstart[row]; jrowstart[row+1]; ++j) - if ((unsigned int)cols->colnums[j] < row) - s -= val[j] * dst(cols->colnums[j]); - - dst(row) = s * om / val[cols->rowstart[row]]; - } -} - -void -dSMatrix::SSOR (dVector& dst, const double om) const -{ - Assert (cols != 0, ExcMatrixNotInitialized()); - Assert (val != 0, ExcMatrixNotInitialized()); - - int p; - const unsigned int n = dst.size(); - unsigned int j; - double s; - - for (unsigned int i=0; irowstart[i]; jrowstart[i+1] ;j++) - { - p = cols->colnums[j]; - if (p>=0) - { - if (i>j) s += val[j] * dst(p); - } - } - dst(i) -= s * om; - dst(i) /= val[cols->rowstart[i]]; - } - - for (int i=n-1; i>=0; i--) // this time, i is signed, but alsways positive! - { - s = 0.; - for (j=cols->rowstart[i]; jrowstart[i+1] ;j++) - { - p = cols->colnums[j]; - if (p>=0) - { - if ((unsigned int)irowstart[i]]; - } -} - - - -const dSMatrixStruct & dSMatrix::get_sparsity_pattern () const { - Assert (cols != 0, ExcMatrixNotInitialized()); - return *cols; -}; - - - -void dSMatrix::print (ostream &out) const { - Assert (cols != 0, ExcMatrixNotInitialized()); - Assert (val != 0, ExcMatrixNotInitialized()); - - for (unsigned int i=0; irows; ++i) - for (unsigned int j=cols->rowstart[i]; jrowstart[i+1]; ++j) - out << "(" << i << "," << cols->colnums[j] << ") " << val[j] << endl; - - AssertThrow (out, ExcIO()); -}; - - - -void dSMatrix::print_formatted (ostream &out, const unsigned int precision) const { - Assert (cols != 0, ExcMatrixNotInitialized()); - Assert (val != 0, ExcMatrixNotInitialized()); - out.precision (precision); - out.setf (ios::scientific, ios::floatfield); // set output format - - for (unsigned int i=0; ioperator()(i,j)] << ' '; - else - out << setw(precision+8) << " "; - out << endl; - }; - AssertThrow (out, ExcIO()); - - out.setf (0, ios::floatfield); // reset output format -}; - diff --git a/deal.II/lac/source/dvector.cc b/deal.II/lac/source/dvector.cc deleted file mode 100644 index 481bf5ffc9..0000000000 --- a/deal.II/lac/source/dvector.cc +++ /dev/null @@ -1,520 +0,0 @@ -// $Id$ - -// This file is part of the DEAL Library -// DEAL is Copyright(1995) by -// Roland Becker, Guido Kanschat, Franz-Theo Suttmeier - -#include -#include -#include - - -static inline double sqr (const double x) { - return x*x; -}; - - - -dVector::dVector () : - dim(0), - maxdim(0), - val(0) -{} - - -dVector::dVector (const unsigned int n) : - dim(0), - maxdim(0), - val(0) -{ - reinit (n, false); -} - - -dVector::dVector (const dVector& v) : - dim(v.size()), - maxdim(v.size()), - val(0) -{ - if (dim) - { - val = new double[maxdim]; - Assert (val != 0, ExcOutOfMemory()); - copy (v.begin(), v.end(), begin()); - } -} - - - -void dVector::reinit (const unsigned int n, const bool fast) { - if (n==0) - { - if (val) delete[] val; - val = 0; - maxdim = dim = 0; - return; - }; - - if (n>maxdim) - { - if (val) delete[] val; - val = new double[n]; - Assert (val != 0, ExcOutOfMemory()); - maxdim = n; - }; - dim = n; - if (fast == false) - clear (); -} - - - -void dVector::reinit (const dVector& v, const bool fast) { - reinit (v.size(), fast); -}; - - - - -dVector::~dVector () -{ - if (val) delete[] val; -} - - - -void dVector::clear () { - if (dim>0) - fill (begin(), end(), 0.); -} - - - -bool dVector::all_zero () const { - Assert (dim!=0, ExcEmptyVector()); - - const_iterator p = begin(), - e = end(); - while (p!=e) - if (*p++ != 0.0) - return false; - return true; -}; - - - -double dVector::operator * (const dVector& v) const -{ - Assert (dim!=0, ExcEmptyVector()); - - if (&v == this) - return norm_sqr(); - - Assert (dim == v.dim, ExcDimensionsDontMatch(dim, v.dim)); - - double sum0 = 0, - sum1 = 0, - sum2 = 0, - sum3 = 0; - - // use modern processors better by - // allowing pipelined commands to be - // executed in parallel - const_iterator ptr = begin(), - vptr = v.begin(), - eptr = ptr + (dim/4)*4; - while (ptr!=eptr) - { - sum0 += (*ptr++ * *vptr++); - sum1 += (*ptr++ * *vptr++); - sum2 += (*ptr++ * *vptr++); - sum3 += (*ptr++ * *vptr++); - }; - // add up remaining elements - while (ptr != end()) - sum0 += *ptr++ * *vptr++; - - return sum0+sum1+sum2+sum3; -} - - - -double dVector::norm_sqr () const -{ - Assert (dim!=0, ExcEmptyVector()); - - double sum0 = 0, - sum1 = 0, - sum2 = 0, - sum3 = 0; - - // use modern processors better by - // allowing pipelined commands to be - // executed in parallel - const_iterator ptr = begin(), - eptr = ptr + (dim/4)*4; - while (ptr!=eptr) - { - sum0 += sqr(*ptr++); - sum1 += sqr(*ptr++); - sum2 += sqr(*ptr++); - sum3 += sqr(*ptr++); - }; - // add up remaining elements - while (ptr != end()) - sum0 += sqr(*ptr++); - - return sum0+sum1+sum2+sum3; -}; - - - -double dVector::mean_value () const -{ - Assert (dim!=0, ExcEmptyVector()); - - double sum0 = 0, - sum1 = 0, - sum2 = 0, - sum3 = 0; - - // use modern processors better by - // allowing pipelined commands to be - // executed in parallel - const_iterator ptr = begin(), - eptr = ptr + (dim/4)*4; - while (ptr!=eptr) - { - sum0 += *ptr++; - sum1 += *ptr++; - sum2 += *ptr++; - sum3 += *ptr++; - }; - // add up remaining elements - while (ptr != end()) - sum0 += *ptr++; - - return (sum0+sum1+sum2+sum3)/size(); -}; - - - -double dVector::l1_norm () const -{ - Assert (dim!=0, ExcEmptyVector()); - - double sum0 = 0, - sum1 = 0, - sum2 = 0, - sum3 = 0; - - // use modern processors better by - // allowing pipelined commands to be - // executed in parallel - const_iterator ptr = begin(), - eptr = ptr + (dim/4)*4; - while (ptr!=eptr) - { - sum0 += fabs(*ptr++); - sum1 += fabs(*ptr++); - sum2 += fabs(*ptr++); - sum3 += fabs(*ptr++); - }; - // add up remaining elements - while (ptr != end()) - sum0 += fabs(*ptr++); - - return sum0+sum1+sum2+sum3; -}; - - - -double dVector::l2_norm () const -{ - return sqrt(norm_sqr()); -}; - - - -double dVector::linfty_norm () const { - Assert (dim!=0, ExcEmptyVector()); - - double max0=0., - max1=0., - max2=0., - max3=0.; - for (unsigned int i=0; i<(dim/4); ++i) - { - if (max0