]> https://gitweb.dealii.org/ - dealii.git/commitdiff
introduce templates for LAPACK functions
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Mon, 8 Nov 2010 00:19:59 +0000 (00:19 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Mon, 8 Nov 2010 00:19:59 +0000 (00:19 +0000)
git-svn-id: https://svn.dealii.org/trunk@22631 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/common/scripts/lapack_templates.pl
deal.II/include/deal.II/lac/full_matrix.templates.h
deal.II/include/deal.II/lac/lapack_full_matrix.h
deal.II/include/deal.II/lac/precondition_block.templates.h
deal.II/include/deal.II/lac/precondition_block_base.h
deal.II/include/deal.II/lac/relaxation_block.templates.h
deal.II/source/lac/lapack_full_matrix.cc

index 2fd0419614f228f73d819fc172a156934c5c104e..886e9f9ee552137882adbbf1bbd47dc27a8760a5 100644 (file)
@@ -105,7 +105,29 @@ while(<>)
        $args0 = $args;
        $args0 =~ s/\*[^,]*,/\*,/g;
        $args0 =~ s/\*[^,]*$/\*/g;
+       # First, do the general template None of these functions is
+       # implemented, but they allow us to link for instance with
+       # long double lapack support
+       my $numbers = 1;
+       my $argst = $args0;
+       my $typet = $type;
+       while ($argst =~ s/double/number$numbers/)
+       {
+           $numbers++;
+       }
+       $typet =~ s/double/number1/g;
        
+       $templates .= "\n\n//--------------------------------------------------------------------------------//\ntemplate<typename number1";
+       for (my $i=2;$i<$numbers;++$i)
+       {
+           $templates .= ", typename number$i";
+       }
+       $templates .= ">\ninline $typet\n$name ($argst)\n";
+       $templates .= "{\n  Assert (false, ExcNotImplemented());\n}";
+
+       # The specialization for double. Note that we always have two
+       # versions, one implemented and calling LAPACK, the other not
+       # implemented
        $templates .= "\n\n#ifdef HAVE_D$capname\_";
        $templates .= "\ninline $type\n$name ($args)\n{\n  d$name\_ ($args2);\n}\n";
        $templates .= "#else\ninline $type\n$name ($args0)\n";
index 7db3267f9f6cf81308d9276889f49b9cfb3ffb46..4095da48260b7635a4d4660a264bfa2a7e6806fc 100644 (file)
@@ -527,37 +527,6 @@ FullMatrix<number>::equ (const number               a,
 
 
 
-namespace internal
-{
-                                  // a namespace into which we import
-                                  // the global definitions of the
-                                  // LAPACK functions getrf for
-                                  // various data types and add
-                                  // general templates for all other
-                                  // types. this allows to call the
-                                  // getrf function for all data
-                                  // types but generated an exception
-                                  // if something is called for a
-                                  // data type not supported by
-                                  // LAPACK.
-  namespace gemm_switch
-  {
-    using ::dealii::gemm;
-
-    template <typename T, typename T2>
-    void
-    gemm (const char*, const char*,  const int*, const int*, const int*,
-         const T*, const T2*, const int*, const T*, const int*,
-         const T*, T2*, const int*)
-    {
-      Assert (false, LAPACKSupport::ExcMissing("dgemm for this data type"));
-    }
-    
-  }
-}
-
-
-
 template <typename number>
 template <typename number2>
 void FullMatrix<number>::mmult (FullMatrix<number2>       &dst,
@@ -616,8 +585,8 @@ void FullMatrix<number>::mmult (FullMatrix<number2>       &dst,
                                       // Use the BLAS function gemm for
                                       // calculating the matrix-matrix
                                       // product.
-      internal::gemm_switch::gemm(notrans, notrans, &m, &n, &k, &alpha, &src(0,0), &m,
-                                 &this->values[0], &k, &beta, &dst(0,0), &m);
+      gemm(notrans, notrans, &m, &n, &k, &alpha, &src(0,0), &m,
+          &this->values[0], &k, &beta, &dst(0,0), &m);
 
       return;
     }
@@ -704,8 +673,8 @@ void FullMatrix<number>::Tmmult (FullMatrix<number2>       &dst,
                                       // Use the BLAS function gemm for
                                       // calculating the matrix-matrix
                                       // product.
-      internal::gemm_switch::gemm(notrans, trans, &m, &n, &k, &alpha, &src(0,0), &m,
-                                 &this->values[0], &n, &beta, &dst(0,0), &m);
+      gemm(notrans, trans, &m, &n, &k, &alpha, &src(0,0), &m,
+          &this->values[0], &n, &beta, &dst(0,0), &m);
 
       return;
     }
@@ -795,8 +764,8 @@ void FullMatrix<number>::mTmult (FullMatrix<number2>       &dst,
                                       // Use the BLAS function gemm for
                                       // calculating the matrix-matrix
                                       // product.
-      internal::gemm_switch::gemm(trans, notrans, &m, &n, &k, &alpha, &src(0,0), &k,
-                                 &this->values[0], &k, &beta, &dst(0,0), &m);
+      gemm(trans, notrans, &m, &n, &k, &alpha, &src(0,0), &k,
+          &this->values[0], &k, &beta, &dst(0,0), &m);
 
       return;
     }
@@ -880,8 +849,8 @@ void FullMatrix<number>::TmTmult (FullMatrix<number2>       &dst,
                                       // Use the BLAS function gemm for
                                       // calculating the matrix-matrix
                                       // product.
-      internal::gemm_switch::gemm(trans, trans, &m, &n, &k, &alpha, &src(0,0), &k,
-                                 &this->values[0], &n, &beta, &dst(0,0), &m);
+      gemm(trans, trans, &m, &n, &k, &alpha, &src(0,0), &k,
+          &this->values[0], &n, &beta, &dst(0,0), &m);
 
       return;
     }
@@ -1714,44 +1683,6 @@ FullMatrix<number>::print_formatted (
 }
 
 
-namespace internal
-{
-                                  // a namespace into which we import
-                                  // the global definitions of the
-                                  // LAPACK functions getrf for
-                                  // various data types and add
-                                  // general templates for all other
-                                  // types. this allows to call the
-                                  // getrf function for all data
-                                  // types but generated an exception
-                                  // if something is called for a
-                                  // data type not supported by
-                                  // LAPACK.
-  namespace getrf_switch
-  {
-    using ::dealii::getrf;
-
-    template <typename T>
-    void
-    getrf (const int*, const int*, T*, const int*, int*, int*)
-    {
-      Assert (false, LAPACKSupport::ExcMissing("dgetrf for this data type"));
-    }
-    
-
-    using ::dealii::getri;
-
-    template <typename T>
-    void
-    getri (const int*, T*, const int*, int*, T*, const int*, int*)
-    {
-      Assert (false, LAPACKSupport::ExcMissing("dgetri for this data type"));
-    }    
-  }
-}
-
-
-
 template <typename number>
 void
 FullMatrix<number>::gauss_jordan ()
@@ -1798,7 +1729,7 @@ FullMatrix<number>::gauss_jordan ()
 
                                       // Use the LAPACK function getrf for 
                                       // calculating the LU factorization.
-      internal::getrf_switch::getrf(&nn, &nn, &this->values[0], &nn, &ipiv[0], &info);
+      getrf(&nn, &nn, &this->values[0], &nn, &ipiv[0], &info);
 
       Assert(info >= 0, ExcInternalError());
       Assert(info == 0, LACExceptions::ExcSingular());
@@ -1809,7 +1740,7 @@ FullMatrix<number>::gauss_jordan ()
                                       // Use the LAPACK function getri for
                                       // calculating the actual inverse using
                                       // the LU factorization.
-      internal::getrf_switch::getri(&nn, &this->values[0], &nn, &ipiv[0], &inv_work[0], &nn, &info);
+      getri(&nn, &this->values[0], &nn, &ipiv[0], &inv_work[0], &nn, &info);
 
       Assert(info >= 0, ExcInternalError());
       Assert(info == 0, LACExceptions::ExcSingular());
index d80c651a2d230427dc8a2385f9ae5d645b4b0e5d..dfafccc881681cd28a5f2d5e0fa06317264168a0 100644 (file)
@@ -161,6 +161,7 @@ class LAPACKFullMatrix : public TransposeTable<number>
               const number factor = 1.,
               const bool transpose = false);
 
+
                                     /**
                                      * Matrix-vector-multiplication.
                                      *
@@ -177,19 +178,30 @@ class LAPACKFullMatrix : public TransposeTable<number>
                                       *
                                       * Source and destination must
                                       * not be the same vector.
+                                     *
+                                     * @note This template only
+                                     * exists for compile-time
+                                     * compatibility with
+                                     * FullMatrix. Implementation is
+                                     * only available for <tt>VECTOR=Vector&lt;number&gt;</tt>
                                      */
-    void vmult (Vector<number>   &w,
-               const Vector<number> &v,
-               const bool            adding=false) const;
+    template <class VECTOR>
+    void vmult(VECTOR& dst, const VECTOR& src, const bool adding = false) const;
                                     /**
                                      * Adding Matrix-vector-multiplication.
                                      *  <i>w += A*v</i>
                                       *
                                       * Source and destination must
                                       * not be the same vector.
+                                     *
+                                     * @note This template only
+                                     * exists for compile-time
+                                     * compatibility with
+                                     * FullMatrix. Implementation is
+                                     * only available for <tt>VECTOR=Vector&lt;number&gt;</tt>
                                      */
-    void vmult_add (Vector<number>       &w,
-                   const Vector<number> &v) const;
+    template <class VECTOR>
+    void vmult_add (VECTOR& w, const VECTOR& v) const;
 
                                     /**
                                      * Transpose
@@ -209,9 +221,15 @@ class LAPACKFullMatrix : public TransposeTable<number>
                                       *
                                       * Source and destination must
                                       * not be the same vector.
-                                     */
-    void Tvmult (Vector<number>       &w,
-                const Vector<number> &v,
+                                     *
+                                     * @note This template only
+                                     * exists for compile-time
+                                     * compatibility with
+                                     * FullMatrix. Implementation is
+                                     * only available for <tt>VECTOR=Vector&lt;number&gt;</tt>
+                                     */
+    template <class VECTOR>
+    void Tvmult (VECTOR& w, const VECTOR& v,
                 const bool            adding=false) const;
 
                                     /**
@@ -221,10 +239,26 @@ class LAPACKFullMatrix : public TransposeTable<number>
                                       *
                                       * Source and destination must
                                       * not be the same vector.
-                                     */
+                                     *
+                                     * @note This template only
+                                     * exists for compile-time
+                                     * compatibility with
+                                     * FullMatrix. Implementation is
+                                     * only available for <tt>VECTOR=Vector&lt;number&gt;</tt>
+                                     */
+    template <class VECTOR>
+    void Tvmult_add (VECTOR& w, const VECTOR& v) const;
+    
+    void vmult (Vector<number>   &w,
+               const Vector<number> &v,
+               const bool            adding=false) const;
+    void vmult_add (Vector<number>       &w,
+                   const Vector<number> &v) const;
+    void Tvmult (Vector<number>       &w,
+                const Vector<number> &v,
+                const bool            adding=false) const;
     void Tvmult_add (Vector<number>       &w,
                     const Vector<number> &v) const;
-
                                     /**
                                      * Compute the LU factorization
                                      * of the matrix using LAPACK
@@ -576,7 +610,7 @@ class PreconditionLU
 
 template <typename number>
 template <class MATRIX>
-void
+inline void
 LAPACKFullMatrix<number>::copy_from (const MATRIX& M)
 {
   this->reinit (M.m(), M.n());
@@ -592,7 +626,7 @@ LAPACKFullMatrix<number>::copy_from (const MATRIX& M)
 
 template <typename number>
 template <class MATRIX>
-void
+inline void
 LAPACKFullMatrix<number>::fill (
   const MATRIX& M,
   const unsigned int dst_offset_i,
@@ -620,7 +654,43 @@ LAPACKFullMatrix<number>::fill (
 
 
 template <typename number>
-std::complex<number>
+template <class VECTOR>
+inline void
+LAPACKFullMatrix<number>::vmult(VECTOR&, const VECTOR&, bool) const
+{
+  Assert(false, ExcNotImplemented());
+}
+
+
+template <typename number>
+template <class VECTOR>
+inline void
+LAPACKFullMatrix<number>::vmult_add(VECTOR&, const VECTOR&) const
+{
+  Assert(false, ExcNotImplemented());
+}
+
+
+template <typename number>
+template <class VECTOR>
+inline void
+LAPACKFullMatrix<number>::Tvmult(VECTOR&, const VECTOR&, bool) const
+{
+  Assert(false, ExcNotImplemented());
+}
+
+
+template <typename number>
+template <class VECTOR>
+inline void
+LAPACKFullMatrix<number>::Tvmult_add(VECTOR&, const VECTOR&) const
+{
+  Assert(false, ExcNotImplemented());
+}
+
+
+template <typename number>
+inline std::complex<number>
 LAPACKFullMatrix<number>::eigenvalue (const unsigned int i) const
 {
   Assert (state & LAPACKSupport::eigenvalues, ExcInvalidState());
@@ -633,7 +703,7 @@ LAPACKFullMatrix<number>::eigenvalue (const unsigned int i) const
 
 
 template <typename number>
-number
+inline number
 LAPACKFullMatrix<number>::singular_value (const unsigned int i) const
 {
   Assert (state == LAPACKSupport::svd || state == LAPACKSupport::inverse_svd, LAPACKSupport::ExcState(state));
index 76eb4cacd3630cf19c34114b049e7809ba2a3549..4e28fa9d23106cdfda5abb66b9426d60792a5df2 100644 (file)
@@ -138,6 +138,10 @@ void PreconditionBlock<MATRIX,inverse_type>::invert_permuted_diagblocks(
          case PreconditionBlockBase<inverse_type>::householder:
                this->inverse_householder(0).initialize(M_cell);
                break;
+         case PreconditionBlockBase<inverse_type>::svd:
+               this->inverse_svd(0) = M_cell;
+               this->inverse_svd(0).compute_inverse_svd(0.);
+               break;
          default:
                Assert(false, ExcNotImplemented());
 
@@ -189,6 +193,10 @@ void PreconditionBlock<MATRIX,inverse_type>::invert_permuted_diagblocks(
              case PreconditionBlockBase<inverse_type>::householder:
                    this->inverse_householder(cell).initialize(M_cell);
                    break;
+         case PreconditionBlockBase<inverse_type>::svd:
+               this->inverse_svd(cell) = M_cell;
+               this->inverse_svd(cell).compute_inverse_svd(0.);
+               break;
              default:
                    Assert(false, ExcNotImplemented());
                    
index d8bc501f1bce0939009965a4ecac20125bc87d3b..0b15a13d436af284411ffc7a6f576256bf2496b4 100644 (file)
@@ -20,6 +20,7 @@
 #include <base/smartpointer.h>
 #include <base/memory_consumption.h>
 #include <lac/householder.h>
+#include <lac/lapack_full_matrix.h>
 
 #include <vector>
 
@@ -68,7 +69,12 @@ class PreconditionBlockBase
                                            * Use QR decomposition of
                                            * the Householder class.
                                            */
-         householder
+         householder,
+                                          /**
+                                           * Use the singular value
+                                           * decomposition of LAPACKFullMatrix.
+                                           */
+         svd
     };
     
                                     /**
@@ -187,6 +193,12 @@ class PreconditionBlockBase
                                      */
     Householder<number>& inverse_householder (unsigned int i);
     
+                                    /**
+                                     * Access to the inverse diagonal
+                                     * blocks if Inversion is #householder.
+                                     */
+    LAPACKFullMatrix<number>& inverse_svd (unsigned int i);
+    
                                     /**
                                      * Access to the inverse diagonal
                                      * blocks.
@@ -243,7 +255,8 @@ class PreconditionBlockBase
                                      * #gauss_jordan is used. Using
                                      * <tt>number=float</tt> saves
                                      * memory in comparison with
-                                     * <tt>number=double</tt>.
+                                     * <tt>number=double</tt>, but
+                                     * may introduce numerical instability.
                                      */
     std::vector<FullMatrix<number> > var_inverse_full;
 
@@ -251,15 +264,30 @@ class PreconditionBlockBase
                                      * Storage of the inverse
                                      * matrices of the diagonal
                                      * blocks matrices as
-                                     * <tt>Householder<number></tt>
+                                     * <tt>Householder</tt>
                                      * matrices if Inversion
                                      * #householder is used. Using
                                      * <tt>number=float</tt> saves
                                      * memory in comparison with
-                                     * <tt>number=double</tt>.
+                                     * <tt>number=double</tt>, but
+                                     * may introduce numerical instability.
                                      */
     std::vector<Householder<number> > var_inverse_householder;
 
+                                    /**
+                                     * Storage of the inverse
+                                     * matrices of the diagonal
+                                     * blocks matrices as
+                                     * <tt>LAPACKFullMatrix</tt>
+                                     * matrices if Inversion
+                                     * #svd is used. Using
+                                     * <tt>number=float</tt> saves
+                                     * memory in comparison with
+                                     * <tt>number=double</tt>, but
+                                     * may introduce numerical instability.
+                                     */
+    std::vector<LAPACKFullMatrix<number> > var_inverse_svd;
+    
                                     /**
                                      * Storage of the original diagonal blocks.
                                      *
@@ -310,8 +338,10 @@ PreconditionBlockBase<number>::clear()
 {
   if (var_inverse_full.size()!=0)
     var_inverse_full.erase(var_inverse_full.begin(), var_inverse_full.end());
-  if (var_inverse_full.size()!=0)
+  if (var_inverse_householder.size()!=0)
     var_inverse_householder.erase(var_inverse_householder.begin(), var_inverse_householder.end());
+  if (var_inverse_svd.size()!=0)
+    var_inverse_svd.erase(var_inverse_svd.begin(), var_inverse_svd.end());
   if (var_diagonal.size()!=0)
     var_diagonal.erase(var_diagonal.begin(), var_diagonal.end());
   var_same_diagonal = false;
@@ -341,6 +371,10 @@ Inversion method)
          case householder:
                var_inverse_householder.resize(1);
                break;
+         case svd:
+               var_inverse_svd.resize(1);
+               var_inverse_svd[0].reinit(b,b);
+               break;
          default:
                Assert(false, ExcNotImplemented());
        }
@@ -380,6 +414,13 @@ Inversion method)
          case householder:
                var_inverse_householder.resize(n);
                break;
+         case svd:
+         {
+           std::vector<LAPACKFullMatrix<number> >
+             tmp(n, LAPACKFullMatrix<number>(b));
+           var_inverse_svd.swap (tmp);
+           break;
+         }
          default:      
                Assert(false, ExcNotImplemented());     
        }
@@ -414,6 +455,10 @@ PreconditionBlockBase<number>::inverse_vmult(
            AssertIndexRange (ii, var_inverse_householder.size());
            var_inverse_householder[ii].vmult(dst, src);
            break;
+      case svd:
+           AssertIndexRange (ii, var_inverse_svd.size());
+           var_inverse_svd[ii].vmult(dst, src);
+           break;
       default:
            Assert(false, ExcNotImplemented());
     }
@@ -439,6 +484,10 @@ PreconditionBlockBase<number>::inverse_Tvmult(
            AssertIndexRange (ii, var_inverse_householder.size());
            var_inverse_householder[ii].Tvmult(dst, src);
            break;
+      case svd:
+           AssertIndexRange (ii, var_inverse_svd.size());
+           var_inverse_svd[ii].Tvmult(dst, src);
+           break;
       default:
            Assert(false, ExcNotImplemented());
     }
@@ -499,6 +548,19 @@ PreconditionBlockBase<number>::inverse_householder(unsigned int i)
 }
 
 
+template <typename number>
+inline
+LAPACKFullMatrix<number>&
+PreconditionBlockBase<number>::inverse_svd(unsigned int i)
+{
+  if (same_diagonal())
+    return var_inverse_svd[0];
+  
+  AssertIndexRange (i, var_inverse_svd.size());
+  return var_inverse_svd[i];
+}
+
+
 template <typename number>
 inline
 FullMatrix<number>&
index 003ea60b6d9062439d9c3b03813c1e16b71667f6..6298b4bfb3687c6720b47feee36b4ce2a7d972a8 100644 (file)
@@ -132,6 +132,11 @@ RelaxationBlock<MATRIX,inverse_type>::invert_diagblocks ()
              case PreconditionBlockBase<inverse_type>::householder:
                    this->inverse_householder(block).initialize(M_cell);
                    break;
+             case PreconditionBlockBase<inverse_type>::svd:
+                   this->inverse_svd(block).reinit(bs, bs);
+                   this->inverse_svd(block) = M_cell;
+                   this->inverse_svd(block).compute_inverse_svd(0.);
+                   break;
              default:
                    Assert(false, ExcNotImplemented());
            }
index 407a63e01495b37c3b8c2a5ffc135a91581d7a57..a8ed9c210d0f9d088e3f1664ca17fc7d0947bd03 100644 (file)
@@ -683,6 +683,14 @@ LAPACKFullMatrix<float>::operator = (const FullMatrix<double>& M);
 template LAPACKFullMatrix<float> &
 LAPACKFullMatrix<float>::operator = (const FullMatrix<float>& M);
 
+template class LAPACKFullMatrix<long double>;
+template LAPACKFullMatrix<long double> &
+LAPACKFullMatrix<long double>::operator = (const FullMatrix<long double>& M);
+template LAPACKFullMatrix<long double> &
+LAPACKFullMatrix<long double>::operator = (const FullMatrix<double>& M);
+template LAPACKFullMatrix<long double> &
+LAPACKFullMatrix<long double>::operator = (const FullMatrix<float>& M);
+
 template class PreconditionLU<double>;
 template class PreconditionLU<float>;
 

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.