From: wolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Sat, 11 Aug 2001 22:03:05 +0000 (+0000)
Subject: Add precondition_Jacobi to FullMatrix.
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=edb33a0e9a01771dd654b8f104e638f6a621a01d;p=dealii-svn.git

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 e7726ec5cd..132b7b6c2f 100644
--- 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 91d92c6dd3..128a0436cb 100644
--- 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 256d706e83..7a7e4894e8 100644
--- 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 15986ce5bd..dc7e460bfd 100644
--- 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 04e8f8f9ab..84a078f7f4 100644
--- 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