]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Make the remaining changes necessary to support complex-valued FullMatrix objects.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Sun, 18 Nov 2007 22:31:00 +0000 (22:31 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Sun, 18 Nov 2007 22:31:00 +0000 (22:31 +0000)
git-svn-id: https://svn.dealii.org/trunk@15518 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/full_matrix.templates.h
deal.II/lac/source/full_matrix.double.cc

index 0d3407c153e54c94c0b237000d7bbab44c67a309..54ab3a673beb5becb0963cf1317b6bae178e3fb7 100644 (file)
@@ -14,6 +14,8 @@
 #define __deal2__full_matrix_templates_h
 
 
+//TODO: this file has a lot of operations between matrices and matrices or matrices and vectors of different precision. we should go through the file and in each case pick the more accurate data type for intermediate results. currently, the choice is pretty much random. this may also allow us some operations where one operand is complex and the other is not
+
 #include <base/config.h>
 #include <lac/vector.h>
 #include <lac/full_matrix.h>
@@ -178,159 +180,29 @@ FullMatrix<number>::vmult (Vector<number2>& dst,
 
   Assert (&src != &dst, ExcSourceEqualsDestination());
 
-  if ((n()==3) && (m()==3))
-  {
-    number2 s;
-    number2 s0,s1,s2;
-    s   = src(0);
-    s0  = s*this->data()[0];
-    s1  = s*this->data()[3];
-    s2  = s*this->data()[6];
-    
-    s   = src(1);
-    s0 += s*this->data()[1];
-    s1 += s*this->data()[4];
-    s2 += s*this->data()[7];
-    
-    s   = src(2);
-    s0 += s*this->data()[2];
-    s1 += s*this->data()[5];
-    s2 += s*this->data()[8];
-
-    if (!adding)
-    {
-      dst(0) = s0;
-      dst(1) = s1;
-      dst(2) = s2;
-    }
-    else
-    {
-      dst(0) += s0;
-      dst(1) += s1;
-      dst(2) += s2;
-    }
-  }
-  else if ((n()==4) && (m()==4))
-  {
-    number2 s;
-    number2 s0,s1,s2,s3;
-    s = src(0);
-    s0  = s*this->data()[0];
-    s1  = s*this->data()[4];
-    s2  = s*this->data()[8];
-    s3  = s*this->data()[12];
-    
-    s = src(1);
-    s0 += s*this->data()[1];
-    s1 += s*this->data()[5];
-    s2 += s*this->data()[9];
-    s3 += s*this->data()[13];
-    
-    s = src(2);
-    s0 += s*this->data()[2];
-    s1 += s*this->data()[6];
-    s2 += s*this->data()[10];
-    s3 += s*this->data()[14];
-    
-    s = src(3);
-    s0 += s*this->data()[3];
-    s1 += s*this->data()[7];
-    s2 += s*this->data()[11];
-    s3 += s*this->data()[15];
-    
-    if (!adding)
-    {
-      dst(0) = s0;
-      dst(1) = s1;
-      dst(2) = s2;
-      dst(3) = s3;
-    }
-    else
-    {
-      dst(0) += s0;
-      dst(1) += s1;
-      dst(2) += s2;
-      dst(3) += s3;
-    }
-  }
-  else if ((n()==8) && (m()==8))
-  {
-    number2 s;
-    number2 s0,s1,s2,s3,s4,s5,s6,s7;
-    s = src(0);
-    s0 = s*this->data()[0]; s1 = s*this->data()[8]; s2 = s*this->data()[16]; s3 = s*this->data()[24];
-    s4 = s*this->data()[32]; s5 = s*this->data()[40]; s6 = s*this->data()[48]; s7 = s*this->data()[56];
-    s = src(1);
-    s0 += s*this->data()[1]; s1 += s*this->data()[9]; s2 += s*this->data()[17]; s3 += s*this->data()[25];
-    s4 += s*this->data()[33]; s5 += s*this->data()[41]; s6 += s*this->data()[49]; s7 += s*this->data()[57];
-    s = src(2);
-    s0 += s*this->data()[2]; s1 += s*this->data()[10]; s2 += s*this->data()[18]; s3 += s*this->data()[26];
-    s4 += s*this->data()[34]; s5 += s*this->data()[42]; s6 += s*this->data()[50]; s7 += s*this->data()[58];
-    s = src(3);
-    s0 += s*this->data()[3]; s1 += s*this->data()[11]; s2 += s*this->data()[19]; s3 += s*this->data()[27];
-    s4 += s*this->data()[35]; s5 += s*this->data()[43]; s6 += s*this->data()[51]; s7 += s*this->data()[59];
-    s = src(4);
-    s0 += s*this->data()[4]; s1 += s*this->data()[12]; s2 += s*this->data()[20]; s3 += s*this->data()[28];
-    s4 += s*this->data()[36]; s5 += s*this->data()[44]; s6 += s*this->data()[52]; s7 += s*this->data()[60];
-    s = src(5);
-    s0 += s*this->data()[5]; s1 += s*this->data()[13]; s2 += s*this->data()[21]; s3 += s*this->data()[29];
-    s4 += s*this->data()[37]; s5 += s*this->data()[45]; s6 += s*this->data()[53]; s7 += s*this->data()[61];
-    s = src(6);
-    s0 += s*this->data()[6]; s1 += s*this->data()[14]; s2 += s*this->data()[22]; s3 += s*this->data()[30];
-    s4 += s*this->data()[38]; s5 += s*this->data()[46]; s6 += s*this->data()[54]; s7 += s*this->data()[62];
-    s = src(7);
-    s0 += s*this->data()[7]; s1 += s*this->data()[15]; s2 += s*this->data()[23]; s3 += s*this->data()[31];
-    s4 += s*this->data()[39]; s5 += s*this->data()[47]; s6 += s*this->data()[55]; s7 += s*this->data()[63];
-    
-    if (!adding)
+  const number* e = this->data();
+  const unsigned int size_m = m(),
+                    size_n = n();
+  if (!adding)
     {
-      dst(0) = s0;
-      dst(1) = s1;
-      dst(2) = s2;
-      dst(3) = s3;
-      dst(4) = s4;
-      dst(5) = s5;
-      dst(6) = s6;
-      dst(7) = s7;
+      for (unsigned int i=0; i<size_m; ++i)
+       {
+         number2 s = 0.;
+         for (unsigned int j=0; j<size_n; ++j)
+           s += number2(src(j)) * number2(*(e++));
+         dst(i) = s;
+       };
     }
-    else
+  else
     {
-      dst(0) += s0;
-      dst(1) += s1;
-      dst(2) += s2;
-      dst(3) += s3;
-      dst(4) += s4;
-      dst(5) += s5;
-      dst(6) += s6;
-      dst(7) += s7;
+      for (unsigned int i=0; i<size_m; ++i)
+       {
+         number2 s = 0.;
+         for (unsigned int j=0; j<size_n; ++j)
+           s += number2(src(j)) * number2(*(e++));
+         dst(i) += s;
+       }
     }
-  }
-  else
-  {    
-    const number* e = this->data();
-    const unsigned int size_m = m(),
-                      size_n = n();
-    if (!adding)
-      {
-       for (unsigned int i=0; i<size_m; ++i)
-         {
-           number s = 0.;
-           for (unsigned int j=0; j<size_n; ++j)
-             s += number(src(j)) * *(e++);
-           dst(i) = s;
-         };
-      }
-    else
-      {
-       for (unsigned int i=0; i<size_m; ++i)
-         {
-           number s = 0.;
-           for (unsigned int j=0; j<size_n; ++j)
-             s += number(src(j)) * *(e++);
-           dst(i) += s;
-         }
-      }
-  }
 }
 
 
index 0361ce0a50d6117e112cf1673c5e6670f2c0b7d8..58bd9a7e845d133f6c8f049c3e8f1274664c5f99 100644 (file)
@@ -118,110 +118,110 @@ FullMatrix<TYPEMAT>::residual<TYPEVEC,TYPERES>(Vector<TYPEVEC>&,
 
 
 
-// #undef TYPEMAT
-
-// #define TYPEMAT std::complex<double>
-
-// template class FullMatrix<TYPEMAT>;
-
-// template void FullMatrix<TYPEMAT>::print(
-//   LogStream&, const unsigned int, const unsigned int) const;
-// template void FullMatrix<TYPEMAT>::print(
-//   std::ostream&, const unsigned int, const unsigned int) const;
-
-// template FullMatrix<TYPEMAT>& FullMatrix<TYPEMAT>::operator =(
-//   const FullMatrix<std::complex<float> >&);
-
-// #undef TYPEMAT2
-// #define TYPEMAT2 std::complex<double>
-
-// template void FullMatrix<TYPEMAT>::fill<TYPEMAT2> (
-//   const FullMatrix<TYPEMAT2>&, unsigned, unsigned, unsigned, unsigned);
-// template void FullMatrix<TYPEMAT>::add<TYPEMAT2> (const TYPEMAT, const FullMatrix<TYPEMAT2>&);
-// template void FullMatrix<TYPEMAT>::add<TYPEMAT2> (const TYPEMAT, const FullMatrix<TYPEMAT2>&,
-//                                               const TYPEMAT, const FullMatrix<TYPEMAT2>&);
-// template void FullMatrix<TYPEMAT>::add<TYPEMAT2> (const TYPEMAT, const FullMatrix<TYPEMAT2>&,
-//                                               const TYPEMAT, const FullMatrix<TYPEMAT2>&,
-//                                               const TYPEMAT, const FullMatrix<TYPEMAT2>&);
-// template void FullMatrix<TYPEMAT>::add<TYPEMAT2> (
-//   const FullMatrix<TYPEMAT2>&, std::complex<double>, unsigned, unsigned, unsigned, unsigned);
-// template void FullMatrix<TYPEMAT>::Tadd<TYPEMAT2> (const TYPEMAT, const FullMatrix<TYPEMAT2>&);
-// template void FullMatrix<TYPEMAT>::Tadd<TYPEMAT2> (
-//   const FullMatrix<TYPEMAT2>&, std::complex<double>, unsigned, unsigned, unsigned, unsigned);
-// template void FullMatrix<TYPEMAT>::equ<TYPEMAT2> (const TYPEMAT, const FullMatrix<TYPEMAT2>&);
-// template void FullMatrix<TYPEMAT>::equ<TYPEMAT2> (const TYPEMAT, const FullMatrix<TYPEMAT2>&,
-//                                               const TYPEMAT, const FullMatrix<TYPEMAT2>&);
-// template void FullMatrix<TYPEMAT>::equ<TYPEMAT2> (const TYPEMAT, const FullMatrix<TYPEMAT2>&,
-//                                               const TYPEMAT, const FullMatrix<TYPEMAT2>&,
-//                                               const TYPEMAT, const FullMatrix<TYPEMAT2>&);
-// template void FullMatrix<TYPEMAT>::mmult<TYPEMAT2> (FullMatrix<TYPEMAT2>&, const FullMatrix<TYPEMAT2>&, const bool) const;
-// template void FullMatrix<TYPEMAT>::Tmmult<TYPEMAT2> (FullMatrix<TYPEMAT2>&, const FullMatrix<TYPEMAT2>&, const bool) const;
-// template void FullMatrix<TYPEMAT>::add_diag<TYPEMAT2> (const TYPEMAT, const FullMatrix<TYPEMAT2>&);
-// template void FullMatrix<TYPEMAT>::invert<TYPEMAT2> (const FullMatrix<TYPEMAT2>&);
-
-
-// #undef TYPEVEC
-// #undef TYPERES
-// #define TYPEVEC std::complex<double>
-// #define TYPERES std::complex<double>
-
-// template void FullMatrix<TYPEMAT>::fill_permutation<TYPEVEC> (
-//   const FullMatrix<TYPEVEC>&,
-//   const std::vector<unsigned int>&,
-//   const std::vector<unsigned int>&);
-// template void FullMatrix<TYPEMAT>::vmult<TYPEVEC>(
-//   Vector<TYPEVEC>&, const Vector<TYPEVEC>&, bool) const;
-// template void FullMatrix<TYPEMAT>::Tvmult<TYPEVEC>(
-//   Vector<TYPEVEC>&, const Vector<TYPEVEC>&, bool) const;
-// template std::complex<double> FullMatrix<TYPEMAT>::residual<TYPEVEC>(
-//   Vector<TYPEVEC>&, const Vector<TYPEVEC>&, const Vector<TYPERES>&) const;
-// template TYPEVEC FullMatrix<TYPEMAT>::matrix_norm_square<TYPEVEC> (
-//   const Vector<TYPEVEC> &) const;
-// template TYPEVEC FullMatrix<TYPEMAT>::matrix_scalar_product<TYPEVEC>(
-//   const Vector<TYPEVEC>&, const Vector<TYPEVEC>&) const;
-// template void FullMatrix<TYPEMAT>::forward<TYPEVEC>(
-//   Vector<TYPEVEC>&, const Vector<TYPEVEC>&) const;
-// template void FullMatrix<TYPEMAT>::backward<TYPEVEC>(
-//   Vector<TYPEVEC>&, const Vector<TYPEVEC>&) const;
-
-// template
-// void FullMatrix<TYPEMAT>::precondition_Jacobi<TYPEVEC> (
-//   Vector<TYPEVEC> &, const Vector<TYPEVEC> &, const TYPEMAT) const;
-
-// #undef TYPEVEC
-// #define TYPEVEC std::complex<float>
-
-// template void FullMatrix<TYPEMAT>::fill_permutation<TYPEVEC> (
-//   const FullMatrix<TYPEVEC>&,
-//   const std::vector<unsigned int>&,
-//   const std::vector<unsigned int>&);
-// template void FullMatrix<TYPEMAT>::vmult<TYPEVEC>(
-//   Vector<TYPEVEC>&, const Vector<TYPEVEC>&, bool) const;
-// template void FullMatrix<TYPEMAT>::Tvmult<TYPEVEC>(
-//   Vector<TYPEVEC>&, const Vector<TYPEVEC>&, bool) const;
-// template std::complex<double> FullMatrix<TYPEMAT>::residual<TYPEVEC>(
-//   Vector<TYPEVEC>&, const Vector<TYPEVEC>&, const Vector<TYPERES>&) const;
-// template TYPEVEC FullMatrix<TYPEMAT>::matrix_norm_square<TYPEVEC> (
-//   const Vector<TYPEVEC> &) const;
-// template TYPEVEC FullMatrix<TYPEMAT>::matrix_scalar_product<TYPEVEC>(
-//   const Vector<TYPEVEC>&, const Vector<TYPEVEC>&) const;
-// template void FullMatrix<TYPEMAT>::forward<TYPEVEC>(
-//   Vector<TYPEVEC>&, const Vector<TYPEVEC>&) const;
-// template void FullMatrix<TYPEMAT>::backward<TYPEVEC>(
-//   Vector<TYPEVEC>&, const Vector<TYPEVEC>&) const;
-// template
-// void FullMatrix<TYPEMAT>::precondition_Jacobi<TYPEVEC> (
-//   Vector<TYPEVEC> &, const Vector<TYPEVEC> &, const TYPEMAT) const;
-
-
-// #undef TYPERES
-// #define TYPERES std::complex<float>
-
-// template
-// std::complex<double>
-// FullMatrix<TYPEMAT>::residual<TYPEVEC,TYPERES>(Vector<TYPEVEC>&,
-//                                            const Vector<TYPEVEC>&,
-//                                            const Vector<TYPERES>&) const;
+#undef TYPEMAT
+
+#define TYPEMAT std::complex<double>
+
+template class FullMatrix<TYPEMAT>;
+
+template void FullMatrix<TYPEMAT>::print(
+  LogStream&, const unsigned int, const unsigned int) const;
+template void FullMatrix<TYPEMAT>::print(
+  std::ostream&, const unsigned int, const unsigned int) const;
+
+template FullMatrix<TYPEMAT>& FullMatrix<TYPEMAT>::operator =(
+  const FullMatrix<std::complex<float> >&);
+
+#undef TYPEMAT2
+#define TYPEMAT2 std::complex<double>
+
+template void FullMatrix<TYPEMAT>::fill<TYPEMAT2> (
+  const FullMatrix<TYPEMAT2>&, unsigned, unsigned, unsigned, unsigned);
+template void FullMatrix<TYPEMAT>::add<TYPEMAT2> (const TYPEMAT, const FullMatrix<TYPEMAT2>&);
+template void FullMatrix<TYPEMAT>::add<TYPEMAT2> (const TYPEMAT, const FullMatrix<TYPEMAT2>&,
+                                                 const TYPEMAT, const FullMatrix<TYPEMAT2>&);
+template void FullMatrix<TYPEMAT>::add<TYPEMAT2> (const TYPEMAT, const FullMatrix<TYPEMAT2>&,
+                                                 const TYPEMAT, const FullMatrix<TYPEMAT2>&,
+                                                 const TYPEMAT, const FullMatrix<TYPEMAT2>&);
+template void FullMatrix<TYPEMAT>::add<TYPEMAT2> (
+  const FullMatrix<TYPEMAT2>&, std::complex<double>, unsigned, unsigned, unsigned, unsigned);
+template void FullMatrix<TYPEMAT>::Tadd<TYPEMAT2> (const TYPEMAT, const FullMatrix<TYPEMAT2>&);
+template void FullMatrix<TYPEMAT>::Tadd<TYPEMAT2> (
+  const FullMatrix<TYPEMAT2>&, std::complex<double>, unsigned, unsigned, unsigned, unsigned);
+template void FullMatrix<TYPEMAT>::equ<TYPEMAT2> (const TYPEMAT, const FullMatrix<TYPEMAT2>&);
+template void FullMatrix<TYPEMAT>::equ<TYPEMAT2> (const TYPEMAT, const FullMatrix<TYPEMAT2>&,
+                                                 const TYPEMAT, const FullMatrix<TYPEMAT2>&);
+template void FullMatrix<TYPEMAT>::equ<TYPEMAT2> (const TYPEMAT, const FullMatrix<TYPEMAT2>&,
+                                                 const TYPEMAT, const FullMatrix<TYPEMAT2>&,
+                                                 const TYPEMAT, const FullMatrix<TYPEMAT2>&);
+template void FullMatrix<TYPEMAT>::mmult<TYPEMAT2> (FullMatrix<TYPEMAT2>&, const FullMatrix<TYPEMAT2>&, const bool) const;
+template void FullMatrix<TYPEMAT>::Tmmult<TYPEMAT2> (FullMatrix<TYPEMAT2>&, const FullMatrix<TYPEMAT2>&, const bool) const;
+template void FullMatrix<TYPEMAT>::add_diag<TYPEMAT2> (const TYPEMAT, const FullMatrix<TYPEMAT2>&);
+template void FullMatrix<TYPEMAT>::invert<TYPEMAT2> (const FullMatrix<TYPEMAT2>&);
+
+
+#undef TYPEVEC
+#undef TYPERES
+#define TYPEVEC std::complex<double>
+#define TYPERES std::complex<double>
+
+template void FullMatrix<TYPEMAT>::fill_permutation<TYPEVEC> (
+  const FullMatrix<TYPEVEC>&,
+  const std::vector<unsigned int>&,
+  const std::vector<unsigned int>&);
+template void FullMatrix<TYPEMAT>::vmult<TYPEVEC>(
+  Vector<TYPEVEC>&, const Vector<TYPEVEC>&, bool) const;
+template void FullMatrix<TYPEMAT>::Tvmult<TYPEVEC>(
+  Vector<TYPEVEC>&, const Vector<TYPEVEC>&, bool) const;
+template std::complex<double> FullMatrix<TYPEMAT>::residual<TYPEVEC>(
+  Vector<TYPEVEC>&, const Vector<TYPEVEC>&, const Vector<TYPERES>&) const;
+template TYPEVEC FullMatrix<TYPEMAT>::matrix_norm_square<TYPEVEC> (
+  const Vector<TYPEVEC> &) const;
+template TYPEVEC FullMatrix<TYPEMAT>::matrix_scalar_product<TYPEVEC>(
+  const Vector<TYPEVEC>&, const Vector<TYPEVEC>&) const;
+template void FullMatrix<TYPEMAT>::forward<TYPEVEC>(
+  Vector<TYPEVEC>&, const Vector<TYPEVEC>&) const;
+template void FullMatrix<TYPEMAT>::backward<TYPEVEC>(
+  Vector<TYPEVEC>&, const Vector<TYPEVEC>&) const;
+
+template
+void FullMatrix<TYPEMAT>::precondition_Jacobi<TYPEVEC> (
+  Vector<TYPEVEC> &, const Vector<TYPEVEC> &, const TYPEMAT) const;
+
+#undef TYPEVEC
+#define TYPEVEC std::complex<float>
+
+template void FullMatrix<TYPEMAT>::fill_permutation<TYPEVEC> (
+  const FullMatrix<TYPEVEC>&,
+  const std::vector<unsigned int>&,
+  const std::vector<unsigned int>&);
+template void FullMatrix<TYPEMAT>::vmult<TYPEVEC>(
+  Vector<TYPEVEC>&, const Vector<TYPEVEC>&, bool) const;
+template void FullMatrix<TYPEMAT>::Tvmult<TYPEVEC>(
+  Vector<TYPEVEC>&, const Vector<TYPEVEC>&, bool) const;
+template std::complex<double> FullMatrix<TYPEMAT>::residual<TYPEVEC>(
+  Vector<TYPEVEC>&, const Vector<TYPEVEC>&, const Vector<TYPERES>&) const;
+template TYPEVEC FullMatrix<TYPEMAT>::matrix_norm_square<TYPEVEC> (
+  const Vector<TYPEVEC> &) const;
+template TYPEVEC FullMatrix<TYPEMAT>::matrix_scalar_product<TYPEVEC>(
+  const Vector<TYPEVEC>&, const Vector<TYPEVEC>&) const;
+template void FullMatrix<TYPEMAT>::forward<TYPEVEC>(
+  Vector<TYPEVEC>&, const Vector<TYPEVEC>&) const;
+template void FullMatrix<TYPEMAT>::backward<TYPEVEC>(
+  Vector<TYPEVEC>&, const Vector<TYPEVEC>&) const;
+template
+void FullMatrix<TYPEMAT>::precondition_Jacobi<TYPEVEC> (
+  Vector<TYPEVEC> &, const Vector<TYPEVEC> &, const TYPEMAT) const;
+
+
+#undef TYPERES
+#define TYPERES std::complex<float>
+
+template
+std::complex<double>
+FullMatrix<TYPEMAT>::residual<TYPEVEC,TYPERES>(Vector<TYPEVEC>&,
+                                              const Vector<TYPEVEC>&,
+                                              const Vector<TYPERES>&) const;
 
 
 DEAL_II_NAMESPACE_CLOSE

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.