From: Guido Kanschat Date: Thu, 29 Jun 2000 20:04:58 +0000 (+0000) Subject: FullMatrix::invert inverts all square matrices X-Git-Tag: v8.0.0~20317 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=947f34a00146ef51aa304f2c24dc6cef3aa98bee;p=dealii.git FullMatrix::invert inverts all square matrices git-svn-id: https://svn.dealii.org/trunk@3106 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/full_matrix.h b/deal.II/lac/include/lac/full_matrix.h index b7e9af2998..dd9813f903 100644 --- a/deal.II/lac/include/lac/full_matrix.h +++ b/deal.II/lac/include/lac/full_matrix.h @@ -365,11 +365,17 @@ class FullMatrix : public Subscriptor number norm2 () const; /** - * A=Inverse(A). Inversion of this by - * Gauss-Jordan-algorithm. Note that this - * is a rather expensive operation, so - * you may not want to use it for - * larger matrices if not necessary. + * A=Inverse(A). Inversion of + * this by Gauss-Jordan-algorithm + * with partial pivoting. This + * process is well-behaved for + * positive definite matrices, + * but be aware of round-off + * errors in the indefinite case. + * + * The numerical effort to invert + * an @p{n x n} matrix is of the + * order @p{n³}. */ void gauss_jordan (); @@ -384,28 +390,16 @@ class FullMatrix : public Subscriptor double determinant () const; /** - * Assign the inverse of the given - * matrix to @p{*this}. This function is - * only implemented (hardcoded) for - * square matrices of dimension one, - * two, three and four, since the - * amount of code needed grows quickly. - * The implementation does not use - * an elimination method like the - * Gauss-Jordan one, but rather sets - * the element directly; their values - * are precomputed symbolically using - * Maple. This way, we can avoid the - * overhead of loops and local variables - * but the number of lines of code - * grows rapidly. - * - * For all other sizes than the ones given - * above, an exception of type - * @p{ExcNotImplemented(dim_range)} is - * thrown, which you can catch and use - * some other method to invert the matrix, - * e.g. the @p{gauss_jordan} function. + * Assign the inverse of the + * given matrix to + * @p{*this}. This function is + * hardcoded for square matrices + * of dimension one to four, + * since the amount of code + * needed grows quickly. For + * larger matrices, the method + * @ref{gauss_jordan} is invoked + * implicitly. */ void invert (const FullMatrix &M); diff --git a/deal.II/lac/include/lac/full_matrix.templates.h b/deal.II/lac/include/lac/full_matrix.templates.h index 9a2f13c493..d423afa1c3 100644 --- a/deal.II/lac/include/lac/full_matrix.templates.h +++ b/deal.II/lac/include/lac/full_matrix.templates.h @@ -1276,7 +1276,8 @@ FullMatrix::invert (const FullMatrix &M) } default: - AssertThrow (false, ExcNotImplemented(dim_range)); + *this = M; + gauss_jordan(); }; }; @@ -1437,5 +1438,4 @@ FullMatrix::least_squares(Vector& dst, Vector& src) return sqrt(sum); } - #endif diff --git a/deal.II/lac/include/lac/precondition_block.templates.h b/deal.II/lac/include/lac/precondition_block.templates.h index b71ed8f3e1..22a3404084 100644 --- a/deal.II/lac/include/lac/precondition_block.templates.h +++ b/deal.II/lac/include/lac/precondition_block.templates.h @@ -15,6 +15,7 @@ #include +#include #include #include #include @@ -112,6 +113,7 @@ void PreconditionBlock::invert_diagblocks() if (same_diagonal) { + deallog << "PreconditionBlock uses only one diagonal block" << endl; // Invert only the first block // This is a copy of the code in the // 'else' part, stripped of the outer loop @@ -121,15 +123,7 @@ void PreconditionBlock::invert_diagblocks() for (unsigned int column_cell=0; column_cell::invert_diagblocks() for (unsigned int column_cell=0, column=cell*blocksize; column_cell