new permuted SOR, must be tested

git-svn-id: https://svn.dealii.org/trunk@7446 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/lac/include/lac/precondition.h b/deal.II/lac/include/lac/precondition.h
index 62e6652285944f1de7fa9b9262a0b972f2f7f35e..2e2d88b86299bc31e873b3fbc1d437b08c27a874 100644 (file)
--- a/deal.II/lac/include/lac/precondition.h
+++ b/deal.II/lac/include/lac/precondition.h
@@ -355,6 +355,87 @@ class PreconditionSSOR : public PreconditionRelaxation<MATRIX>
 };
 
 
+/**
+ * Permuted SOR preconditioner using matrix built-in function.  The MATRIX
+ * class used is required to have functions
+ * @p{PSOR(VECTOR&, const VECTOR&, double)} and
+ * @p{TPSOR(VECTOR&, const VECTOR&, double)}.
+ *
+ *
+ * @sect2{Usage example}
+ * @begin{itemize}
+ *     // Declare related objects
+ *
+ * SparseMatrix<double> A;
+ * Vector<double> x;
+ * Vector<double> b;
+ * SolverCG<> solver(...);
+ *
+ * //...initialize and build A
+ *
+ * std::vector<unsigned int> permutation(x.size());
+ *
+ * //...fill permutation with reasonable values
+ *
+ *     // Define and initialize preconditioner
+ *
+ * PreconditionPSOR<SparseMatrix<double> > precondition;
+ * precondition.initialize (A, permutation, .6);
+ *
+ * solver.solve (A, x, b, precondition);
+ * @end{itemize}
+ *
+ * @author Guido Kanschat, 2003
+ */
+template <class MATRIX = SparseMatrix<double> >
+class PreconditionPSOR : public PreconditionRelaxation<MATRIX>
+{
+  public:
+                                    /**
+                                     * Initialize matrix and
+                                     * relaxation parameter. The
+                                     * matrix is just stored in the
+                                     * preconditioner object.
+                                     *
+                                     * The permutation vector is
+                                     * stored as a
+                                     * pointer. Therefore, it has to
+                                     * be assured that the lifetime
+                                     * of the vector exceeds the
+                                     * lifetime of the
+                                     * preconditioner.
+                                     *
+                                     * The relaxation parameter
+                                     * should be larger than zero and
+                                     * smaller than 2 for numerical
+                                     * reasons. It defaults to 1.
+                                     */
+    void initialize (const MATRIX& A,
+                    const std::vector<unsigned int>& permutation,
+                    typename PreconditionRelaxation<MATRIX>::AdditionalData
+                    parameters = AdditionalData());
+    
+                                    /**
+                                     * Apply preconditioner.
+                                     */
+    template<class VECTOR>
+    void vmult (VECTOR&, const VECTOR&) const;
+
+                                    /**
+                                     * Apply transpose
+                                     * preconditioner.
+                                     */
+    template<class VECTOR>
+    void Tvmult (VECTOR&, const VECTOR&) const;
+  private:
+                                    /**
+                                     * Storage for the permutation vector.
+                                     */
+    const std::vector<unsigned int>* permutation;
+};
+
+
+
 /**
  * Preconditioner using an iterative solver.  This preconditioner uses
  * a fully initialized LAC iterative solver for the approximate
@@ -624,6 +705,43 @@ PreconditionSSOR<MATRIX>::Tvmult (VECTOR& dst, const VECTOR& src) const
 }
 
 
+//----------------------------------------------------------------------//
+
+template <class MATRIX>
+inline void
+PreconditionPSOR<MATRIX>::initialize (
+  const MATRIX &rA,
+  const std::vector<unsigned int>& p,
+  typename PreconditionRelaxation<MATRIX>::AdditionalData parameters)
+{
+  permutation = &p;
+  PreconditionRelaxation<MATRIX>::initalize(rA, parameters);
+}
+
+
+template <class MATRIX>
+template<class VECTOR>
+inline void
+PreconditionPSOR<MATRIX>::vmult (VECTOR& dst, const VECTOR& src) const
+{
+  Assert (this->A!=0, ExcNotInitialized());
+  dst = src;
+  this->A->PSOR (dst, *permutation, this->relaxation);
+}
+
+
+
+template <class MATRIX>
+template<class VECTOR>
+inline void
+PreconditionPSOR<MATRIX>::Tvmult (VECTOR& dst, const VECTOR& src) const
+{
+  Assert (this->A!=0, ExcNotInitialized());
+  dst = src;
+  this->A->TPSOR (dst, *permutation, this->relaxation);
+}
+
+
 //----------------------------------------------------------------------//
 
 

diff --git a/deal.II/lac/include/lac/sparse_matrix.h b/deal.II/lac/include/lac/sparse_matrix.h
index d5d6680460739f1956bd22e4c82670f681d36837..e69fc38545b9000cd6d5abcfef805826049fa92a 100644 (file)
--- a/deal.II/lac/include/lac/sparse_matrix.h
+++ b/deal.II/lac/include/lac/sparse_matrix.h
@@ -38,7 +38,7 @@ template<typename number> class FullMatrix;
  * corresponding ``.templates.h'' file and instantiate the respective
  * class yourself.
  *
- * @author Original version by Roland Becker, Guido Kanschat, Franz-Theo Suttmeier; lots of enhancements, reorganisation and documentation by Wolfgang Bangerth 1998
+ * @author several, 1994-2003
  */
 template <typename number>
 class SparseMatrix : public Subscriptor
@@ -802,6 +802,46 @@ class SparseMatrix : public Subscriptor
                                      */
     template <typename somenumber>
     void TSOR (Vector<somenumber> &v,
+             const number        om = 1.) const;
+
+                                    /**
+                                     * Perform a permuted SOR
+                                     * preconditioning in-place.
+                                     *
+                                     * The standard SOR method is
+                                     * applied in the order
+                                     * prescribed by @p{permutation},
+                                     * that is, first the row
+                                     * @p{permutation[0]}, then
+                                     * @p{permutation[1]} and so
+                                     * on.
+                                     *
+                                     * @p{omega} is the relaxation
+                                     * parameter.
+                                     */
+    template <typename somenumber>
+    void PSOR (Vector<somenumber> &v,
+             const std::vector<unsigned int> permutation,
+             const number        om = 1.) const;
+
+                                    /**
+                                     * Perform a transposed permuted SOR
+                                     * preconditioning in-place.
+                                     *
+                                     * The transposed SOR method is
+                                     * applied in the order
+                                     * prescribed by @p{permutation},
+                                     * that is, first the row
+                                     * @p{permutation[m()-1]}, then
+                                     * @p{permutation[m()-2]} and so
+                                     * on.
+                                     *
+                                     * @p{omega} is the relaxation
+                                     * parameter.
+                                     */
+    template <typename somenumber>
+    void TPSOR (Vector<somenumber> &v,
+             const std::vector<unsigned int> permutation,
              const number        om = 1.) const;
 
                                     /**

diff --git a/deal.II/lac/include/lac/sparse_matrix.templates.h b/deal.II/lac/include/lac/sparse_matrix.templates.h
index 3464c7da3e418160be1d837a12b374572ae363bd..f5cf298da69d515c92800ed883c12503866eff1a 100644 (file)
--- a/deal.II/lac/include/lac/sparse_matrix.templates.h
+++ b/deal.II/lac/include/lac/sparse_matrix.templates.h
@@ -993,6 +993,69 @@ SparseMatrix<number>::TSOR (Vector<somenumber>& dst, const number om) const
 }
 
 
+template <typename number>
+template <typename somenumber>
+void
+SparseMatrix<number>::PSOR (
+  Vector<somenumber>& dst,
+  const std::vector<unsigned int> permutation,
+  const number om) const
+{
+  Assert (cols != 0, ExcMatrixNotInitialized());
+  Assert (val != 0, ExcMatrixNotInitialized());
+  Assert (m() == n(), ExcMatrixNotSquare());
+  Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
+  Assert (m() == permutation.size(),
+         ExcDimensionMismatch(m(), permutation.size()));
+
+  for (unsigned int urow=0; urow<m(); ++urow)
+    {
+      const unsigned int row = permutation[urow];
+      somenumber s = dst(row);
+      for (unsigned int j=cols->rowstart[row]; j<cols->rowstart[row+1]; ++j)
+       {
+         const unsigned int col = permutation[cols->colnums[j]];
+         if (col < row)
+           s -= val[j] * dst(col);
+       }
+
+      dst(row) = s * om / val[cols->rowstart[row]];
+    }
+}
+
+
+template <typename number>
+template <typename somenumber>
+void
+SparseMatrix<number>::TPSOR (
+  Vector<somenumber>& dst,
+  const std::vector<unsigned int> permutation,
+  const number om) const
+{
+  Assert (cols != 0, ExcMatrixNotInitialized());
+  Assert (val != 0, ExcMatrixNotInitialized());
+  Assert (m() == n(), ExcMatrixNotSquare());
+  Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
+  Assert (m() == permutation.size(),
+         ExcDimensionMismatch(m(), permutation.size()));
+
+  for (unsigned int urow=m(); urow != 0;)
+    {
+      --urow;
+      const unsigned int row = permutation[urow];
+      somenumber s = dst(row);
+      for (unsigned int j=cols->rowstart[row]; j<cols->rowstart[row+1]; ++j)
+       {
+         const unsigned int col = permutation[cols->colnums[j]];
+         if (col > row)
+           s -= val[j] * dst(col);
+       }
+
+      dst(row) = s * om / val[cols->rowstart[row]];
+    }
+}
+
+
 template <typename number>
 template <typename somenumber>
 void