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;
+
--- /dev/null
+// $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
--- /dev/null
+// $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);
+ }
+ }
+}
+
#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
#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
#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