parameterized reinit functions and some constructor fixes

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

diff --git a/deal.II/lac/include/lac/block_vector.h b/deal.II/lac/include/lac/block_vector.h
index 5aeff86099ed43042900d54218645add6c73db3d..a9fd7c994c5d36723106d4db7c9c215f64f3bace 100644 (file)
--- a/deal.II/lac/include/lac/block_vector.h
+++ b/deal.II/lac/include/lac/block_vector.h
@@ -639,8 +639,8 @@ class BlockVector
                                      *  use blocks of different
                                      *  sizes.
                                      */
-    BlockVector (unsigned int num_blocks = 0,
-                unsigned int block_size = 0);
+    explicit BlockVector (unsigned int num_blocks = 0,
+                         unsigned int block_size = 0);
     
                                     /**
                                      * Copy-Constructor. Dimension set to
@@ -748,7 +748,8 @@ class BlockVector
                                      * this function is the same as calling
                                      * @p{reinit (V.size(), fast)}.
                                      */
-    void reinit (const BlockVector<Number> &V,
+    template <typename Number2>
+    void reinit (const BlockVector<Number2> &V,
                 const bool                 fast=false);
     
                                     /**
@@ -1129,6 +1130,8 @@ class BlockVector
                                      */
     friend class iterator;
     friend class const_iterator;
+
+    template <typename Number2> friend class BlockVector;
 };
 
 

diff --git a/deal.II/lac/include/lac/block_vector.templates.h b/deal.II/lac/include/lac/block_vector.templates.h
index 59d31a3a93e8565b2b31a621f8e17bd547e9bbac..2b53f56671751ea8a86208ceeced82b321633651 100644 (file)
--- a/deal.II/lac/include/lac/block_vector.templates.h
+++ b/deal.II/lac/include/lac/block_vector.templates.h
@@ -96,7 +96,8 @@ void BlockVector<Number>::reinit (const std::vector<unsigned int> &n,
 
 
 template <typename Number>
-void BlockVector<Number>::reinit (const BlockVector<Number>& v,
+template <typename Number2>
+void BlockVector<Number>::reinit (const BlockVector<Number2>& v,
                                           const bool fast)
 {
   block_indices = v.block_indices;
@@ -493,7 +494,7 @@ BlockVector<Number>::operator = (const BlockVector<Number>& v)
 template <typename Number>
 template<typename Number2>
 BlockVector<Number>&
-BlockVector<Number>::operator = (const BlockVector< Number2>& v)
+BlockVector<Number>::operator = (const BlockVector<Number2>& v)
 {
   Assert (num_blocks == v.num_blocks,
          ExcDimensionMismatch(num_blocks, v.num_blocks));

diff --git a/deal.II/lac/include/lac/solver_richardson.h b/deal.II/lac/include/lac/solver_richardson.h
index 7a825b8078ce9f340f6238b5e0e08824b0122642..a59da8dc8af6ed6734d61c257e2f8f89984f53a5 100644 (file)
--- a/deal.II/lac/include/lac/solver_richardson.h
+++ b/deal.II/lac/include/lac/solver_richardson.h
@@ -57,14 +57,21 @@ class SolverRichardson : public Solver<VECTOR>
                                          * set the damping parameter
                                          * to one.
                                          */
-       AdditionalData(double omega=1):
-                       omega(omega)
+       AdditionalData(double omega=1,
+                      bool use_preconditioned_residual = true):
+                       omega(omega),
+                       use_preconditioned_residual(use_preconditioned_residual)
          {};
        
                                         /**
                                          * Relaxation parameter.
                                          */
        double omega;
+                                        /**
+                                         * Parameter for stopping criterion.
+                                         */
+       bool use_preconditioned_residual;
+         
     };
        
                                     /**
@@ -125,7 +132,7 @@ class SolverRichardson : public Solver<VECTOR>
                                      * Implementation of the computation of
                                      * the norm of the residual.
                                      */
-    virtual double criterion();
+    virtual typename VECTOR::value_type criterion();
     
                                     /**
                                      * Temporary vectors, allocated through
@@ -151,7 +158,7 @@ class SolverRichardson : public Solver<VECTOR>
                                      * norm of the residual vector and thus
                                      * the square root of the @p{res2} value.
                                      */
-    double res;
+    typename VECTOR::value_type res;
 };
 
 
@@ -196,13 +203,31 @@ SolverRichardson<VECTOR>::solve (const MATRIX         &A,
                                       // but do it in 2 steps
       A.vmult(r,x);
       r.sadd(-1.,1.,b);
-      res=sqrt(r*r);
-
-      conv = control().check (iter, criterion());
-      if (conv != SolverControl::iterate)
-       break;
 
+      if (!additional_data.use_preconditioned_residual)
+       {
+         res = r*r;
+//      cout << '[' << res << ' ';
+         res=sqrt(res);
+//      cout << res << ' ' << r.l1_norm() << ']';
+//      r.print(cout);
+         conv = control().check (iter, criterion());
+         if (conv != SolverControl::iterate)
+           break;
+       }
+      
       precondition.vmult(d,r);
+      if (additional_data.use_preconditioned_residual)
+       {
+         res = d*d;
+//      cout << '[' << res << ' ';
+         res=sqrt(res);
+//      cout << res << ' ' << r.l1_norm() << ']';
+//      r.print(cout);
+         conv = control().check (iter, criterion());
+         if (conv != SolverControl::iterate)
+           break;
+       }
       x.add(additional_data.omega,d);
       print_vectors(iter,x,r,d);
     }
@@ -281,7 +306,7 @@ SolverRichardson<VECTOR>::print_vectors(const unsigned int,
 
 
 template<class VECTOR>
-inline double
+inline typename VECTOR::value_type
 SolverRichardson<VECTOR>::criterion()
 {
   return res;

diff --git a/deal.II/lac/include/lac/vector.h b/deal.II/lac/include/lac/vector.h
index e4513eb27ed8166c4c3040f8bd6773189809c583..83e1363e1a10ba48d91d34efa34dc29ebb2215fa 100644 (file)
--- a/deal.II/lac/include/lac/vector.h
+++ b/deal.II/lac/include/lac/vector.h
@@ -76,7 +76,11 @@ class Vector
                                     /**
                                      * Copy-Constructor. Dimension set to
                                      * that of V, all components are copied
-                                     * from V
+                                     * from V.
+                                     *
+                                     * We would like to make this
+                                     * constructor explicit, but STL
+                                     * insists on using it implicitly.
                                      */
     Vector (const Vector<Number>& V);
 
@@ -102,7 +106,7 @@ class Vector
                                      * Constructor. Set dimension to @p{n} and
                                      * initialize all elements with zero.
                                      */
-    Vector (const unsigned int n);
+    explicit Vector (const unsigned int n);
 
                                     /**
                                      * Initialize the vector with a
@@ -160,7 +164,8 @@ class Vector
                                      * this function is the same as calling
                                      * @p{reinit (V.size(), fast)}.
                                      */
-    void reinit (const Vector<Number> &V,
+    template <typename Number2>
+    void reinit (const Vector<Number2> &V,
                 const bool            fast=false);
 
                                     /**
@@ -535,6 +540,8 @@ class Vector
                                      * Pointer to the array of components.
                                      */
     Number *val;
+
+    template <typename Number2> friend class Vector;    
 };
 
 

diff --git a/deal.II/lac/include/lac/vector.templates.h b/deal.II/lac/include/lac/vector.templates.h
index 842fdf818fa9dacab5addbd8df88e0ce34ded7cc..a1850e4fd7c9f716d795ca0bc1b07ff28c673bbd 100644 (file)
--- a/deal.II/lac/include/lac/vector.templates.h
+++ b/deal.II/lac/include/lac/vector.templates.h
@@ -71,7 +71,8 @@ Vector<Number>::Vector (const Vector<Number>& v) :
 
 
 template <typename Number>
-void Vector<Number>::reinit (const Vector<Number>& v, const bool fast)
+template <typename Number2>
+void Vector<Number>::reinit (const Vector<Number2>& v, const bool fast)
 {
   reinit (v.size(), fast);
 };

diff --git a/deal.II/lac/source/block_vector.cc b/deal.II/lac/source/block_vector.cc
index 601deb685b513db3b8e968e5255c460a26ee605e..b4cf91a8e9a4b7e49a33af7be20f2786ce0d33d1 100644 (file)
--- a/deal.II/lac/source/block_vector.cc
+++ b/deal.II/lac/source/block_vector.cc
@@ -15,7 +15,16 @@
 
 // explicit instantiations
 template class BlockVector<double>;
+template BlockVector<double>& BlockVector<double>::operator=(
+  const BlockVector<float>&);
+template void BlockVector<double>::reinit(const BlockVector<double>&, bool);
+template void BlockVector<double>::reinit(const BlockVector<float>&, bool);
+
 template class BlockVector<float>;
+template BlockVector<float>& BlockVector<float>::operator=(
+  const BlockVector<double>&);
+template void BlockVector<float>::reinit(const BlockVector<double>&, bool);
+template void BlockVector<float>::reinit(const BlockVector<float>&, bool);
 
 namespace BlockVectorIterators
 {

diff --git a/deal.II/lac/source/vector.cc b/deal.II/lac/source/vector.cc
index d1e40da33be40f3177fda411fd33bd8d6caedf2f..73a210a9be24d298f7843aa498d3b509dd655fce 100644 (file)
--- a/deal.II/lac/source/vector.cc
+++ b/deal.II/lac/source/vector.cc
@@ -19,11 +19,15 @@ template class Vector<double>;
 template Vector<double>& Vector<double>::operator=(const Vector<float>&);
 template double Vector<double>::operator*(const Vector<float>&) const;
 template double Vector<double>::operator*(const Vector<double>&) const;
+template void Vector<double>::reinit(const Vector<double>&, bool);
+template void Vector<double>::reinit(const Vector<float>&, bool);
 
 template class Vector<float>;
 template Vector<float>& Vector<float>::operator=(const Vector<double>&);
 template float Vector<float>::operator*(const Vector<float>&) const;
 template float Vector<float>::operator*(const Vector<double>&) const;
+template void Vector<float>::reinit(const Vector<double>&, bool);
+template void Vector<float>::reinit(const Vector<float>&, bool);
 
 // see the .h file for why these functions are disabled.
 // template Vector<float>::Vector (const Vector<double>& v);