From 0f4064cf72517566bd30101a88468b6cbd649d10 Mon Sep 17 00:00:00 2001 From: hartmann Date: Tue, 9 Feb 1999 13:03:06 +0000 Subject: [PATCH] implementation of householder and gauss-jordan, least-square git-svn-id: https://svn.dealii.org/trunk@773 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/dfmatrix.h | 20 ++- deal.II/lac/source/dfmatrix.cc | 218 ++++++++++++++++++++++++++++- 2 files changed, 235 insertions(+), 3 deletions(-) diff --git a/deal.II/lac/include/lac/dfmatrix.h b/deal.II/lac/include/lac/dfmatrix.h index 34e85950a1..1cbdff75d7 100644 --- a/deal.II/lac/include/lac/dfmatrix.h +++ b/deal.II/lac/include/lac/dfmatrix.h @@ -277,7 +277,7 @@ class dFMatrix double matrix_scalar_product (const dVector &u, const dVector &v) const; /** - * A=Inverse(A). Inversion of this by + * A=Inverse(A). Inversion of (*this) by * Gauss-Jordan-algorithm */ void gauss_jordan (); @@ -381,11 +381,23 @@ class dFMatrix * After execution of householder, the upper * triangle contains the resulting matrix R,

* the lower the incomplete factorization matrices. + * + * #householder(src); backward(dst, src);# gives + * the solution #dst# of the linear system + * #(*this)dst=src#. + * + * Note that #src# and #(*this)# (i.e. the + * matrix itself) is changed in + * the process of the #householder(src)# function!! */ - void householder (dVector& y); + void householder (dVector& src); /** * Least - Squares - Approximation by QR-factorization. + * + * Note that #src# and #(*this)# (i.e. the + * matrix itself) is changed in + * the process of this function!! */ double least_squares (dVector& dst, dVector& src); @@ -463,6 +475,10 @@ class dFMatrix /** * Exception */ + DeclException0 (ExcNotRegular); + /** + * Exception + */ DeclException0 (ExcInternalError); /** * Exception diff --git a/deal.II/lac/source/dfmatrix.cc b/deal.II/lac/source/dfmatrix.cc index fb00b87441..040141a95c 100644 --- a/deal.II/lac/source/dfmatrix.cc +++ b/deal.II/lac/source/dfmatrix.cc @@ -1047,7 +1047,7 @@ void dFMatrix::clear () { void dFMatrix::invert (const dFMatrix &M) { Assert (dim_range == dim_image, ExcNotQuadratic()); - Assert ((dim_range>=1) && (dim_range<=3), ExcNotImplemented(dim_range)); + 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, @@ -1091,6 +1091,105 @@ void dFMatrix::invert (const dFMatrix &M) { 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; + } }; }; @@ -1115,3 +1214,120 @@ void dFMatrix::print_formatted (ostream &out, const unsigned int precision) cons out.setf (0, ios::floatfield); // reset output format }; + + +// Gauss-Jordan-Algorithmus +// cf. Stoer I (4th Edition) p. 153 + +void dFMatrix::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 n, src.n() = m + Assert (dim_range <= dim_image, ExcDimensionMismatch(dim_range, dim_image)); + Assert (src.size() == dim_range, ExcDimensionMismatch(src.size(), dim_range)); + + for (unsigned int j=0 ; j n, m = src.n, n = dst.n + + householder(src); + backward(dst, src); + + double sum = 0.; + for (unsigned int i=n() ; i