After half a year, finally check in some more of the conversion that allows to use...

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

diff --git a/deal.II/lac/include/lac/sparse_matrix.h b/deal.II/lac/include/lac/sparse_matrix.h
index f3948e2ee0042c7fe3d7976e2e0f8b5d4b856c5e..b6cf3d02ba7226b8a5b077c54617af25104ce2df 100644 (file)
--- a/deal.II/lac/include/lac/sparse_matrix.h
+++ b/deal.II/lac/include/lac/sparse_matrix.h
@@ -441,7 +441,7 @@ namespace internals
  * @<double@></tt>; others can be generated in application programs (see the
  * section on @ref Instantiations in the manual).
  *
- * @author Essentially everyone who has ever worked on deal.II, 1994-2004
+ * @author Essentially everyone who has ever worked on deal.II, 1994-2007
  */
 template <typename number>
 class SparseMatrix : public virtual Subscriptor
@@ -453,6 +453,27 @@ class SparseMatrix : public virtual Subscriptor
                                      */
     typedef number value_type;
 
+                                    /**
+                                     * Declare a type that has holds
+                                     * real-valued numbers with the
+                                     * same precision as the template
+                                     * argument to this class. If the
+                                     * template argument of this
+                                     * class is a real data type,
+                                     * then real_type equals the
+                                     * template argument. If the
+                                     * template argument is a
+                                     * std::complex type then
+                                     * real_type equals the type
+                                     * underlying the complex
+                                     * numbers.
+                                     *
+                                     * This typedef is used to
+                                     * represent the return type of
+                                     * norms.
+                                     */
+    typedef typename numbers::NumberTraits<number>::real_type real_type;
+    
                                      /**
                                       * Typedef of an STL conforming iterator
                                       * class walking over all the nonzero
@@ -1376,8 +1397,17 @@ class SparseMatrix : public virtual Subscriptor
                                      * nodal values of the finite
                                      * element function.
                                      *
-                                     * Obviously, the matrix needs to
-                                     * be quadratic for this operation.
+                                     * Obviously, the matrix needs to be
+                                     * quadratic for this operation, and for
+                                     * the result to actually be a norm it
+                                     * also needs to be either real symmetric
+                                     * or complex hermitian.
+                                     *
+                                     * The underlying template types of both
+                                     * this matrix and the given vector
+                                     * should either both be real or
+                                     * complex-valued, but not mixed, for
+                                     * this function to make sense.
                                      */
     template <typename somenumber>
     somenumber matrix_norm_square (const Vector<somenumber> &v) const;
@@ -1389,6 +1419,7 @@ class SparseMatrix : public virtual Subscriptor
     template <typename somenumber>
     somenumber matrix_scalar_product (const Vector<somenumber> &u,
                                      const Vector<somenumber> &v) const;
+    
                                     /**
                                      * Compute the residual of an
                                      * equation <i>Mx=b</i>, where
@@ -1426,7 +1457,7 @@ class SparseMatrix : public virtual Subscriptor
                                      * $|Mv|_1\leq |M|_1 |v|_1$.
                                      * (cf. Haemmerlin-Hoffmann : Numerische Mathematik)
                                      */
-    number l1_norm () const;
+    real_type l1_norm () const;
 
                                     /**
                                      * Return the linfty-norm of the
@@ -1440,7 +1471,7 @@ class SparseMatrix : public virtual Subscriptor
                                      * $|Mv|_infty \leq |M|_infty |v|_infty$.
                                      * (cf. Haemmerlin-Hoffmann : Numerische Mathematik)
                                      */
-    number linfty_norm () const;
+    real_type linfty_norm () const;
 
                                      /**
                                       * Return the frobenius norm of the
@@ -1448,7 +1479,7 @@ class SparseMatrix : public virtual Subscriptor
                                       * sum of squares of all entries in the
                                       * matrix.
                                       */
-    number frobenius_norm () const;
+    real_type frobenius_norm () const;
 //@}
 /**
  * @name Preconditioning methods

diff --git a/deal.II/lac/include/lac/sparse_matrix.templates.h b/deal.II/lac/include/lac/sparse_matrix.templates.h
index dc3b73c97c7da4dfdf1a4371e24a24e073cab11d..3fe59f41ffe5c94d807b5c20cfc51c1dfb05c297 100644 (file)
--- a/deal.II/lac/include/lac/sparse_matrix.templates.h
+++ b/deal.II/lac/include/lac/sparse_matrix.templates.h
@@ -41,6 +41,7 @@
 #include <base/thread_management.h>
 #include <base/multithread_info.h>
 
+
 DEAL_II_NAMESPACE_OPEN
 
 
@@ -778,11 +779,11 @@ SparseMatrix<number>::matrix_norm_square (const Vector<somenumber>& v) const
          while (val_ptr != val_end_of_row)
            s += *val_ptr++ * v(*colnum_ptr++);
          
-         sum += s* v(row);
-       };
+         sum += v(row) * numbers::NumberTraits<somenumber>::conjugate(s);
+       }
       
       return sum;
-    };
+    }
 }
 
 
@@ -810,8 +811,9 @@ threaded_matrix_norm_square (const Vector<somenumber> &v,
       while (val_ptr != val_end_of_row)
        s += *val_ptr++ * v(*colnum_ptr++);
 
-      sum += s* v(row);
-    };
+      sum += v(row) * numbers::NumberTraits<somenumber>::conjugate(s);
+    }
+  
   *partial_sum = sum;
 }
 
@@ -891,11 +893,11 @@ SparseMatrix<number>::matrix_scalar_product (const Vector<somenumber>& u,
          while (val_ptr != val_end_of_row)
            s += *val_ptr++ * v(*colnum_ptr++);
          
-         sum += s* u(row);
-       };
+         sum += u(row) * numbers::NumberTraits<somenumber>::conjugate(s);
+       }
       
       return sum;
-    };
+    }
 }
 
 
@@ -924,8 +926,9 @@ threaded_matrix_scalar_product (const Vector<somenumber> &u,
       while (val_ptr != val_end_of_row)
        s += *val_ptr++ * v(*colnum_ptr++);
 
-      sum += s* u(row);
-    };
+      sum += u(row) * numbers::NumberTraits<somenumber>::conjugate(s);
+    }
+  
   *partial_sum = sum;
 }
 
@@ -934,16 +937,17 @@ threaded_matrix_scalar_product (const Vector<somenumber> &u,
 
 
 template <typename number>
-number SparseMatrix<number>::l1_norm () const
+typename SparseMatrix<number>::real_type
+SparseMatrix<number>::l1_norm () const
 {
   Assert (cols != 0, ExcNotInitialized());
   Assert (val != 0, ExcNotInitialized());
 
-  Vector<number> column_sums(n());
+  Vector<real_type> column_sums(n());
   const unsigned int n_rows = m();
   for (unsigned int row=0; row<n_rows; ++row)
     for (unsigned int j=cols->rowstart[row]; j<cols->rowstart[row+1] ; ++j)
-      column_sums(cols->colnums[j]) += std::fabs(val[j]);
+      column_sums(cols->colnums[j]) += numbers::NumberTraits<number>::abs(val[j]);
 
   return column_sums.linfty_norm();
 }
@@ -951,21 +955,22 @@ number SparseMatrix<number>::l1_norm () const
 
 
 template <typename number>
-number SparseMatrix<number>::linfty_norm () const
+typename SparseMatrix<number>::real_type
+SparseMatrix<number>::linfty_norm () const
 {
   Assert (cols != 0, ExcNotInitialized());
   Assert (val != 0, ExcNotInitialized());
 
   const number *val_ptr = &val[cols->rowstart[0]];
 
-  number sum, max=0;
+  real_type max=0;
   const unsigned int n_rows = m();
   for (unsigned int row=0; row<n_rows; ++row)
     {
-      sum=0;
+      real_type sum = 0;
       const number *const val_end_of_row = &val[cols->rowstart[row+1]];
       while (val_ptr != val_end_of_row)
-       sum += std::fabs(*val_ptr++);
+       sum += numbers::NumberTraits<number>::abs(*val_ptr++);
       if (sum > max)
        max = sum;
     }
@@ -975,16 +980,17 @@ number SparseMatrix<number>::linfty_norm () const
 
 
 template <typename number>
-number SparseMatrix<number>::frobenius_norm () const
+typename SparseMatrix<number>::real_type
+SparseMatrix<number>::frobenius_norm () const
 {
                                    // simply add up all entries in the
                                    // sparsity pattern, without taking any
                                    // reference to rows or columns
-  number norm_sqr = 0;
+  real_type norm_sqr = 0;
   const unsigned int n_rows = m();
   for (const number *ptr = &val[0];
        ptr != &val[cols->rowstart[n_rows]]; ++ptr)
-    norm_sqr += *ptr * *ptr;
+    norm_sqr +=  numbers::NumberTraits<number>::abs_square(*ptr);
 
   return std::sqrt (norm_sqr);
 }

diff --git a/deal.II/lac/source/sparse_matrix.inst.in b/deal.II/lac/source/sparse_matrix.inst.in
index 0edd36c17dd626355e84f01bef96b07369da5c5a..a8dafdb7cbcf45252f886847d2e47d161bf22a40 100644 (file)
--- a/deal.II/lac/source/sparse_matrix.inst.in
+++ b/deal.II/lac/source/sparse_matrix.inst.in
@@ -12,6 +12,8 @@
 //----------------------------  sparse_matrix.inst.in  ---------------------------
 
 
+// real instantiations
+
 for (S : REAL_SCALARS)
   {
     template class SparseMatrix<S>;
@@ -128,3 +130,110 @@ for (S1, S2, S3 : REAL_SCALARS;
     template void SparseMatrix<S1>::
       Tvmult_add (V1<S2> &, const V2<S3> &) const;
   }
+
+
+
+// complex instantiations
+
+// for (S : COMPLEX_SCALARS)
+//   {
+//     template class SparseMatrix<S>;
+//   }
+
+
+
+// for (S1, S2 : COMPLEX_SCALARS)
+//   {
+//     template SparseMatrix<S1> &
+//       SparseMatrix<S1>::copy_from<S2> (const SparseMatrix<S2> &);
+
+//     template 
+//       void SparseMatrix<S1>::copy_from<S2> (const FullMatrix<S2> &);
+
+//     template void SparseMatrix<S1>::add<S2> (const S1,
+//                                          const SparseMatrix<S2> &);
+//   }
+
+
+// for (S1, S2 : COMPLEX_SCALARS)
+//   {
+//     template S2
+//       SparseMatrix<S1>::
+//       matrix_norm_square<S2> (const Vector<S2> &) const;
+
+//     template S2
+//       SparseMatrix<S1>::
+//       matrix_scalar_product<S2> (const Vector<S2> &,
+//                              const Vector<S2> &) const;
+
+//     template S2 SparseMatrix<S1>::
+//       residual<S2> (Vector<S2> &,
+//                 const Vector<S2> &,
+//                 const Vector<S2> &) const;
+
+//     template void SparseMatrix<S1>::
+//       precondition_SSOR<S2> (Vector<S2> &,
+//                          const Vector<S2> &,
+//                          const S1) const;
+
+//     template void SparseMatrix<S1>::
+//       precondition_SOR<S2> (Vector<S2> &,
+//                         const Vector<S2> &,
+//                         const S1) const;
+
+//     template void SparseMatrix<S1>::
+//       precondition_TSOR<S2> (Vector<S2> &,
+//                          const Vector<S2> &,
+//                          const S1) const;
+
+//     template void SparseMatrix<S1>::
+//       precondition_Jacobi<S2> (Vector<S2> &,
+//                            const Vector<S2> &,
+//                            const S1) const;
+
+//     template void SparseMatrix<S1>::
+//       SOR<S2> (Vector<S2> &,
+//            const S1) const;
+//     template void SparseMatrix<S1>::
+//       TSOR<S2> (Vector<S2> &,
+//             const S1) const;
+//     template void SparseMatrix<S1>::
+//       SSOR<S2> (Vector<S2> &,
+//             const S1) const;
+//     template void SparseMatrix<S1>::
+//       PSOR<S2> (Vector<S2> &,
+//             const std::vector<unsigned int>&,
+//             const std::vector<unsigned int>&,
+//             const S1) const;
+//     template void SparseMatrix<S1>::
+//       TPSOR<S2> (Vector<S2> &,
+//              const std::vector<unsigned int>&,
+//              const std::vector<unsigned int>&,
+//              const S1) const;
+//     template void SparseMatrix<S1>::
+//       SOR_step<S2> (Vector<S2> &,
+//                 const Vector<S2> &,
+//                 const S1) const;
+//     template void SparseMatrix<S1>::
+//       TSOR_step<S2> (Vector<S2> &,
+//                  const Vector<S2> &,
+//                  const S1) const;
+//     template void SparseMatrix<S1>::
+//       SSOR_step<S2> (Vector<S2> &,
+//                  const Vector<S2> &, 
+//                  const S1) const;
+//   }
+
+
+// for (S1, S2, S3 : COMPLEX_SCALARS;
+//      V1, V2     : DEAL_II_VEC_TEMPLATES)
+//   {
+//     template void SparseMatrix<S1>::
+//       vmult (V1<S2> &, const V2<S3> &) const;
+//     template void SparseMatrix<S1>::
+//       Tvmult (V1<S2> &, const V2<S3> &) const;
+//     template void SparseMatrix<S1>::
+//       vmult_add (V1<S2> &, const V2<S3> &) const;
+//     template void SparseMatrix<S1>::
+//       Tvmult_add (V1<S2> &, const V2<S3> &) const;
+//   }