]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
PSOR
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 25 Apr 2003 12:52:33 +0000 (12:52 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 25 Apr 2003 12:52:33 +0000 (12:52 +0000)
git-svn-id: https://svn.dealii.org/trunk@7478 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/precondition.h
deal.II/lac/include/lac/sparse_matrix.2.templates
deal.II/lac/include/lac/sparse_matrix.h
deal.II/lac/include/lac/sparse_matrix.templates.h
tests/lac/solver.cc
tests/results/i686-pc-linux-gnu+gcc3.2/lac/solver.output

index b30f2d119f0b78915eca61ee3f5a53f9d2843f5d..7259f27fce28e484316b967cb404f6e3c962efd8 100644 (file)
@@ -374,13 +374,14 @@ class PreconditionSSOR : public PreconditionRelaxation<MATRIX>
  * //...initialize and build A
  *
  * std::vector<unsigned int> permutation(x.size());
+ * std::vector<unsigned int> inverse_permutation(x.size());
  *
- * //...fill permutation with reasonable values
+ * //...fill permutation and its inverse with reasonable values
  *
  *     // Define and initialize preconditioner
  *
  * PreconditionPSOR<SparseMatrix<double> > precondition;
- * precondition.initialize (A, permutation, .6);
+ * precondition.initialize (A, permutation, inverse_permutation, .6);
  *
  * solver.solve (A, x, b, precondition);
  * @end{itemize}
@@ -412,6 +413,7 @@ class PreconditionPSOR : public PreconditionRelaxation<MATRIX>
                                      */
     void initialize (const MATRIX& A,
                     const std::vector<unsigned int>& permutation,
+                    const std::vector<unsigned int>& inverse_permutation,
                     typename PreconditionRelaxation<MATRIX>::AdditionalData
                     parameters = typename PreconditionRelaxation<MATRIX>::AdditionalData());
     
@@ -432,6 +434,11 @@ class PreconditionPSOR : public PreconditionRelaxation<MATRIX>
                                      * Storage for the permutation vector.
                                      */
     const std::vector<unsigned int>* permutation;
+                                    /**
+                                     * Storage for the inverse
+                                     * permutation vector.
+                                     */
+    const std::vector<unsigned int>* inverse_permutation;
 };
 
 
@@ -712,9 +719,11 @@ inline void
 PreconditionPSOR<MATRIX>::initialize (
   const MATRIX &rA,
   const std::vector<unsigned int>& p,
+  const std::vector<unsigned int>& ip,
   typename PreconditionRelaxation<MATRIX>::AdditionalData parameters)
 {
   permutation = &p;
+  inverse_permutation = &ip;
   PreconditionRelaxation<MATRIX>::initialize(rA, parameters);
 }
 
@@ -726,7 +735,7 @@ PreconditionPSOR<MATRIX>::vmult (VECTOR& dst, const VECTOR& src) const
 {
   Assert (this->A!=0, ExcNotInitialized());
   dst = src;
-  this->A->PSOR (dst, *permutation, this->relaxation);
+  this->A->PSOR (dst, *permutation, *inverse_permutation, this->relaxation);
 }
 
 
@@ -738,7 +747,7 @@ PreconditionPSOR<MATRIX>::Tvmult (VECTOR& dst, const VECTOR& src) const
 {
   Assert (this->A!=0, ExcNotInitialized());
   dst = src;
-  this->A->TPSOR (dst, *permutation, this->relaxation);
+  this->A->TPSOR (dst, *permutation, *inverse_permutation, this->relaxation);
 }
 
 
index caa15d8dcbf83437b1e570967ca46a53db050646..0f664b8a34ac0455095cf689d4ebecbc143284e4 100644 (file)
@@ -67,6 +67,16 @@ template void SparseMatrix<TYPEMAT>::TSOR<TYPE2> (Vector<TYPE2> &,
                                                  const TYPEMAT) const;
 template void SparseMatrix<TYPEMAT>::SSOR<TYPE2> (Vector<TYPE2> &,
                                                  const TYPEMAT) const;
+template void SparseMatrix<TYPEMAT>::PSOR<TYPE2> (
+  Vector<TYPE2> &,
+  const std::vector<unsigned int>&,
+  const std::vector<unsigned int>&,
+  const TYPEMAT) const;
+template void SparseMatrix<TYPEMAT>::TPSOR<TYPE2> (
+  Vector<TYPE2> &,
+  const std::vector<unsigned int>&,
+  const std::vector<unsigned int>&,
+  const TYPEMAT) const;
 template void SparseMatrix<TYPEMAT>::SOR_step<TYPE2> (Vector<TYPE2> &, const Vector<TYPE2> &,
                                                      const TYPEMAT) const;
 template void SparseMatrix<TYPEMAT>::TSOR_step<TYPE2> (Vector<TYPE2> &,
index e69fc38545b9000cd6d5abcfef805826049fa92a..b35345df640b41d83f21f92345915c01b4a69f74 100644 (file)
@@ -814,14 +814,17 @@ class SparseMatrix : public Subscriptor
                                      * that is, first the row
                                      * @p{permutation[0]}, then
                                      * @p{permutation[1]} and so
-                                     * on.
+                                     * on. For efficiency reasons,
+                                     * the permutation as well as its
+                                     * inverse are required.
                                      *
                                      * @p{omega} is the relaxation
                                      * parameter.
                                      */
     template <typename somenumber>
     void PSOR (Vector<somenumber> &v,
-             const std::vector<unsigned int> permutation,
+             const std::vector<unsigned int>& permutation,
+             const std::vector<unsigned int>& inverse_permutation,
              const number        om = 1.) const;
 
                                     /**
@@ -834,14 +837,17 @@ class SparseMatrix : public Subscriptor
                                      * that is, first the row
                                      * @p{permutation[m()-1]}, then
                                      * @p{permutation[m()-2]} and so
-                                     * on.
+                                     * on. For efficiency reasons,
+                                     * the permutation as well as its
+                                     * inverse are required.
                                      *
                                      * @p{omega} is the relaxation
                                      * parameter.
                                      */
     template <typename somenumber>
     void TPSOR (Vector<somenumber> &v,
-             const std::vector<unsigned int> permutation,
+             const std::vector<unsigned int>& permutation,
+             const std::vector<unsigned int>& inverse_permutation,
              const number        om = 1.) const;
 
                                     /**
index f5cf298da69d515c92800ed883c12503866eff1a..3a9659c2dbab2371a8901cca7569b607cf85722a 100644 (file)
@@ -962,9 +962,12 @@ SparseMatrix<number>::SOR (Vector<somenumber>& dst, const number om) const
     {
       somenumber s = dst(row);
       for (unsigned int j=cols->rowstart[row]; j<cols->rowstart[row+1]; ++j)
-       if (cols->colnums[j] < row)
-         s -= val[j] * dst(cols->colnums[j]);
-
+       {
+         const unsigned int col = cols->colnums[j];
+         if (col < row)
+         s -= val[j] * dst(col);
+       }
+      
       dst(row) = s * om / val[cols->rowstart[row]];
     }
 }
@@ -998,7 +1001,8 @@ template <typename somenumber>
 void
 SparseMatrix<number>::PSOR (
   Vector<somenumber>& dst,
-  const std::vector<unsigned int> permutation,
+  const std::vector<unsigned int>& permutation,
+  const std::vector<unsigned int>& inverse_permutation,
   const number om) const
 {
   Assert (cols != 0, ExcMatrixNotInitialized());
@@ -1007,6 +1011,8 @@ SparseMatrix<number>::PSOR (
   Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
   Assert (m() == permutation.size(),
          ExcDimensionMismatch(m(), permutation.size()));
+  Assert (m() == inverse_permutation.size(),
+         ExcDimensionMismatch(m(), inverse_permutation.size()));
 
   for (unsigned int urow=0; urow<m(); ++urow)
     {
@@ -1014,8 +1020,8 @@ SparseMatrix<number>::PSOR (
       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)
+         const unsigned int col = cols->colnums[j];
+         if (inverse_permutation[col] < urow)
            s -= val[j] * dst(col);
        }
 
@@ -1029,7 +1035,8 @@ template <typename somenumber>
 void
 SparseMatrix<number>::TPSOR (
   Vector<somenumber>& dst,
-  const std::vector<unsigned int> permutation,
+  const std::vector<unsigned int>& permutation,
+  const std::vector<unsigned int>& inverse_permutation,
   const number om) const
 {
   Assert (cols != 0, ExcMatrixNotInitialized());
@@ -1038,6 +1045,8 @@ SparseMatrix<number>::TPSOR (
   Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
   Assert (m() == permutation.size(),
          ExcDimensionMismatch(m(), permutation.size()));
+  Assert (m() == inverse_permutation.size(),
+         ExcDimensionMismatch(m(), inverse_permutation.size()));
 
   for (unsigned int urow=m(); urow != 0;)
     {
@@ -1046,8 +1055,8 @@ SparseMatrix<number>::TPSOR (
       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)
+         const unsigned int col = cols->colnums[j];
+         if (inverse_permutation[col] > urow)
            s -= val[j] * dst(col);
        }
 
index dc54ef4c80f3c37794cd127d20fbc10063e081c7..538d9c71b610757fa9dbdf5cc6163be3d8a8a1cc 100644 (file)
@@ -31,7 +31,7 @@
 
 template<class SOLVER, class MATRIX, class VECTOR, class PRECONDITION>
 void
-check_method( SOLVER& solver, const MATRIX& A,
+check_solve( SOLVER& solver, const MATRIX& A,
             VECTOR& u, VECTOR& f, const PRECONDITION& P)
 {
   u = 0.;
@@ -46,6 +46,23 @@ check_method( SOLVER& solver, const MATRIX& A,
     }  
 }
 
+template<class SOLVER, class MATRIX, class VECTOR, class PRECONDITION>
+void
+check_Tsolve(SOLVER& solver, const MATRIX& A,
+            VECTOR& u, VECTOR& f, const PRECONDITION& P)
+{
+  u = 0.;
+  f = 1.;
+  try 
+    {
+      solver.Tsolve(A,u,f,P);
+    }
+  catch (std::exception& e)
+    {
+      deallog << e.what() << std::endl;
+    }  
+}
+
 int main()
 {
   std::ofstream logfile("solver.output");
@@ -59,6 +76,8 @@ int main()
   SolverControl verbose_control(100, 1.e-3, true);
   SolverCG<> cg(control, mem);
   SolverGMRES<> gmres(control, mem, 8);
+  SolverGMRES<>::AdditionalData(8, true);
+  SolverGMRES<> gmresright(control, mem, 8);  
   SolverBicgstab<> bicgstab(control, mem);
   SolverRichardson<> rich(control, mem);
   SolverQMRS<> qmrs(control, mem);
@@ -84,11 +103,17 @@ int main()
       prec_ssor.initialize(A, 1.2);
 
       std::vector<unsigned int> permutation(dim);
+      std::vector<unsigned int> inverse_permutation(dim);
       for (unsigned int i=0;i<dim;++i)
-       permutation(i) = dim-i-1;
+       permutation[i] = dim-i-1;
+      for (unsigned int i=0;i<dim;++i)
+       inverse_permutation[permutation[i]] = i;
+
+      for (unsigned int i=0;i<dim;++i)
+       permutation[i] = dim-i-1;
       
       PreconditionPSOR<> prec_psor;
-      prec_psor.initialize(A, permutation, 1.2);
+      prec_psor.initialize(A, permutation, inverse_permutation, 1.2);
       
       Vector<double>  f(dim);
       Vector<double>  u(dim);
@@ -110,46 +135,56 @@ int main()
          deallog.push("no-fail");
 
          control.set_max_steps(10);
-         check_method(cg,A,u,f,prec_no);
-         check_method(bicgstab,A,u,f,prec_no);
-         check_method(gmres,A,u,f,prec_no);
-         check_method(qmrs,A,u,f,prec_no);
+         check_solve(cg,A,u,f,prec_no);
+         check_solve(bicgstab,A,u,f,prec_no);
+         check_solve(gmres,A,u,f,prec_no);
+         check_solve(gmresright,A,u,f,prec_no);
+         check_solve(qmrs,A,u,f,prec_no);
          control.set_max_steps(100);
          
          deallog.pop();
          
          deallog.push("no");
          
-         check_method(cg,A,u,f,prec_no);
-         check_method(bicgstab,A,u,f,prec_no);
-         check_method(gmres,A,u,f,prec_no);
-         check_method(qmrs,A,u,f,prec_no);
+         check_solve(cg,A,u,f,prec_no);
+         check_solve(bicgstab,A,u,f,prec_no);
+         check_solve(gmres,A,u,f,prec_no);
+         check_solve(gmresright,A,u,f,prec_no);
+         check_solve(qmrs,A,u,f,prec_no);
          
          deallog.pop();
          
          deallog.push("ssor");
          
-         check_method(rich,A,u,f,prec_ssor);
-         check_method(cg,A,u,f,prec_ssor);
-         check_method(bicgstab,A,u,f,prec_ssor);
-         check_method(gmres,A,u,f,prec_ssor);
-         check_method(qmrs,A,u,f,prec_ssor);
+         check_Tsolve(rich,A,u,f,prec_ssor);
+         check_solve(rich,A,u,f,prec_ssor);
+         check_solve(cg,A,u,f,prec_ssor);
+         check_solve(bicgstab,A,u,f,prec_ssor);
+         check_solve(gmres,A,u,f,prec_ssor);
+         check_solve(gmresright,A,u,f,prec_ssor);
+         check_solve(qmrs,A,u,f,prec_ssor);
          
          deallog.pop();
          
          deallog.push("sor");
          
-         check_method(rich,A,u,f,prec_sor);
-         check_method(cg,A,u,f,prec_sor);
-         check_method(bicgstab,A,u,f,prec_sor);
-         check_method(gmres,A,u,f,prec_sor);
+         check_Tsolve(rich,A,u,f,prec_sor);
+         check_solve(rich,A,u,f,prec_sor);
+         check_solve(cg,A,u,f,prec_sor);
+         check_solve(bicgstab,A,u,f,prec_sor);
+         check_solve(gmres,A,u,f,prec_sor);
+         check_solve(gmresright,A,u,f,prec_sor);
+         
+         deallog.pop();
          
          deallog.push("psor");
          
-         check_method(rich,A,u,f,prec_psor);
-         check_method(cg,A,u,f,prec_psor);
-         check_method(bicgstab,A,u,f,prec_psor);
-         check_method(gmres,A,u,f,prec_psor);
+         check_Tsolve(rich,A,u,f,prec_psor);
+         check_solve(rich,A,u,f,prec_psor);
+         check_solve(cg,A,u,f,prec_psor);
+         check_solve(bicgstab,A,u,f,prec_psor);
+         check_solve(gmres,A,u,f,prec_psor);
+         check_solve(gmresright,A,u,f,prec_psor);
          
          deallog.pop();
        }
index 4513f30b2e470a1af066187d2aa6044df4a2618c..a964277c44e5b4a380d3f73ca2ab4b580f78de7b 100644 (file)
@@ -7,6 +7,8 @@ DEAL:no-fail:Bicgstab::Starting value 3.000
 DEAL:no-fail:Bicgstab::Convergence step 3 value 7.850e-16
 DEAL:no-fail:GMRES::Starting value 3.000
 DEAL:no-fail:GMRES::Convergence step 3 value 2.266e-17
+DEAL:no-fail:GMRES::Starting value 3.000
+DEAL:no-fail:GMRES::Convergence step 3 value 2.266e-17
 DEAL:no-fail:QMRS::Starting value 3.000
 DEAL:no-fail:QMRS::Convergence step 3 value 1.183e-16
 DEAL:no:cg::Starting value 3.000
@@ -15,8 +17,12 @@ DEAL:no:Bicgstab::Starting value 3.000
 DEAL:no:Bicgstab::Convergence step 3 value 7.850e-16
 DEAL:no:GMRES::Starting value 3.000
 DEAL:no:GMRES::Convergence step 3 value 2.266e-17
+DEAL:no:GMRES::Starting value 3.000
+DEAL:no:GMRES::Convergence step 3 value 2.266e-17
 DEAL:no:QMRS::Starting value 3.000
 DEAL:no:QMRS::Convergence step 3 value 1.183e-16
+DEAL:ssor:RichardsonT::Starting value 3.000
+DEAL:ssor:RichardsonT::Convergence step 10 value 0.0004418
 DEAL:ssor:Richardson::Starting value 3.000
 DEAL:ssor:Richardson::Convergence step 10 value 0.0004418
 DEAL:ssor:cg::Starting value 3.000
@@ -25,8 +31,12 @@ DEAL:ssor:Bicgstab::Starting value 3.000
 DEAL:ssor:Bicgstab::Convergence step 2 value 0.0003670
 DEAL:ssor:GMRES::Starting value 1.600
 DEAL:ssor:GMRES::Convergence step 3 value 0.0002978
+DEAL:ssor:GMRES::Starting value 1.600
+DEAL:ssor:GMRES::Convergence step 3 value 0.0002978
 DEAL:ssor:QMRS::Starting value 3.000
 DEAL:ssor:QMRS::Convergence step 4 value 0.0001345
+DEAL:sor:RichardsonT::Starting value 3.000
+DEAL:sor:RichardsonT::Convergence step 7 value 0.0004339
 DEAL:sor:Richardson::Starting value 3.000
 DEAL:sor:Richardson::Convergence step 7 value 0.0004339
 DEAL:sor:cg::Starting value 3.000
@@ -36,6 +46,21 @@ DEAL:sor:Bicgstab::Starting value 3.000
 DEAL:sor:Bicgstab::Convergence step 5 value 1.632e-15
 DEAL:sor:GMRES::Starting value 1.462
 DEAL:sor:GMRES::Convergence step 5 value 1.278e-18
+DEAL:sor:GMRES::Starting value 1.462
+DEAL:sor:GMRES::Convergence step 5 value 1.278e-18
+DEAL:psor:RichardsonT::Starting value 3.000
+DEAL:psor:RichardsonT::Convergence step 7 value 0.0004339
+DEAL:psor:Richardson::Starting value 3.000
+DEAL:psor:Richardson::Convergence step 7 value 0.0004339
+DEAL:psor:cg::Starting value 3.000
+DEAL:psor:cg::Failure step 100 value 0.2585
+DEAL:psor::Iterative method reported convergence failure in step 100 with residual 0.258509
+DEAL:psor:Bicgstab::Starting value 3.000
+DEAL:psor:Bicgstab::Convergence step 5 value 1.797e-15
+DEAL:psor:GMRES::Starting value 1.462
+DEAL:psor:GMRES::Convergence step 5 value 1.208e-18
+DEAL:psor:GMRES::Starting value 1.462
+DEAL:psor:GMRES::Convergence step 5 value 1.208e-18
 DEAL::Size 12 Unknowns 121
 DEAL::SOR-diff:0.000
 DEAL:no-fail:cg::Starting value 11.00
@@ -47,6 +72,9 @@ DEAL:no-fail::Iterative method reported convergence failure in step 10 with resi
 DEAL:no-fail:GMRES::Starting value 11.00
 DEAL:no-fail:GMRES::Failure step 10 value 0.8414
 DEAL:no-fail::Iterative method reported convergence failure in step 10 with residual 0.841354
+DEAL:no-fail:GMRES::Starting value 11.00
+DEAL:no-fail:GMRES::Failure step 10 value 0.8414
+DEAL:no-fail::Iterative method reported convergence failure in step 10 with residual 0.841354
 DEAL:no-fail:QMRS::Starting value 11.00
 DEAL:no-fail:QMRS::Failure step 10 value 0.4215
 DEAL:no-fail::Iterative method reported convergence failure in step 10 with residual 0.421504
@@ -56,8 +84,12 @@ DEAL:no:Bicgstab::Starting value 11.00
 DEAL:no:Bicgstab::Convergence step 11 value 0.0003990
 DEAL:no:GMRES::Starting value 11.00
 DEAL:no:GMRES::Convergence step 43 value 0.0009788
+DEAL:no:GMRES::Starting value 11.00
+DEAL:no:GMRES::Convergence step 43 value 0.0009788
 DEAL:no:QMRS::Starting value 11.00
 DEAL:no:QMRS::Convergence step 16 value 0.0002583
+DEAL:ssor:RichardsonT::Starting value 11.00
+DEAL:ssor:RichardsonT::Convergence step 59 value 0.0008892
 DEAL:ssor:Richardson::Starting value 11.00
 DEAL:ssor:Richardson::Convergence step 59 value 0.0008892
 DEAL:ssor:cg::Starting value 11.00
@@ -66,8 +98,12 @@ DEAL:ssor:Bicgstab::Starting value 11.00
 DEAL:ssor:Bicgstab::Convergence step 5 value 9.315e-05
 DEAL:ssor:GMRES::Starting value 11.05
 DEAL:ssor:GMRES::Convergence step 7 value 0.0006838
+DEAL:ssor:GMRES::Starting value 11.05
+DEAL:ssor:GMRES::Convergence step 7 value 0.0006838
 DEAL:ssor:QMRS::Starting value 11.00
 DEAL:ssor:QMRS::Convergence step 9 value 0.0004140
+DEAL:sor:RichardsonT::Starting value 11.00
+DEAL:sor:RichardsonT::Convergence step 88 value 0.0009636
 DEAL:sor:Richardson::Starting value 11.00
 DEAL:sor:Richardson::Convergence step 88 value 0.0009636
 DEAL:sor:cg::Starting value 11.00
@@ -77,5 +113,20 @@ DEAL:sor:Bicgstab::Starting value 11.00
 DEAL:sor:Bicgstab::Convergence step 14 value 0.0009987
 DEAL:sor:GMRES::Starting value 7.322
 DEAL:sor:GMRES::Convergence step 23 value 0.0006981
-DEAL::GrowingVectorMemory:Overall allocated vectors: 188
+DEAL:sor:GMRES::Starting value 7.322
+DEAL:sor:GMRES::Convergence step 23 value 0.0006981
+DEAL:psor:RichardsonT::Starting value 11.00
+DEAL:psor:RichardsonT::Convergence step 88 value 0.0009636
+DEAL:psor:Richardson::Starting value 11.00
+DEAL:psor:Richardson::Convergence step 88 value 0.0009636
+DEAL:psor:cg::Starting value 11.00
+DEAL:psor:cg::Failure step 100 value 5.643
+DEAL:psor::Iterative method reported convergence failure in step 100 with residual 5.64329
+DEAL:psor:Bicgstab::Starting value 11.00
+DEAL:psor:Bicgstab::Convergence step 14 value 0.0009987
+DEAL:psor:GMRES::Starting value 7.322
+DEAL:psor:GMRES::Convergence step 23 value 0.0006981
+DEAL:psor:GMRES::Starting value 7.322
+DEAL:psor:GMRES::Convergence step 23 value 0.0006981
+DEAL::GrowingVectorMemory:Overall allocated vectors: 312
 DEAL::GrowingVectorMemory:Maximum allocated vectors: 8

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.