From: Wolfgang Bangerth Date: Fri, 6 Mar 1998 15:04:07 +0000 (+0000) Subject: Add two functions which clear the matrix and compute its inverse. X-Git-Tag: v8.0.0~23217 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d63a4c0307e083c4f668f5fa17ccdedeabb4362a;p=dealii.git Add two functions which clear the matrix and compute its inverse. git-svn-id: https://svn.dealii.org/trunk@34 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/dfmatrix.h b/deal.II/lac/include/lac/dfmatrix.h index 04943233d0..0f355658d8 100644 --- a/deal.II/lac/include/lac/dfmatrix.h +++ b/deal.II/lac/include/lac/dfmatrix.h @@ -164,5 +164,20 @@ class dFMatrix * for this function. */ double determinant () const; + + /** + * Assign the inverse of the given + * matrix to #*this#. This function is + * only implemented (hardcoded) for + * square matrices of dimension one, + * two and three. + */ + void invert (const dFMatrix &M); + + /** + * Set all entries in the matrix to + * zero. + */ + void clear (); }; #endif diff --git a/deal.II/lac/source/dfmatrix.cc b/deal.II/lac/source/dfmatrix.cc index 9e8c807be7..6f583a99f0 100644 --- a/deal.II/lac/source/dfmatrix.cc +++ b/deal.II/lac/source/dfmatrix.cc @@ -21,7 +21,7 @@ void dFMatrix::init(int mm, int nn) val = new double[val_size]; dim_range = nn; dim_image = mm; - for (int i = nn*mm-1; i>=0 ; i--) val[i] = 0; + clear (); } dFMatrix::~dFMatrix() @@ -863,8 +863,8 @@ dFMatrix::operator == (const dFMatrix &m) const double dFMatrix::determinant () const { - THROW2 (dim_range != dim_image, ); - THROW2 ((dim_range<1) || (dim_range>3), IntError::Range(dim_range)); +// THROW2 (dim_range != dim_image, ); +// THROW2 ((dim_range<1) || (dim_range>3), IntError::Range(dim_range)); switch (dim_range) { @@ -883,3 +883,54 @@ double dFMatrix::determinant () const { return 0; }; }; + + + +void dFMatrix::clear () { + for (unsigned int i=0; i3), IntError::Range(dim_range)); +// THROW2 ((dim_range != M.dim_range) || (dim_image != M.dim_image), +// ExcXXX()); + 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; + }; +); +