]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Vanka compiles
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 2 Jul 1999 12:36:14 +0000 (12:36 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 2 Jul 1999 12:36:14 +0000 (12:36 +0000)
git-svn-id: https://svn.dealii.org/trunk@1530 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/sparse_matrix.2.templates
deal.II/lac/include/lac/sparsevanka.h [new file with mode: 0644]
deal.II/lac/include/lac/sparsevanka.templates.h [new file with mode: 0644]
deal.II/lac/source/sparse_matrix.double.cc
deal.II/lac/source/sparse_matrix.float.cc
deal.II/lac/source/sparse_matrix.long_double.cc

index 9dbccff468481268320d974db8741a1ca192248e..177b86395dce1f44eda81d88fa1b0b2530b6de0d 100644 (file)
@@ -42,3 +42,6 @@ template void SparseMatrix<TYPEMAT>::SSOR (Vector<TYPE2> &, TYPEMAT) const;
 template void SparseMatrix<TYPEMAT>::SOR_step (Vector<TYPE2> &, const Vector<TYPE2> &, TYPEMAT) const;
 template void SparseMatrix<TYPEMAT>::TSOR_step (Vector<TYPE2> &, const Vector<TYPE2> &, TYPEMAT) const;
 template void SparseMatrix<TYPEMAT>::SSOR_step (Vector<TYPE2> &, const Vector<TYPE2> &, TYPEMAT) const;
+
+template void SparseVanka<TYPEMAT>::apply(Vector<TYPE2>& dst) const;
+
diff --git a/deal.II/lac/include/lac/sparsevanka.h b/deal.II/lac/include/lac/sparsevanka.h
new file mode 100644 (file)
index 0000000..0a06642
--- /dev/null
@@ -0,0 +1,89 @@
+// $Id$
+// Copyright Guido Kanschat, 1999
+
+#ifndef __lac_sparsematrix_H
+#define __lac_sparsematrix_H
+
+#include <base/smartpointer.h>
+#include <lac/forward-declarations.h>
+
+#include <vector>
+
+/**
+ * Point-wise Vanka preconditioning.
+ * This class does Vanka preconditioning  on a point-wise base.
+ * Vanka preconditioners are used for saddle point problems. There the
+ * application of Jacobi or Gauß-Seidel methods is impossible, because
+ * the diagonal elements are zero in the rows of the Lagrange multiplier.
+ *
+ * It is constructed initializing a vector of indices to the degrees of
+ * freedom of the Lagrange multiplier.
+ *
+ * In the actual preconditioning method, these rows are traversed in
+ * original order. Since this is a Gauß-Seidel like procedure,
+ * remember to have a good ordering in advance.
+ *
+ * For each row, a local system of equations is built by the degree of
+ * freedom itself and all other values coupling immediately. The right
+ * hand side is augmented by all further couplings.
+ *
+ * This local system is solved and the values are updated into the
+ * destination vector.
+ * @author Guido Kanschat
+ */
+template<typename number>
+class SparseVanka
+{
+  public:
+                                    /**
+                                     * Constructor.
+                                     * Take a vector of the indices
+                                     * of the Lagrange multiplier as
+                                     * argument. A reference to this
+                                     * vector will be stored, so it
+                                     * must persist longer than the
+                                     * Vanka object. The same is true
+                                     * for the matrix.
+                                     */
+    SparseVanka(const SparseMatrix<number>& M,
+               const vector<int>& indices);
+                                    /**
+                                     * Do the preconditioning.
+                                     */
+    template<typename number2>
+    void operator() (Vector<number2>& dst,
+                    const Vector<number2>& src) const;
+
+                                    /**
+                                     * In-place application of the
+                                     * method.
+                                     */
+    template<typename number2>
+    void apply(Vector<number2>& dst) const;
+    
+  private:
+                                    /**
+                                     * Pointer to the matrix.
+                                     */
+    SmartPointer<const SparseMatrix<number> > matrix;
+    
+                                    /**
+                                     * Indices of Lagrange
+                                     * multipliers.
+                                     */
+    const vector<int>& indices;
+};
+
+template<typename number>
+template<typename number2>
+inline
+void
+SparseVanka<number>::operator() (Vector<number2>& dst,
+                                const Vector<number2>& src) const
+{
+  dst = src;
+  apply(dst);
+}
+
+
+#endif
diff --git a/deal.II/lac/include/lac/sparsevanka.templates.h b/deal.II/lac/include/lac/sparsevanka.templates.h
new file mode 100644 (file)
index 0000000..19e87b7
--- /dev/null
@@ -0,0 +1,77 @@
+// $Id$
+// Copyright Guido Kanschat, 1999
+
+#include <lac/sparsevanka.h>
+#include <lac/fullmatrix.h>
+
+#include <map>
+
+template<typename number>
+SparseVanka<number>::SparseVanka(const SparseMatrix<number>& M,
+                                const vector<int>& indices)
+               :
+               matrix(&M), indices(indices)
+{}
+
+template<typename number>
+template<typename number2>
+void
+SparseVanka<number>::apply(Vector<number2>& dst) const
+{
+  for (unsigned int global_i=0; global_i<indices.size() ; ++global_i)
+    {
+      unsigned int row = indices[global_i];
+      const SparseMatrixStruct& structure = matrix->get_sparsity_pattern();
+      unsigned int n = structure.row_length(row);
+      
+      FullMatrix<number> A(n);
+      Vector<number> b(n);
+      Vector<number> x(n);
+      
+      map<unsigned int, unsigned int> local_index;
+
+                                      // Build local index
+
+      for (unsigned int i=0;i<n;++i)
+       local_index.insert(pair<unsigned int, unsigned int>
+                          (structure.column_number(row, i), i));
+
+                                      // Build local matrix
+
+      for (map<unsigned int, unsigned int>::iterator is=local_index.begin();
+          is!=local_index.end();++is)
+       {
+         unsigned int irow = is->first;
+         unsigned int i = is->second;
+         unsigned int n = structure.row_length(irow);
+         
+         b(i) = dst(irow);
+         
+         for (unsigned int j=0;j<n;++j)
+           {
+             unsigned int col = structure.column_number(irow, j);
+             map<unsigned int, unsigned int>::iterator js
+               = local_index.find(col);
+             if (js == local_index.end())
+               {
+                 b(i) -= matrix->raw_entry(irow,col) * dst(col);
+               } else {
+                 A(i,j) = matrix->raw_entry(irow,col);
+               }
+           }
+       }
+                                      // Compute new values
+      A.gauss_jordan();
+      A.vmult(x,b);
+      
+                                      // Distribute new values
+      for (map<unsigned int, unsigned int>::iterator is=local_index.begin();
+          is!=local_index.end();++is)
+       {
+         unsigned int irow = is->first;
+         unsigned int i = is->second;
+         dst(irow) = x(i);
+       }     
+    }
+}
+
index d129b74e99d8fb0c7ff1e2fc6a886fcecd4743b2..d677b3275fe48bc36442e086abfe4f23b0fc293b 100644 (file)
 
 #include <cmath>
 #include <lac/sparsematrix.templates.h>
+#include <lac/sparsevanka.templates.h>
 
 #define TYPEMAT double
 
 template class SparseMatrix<TYPEMAT>;
+template class SparseVanka<TYPEMAT>;
 
 #define TYPE2 float
 
index f7d2f03495e6493c13a1df1b2aa52bed7cbd6b0e..32b2e0324e50d96bd413bc5618b035891911ef90 100644 (file)
 
 #include <cmath>
 #include <lac/sparsematrix.templates.h>
+#include <lac/sparsevanka.templates.h>
 
 #define TYPEMAT float
 
 template class SparseMatrix<TYPEMAT>;
+template class SparseVanka<TYPEMAT>;
 
 #define TYPE2 float
 
index dbd0e62aa0b63f893caa890d787cefe8849d287a..9afdba45309ae5a6d9334d1146488b7268734938 100644 (file)
 
 #include <cmath>
 #include <lac/sparsematrix.templates.h>
+#include <lac/sparsevanka.templates.h>
 
 #define TYPEMAT long double
 
 template class SparseMatrix<TYPEMAT>;
+template class SparseVanka<TYPEMAT>;
 
 #define TYPE2 float
 

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.