Add precondition_Jacobi to FullMatrix.

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

diff --git a/deal.II/doc/news/2001/c-3-1.html b/deal.II/doc/news/2001/c-3-1.html
index e7726ec5cd4d93bc9c7fd7c4905f8b08ee093e89..132b7b6c2f51478aeb3aae9393d567eaa50a8bd7 100644 (file)
--- a/deal.II/doc/news/2001/c-3-1.html
+++ b/deal.II/doc/news/2001/c-3-1.html
@@ -333,7 +333,16 @@ documentation, etc</a>.
 
 <ol>
   <li> <p>
-       Changed: The <code class="class">Vector</code> and
+       New: There is now a function 
+       <code class="member">FullMatrix::precondition_Jacobi</code>. The 
+       <code class="class">PreconditionJacobi</code> class is
+       therefore now also applicable with the full matrix class.
+       <br>
+       (WB 2001/08/11)
+       </p>
+
+  <li> <p>
+       New: The <code class="class">Vector</code> and
        <code class="class">BlockVector</code> classes can now be
        initialized using a new constructor that takes two iterators
        that denote a range of elements which are to be copied.

diff --git a/deal.II/lac/include/lac/full_matrix.h b/deal.II/lac/include/lac/full_matrix.h
index 91d92c6dd30a41919c5ac24261a37e72161436a5..128a0436cbd13d1e125246d969ab2229930a89d5 100644 (file)
--- a/deal.II/lac/include/lac/full_matrix.h
+++ b/deal.II/lac/include/lac/full_matrix.h
@@ -452,6 +452,21 @@ class FullMatrix : public vector2d<number>
                                      */
     void invert (const FullMatrix<number> &M);
 
+                                    /**
+                                     * Apply the Jacobi
+                                     * preconditioner, which
+                                     * multiplies every element of
+                                     * the @p{src} vector by the
+                                     * inverse of the respective
+                                     * diagonal element and
+                                     * multiplies the result with the
+                                     * damping factor @p{omega}.
+                                     */
+    template <typename somenumber>
+    void precondition_Jacobi (Vector<somenumber>       &dst,
+                             const Vector<somenumber> &src,
+                             const number              omega = 1.) const;
+
                                     /**
                                      * $A(i,1-n)+=s*A(j,1-n)$.
                                      * Simple addition of rows of this

diff --git a/deal.II/lac/include/lac/full_matrix.templates.h b/deal.II/lac/include/lac/full_matrix.templates.h
index 256d706e83a8a41f0094509db95333def4375a54..7a7e4894e8507307088b7630bb6d08fd9935249f 100644 (file)
--- a/deal.II/lac/include/lac/full_matrix.templates.h
+++ b/deal.II/lac/include/lac/full_matrix.templates.h
@@ -1251,12 +1251,37 @@ FullMatrix<number>::invert (const FullMatrix<number> &M)
       }
 
       default:
+                                            // if no inversion is
+                                            // hardcoded, fall back
+                                            // to use the
+                                            // Gauss-Jordan algorithm
            *this = M;
            gauss_jordan();
     };    
 };
 
 
+template <typename number>
+template <typename somenumber>
+void
+FullMatrix<number>::precondition_Jacobi (Vector<somenumber>       &dst,
+                                        const Vector<somenumber> &src,
+                                        const number              om) const
+{
+  Assert (m() == n(), ExcNotQuadratic());
+  Assert (dst.size() == n(), ExcDimensionMismatch (dst.size(), n()));
+  Assert (src.size() == n(), ExcDimensionMismatch (src.size(), n()));
+
+  const unsigned int n = src.size();
+  somenumber              *dst_ptr = dst.begin();
+  const somenumber        *src_ptr = src.begin();
+  
+  for (unsigned int i=0; i<n; ++i, ++dst_ptr, ++src_ptr)
+    *dst_ptr = om * *src_ptr / el(i,i);
+};
+
+
+
 template <typename number>
 void
 FullMatrix<number>::print_formatted (std::ostream       &out,

diff --git a/deal.II/lac/source/full_matrix.double.cc b/deal.II/lac/source/full_matrix.double.cc
index 15986ce5bd920c5e3206a64d8f73641d94b8a562..dc7e460bfd1e3c937a94de5de606090292e9cc3c 100644 (file)
--- a/deal.II/lac/source/full_matrix.double.cc
+++ b/deal.II/lac/source/full_matrix.double.cc
@@ -46,6 +46,12 @@ template void FullMatrix<TYPEMAT>::backward(Vector<TYPEVEC>&, const Vector<TYPEV
 template void FullMatrix<TYPEMAT>::householder(Vector<TYPEVEC>&);
 template double FullMatrix<TYPEMAT>::least_squares(Vector<TYPEVEC>&, Vector<TYPEVEC>&);
 
+template
+void FullMatrix<TYPEMAT>::precondition_Jacobi (Vector<TYPEVEC> &,
+                                              const Vector<TYPEVEC> &,
+                                              const TYPEMAT) const;
+
+
 #undef TYPEVEC
 #define TYPEVEC float
 
@@ -62,6 +68,11 @@ template void FullMatrix<TYPEMAT>::backward(Vector<TYPEVEC>&, const Vector<TYPEV
 template void FullMatrix<TYPEMAT>::householder(Vector<TYPEVEC>&);
 template double FullMatrix<TYPEMAT>::least_squares(Vector<TYPEVEC>&, Vector<TYPEVEC>&);
 
+template
+void FullMatrix<TYPEMAT>::precondition_Jacobi (Vector<TYPEVEC> &,
+                                              const Vector<TYPEVEC> &,
+                                              const TYPEMAT) const;
+
 #undef TYPERES
 #define TYPERES float
 

diff --git a/deal.II/lac/source/full_matrix.float.cc b/deal.II/lac/source/full_matrix.float.cc
index 04e8f8f9aba84bea47db4bea71dfaf97e2e40394..84a078f7f4307a2c562572650eb03e92711bd415 100644 (file)
--- a/deal.II/lac/source/full_matrix.float.cc
+++ b/deal.II/lac/source/full_matrix.float.cc
@@ -44,6 +44,10 @@ template void FullMatrix<TYPEMAT>::forward(Vector<TYPEVEC>&, const Vector<TYPEVE
 template void FullMatrix<TYPEMAT>::backward(Vector<TYPEVEC>&, const Vector<TYPEVEC>&) const;
 template void FullMatrix<TYPEMAT>::householder(Vector<TYPEVEC>&);
 template double FullMatrix<TYPEMAT>::least_squares(Vector<TYPEVEC>&, Vector<TYPEVEC>&);
+template
+void FullMatrix<TYPEMAT>::precondition_Jacobi (Vector<TYPEVEC> &,
+                                              const Vector<TYPEVEC> &,
+                                              const TYPEMAT) const;
 
 #undef TYPEVEC
 #define TYPEVEC float
@@ -60,6 +64,10 @@ template void FullMatrix<TYPEMAT>::forward(Vector<TYPEVEC>&, const Vector<TYPEVE
 template void FullMatrix<TYPEMAT>::backward(Vector<TYPEVEC>&, const Vector<TYPEVEC>&) const;
 template void FullMatrix<TYPEMAT>::householder(Vector<TYPEVEC>&);
 template double FullMatrix<TYPEMAT>::least_squares(Vector<TYPEVEC>&, Vector<TYPEVEC>&);
+template
+void FullMatrix<TYPEMAT>::precondition_Jacobi (Vector<TYPEVEC> &,
+                                              const Vector<TYPEVEC> &,
+                                              const TYPEMAT) const;
 
 
 #undef TYPERES