]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
add function for Schur complement; needs testing
authorkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 28 Jul 2011 17:48:13 +0000 (17:48 +0000)
committerkanschat <kanschat@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 28 Jul 2011 17:48:13 +0000 (17:48 +0000)
git-svn-id: https://svn.dealii.org/trunk@23981 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/lac/full_matrix.h
deal.II/include/deal.II/lac/full_matrix.templates.h

index 7d5031207de01fd532c3a8d241542c27c8d8178b..94f3f70e19b9cff43e6bd1dbbb9d32ed49bcc950 100644 (file)
@@ -1,7 +1,7 @@
 //---------------------------------------------------------------------------
 //    $Id$
 //
-//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 by the deal.II authors
+//    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -1156,7 +1156,23 @@ class FullMatrix : public Table<2,number>
     void TmTmult (FullMatrix<number2>       &C,
                  const FullMatrix<number2> &B,
                  const bool                 adding=false) const;
-
+    
+                                    /**
+                                     * Add to the current matrix the
+                                     * Schur complement <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);
+    
                                     /**
                                      * Matrix-vector-multiplication.
                                      *
index 7e8337f1c540029b1bec014976d48326650c875e..5b756750dd948b29f3c338fddd1b1cd0d6100ccc 100644 (file)
@@ -880,6 +880,67 @@ void FullMatrix<number>::TmTmult (FullMatrix<number2>       &dst,
 }
 
 
+template <typename number>
+void
+FullMatrix<number>::schur_complement(
+  const FullMatrix<number>& Ainverse,
+  const FullMatrix<number>& B,
+  const FullMatrix<number>& D,
+  const bool transpose_B,
+  const bool transpose_D)
+{  
+  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);
+    }
+  else
+    {
+      AssertDimension(B.n(), M);
+      AssertDimension(B.m(), N);
+    }
+  if (transpose_D)
+    {
+      AssertDimension(D.n(), M);
+      AssertDimension(D.m(), N);
+    }
+  else
+    {
+      AssertDimension(D.m(), M);
+      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)
+      {
+                                        // Compute the entry
+       number ADij = 0.;
+       if (transpose_D)
+         for (unsigned int k=0;k<M;++k)
+           ADij += Ainverse(i,k)*D(j,k);
+       else
+         for (unsigned int k=0;k<M;++k)
+           ADij += Ainverse(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);
+      }
+}
+    
 
 template <typename number>
 template <typename number2>

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.