]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
improve and rename triple_product
authorkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 5 Aug 2011 14:46:28 +0000 (14:46 +0000)
committerkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 5 Aug 2011 14:46:28 +0000 (14:46 +0000)
git-svn-id: https://svn.dealii.org/trunk@24018 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/include/deal.II/lac/full_matrix.h
deal.II/include/deal.II/lac/full_matrix.templates.h

index 73f578c4f2dbc60f5ce971d3947b377339d44715..f27eaee2dc60e4c5d5ea815e47f70c537294ef39 100644 (file)
@@ -240,7 +240,8 @@ should be fixed now.
 <br>
 (Wolfgang Bangerth, 2011/01/18)
 
-<li> Extended: Several missing instantiations of functions for triangulations and DoF handlers embedded in higher dimensional space have been added.
+<li> Extended: Several missing instantiations of functions for triangulations
+and DoF handlers embedded in higher dimensional space have been added.
 <br>
 (Wolfgang Bangerth, 2011/01/15)
 </ol>
@@ -251,6 +252,12 @@ should be fixed now.
 <h3>Specific improvements</h3>
 
 <ol>
+
+<li> New: The function FullMatrix::triple_product() adds triple products
+like Schur complements to existing matrices.
+<br>
+(Guido Kanschat, 2011/08/05)
+
 <li> Improved: The PETScWrapper::VectorBase class was rather generous in
 calling the PETSc <code>VecAssembleBegin/End</code> functions that incur
 communication in the %parallel case and are therefore causes of potential
index 1f61a511a8f8eb1d753e40675af4256b28700602..74477a446b3f41f128152823c3ca7619d1482e5e 100644 (file)
@@ -1159,19 +1159,31 @@ class FullMatrix : public Table<2,number>
     
                                     /**
                                      * Add to the current matrix the
-                                     * Schur complement <b>B
+                                     * triple product <b>B
                                      * A<sup>-1</sup>
                                      * D</b>. Optionally, use the
                                      * transposes of the matrices
-                                     * <b>B</b> and <b>D</b>. Note
-                                     * that the argument for <b>A</b>
-                                     * is already the inverse.
-                                     */
-    void schur_complement(const FullMatrix<number>& Ainverse,
-                         const FullMatrix<number>& B,
-                         const FullMatrix<number>& D,
-                         const bool transpose_B = false,
-                         const bool transpose_D = false);
+                                     * <b>B</b> and <b>D</b>. The
+                                     * scaling factor scales the
+                                     * whole product, which is
+                                     * helpful when adding a multiple
+                                     * of the triple product to the
+                                     * matrix.
+                                     *
+                                     * This product was written with
+                                     * the Schur complement
+                                     * <b>B<sup>T</sup>
+                                     * A<sup>-1</sup> D</b> in mind.
+                                     * Note that in this case the
+                                     * argument for <tt>A<tt> must be
+                                     * the inverse of the matrix <b>A</b>.
+                                     */
+    void triple_product(const FullMatrix<number>& A,
+                       const FullMatrix<number>& B,
+                       const FullMatrix<number>& D,
+                       const bool transpose_B = false,
+                       const bool transpose_D = false,
+                       const number scaling = number(1.));
     
                                     /**
                                      * Matrix-vector-multiplication.
index 5b756750dd948b29f3c338fddd1b1cd0d6100ccc..a912f598114ebe0232a1762e8319788039b77025 100644 (file)
@@ -882,62 +882,57 @@ void FullMatrix<number>::TmTmult (FullMatrix<number2>       &dst,
 
 template <typename number>
 void
-FullMatrix<number>::schur_complement(
-  const FullMatrix<number>& Ainverse,
+FullMatrix<number>::triple_product(
+  const FullMatrix<number>& A,
   const FullMatrix<number>& B,
   const FullMatrix<number>& D,
   const bool transpose_B,
-  const bool transpose_D)
+  const bool transpose_D,
+  const number scaling)
 {  
-  AssertDimension (m(), n());
-  AssertDimension (Ainverse.m(), Ainverse.n());
-
-  const unsigned int N = n();
-  const unsigned int M = Ainverse.n();
-  
   if (transpose_B)
     {
-      AssertDimension(B.m(), M);
-      AssertDimension(B.n(), N);
+      AssertDimension(B.m(), A.m());
+      AssertDimension(B.n(), m());
     }
   else
     {
-      AssertDimension(B.n(), M);
-      AssertDimension(B.m(), N);
+      AssertDimension(B.n(), A.m());
+      AssertDimension(B.m(), m());
     }
   if (transpose_D)
     {
-      AssertDimension(D.n(), M);
-      AssertDimension(D.m(), N);
+      AssertDimension(D.n(), A.n());
+      AssertDimension(D.m(), n());
     }
   else
     {
-      AssertDimension(D.m(), M);
-      AssertDimension(D.n(), N);
+      AssertDimension(D.m(), A.n());
+      AssertDimension(D.n(), n());
     }
 
                                   // For all entries of the product
                                   // AD
-  for (unsigned int i=0; i<M;++i)
-    for (unsigned int j=0; j<N;++j)
+  for (unsigned int i=0; i<A.m();++i)
+    for (unsigned int j=0; j<n();++j)
       {
                                         // Compute the entry
        number ADij = 0.;
        if (transpose_D)
-         for (unsigned int k=0;k<M;++k)
-           ADij += Ainverse(i,k)*D(j,k);
+         for (unsigned int k=0;k<A.n();++k)
+           ADij += A(i,k)*D(j,k);
        else
-         for (unsigned int k=0;k<M;++k)
-           ADij += Ainverse(i,k)*D(k,j);
+         for (unsigned int k=0;k<A.n();++k)
+           ADij += A(i,k)*D(k,j);
                                   // And add it to this after
                                   // multiplying with the right
                                   // factor from B
-  if (transpose_B)
-    for (unsigned int k=0;k<N;++k)
-      this->operator()(k,j) += ADij*B(i,k);
-  else
-    for (unsigned int k=0;k<N;++k)
-      this->operator()(k,j) += ADij*B(k,i);
+       if (transpose_B)
+         for (unsigned int k=0;k<m();++k)
+           this->operator()(k,j) += scaling * ADij * B(i,k);
+       else
+         for (unsigned int k=0;k<m();++k)
+           this->operator()(k,j) += scaling * ADij * B(k,i);
       }
 }
     

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.