From 046cb4a4c093716d8ac73bf372222aa998a92782 Mon Sep 17 00:00:00 2001 From: guido Date: Tue, 23 Feb 1999 14:45:11 +0000 Subject: [PATCH] FullMatrix template implemented git-svn-id: https://svn.dealii.org/trunk@880 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/full_matrix.h | 581 +++++++ .../lac/include/lac/full_matrix.templates.h | 1414 +++++++++++++++++ deal.II/lac/include/lac/fullmatrix.h~ | 577 +++++++ .../lac/include/lac/fullmatrix.templates.h~ | 1367 ++++++++++++++++ deal.II/lac/source/full_matrix.double.cc | 58 + deal.II/lac/source/full_matrix.float.cc | 58 + 6 files changed, 4055 insertions(+) create mode 100644 deal.II/lac/include/lac/full_matrix.h create mode 100644 deal.II/lac/include/lac/full_matrix.templates.h create mode 100644 deal.II/lac/include/lac/fullmatrix.h~ create mode 100644 deal.II/lac/include/lac/fullmatrix.templates.h~ create mode 100644 deal.II/lac/source/full_matrix.double.cc create mode 100644 deal.II/lac/source/full_matrix.float.cc diff --git a/deal.II/lac/include/lac/full_matrix.h b/deal.II/lac/include/lac/full_matrix.h new file mode 100644 index 0000000000..a857256418 --- /dev/null +++ b/deal.II/lac/include/lac/full_matrix.h @@ -0,0 +1,581 @@ +/*---------------------------- 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 + +template class Vector; + +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 FullMatrix (const unsigned int n = 1); + + /** + * Constructor. Initialize the matrix as + * a rectangular #m# times #n# matrix. + */ + FullMatrix (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 FullMatrix (const FullMatrix&); + + /** + * Destructor. Release all memory. + */ + ~FullMatrix(); + + /** + * 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 FullMatrix &) const; + + /** + * A = B . Copy all elements + */ + template + FullMatrix& operator = (const FullMatrix& B); + + + /** + * U(0-m,0-n) = s . Fill all elements + */ + template + void fill (const FullMatrix& 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 ) + */ + template + void reinit (const FullMatrix &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 FullMatrix& B); + + /** + * A+=Transp(B). + * Simple addition of the transpose of B to this + */ + template + void Tadd (const number s, const FullMatrix& B); + + /** + * C=A*B. + * Matrix-matrix-multiplication + */ + + template + void mmult (FullMatrix& C, const FullMatrix& B) const; + + /** + * C=Transp(A)*B. + * Matrix-matrix-multiplication using + * transpose of this + */ + template + void Tmmult (FullMatrix& C, const FullMatrix& 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 FullMatrix &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 + double 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 FullMatrix& 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 (ExcNotRegular); + /** + * 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 & +FullMatrix::el (const unsigned int i, const unsigned int j) +{ + return val[i*dim_range+j]; +}; + + +template +inline number +FullMatrix::el (const unsigned int i, const unsigned int j) const +{ + return val[i*dim_range+j]; +}; + + +template +inline unsigned int +FullMatrix::m() const +{ + return dim_image; +}; + + +template +inline unsigned int +FullMatrix::n() const +{ + return dim_range; +}; + + +template +inline number +FullMatrix::el (const unsigned int i) const +{ + return val[i]; +}; + + +template +inline number +FullMatrix::operator() (const unsigned int i, const unsigned int j) const +{ + Assert (i +inline number & +FullMatrix::operator() (const unsigned int i, const unsigned int j) +{ + Assert (i +#include +#include + +#include +#include +#include +#include + + +template +FullMatrix::FullMatrix (const unsigned int n) +{ + init (n,n); +}; + + +template +FullMatrix::FullMatrix (const unsigned int m, const unsigned int n) +{ + init (m,n); +}; + + +template +FullMatrix::FullMatrix (const FullMatrix &m) +{ + init (m.dim_image, m.dim_range); + number * p = &val[0]; + const number * vp = &m.val[0]; + const number * const e = &val[dim_image*dim_range]; + + while (p!=e) + *p++ = *vp++; +}; + + +template +void +FullMatrix::init (const unsigned int mm, const unsigned int nn) +{ + val_size = nn*mm; + val = new number[val_size]; + dim_range = nn; + dim_image = mm; + clear (); +}; + + +template +FullMatrix::~FullMatrix () +{ + delete[] val; +}; + + +template +bool +FullMatrix::all_zero () const +{ + const number *p = &val[0], + *e = &val[n()*m()]; + while (p!=e) + if (*p++ != 0.0) + return false; + + return true; +}; + + +template +void +FullMatrix::reinit (const unsigned int mm, const unsigned int nn) +{ + if (val_size +void +FullMatrix::reinit (const unsigned int n) +{ + reinit (n, n); +}; + + +template +template +void +FullMatrix::reinit (const FullMatrix &B) +{ + reinit (B.m(), B.n()); +}; + + +template +template +void +FullMatrix::vmult (Vector& dst, + const Vector& 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 + { + number* e = val; + const unsigned int size_m = m(), + size_n = n(); + for (unsigned int i=0; i +template +void FullMatrix::gsmult (Vector& dst, const Vector& 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) +template +void FullMatrix::Tvmult (Vector& dst, const Vector& 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 +template +double FullMatrix::residual (Vector& dst, const Vector& src, + const Vector& 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 +template +void FullMatrix::forward (Vector& dst, const Vector& 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() +template +void FullMatrix::backward (Vector& dst, const Vector& 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 +template +FullMatrix& +FullMatrix::operator = (const FullMatrix& m) +{ + reinit(m); + + number * p = &val[0]; + const number2 * vp = &m.val[0]; + const number * const e = &val[dim_image*dim_range]; + + while (p!=e) + *p++ = *vp++; + + return *this; +} + +template +template +void FullMatrix::fill (const FullMatrix& 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 FullMatrix::add_row (const unsigned int i, + const number s, const unsigned int j) +{ + for (unsigned int k=0; k +void FullMatrix::add_row (const unsigned int i, const number s, + const unsigned int j, const number t, + const unsigned int k) +{ + const unsigned int size_m = m(); + for (unsigned l=0; l +void FullMatrix::add_col (const unsigned int i, const number s, + const unsigned int j) +{ + for (unsigned int k=0; k +void FullMatrix::add_col (const unsigned int i, const number s, + const unsigned int j, const number t, + const unsigned int k) +{ + for (unsigned int l=0; l +void FullMatrix::swap_row (const unsigned int i, const unsigned int j) +{ + number s; + for (unsigned int k=0; k +void FullMatrix::swap_col (const unsigned int i, const unsigned int j) +{ + number s; + for (unsigned int k=0; k +void FullMatrix::diagadd (const number src) +{ + Assert (m() == n(), ExcDimensionMismatch(m(),n())); + for (unsigned int i=0; i +template +void FullMatrix::mmult (FullMatrix& dst, const FullMatrix& 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::mmult (FullMatrix& dst, const FullMatrix& src) const +{ + Assert (m() == src.n(), ExcDimensionMismatch(m(), src.n())); + + unsigned int i,j,k; + double s = 1.; + + dst.reinit(n(), src.m()); + + for (i=0;i +template +void FullMatrix::Tmmult (FullMatrix& dst, const FullMatrix& 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::Tmmult(FullMatrix& dst, const FullMatrix& src) const +{ + Assert (m() == src.n(), ExcDimensionMismatch(m(), src.n())); + + unsigned int i,j,k; + double s = 1.; + + dst.reinit(n(), src.m()); + + for (i=0;i +template +double FullMatrix::matrix_norm (const Vector &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 number *val_ptr = &val[0]; + const number2 *v_ptr; + + for (unsigned int row=0; row +template +double FullMatrix::matrix_scalar_product (const Vector &u, const Vector &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 number *val_ptr = &val[0]; + const number2 *v_ptr; + + for (unsigned int row=0; row +void +FullMatrix::print (ostream& s, int w, int p) const +{ + unsigned int i,j; + for (i=0;i +template +void +FullMatrix::add (const number s,const FullMatrix& 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 +template +void +FullMatrix::add_diag (const number s, const FullMatrix& 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 +template +void +FullMatrix::Tadd (const number s, const FullMatrix& 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, ExcNotImplemented(n())); +} + + +template +bool +FullMatrix::operator == (const FullMatrix &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 +FullMatrix::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 +FullMatrix::norm2 () const +{ + double s = 0.; + for (unsigned int i=0;i +void FullMatrix::clear () +{ + number *val_ptr = &val[0]; + const number *end_ptr = &val[n()*m()]; + while (val_ptr != end_ptr) + *val_ptr++ = 0.; +}; + + + +template +void +FullMatrix::invert (const FullMatrix &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 +FullMatrix::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 +void +FullMatrix::gauss_jordan() +{ + Assert (dim_range == dim_image, ExcNotQuadratic()); + iVector p(n()); + + unsigned int i,j,k,r; + double max, hr; + + for (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 hv(n()); + for (i=0; i +template +void +FullMatrix::householder(Vector& 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 +template +double +FullMatrix::least_squares(Vector& dst, Vector& src) +{ + // m > n, m = src.n, n = dst.n + + householder(src); + backward(dst, src); + + double sum = 0.; + for (unsigned int i=n() ; i + + +// 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 + +#define TYPEMAT double + +template class FullMatrix; + +#define TYPEMAT2 double + +//template FullMatrix& FullMatrix::operator =(const FullMatrix&); +template void FullMatrix::fill (const FullMatrix&, const unsigned, const unsigned); +template void FullMatrix::reinit (const FullMatrix&); +template void FullMatrix::add (const TYPEMAT, const FullMatrix&); +template void FullMatrix::Tadd (const TYPEMAT, const FullMatrix&); +template void FullMatrix::mmult (FullMatrix&, const FullMatrix&) const; +template void FullMatrix::Tmmult (FullMatrix&, const FullMatrix&) const; +template void FullMatrix::add_diag (const TYPEMAT, const FullMatrix&); + + +#define TYPEVEC double +#define TYPERES double + +template void FullMatrix::vmult(Vector&, const Vector&, const bool) const; +template void FullMatrix::Tvmult(Vector&, const Vector&, const bool) const; +template double FullMatrix::residual(Vector&, const Vector&, const Vector&) const; +template double FullMatrix::matrix_norm (const Vector &) const; +template double FullMatrix::matrix_scalar_product(const Vector&, const Vector&) const; +template void FullMatrix::forward(Vector&, const Vector&) const; +template void FullMatrix::backward(Vector&, const Vector&) const; +template void FullMatrix::householder(Vector&); +template double FullMatrix::least_squares(Vector&, Vector&); +template void FullMatrix::gsmult(Vector&, const Vector&, const iVector&) const; + +#undef TYPEVEC +#define TYPEVEC float + +template void FullMatrix::vmult(Vector&, const Vector&, const bool) const; +template void FullMatrix::Tvmult(Vector&, const Vector&, const bool) const; +template double FullMatrix::residual(Vector&, const Vector&, const Vector&) const; +template double FullMatrix::matrix_norm (const Vector &) const; +template double FullMatrix::matrix_scalar_product(const Vector&, const Vector&) const; +template void FullMatrix::forward(Vector&, const Vector&) const; +template void FullMatrix::backward(Vector&, const Vector&) const; +template void FullMatrix::householder(Vector&); +template double FullMatrix::least_squares(Vector&, Vector&); +template void FullMatrix::gsmult(Vector&, const Vector&, const iVector&) const; + diff --git a/deal.II/lac/source/full_matrix.float.cc b/deal.II/lac/source/full_matrix.float.cc new file mode 100644 index 0000000000..4126b44ed6 --- /dev/null +++ b/deal.II/lac/source/full_matrix.float.cc @@ -0,0 +1,58 @@ +// $Id$ + +// Driver for FullMatrix template instantiation. + +/* Instantiation is controlled by preprocessor symbols: + * + * 1. TYPEMAT : numerical type used in the matrix + * 2. TYPEVEC : numerical type for vector entries + * 3. TYPERES : numerical type for entries in the right hand side vector + * 4. TYPEMAT2: numerical type for the second matrix + */ + +#include + +#define TYPEMAT float + +template class FullMatrix; + +#define TYPEMAT2 float + +//template FullMatrix& FullMatrix::operator =(const FullMatrix&); +template void FullMatrix::fill (const FullMatrix&, const unsigned, const unsigned); +template void FullMatrix::reinit (const FullMatrix&); +template void FullMatrix::add (const TYPEMAT, const FullMatrix&); +template void FullMatrix::Tadd (const TYPEMAT, const FullMatrix&); +template void FullMatrix::mmult (FullMatrix&, const FullMatrix&) const; +template void FullMatrix::Tmmult (FullMatrix&, const FullMatrix&) const; +template void FullMatrix::add_diag (const TYPEMAT, const FullMatrix&); + + +#define TYPEVEC double +#define TYPERES double + +template void FullMatrix::vmult(Vector&, const Vector&, const bool) const; +template void FullMatrix::Tvmult(Vector&, const Vector&, const bool) const; +template double FullMatrix::residual(Vector&, const Vector&, const Vector&) const; +template double FullMatrix::matrix_norm (const Vector &) const; +template double FullMatrix::matrix_scalar_product(const Vector&, const Vector&) const; +template void FullMatrix::forward(Vector&, const Vector&) const; +template void FullMatrix::backward(Vector&, const Vector&) const; +template void FullMatrix::householder(Vector&); +template double FullMatrix::least_squares(Vector&, Vector&); +template void FullMatrix::gsmult(Vector&, const Vector&, const iVector&) const; + +#undef TYPEVEC +#define TYPEVEC float + +template void FullMatrix::vmult(Vector&, const Vector&, const bool) const; +template void FullMatrix::Tvmult(Vector&, const Vector&, const bool) const; +template double FullMatrix::residual(Vector&, const Vector&, const Vector&) const; +template double FullMatrix::matrix_norm (const Vector &) const; +template double FullMatrix::matrix_scalar_product(const Vector&, const Vector&) const; +template void FullMatrix::forward(Vector&, const Vector&) const; +template void FullMatrix::backward(Vector&, const Vector&) const; +template void FullMatrix::householder(Vector&); +template double FullMatrix::least_squares(Vector&, Vector&); +template void FullMatrix::gsmult(Vector&, const Vector&, const iVector&) const; + -- 2.39.5