From 567040827945d6cb317b31793bf035efb155febc Mon Sep 17 00:00:00 2001 From: guido Date: Tue, 23 Feb 1999 14:46:23 +0000 Subject: [PATCH] tildas removed git-svn-id: https://svn.dealii.org/trunk@881 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/fullmatrix.h~ | 577 ------- .../lac/include/lac/fullmatrix.templates.h~ | 1367 ----------------- 2 files changed, 1944 deletions(-) delete mode 100644 deal.II/lac/include/lac/fullmatrix.h~ delete mode 100644 deal.II/lac/include/lac/fullmatrix.templates.h~ diff --git a/deal.II/lac/include/lac/fullmatrix.h~ b/deal.II/lac/include/lac/fullmatrix.h~ deleted file mode 100644 index 6d2f44ff11..0000000000 --- a/deal.II/lac/include/lac/fullmatrix.h~ +++ /dev/null @@ -1,577 +0,0 @@ -/*---------------------------- fmatrix.h ---------------------------*/ -/* $Id$ */ -#ifndef __lac_fullmatrix_H -#define __lac_fullmatrix_H -/*---------------------------- fmatrix.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 iVector; - - - -/** - * Rectangular/quadratic 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 - * - */ -template -class FullMatrix -{ - private: - /** - * Component-array. - */ - number* 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. - */ - number& el (const unsigned int i, const unsigned int j); - - /** - * Return the value of the element #(i,j)#. - * - * This function does no bounds checking. - */ - number el (const unsigned int i, const unsigned int j) const; - - - public: - /** - * Constructor. Initialize the matrix as - * a square matrix with dimension #n#. - */ - explicit FMatrix (const unsigned int n = 1); - - /** - * Constructor. Initialize the matrix as - * a rectangular #m# times #n# matrix. - */ - FMatrix (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!! - */ - explicit FMatrix (const FMatrix&); - - /** - * Destructor. Release all memory. - */ - ~FMatrix(); - - /** - * 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 FMatrix &) const; - - /** - * A = B . Copy all elements - */ - template - FMatrix& operator = (const FMatrix& B); - - - /** - * U(0-m,0-n) = s . Fill all elements - */ - template - void fill (const FMatrix& 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 FMatrix &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) ) - */ - number 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. - */ - 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 #el(i,j)# function - * but does bounds checking. - */ - number& 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 - */ - template - void add (const number s, const FMatrix& B); - - /** - * A+=Transp(B). - * Simple addition of the transpose of B to this - */ - template - void Tadd (const number s, const FMatrix& B); - - /** - * C=A*B. - * Matrix-matrix-multiplication - */ - - template - void mmult (FMatrix& C, const FMatrix& B) const; - - /** - * C=Transp(A)*B. - * Matrix-matrix-multiplication using - * transpose of this - */ - template - void Tmmult (FMatrix& C, const FMatrix& B) const; - - /** - * w (+)= A*v. - * Matrix-vector-multiplication ;

- * ( application of this to a vector v ) - * flag adding=true : w+=A*v - */ - template - void vmult (Vector& w, const Vector& 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 - */ - template - void Tvmult (Vector& w, const Vector& 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. - */ - template - double matrix_norm (const Vector &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. - */ - template - double matrix_scalar_product (const Vector &u, const Vector &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 FMatrix &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 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 - */ - 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 - */ - 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 - */ - void add_col (const unsigned int i, - const number s, const unsigned int j, - const number 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| - */ - template - double residual (Vector& w, const Vector& v, const Vector& b) const; - - /** - * Inversion of lower triangle . - */ - template - void forward (Vector& dst, const Vector& src) const; - - /** - * Inversion of upper triangle . - */ - template - void backward (Vector& dst, const Vector& 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. - */ - template - void householder (Vector& y); - - /** - * Least - Squares - Approximation by QR-factorization. - */ - template - number least_squares (Vector& dst, Vector& src); - - /** - * A(i,i)+=B(i,1-n). Addition of complete - * rows of B to diagonal-elements of this ;

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

- * (used elements of v determined by x) - */ - template - void gsmult (Vector& w, const Vector& 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 (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 -------------------------------*/ - -template -inline number & -FMatrix::el (const unsigned int i, const unsigned int j) -{ - return val[i*dim_range+j]; -}; - - -template -inline number -FMatrix::el (const unsigned int i, const unsigned int j) const -{ - return val[i*dim_range+j]; -}; - - -template -inline unsigned int -FMatrix::m() const -{ - return dim_image; -}; - - -template -inline unsigned int -FMatrix::n() const -{ - return dim_range; -}; - - -template -inline number -FMatrix::el (const unsigned int i) const -{ - return val[i]; -}; - - -template -inline number -FMatrix::operator() (const unsigned int i, const unsigned int j) const -{ - Assert (i -inline number & -FMatrix::operator() (const unsigned int i, const unsigned int j) -{ - Assert (i -#include -#include - -#include -#include -#include -#include - - -template -dFMatrix::dFMatrix (const unsigned int n) { - init (n,n); -}; - - -template -dFMatrix::dFMatrix (const unsigned int m, const unsigned int n) { - init (m,n); -}; - - -template -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++; -}; - - -template -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 (); -}; - - -template -dFMatrix::~dFMatrix () { - delete[] val; -}; - - -template -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; -}; - - -template -void dFMatrix::reinit (const unsigned int mm, const unsigned int nn) -{ - if (val_size -void dFMatrix::reinit (const unsigned int n) { - reinit (n, n); -}; - - -template -void dFMatrix::reinit (const dFMatrix &B) { - reinit (B.m(), B.n()); -}; - - -template -void dFMatrix::vmult (dVector& dst, const dVector& src, - const bool adding) const -{ - Assert(dst.size() == m(), ExcDimensionMismatch(dst.size(), m())); - Assert(src.size() == n(), ExcDimensionMismatch(src.size(), n())); - - double s; - if ((n()==3) && (m()==3)) - { - double s0,s1,s2; - s = src(0); - s0 = s*val[0]; s1 = s*val[3]; s2 = s*val[6]; - s = src(1); - s0 += s*val[1]; s1 += s*val[4]; s2 += s*val[7]; - s = src(2); - s0 += s*val[2]; s1 += s*val[5]; s2 += s*val[8]; - - if (!adding) - { - dst(0) = s0; - dst(1) = s1; - dst(2) = s2; - } - else - { - dst(0) += s0; - dst(1) += s1; - dst(2) += s2; - } - } - else if ((n()==4) && (m()==4)) - { - double s0,s1,s2,s3; - s = src(0); - s0 = s*val[0]; s1 = s*val[4]; s2 = s*val[8]; s3 = s*val[12]; - s = src(1); - s0 += s*val[1]; s1 += s*val[5]; s2 += s*val[9]; s3 += s*val[13]; - s = src(2); - s0 += s*val[2]; s1 += s*val[6]; s2 += s*val[10]; s3 += s*val[14]; - s = src(3); - s0 += s*val[3]; s1 += s*val[7]; s2 += s*val[11]; s3 += s*val[15]; - - if (!adding) - { - dst(0) = s0; - dst(1) = s1; - dst(2) = s2; - dst(3) = s3; - } - else - { - dst(0) += s0; - dst(1) += s1; - dst(2) += s2; - dst(3) += s3; - } - } - else if ((n()==8) && (m()==8)) - { - double s0,s1,s2,s3,s4,s5,s6,s7; - s = src(0); - s0 = s*val[0]; s1 = s*val[8]; s2 = s*val[16]; s3 = s*val[24]; - s4 = s*val[32]; s5 = s*val[40]; s6 = s*val[48]; s7 = s*val[56]; - s = src(1); - s0 += s*val[1]; s1 += s*val[9]; s2 += s*val[17]; s3 += s*val[25]; - s4 += s*val[33]; s5 += s*val[41]; s6 += s*val[49]; s7 += s*val[57]; - s = src(2); - s0 += s*val[2]; s1 += s*val[10]; s2 += s*val[18]; s3 += s*val[26]; - s4 += s*val[34]; s5 += s*val[42]; s6 += s*val[50]; s7 += s*val[58]; - s = src(3); - s0 += s*val[3]; s1 += s*val[11]; s2 += s*val[19]; s3 += s*val[27]; - s4 += s*val[35]; s5 += s*val[43]; s6 += s*val[51]; s7 += s*val[59]; - s = src(4); - s0 += s*val[4]; s1 += s*val[12]; s2 += s*val[20]; s3 += s*val[28]; - s4 += s*val[36]; s5 += s*val[44]; s6 += s*val[52]; s7 += s*val[60]; - s = src(5); - s0 += s*val[5]; s1 += s*val[13]; s2 += s*val[21]; s3 += s*val[29]; - s4 += s*val[37]; s5 += s*val[45]; s6 += s*val[53]; s7 += s*val[61]; - s = src(6); - s0 += s*val[6]; s1 += s*val[14]; s2 += s*val[22]; s3 += s*val[30]; - s4 += s*val[38]; s5 += s*val[46]; s6 += s*val[54]; s7 += s*val[62]; - s = src(7); - s0 += s*val[7]; s1 += s*val[15]; s2 += s*val[23]; s3 += s*val[31]; - s4 += s*val[39]; s5 += s*val[47]; s6 += s*val[55]; s7 += s*val[63]; - - if (!adding) - { - dst(0) = s0; - dst(1) = s1; - dst(2) = s2; - dst(3) = s3; - dst(4) = s4; - dst(5) = s5; - dst(6) = s6; - dst(7) = s7; - } - else - { - dst(0) += s0; - dst(1) += s1; - dst(2) += s2; - dst(3) += s3; - dst(4) += s4; - dst(5) += s5; - dst(6) += s6; - dst(7) += s7; - } - } - else - { - double* e = val; - const unsigned int size_m = m(), - size_n = n(); - for (unsigned int i=0; i -void dFMatrix::gsmult (dVector& dst, const dVector& src, const iVector& gl) const -{ - Assert(n() == m(), ExcNotQuadratic()); - Assert(dst.size() == n(), ExcDimensionMismatch(dst.size(), n())); - Assert(src.size() == n(), ExcDimensionMismatch(src.size(), n())); - Assert(gl.n() == n(), ExcDimensionMismatch(gl.n(), n())); - - double s; - if ((n()==3) && (m()==3)) - { - double s0=0.,s1=0.,s2=0.; - s = src(0); - if(gl(1) -void dFMatrix::Tvmult (dVector& dst, const dVector& src, const bool adding) const -{ - Assert(dst.size() == n(), ExcDimensionMismatch(dst.size(), n())); - Assert(src.size() == m(), ExcDimensionMismatch(src.size(), m())); - - unsigned int i,j; - double s; - const unsigned int size_m = m(), - size_n = n(); - for (i=0; i -double dFMatrix::residual (dVector& dst, const dVector& src, - const dVector& right) const -{ - Assert(dst.size() == m(), ExcDimensionMismatch(dst.size(), m())); - Assert(src.size() == n(), ExcDimensionMismatch(src.size(), n())); - Assert(right.size() == m(), ExcDimensionMismatch(right.size(), m())); - - unsigned int i,j; - double s, res = 0.; - const unsigned int size_m = m(), - size_n = n(); - for (i=0; i -void dFMatrix::forward (dVector& dst, const dVector& src) const -{ - Assert(n() == m(), ExcNotQuadratic()); - Assert(dst.size() == n(), ExcDimensionMismatch(dst.size(), n())); - Assert(src.size() == n(), ExcDimensionMismatch(src.size(), n())); - - unsigned int i,j; - unsigned int nu = (m() -void dFMatrix::backward (dVector& dst, const dVector& src) const -{ - Assert(n() == m(), ExcNotQuadratic()); - Assert(dst.size() == n(), ExcDimensionMismatch(dst.size(), n())); - Assert(src.size() == n(), ExcDimensionMismatch(src.size(), n())); - - unsigned int j; - unsigned int nu = (m()=0; --i) - { - s = src(i); - for (j=i+1; j -dFMatrix& -dFMatrix::operator = (const dFMatrix& m) { - reinit(m); - - double * p = &val[0]; - const double * vp = &m.val[0]; - const double * const e = &val[dim_image*dim_range]; - - while (p!=e) - *p++ = *vp++; - - return *this; -} - -template -void dFMatrix::fill (const dFMatrix& src, - const unsigned int i, const unsigned int j) -{ - Assert (n() >= src.n() + j, ExcInvalidDestination(n(), src.n(), j)); - Assert (m() >= src.m() + i, ExcInvalidDestination(m(), src.m(), i)); - - for (unsigned int ii=0; ii -void dFMatrix::add_row (const unsigned int i, - const double s, const unsigned int j) -{ - for (unsigned int k=0; k -void dFMatrix::add_row (const unsigned int i, const double s, - const unsigned int j, const double t, - const unsigned int k) -{ - const unsigned int size_m = m(); - for (unsigned l=0; l -void dFMatrix::add_col (const unsigned int i, const double s, - const unsigned int j) -{ - for (unsigned int k=0; k -void dFMatrix::add_col (const unsigned int i, const double s, - const unsigned int j, const double t, - const unsigned int k) -{ - for (unsigned int l=0; l -void dFMatrix::swap_row (const unsigned int i, const unsigned int j) -{ - double s; - for (unsigned int k=0; k -void dFMatrix::swap_col (const unsigned int i, const unsigned int j) -{ - double s; - for (unsigned int k=0; k -void dFMatrix::diagadd (const double& src) -{ - Assert (m() == n(), ExcDimensionMismatch(m(),n())); - for (unsigned int i=0; i -void dFMatrix::mmult (dFMatrix& dst, const dFMatrix& src) const -{ - Assert (n() == src.m(), ExcDimensionMismatch(n(), src.m())); - unsigned int i,j,k; - double s = 1.; - dst.reinit(m(), src.n()); - - for (i=0;i -void dFMatrix::Tmmult (dFMatrix& dst, const dFMatrix& src) const -{ - Assert (n() == src.n(), ExcDimensionMismatch(n(), src.n())); - - unsigned int i,j,k; - double s = 1.; - dst.reinit(m(), src.m()); - - for (i=0;i -double dFMatrix::matrix_norm (const dVector &v) const { - Assert(m() == v.size(), ExcDimensionMismatch(m(),v.size())); - Assert(n() == v.size(), ExcDimensionMismatch(n(),v.size())); - - double sum = 0.; - const unsigned int n_rows = m(); - const double *val_ptr = &val[0]; - const double *v_ptr; - - for (unsigned int row=0; row -double dFMatrix::matrix_scalar_product (const dVector &u, const dVector &v) const { - Assert(m() == u.size(), ExcDimensionMismatch(m(),v.size())); - Assert(n() == v.size(), ExcDimensionMismatch(n(),v.size())); - - double sum = 0.; - const unsigned int n_rows = m(); - const unsigned int n_cols = n(); - const double *val_ptr = &val[0]; - const double *v_ptr; - - for (unsigned int row=0; row -void dFMatrix::print (ostream& s, int w, int p) const -{ - unsigned int i,j; - for (i=0;i -void dFMatrix::add (const double s,const dFMatrix& src) -{ - Assert (m() == src.m(), ExcDimensionMismatch(m(), src.m())); - Assert (n() == src.n(), ExcDimensionMismatch(n(), src.n())); - if ((n()==3) && (m()==3)) - { - val[0] += s * src.el(0); - val[1] += s * src.el(1); - val[2] += s * src.el(2); - val[3] += s * src.el(3); - val[4] += s * src.el(4); - val[5] += s * src.el(5); - val[6] += s * src.el(6); - val[7] += s * src.el(7); - val[8] += s * src.el(8); - } - else if ((n()==4) && (m()==4)) - { - val[0] += s * src.el(0); - val[1] += s * src.el(1); - val[2] += s * src.el(2); - val[3] += s * src.el(3); - val[4] += s * src.el(4); - val[5] += s * src.el(5); - val[6] += s * src.el(6); - val[7] += s * src.el(7); - val[8] += s * src.el(8); - val[9] += s * src.el(9); - val[10] += s * src.el(10); - val[11] += s * src.el(11); - val[12] += s * src.el(12); - val[13] += s * src.el(13); - val[14] += s * src.el(14); - val[15] += s * src.el(15); - } - else if ((n()==8) && (m()==8)) - { - val[0] += s * src.el(0); - val[1] += s * src.el(1); - val[2] += s * src.el(2); - val[3] += s * src.el(3); - val[4] += s * src.el(4); - val[5] += s * src.el(5); - val[6] += s * src.el(6); - val[7] += s * src.el(7); - val[8] += s * src.el(8); - val[9] += s * src.el(9); - val[10] += s * src.el(10); - val[11] += s * src.el(11); - val[12] += s * src.el(12); - val[13] += s * src.el(13); - val[14] += s * src.el(14); - val[15] += s * src.el(15); - val[16] += s * src.el(16); - val[17] += s * src.el(17); - val[18] += s * src.el(18); - val[19] += s * src.el(19); - - val[20] += s * src.el(20); - val[21] += s * src.el(21); - val[22] += s * src.el(22); - val[23] += s * src.el(23); - val[24] += s * src.el(24); - val[25] += s * src.el(25); - val[26] += s * src.el(26); - val[27] += s * src.el(27); - val[28] += s * src.el(28); - val[29] += s * src.el(29); - - val[30] += s * src.el(30); - val[31] += s * src.el(31); - val[32] += s * src.el(32); - val[33] += s * src.el(33); - val[34] += s * src.el(34); - val[35] += s * src.el(35); - val[36] += s * src.el(36); - val[37] += s * src.el(37); - val[38] += s * src.el(38); - val[39] += s * src.el(39); - - val[40] += s * src.el(40); - val[41] += s * src.el(41); - val[42] += s * src.el(42); - val[43] += s * src.el(43); - val[44] += s * src.el(44); - val[45] += s * src.el(45); - val[46] += s * src.el(46); - val[47] += s * src.el(47); - val[48] += s * src.el(48); - val[49] += s * src.el(49); - - val[50] += s * src.el(50); - val[51] += s * src.el(51); - val[52] += s * src.el(52); - val[53] += s * src.el(53); - val[54] += s * src.el(54); - val[55] += s * src.el(55); - val[56] += s * src.el(56); - val[57] += s * src.el(57); - val[58] += s * src.el(58); - val[59] += s * src.el(59); - - val[60] += s * src.el(60); - val[61] += s * src.el(61); - val[62] += s * src.el(62); - val[63] += s * src.el(63); - } - else - { - const unsigned int size = n()*m(); - for (unsigned int i=0; i -void dFMatrix::add_diag (const double s, const dFMatrix& src) -{ - Assert (m() == src.m(), ExcDimensionMismatch(m(), src.m())); - Assert (n() == src.n(), ExcDimensionMismatch(n(), src.n())); - - if ((n()==3) && (m()==3)) - { - val[0] += s * src.el(0); - val[0] += s * src.el(1); - val[0] += s * src.el(2); - val[3] += s * src.el(3); - val[3] += s * src.el(4); - val[3] += s * src.el(5); - val[6] += s * src.el(6); - val[6] += s * src.el(7); - val[6] += s * src.el(8); - } - else if ((n()==4) && (m()==4)) - { - val[0] += s * src.el(0); - val[0] += s * src.el(1); - val[0] += s * src.el(2); - val[0] += s * src.el(3); - val[4] += s * src.el(4); - val[4] += s * src.el(5); - val[4] += s * src.el(6); - val[4] += s * src.el(7); - val[8] += s * src.el(8); - val[8] += s * src.el(9); - val[8] += s * src.el(10); - val[8] += s * src.el(11); - val[12] += s * src.el(12); - val[12] += s * src.el(13); - val[12] += s * src.el(14); - val[12] += s * src.el(15); - } - else if ((n()==8) && (m()==8)) - { - val[0] += s * src.el(0); - val[0] += s * src.el(1); - val[0] += s * src.el(2); - val[0] += s * src.el(3); - val[0] += s * src.el(4); - val[0] += s * src.el(5); - val[0] += s * src.el(6); - val[0] += s * src.el(7); - val[8] += s * src.el(8); - val[8] += s * src.el(9); - val[8] += s * src.el(10); - val[8] += s * src.el(11); - val[8] += s * src.el(12); - val[8] += s * src.el(13); - val[8] += s * src.el(14); - val[8] += s * src.el(15); - val[16] += s * src.el(16); - val[16] += s * src.el(17); - val[16] += s * src.el(18); - val[16] += s * src.el(19); - - val[16] += s * src.el(20); - val[16] += s * src.el(21); - val[16] += s * src.el(22); - val[16] += s * src.el(23); - val[24] += s * src.el(24); - val[24] += s * src.el(25); - val[24] += s * src.el(26); - val[24] += s * src.el(27); - val[24] += s * src.el(28); - val[24] += s * src.el(29); - - val[24] += s * src.el(30); - val[24] += s * src.el(31); - val[32] += s * src.el(32); - val[32] += s * src.el(33); - val[32] += s * src.el(34); - val[32] += s * src.el(35); - val[32] += s * src.el(36); - val[32] += s * src.el(37); - val[32] += s * src.el(38); - val[32] += s * src.el(39); - - val[40] += s * src.el(40); - val[40] += s * src.el(41); - val[40] += s * src.el(42); - val[40] += s * src.el(43); - val[40] += s * src.el(44); - val[40] += s * src.el(45); - val[40] += s * src.el(46); - val[40] += s * src.el(47); - val[48] += s * src.el(48); - val[48] += s * src.el(49); - - val[48] += s * src.el(50); - val[48] += s * src.el(51); - val[48] += s * src.el(52); - val[48] += s * src.el(53); - val[48] += s * src.el(54); - val[48] += s * src.el(55); - val[56] += s * src.el(56); - val[56] += s * src.el(57); - val[56] += s * src.el(58); - val[56] += s * src.el(59); - - val[56] += s * src.el(60); - val[56] += s * src.el(61); - val[56] += s * src.el(62); - val[56] += s * src.el(63); - } - else - { - const unsigned int size = n()*m(); - for (unsigned int i=0; i -void dFMatrix::Tadd (const double s, const dFMatrix& src) -{ - Assert (m() == n(), ExcNotQuadratic()); - Assert (m() == src.m(), ExcDimensionMismatch(m(), src.m())); - Assert (n() == src.n(), ExcDimensionMismatch(n(), src.n())); - - if ((n()==3) && (m()==3)) - { - val[0] += s * src.el(0); - val[1] += s * src.el(3); - val[2] += s * src.el(6); - - val[3] += s * src.el(1); - val[4] += s * src.el(4); - val[5] += s * src.el(7); - - val[6] += s * src.el(2); - val[7] += s * src.el(5); - val[8] += s * src.el(8); - } - else if ((n()==4) && (m()==4)) - { - val[0] += s * src.el(0); - val[1] += s * src.el(4); - val[2] += s * src.el(8); - val[3] += s * src.el(12); - - val[4] += s * src.el(1); - val[5] += s * src.el(5); - val[6] += s * src.el(9); - val[7] += s * src.el(13); - - val[8] += s * src.el(2); - val[9] += s * src.el(6); - val[10] += s * src.el(10); - val[11] += s * src.el(14); - - val[12] += s * src.el(3); - val[13] += s * src.el(7); - val[14] += s * src.el(11); - val[15] += s * src.el(15); - } - else if ((n()==8) && (m()==8)) - { - val[0] += s * src.el(0); - val[1] += s * src.el(8); - val[2] += s * src.el(16); - val[3] += s * src.el(24); - val[4] += s * src.el(32); - val[5] += s * src.el(40); - val[6] += s * src.el(48); - val[7] += s * src.el(56); - - val[8] += s * src.el(1); - val[9] += s * src.el(9); - val[10] += s * src.el(17); - val[11] += s * src.el(25); - val[12] += s * src.el(33); - val[13] += s * src.el(41); - val[14] += s * src.el(49); - val[15] += s * src.el(57); - - val[16] += s * src.el(2); - val[17] += s * src.el(10); - val[18] += s * src.el(18); - val[19] += s * src.el(26); - val[20] += s * src.el(34); - val[21] += s * src.el(42); - val[22] += s * src.el(50); - val[23] += s * src.el(58); - - val[24] += s * src.el(3); - val[25] += s * src.el(11); - val[26] += s * src.el(19); - val[27] += s * src.el(27); - val[28] += s * src.el(35); - val[29] += s * src.el(43); - val[30] += s * src.el(51); - val[31] += s * src.el(59); - - val[32] += s * src.el(4); - val[33] += s * src.el(12); - val[34] += s * src.el(20); - val[35] += s * src.el(28); - val[36] += s * src.el(36); - val[37] += s * src.el(44); - val[38] += s * src.el(52); - val[39] += s * src.el(60); - - val[40] += s * src.el(5); - val[41] += s * src.el(13); - val[42] += s * src.el(21); - val[43] += s * src.el(29); - val[44] += s * src.el(37); - val[45] += s * src.el(45); - val[46] += s * src.el(53); - val[47] += s * src.el(61); - - val[48] += s * src.el(6); - val[49] += s * src.el(14); - val[50] += s * src.el(22); - val[51] += s * src.el(30); - val[52] += s * src.el(38); - val[53] += s * src.el(46); - val[54] += s * src.el(54); - val[55] += s * src.el(62); - - val[56] += s * src.el(7); - val[57] += s * src.el(15); - val[58] += s * src.el(23); - val[59] += s * src.el(31); - val[60] += s * src.el(39); - val[61] += s * src.el(47); - val[62] += s * src.el(55); - val[63] += s * src.el(63); - } - else - Assert (false, ExcInternalError()); -} - - -template -bool -dFMatrix::operator == (const dFMatrix &m) const -{ - bool q = (dim_range==m.dim_range) && (dim_image==m.dim_image); - if (!q) return false; - - for (unsigned int i=0; i -double dFMatrix::determinant () const { - Assert (dim_range == dim_image, - ExcDimensionMismatch(dim_range, dim_image)); - Assert ((dim_range>=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; - }; -}; - -template -double dFMatrix::norm2 () const -{ - double s = 0.; - for (unsigned int i=0;i -void dFMatrix::clear () { - double *val_ptr = &val[0]; - const double *end_ptr = &val[n()*m()]; - while (val_ptr != end_ptr) - *val_ptr++ = 0.; -}; - - - -template -void dFMatrix::invert (const dFMatrix &M) { - Assert (dim_range == dim_image, ExcNotQuadratic()); - Assert ((dim_range>=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; - } - }; -}; - - - -template -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 -void -dFMatrix::householder(dVector& src) -{ - // m > 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 -double -dFMatrix::least_squares(dVector& dst, dVector& src) -{ - // m > n, m = src.n, n = dst.n - - householder(src); - backward(dst, src); - - double sum = 0.; - for (unsigned int i=n() ; i