]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Added cholesky and outer_product to FullMatrix class, various typo edits
authorlinhart <linhart@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 28 Jul 2009 15:51:18 +0000 (15:51 +0000)
committerlinhart <linhart@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 28 Jul 2009 15:51:18 +0000 (15:51 +0000)
git-svn-id: https://svn.dealii.org/trunk@19128 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/include/base/exceptions.h
deal.II/doc/news/changes.h
deal.II/lac/include/lac/full_matrix.h
deal.II/lac/include/lac/full_matrix.templates.h
deal.II/lac/source/full_matrix.inst.in

index 17fc4dfc54a559929660b4e8d395cd2ed7a72ed6..e03628133f9eeecf4c96cb24e280c35aca2a22a9 100644 (file)
@@ -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 <tt>char*</tt>
index 997bd38ffa59435b52d01173f03691109afc11ee..bc73458fadc48321ddad539cd183fb7ef21b6b98 100644 (file)
@@ -193,6 +193,28 @@ inconvenience this causes.
 <h3>lac</h3>
 
 <ol>
+  <li>
+  <p>
+  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 <tt>*this</tt> $= VW^T$ where $V$ 
+  and $W$ are vectors.
+  <br>
+  (Jean Marie Linhart 2009/07/27)
+  </p>
+  </li>
+
+  <li>
+  <p>
+  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.
+  <br>
+  (WB 2009/06/29)
+  </p>
+  </li>
+
   <li>
   <p>
   New: Based on work by Francisco Alvaro, the existing SLEPcWrappers now 
@@ -219,16 +241,6 @@ inconvenience this causes.
   </p> 
   </li>
 
-  <li>
-  <p>
-  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.
-  <br>
-  (WB 2009/06/29)
-  </p>
-  </li>
-
   <li>
   <p>
   Fixed: The TrilinosWrappers::BlockVector class declares an assignment
index ee6df736c7305390812ea030d4c122c117bb23fc..89e89e07d0ade663185777d45b247f23ce7c5069 100644 (file)
@@ -772,7 +772,7 @@ class FullMatrix : public Table<2,number>
               const FullMatrix<number2> &B);
 
                                     /**
-                                     * Add transose of a rectangular block.
+                                     * Add transpose of a rectangular block.
                                      *
                                      * A rectangular block of the
                                      * matrix <tt>src</tt> is
@@ -940,6 +940,27 @@ class FullMatrix : public Table<2,number>
     template <typename number2>
     void invert (const FullMatrix<number2> &M);
 
+                                    /**
+                                     * Assign the Cholesky decomposition
+                                     * of the given matrix to <tt>*this</tt>. 
+                                     * The given matrix must be symmetric
+                                     * positive definite.
+                                     *
+                                     * ExcMatrixNotPositiveDefinite
+                                     * will be thrown in the case that the
+                                     * matrix is not positive definite.
+                                     */
+    template <typename number2>
+    void cholesky (const FullMatrix<number2> &A);
+
+                                    /**
+                                     * <tt>*this</tt> = $V W^T$ where $V,W$
+                                     * are vectors of the same length.
+                                     */
+    template <typename number2>
+    void outer_product (const Vector<number2> &V,
+                       const Vector<number2> &W);
+    
                                     /**
                                      * Assign the left_inverse of the given matrix
                                      * to <tt>*this</tt>. The calculation being 
@@ -1234,6 +1255,10 @@ class FullMatrix : public Table<2,number>
                                       * Exception
                                       */
     DeclException0 (ExcSourceEqualsDestination);
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcMatrixNotPositiveDefinite);
                                     //@}
     
     friend class Accessor;
index f56acaf1c44a502c87cf4fa5fbb0ed45134b0b8b..1c32461da36c426ce71618d4292212a10fcc4dc1 100644 (file)
@@ -1281,6 +1281,7 @@ FullMatrix<number>::invert (const FullMatrix<number2> &M)
        
        break;
       }
+      
 
       default:
                                             // if no inversion is
@@ -1292,6 +1293,70 @@ FullMatrix<number>::invert (const FullMatrix<number2> &M)
     };    
 }
 
+
+template <typename number>
+template <typename number2>
+void
+FullMatrix<number>::cholesky (const FullMatrix<number2> &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<number2> 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<j; k++)
+           {
+             SLikLjk = (*this)(i,k)*(*this)(j,k);
+           };
+         (*this)(i,j) = (1./(*this)(j,j))*(A(i,j) - SLikLjk);
+       }
+       AssertThrow (A(i,i) - SLik2 >= 0,
+                    ExcMatrixNotPositiveDefinite());
+       
+       (*this)(i,i) = std::sqrt(A(i,i) - SLik2);
+      }
+    }
+  }
+
+
+template <typename number>
+template <typename number2>
+void
+FullMatrix<number>::outer_product (const Vector<number2> &V,
+                                  const Vector<number2> &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; i<this->n(); i++)
+    {
+      for(unsigned int j = 0; j< this->n(); j++)
+       {
+         (*this)(i,j) = V(i)*W(j);
+       }
+    }
+}
+
+
 template <typename number>
 template <typename number2>
 void
index 4b9c727460a27309a6a2b9077e47ab687e6e434f..2d7537b2abdb03277d5ff0a8a862a04732ea4dce 100644 (file)
@@ -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<S>::copy_from<1>(
       Tensor<2,1>&, const unsigned, const unsigned, const unsigned, const unsigned, const unsigned, const unsigned);
-                                               
+
     template void FullMatrix<S>::copy_from<2>(
       Tensor<2,2>&, const unsigned, const unsigned, const unsigned, const unsigned, const unsigned, const unsigned);
-                                               
+
     template void FullMatrix<S>::copy_from<3>(
       Tensor<2,3>&, const unsigned, const unsigned, const unsigned, const unsigned, const unsigned, const unsigned);
-                                               
+
     template void FullMatrix<S>::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<S1>::precondition_Jacobi<S2> (
        Vector<S2> &, const Vector<S2> &, const S1) const;
+
+    template
+      void FullMatrix<S1>::cholesky<S2> (const FullMatrix<S2>&);
+
+    template
+      void FullMatrix<S1>::outer_product<S2> (const Vector<S2>&,
+                                              const Vector<S2>&);
   }
 
 for (S1, S2, S3 : REAL_SCALARS)

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.