From c139075c4981488b08f729837f0b90f80db0ce38 Mon Sep 17 00:00:00 2001 From: linhart Date: Tue, 28 Jul 2009 15:51:18 +0000 Subject: [PATCH] Added cholesky and outer_product to FullMatrix class, various typo edits git-svn-id: https://svn.dealii.org/trunk@19128 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/exceptions.h | 8 +-- deal.II/doc/news/changes.h | 32 ++++++--- deal.II/lac/include/lac/full_matrix.h | 27 +++++++- .../lac/include/lac/full_matrix.templates.h | 65 +++++++++++++++++++ deal.II/lac/source/full_matrix.inst.in | 15 +++-- 5 files changed, 128 insertions(+), 19 deletions(-) diff --git a/deal.II/base/include/base/exceptions.h b/deal.II/base/include/base/exceptions.h index 17fc4dfc54..e03628133f 100644 --- a/deal.II/base/include/base/exceptions.h +++ b/deal.II/base/include/base/exceptions.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -590,20 +590,20 @@ namespace StandardExceptions * where all allocated objects * should have been released. Since * this exception is thrown, some - * where still allocated. + * were still allocated. */ DeclException1 (ExcMemoryLeak, int, << "Destroying memory handler while " << arg1 << " objects are still allocated"); /** - * An error occured reading or + * An error occurred reading or * writing a file. */ DeclException0 (ExcIO); /** - * An error occured opening the named file. + * An error occurred opening the named file. * * The constructor takes a single * argument of type char* diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 997bd38ffa..bc73458fad 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -193,6 +193,28 @@ inconvenience this causes.

lac

    +
  1. +

    + New: There are new functions FullMatrix::cholesky and + FullMatrix::outer_product. FullMatrix::cholesky finds the Cholesky + decomposition of a matrix in lower triangular form. + FullMatrix::outer_product calculates *this $= VW^T$ where $V$ + and $W$ are vectors. +
    + (Jean Marie Linhart 2009/07/27) +

    +
  2. + +
  3. +

    + Fixed: The TrilinosWrappers::MPI::BlockVector class declares an assignment + operator from the non-Trilinos BlockVector class but it could not be + compiled due to an oversight. This is now fixed. +
    + (WB 2009/06/29) +

    +
  4. +
  5. New: Based on work by Francisco Alvaro, the existing SLEPcWrappers now @@ -219,16 +241,6 @@ inconvenience this causes.

  6. -
  7. -

    - Fixed: The TrilinosWrappers::MPI::BlockVector class declares an assignment - operator from the non-Trilinos BlockVector class but it could not be - compiled due to an oversight. This is now fixed. -
    - (WB 2009/06/29) -

    -
  8. -
  9. Fixed: The TrilinosWrappers::BlockVector class declares an assignment diff --git a/deal.II/lac/include/lac/full_matrix.h b/deal.II/lac/include/lac/full_matrix.h index ee6df736c7..89e89e07d0 100644 --- a/deal.II/lac/include/lac/full_matrix.h +++ b/deal.II/lac/include/lac/full_matrix.h @@ -772,7 +772,7 @@ class FullMatrix : public Table<2,number> const FullMatrix &B); /** - * Add transose of a rectangular block. + * Add transpose of a rectangular block. * * A rectangular block of the * matrix src is @@ -940,6 +940,27 @@ class FullMatrix : public Table<2,number> template void invert (const FullMatrix &M); + /** + * Assign the Cholesky decomposition + * of the given matrix to *this. + * The given matrix must be symmetric + * positive definite. + * + * ExcMatrixNotPositiveDefinite + * will be thrown in the case that the + * matrix is not positive definite. + */ + template + void cholesky (const FullMatrix &A); + + /** + * *this = $V W^T$ where $V,W$ + * are vectors of the same length. + */ + template + void outer_product (const Vector &V, + const Vector &W); + /** * Assign the left_inverse of the given matrix * to *this. The calculation being @@ -1234,6 +1255,10 @@ class FullMatrix : public Table<2,number> * Exception */ DeclException0 (ExcSourceEqualsDestination); + /** + * Exception + */ + DeclException0 (ExcMatrixNotPositiveDefinite); //@} friend class Accessor; diff --git a/deal.II/lac/include/lac/full_matrix.templates.h b/deal.II/lac/include/lac/full_matrix.templates.h index f56acaf1c4..1c32461da3 100644 --- a/deal.II/lac/include/lac/full_matrix.templates.h +++ b/deal.II/lac/include/lac/full_matrix.templates.h @@ -1281,6 +1281,7 @@ FullMatrix::invert (const FullMatrix &M) break; } + default: // if no inversion is @@ -1292,6 +1293,70 @@ FullMatrix::invert (const FullMatrix &M) }; } + +template +template +void +FullMatrix::cholesky (const FullMatrix &A) +{ + Assert (!A.empty(), ExcEmptyMatrix()); + Assert (A.n() == A.m(), + ExcNotQuadratic()); + // Matrix must be symmetric. + Assert(A.relative_symmetry_norm2() < 1.0e-10, ExcMessage("A must be symmetric.")); + + if (PointerComparison::equal(&A, this)) + { + // avoid overwriting source + // by destination matrix: + FullMatrix A2 = A; + cholesky(A2); + } + else + { + /* reinit *this to 0 */ + this->reinit(A.m(), A.n()); + + double SLik2 = 0.0, SLikLjk = 0.0; + for (unsigned int i=0; i< this->n_cols(); i++){ + SLik2 = 0.0; + for (unsigned int j = 0; j < i; j++){ + SLik2 += (*this)(i,j)*(*this)(i,j); + SLikLjk = 0.0; + for (unsigned int k =0; k= 0, + ExcMatrixNotPositiveDefinite()); + + (*this)(i,i) = std::sqrt(A(i,i) - SLik2); + } + } + } + + +template +template +void +FullMatrix::outer_product (const Vector &V, + const Vector &W) +{ + Assert (V.size() == W.size(), ExcMessage("Vectors V, W must be the same size.")); + this->reinit(V.size(), V.size()); + + for(unsigned int i = 0; in(); i++) + { + for(unsigned int j = 0; j< this->n(); j++) + { + (*this)(i,j) = V(i)*W(j); + } + } +} + + template template void diff --git a/deal.II/lac/source/full_matrix.inst.in b/deal.II/lac/source/full_matrix.inst.in index 4b9c727460..2d7537b2ab 100644 --- a/deal.II/lac/source/full_matrix.inst.in +++ b/deal.II/lac/source/full_matrix.inst.in @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors +// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -23,13 +23,13 @@ for (S : REAL_SCALARS) template void FullMatrix::copy_from<1>( Tensor<2,1>&, const unsigned, const unsigned, const unsigned, const unsigned, const unsigned, const unsigned); - + template void FullMatrix::copy_from<2>( Tensor<2,2>&, const unsigned, const unsigned, const unsigned, const unsigned, const unsigned, const unsigned); - + template void FullMatrix::copy_from<3>( Tensor<2,3>&, const unsigned, const unsigned, const unsigned, const unsigned, const unsigned, const unsigned); - + template void FullMatrix::copy_to<1>( Tensor<2,1>&, const unsigned, const unsigned, const unsigned, const unsigned, const unsigned, const unsigned); @@ -116,6 +116,13 @@ for (S1, S2 : REAL_SCALARS) template void FullMatrix::precondition_Jacobi ( Vector &, const Vector &, const S1) const; + + template + void FullMatrix::cholesky (const FullMatrix&); + + template + void FullMatrix::outer_product (const Vector&, + const Vector&); } for (S1, S2, S3 : REAL_SCALARS) -- 2.39.5