From: bangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Tue, 21 Apr 2009 13:44:29 +0000 (+0000)
Subject: After half a year, finally check in some more of the conversion that allows to use... 
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=847d7ec66f1e3b61f84805e0c9f5b50266f59af4;p=dealii-svn.git

After half a year, finally check in some more of the conversion that allows to use the SparseMatrix also for complex-valued template arguments.

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 f3948e2ee0..b6cf3d02ba 100644
--- 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 dc3b73c97c..3fe59f41ff 100644
--- 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 0edd36c17d..a8dafdb7cb 100644
--- 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;
+//   }