]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Always insert something into the diagonal with distribute_local_to_global. Use BLAS...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 31 Mar 2009 15:34:35 +0000 (15:34 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 31 Mar 2009 15:34:35 +0000 (15:34 +0000)
git-svn-id: https://svn.dealii.org/trunk@18539 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/base/include/base/config.h.in
deal.II/configure
deal.II/configure.in
deal.II/lac/include/lac/constraint_matrix.templates.h
deal.II/lac/include/lac/full_matrix.templates.h

index 97b4bb871d5217746d35ba561961b959a2f8f71f..6263f638740c05f2df6df30686e274a8a3f097e0 100644 (file)
 /* Define to 1 if you have the `dgeev_' function. */
 #undef HAVE_DGEEV_
 
+/* Define to 1 if you have the `dgemm_' function. */
+#undef HAVE_DGEMM_
+
 /* Define to 1 if you have the `dgemv_' function. */
 #undef HAVE_DGEMV_
 
 /* Define to 1 if you have the `sgeev_' function. */
 #undef HAVE_SGEEV_
 
+/* Define to 1 if you have the `sgemm_' function. */
+#undef HAVE_SGEMM_
+
 /* Define to 1 if you have the `sgemv_' function. */
 #undef HAVE_SGEMV_
 
index e397ac630777f2adaf7befc1945a2966c02eb285..783843e852a4127eb6068872b212d821b961f62b 100755 (executable)
@@ -16109,7 +16109,9 @@ done
 
 
 
-for ac_func in dgesvd_ sgesvd_ dgetrf_ sgetrf_ dgetri_ sgetri_
+
+
+for ac_func in dgesvd_ sgesvd_ dgemm_ sgemm_ dgetrf_ sgetrf_ dgetri_ sgetri_
 do
 as_ac_var=`echo "ac_cv_func_$ac_func" | $as_tr_sh`
 { echo "$as_me:$LINENO: checking for $ac_func" >&5
index 59135b6ef58a3786c59866653548d4ac90e88e7a..aa9dccfa0f4c7b9aa74a78b28996a2fd60f1a604 100644 (file)
@@ -551,7 +551,7 @@ if test "x$with_blas" != "x" -a "x$with_blas" != "xno" ; then
 fi
 
 AC_CHECK_FUNCS([daxpy_ saxpy_ dgemv_ sgemv_ dgeev_ sgeev_ dgeevx_ sgeevx_])
-AC_CHECK_FUNCS([dgesvd_ sgesvd_ dgetrf_ sgetrf_ dgetri_ sgetri_])
+AC_CHECK_FUNCS([dgesvd_ sgesvd_ dgemm_ sgemm_ dgetrf_ sgetrf_ dgetri_ sgetri_])
 AC_CHECK_FUNCS([dgetrs_ sgetrs_ dstev_ sstev_])
 
 dnl -------------------------------------------------------------
index 5fb76095ed73477a92e636f8a42cbf27b4c6815c..56de0b588bf9d0397901b199f860de228ada1348 100644 (file)
@@ -1048,7 +1048,7 @@ namespace internals
   list_shellsort (std::vector<std::pair<unsigned int,distributing> > &my_indices)
   {
                                   // now sort the actual dofs using a shell
-                                  // sort (which is very fast in case the
+                                  // sort which is very fast in case the
                                   // indices are already sorted (which is
                                   // the usual case with DG elements)
     unsigned int i, j, j2, temp, templ, istep;
@@ -1121,6 +1121,11 @@ distribute_local_to_global (const FullMatrix<double>        &local_matrix,
 
   const unsigned int n_local_dofs = local_dof_indices.size();
 
+  double average_diagonal = 0;
+  for (unsigned int i=0; i<n_local_dofs; ++i)
+    average_diagonal += std::fabs (local_matrix(i,i));
+  average_diagonal /= n_local_dofs;
+
                                   // when distributing the local data to
                                   // the global matrix, we can quite
                                   // cheaply sort the indices (obviously,
@@ -1234,9 +1239,10 @@ distribute_local_to_global (const FullMatrix<double>        &local_matrix,
                                   // hanging nodes in 3d). however, in the
                                   // line below, we do actually do
                                   // something with this dof
+      const double new_diagonal = std::fabs(local_matrix(local_row,local_row)) != 0 ?
+       std::fabs(local_matrix(local_row,local_row)) : average_diagonal;
       Threads::ThreadMutex::ScopedLock lock(mutex);
-      global_matrix.add(global_row,global_row,
-                       std::fabs(local_matrix(local_row,local_row)));
+      global_matrix.add(global_row, global_row, new_diagonal);
     }
 
   const unsigned int n_actual_dofs = my_indices.size();
@@ -1460,6 +1466,11 @@ distribute_local_to_global (const FullMatrix<double>        &local_matrix,
   const unsigned int n_local_dofs = local_dof_indices.size();
   const unsigned int num_blocks   = global_matrix.n_block_rows();
 
+  double average_diagonal = 0;
+  for (unsigned int i=0; i<n_local_dofs; ++i)
+    average_diagonal += std::fabs (local_matrix(i,i));
+  average_diagonal /= n_local_dofs;
+
   std::vector<std::pair<unsigned int,internals::distributing> > my_indices (n_local_dofs);
   std::vector<std::pair<unsigned int, const ConstraintLine *> > constraint_lines;
   constraint_lines.reserve(n_local_dofs);
@@ -1512,9 +1523,10 @@ distribute_local_to_global (const FullMatrix<double>        &local_matrix,
                                  (local_row, position->entries[q].second));
        }
 
+      const double new_diagonal = std::fabs(local_matrix(local_row,local_row)) != 0 ?
+       std::fabs(local_matrix(local_row,local_row)) : average_diagonal;
       Threads::ThreadMutex::ScopedLock lock(mutex);
-      global_matrix.add(global_row,global_row,
-                       std::fabs(local_matrix(local_row,local_row)));
+      global_matrix.add(global_row, global_row, new_diagonal);
     }
 
   const unsigned int n_actual_dofs = my_indices.size();
index 9598f0fbafb7f593f998e4e1a8b02dd59a9c8100..9efb96c49cb48f4a4c728a8b0f38d57263deb941 100644 (file)
@@ -545,6 +545,37 @@ 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 ::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,
@@ -556,6 +587,60 @@ void FullMatrix<number>::mmult (FullMatrix<number2>       &dst,
   Assert (dst.n() == src.n(), ExcDimensionMismatch(dst.n(), src.n()));
   Assert (dst.m() == m(), ExcDimensionMismatch(m(), dst.m()));
 
+                                  // see if we can use Lapack algorithms
+                                  // for this and if the type for 'number'
+                                  // works for us (it is usually not
+                                  // efficient to use Lapack for very small
+                                  // matrices):
+#if defined(HAVE_DGEMM_) && defined (HAVE_SGEMM_)
+  if ((types_are_equal<number,double>::value
+       ||
+       types_are_equal<number,float>::value)
+      &&
+      types_are_equal<number,number2>::value)
+    if (this->n_cols() > 25)
+    {
+                                      // In case we have the LAPACK
+                                      // function gemm detected at
+                                      // configure, we use that algorithms
+                                      // for matrix-matrix multiplication
+                                      // since it provides better
+                                      // performance than the deal.II
+                                      // native functions (it uses cache
+                                      // and register blocking in order to
+                                      // access local data).
+                                      //
+                                      // Note that BLAS/LAPACK stores
+                                      // matrix elements column-wise (i.e.,
+                                      // all values in one column, then all
+                                      // in the next, etc.), whereas the
+                                      // FullMatrix stores them row-wise.
+                                      // We ignore that difference, and
+                                      // give our row-wise data to LAPACK,
+                                      // let LAPACK build the product of
+                                      // transpose matrices, and read the
+                                      // result as if it were row-wise
+                                      // again. In other words, we calculate 
+                                      // (B^T A^T)^T, which is AB.
+
+      const int m = src.n();
+      const int n = this->m();
+      const int k = this->n();
+      const char *notrans = "n";
+
+      const number alpha = 1.;
+      const number beta = (adding == true) ? 1. : 0.;
+
+                                      // Use the LAPACK function getrf for 
+                                      // calculating the LU factorization.
+      internal::gemm_switch::gemm(notrans, notrans, &m, &n, &k, &alpha, &src(0,0), &m,
+                                 this->val, &k, &beta, &dst(0,0), &m);
+
+      return;
+    }
+
+#endif
+
   if (!adding)
     for (unsigned int i=0; i<m(); i++)
       for (unsigned int j=0; j<src.n(); j++)
@@ -589,6 +674,63 @@ void FullMatrix<number>::Tmmult (FullMatrix<number2>       &dst,
   Assert (n() == dst.m(), ExcDimensionMismatch(n(), dst.m()));
   Assert (src.n() == dst.n(), ExcDimensionMismatch(src.n(), dst.n()));
 
+
+                                  // see if we can use Lapack algorithms
+                                  // for this and if the type for 'number'
+                                  // works for us (it is usually not
+                                  // efficient to use Lapack for very small
+                                  // matrices):
+#if defined(HAVE_DGEMM_) && defined (HAVE_SGEMM_)
+  if ((types_are_equal<number,double>::value
+       ||
+       types_are_equal<number,float>::value)
+      &&
+      types_are_equal<number,number2>::value)
+    if (this->n_cols() > 25)
+    {
+                                      // In case we have the LAPACK
+                                      // function gemm detected at
+                                      // configure, we use that algorithms
+                                      // for matrix-matrix multiplication
+                                      // since it provides better
+                                      // performance than the deal.II
+                                      // native functions (it uses cache
+                                      // and register blocking in order to
+                                      // access local data).
+                                      //
+                                      // Note that BLAS/LAPACK stores
+                                      // matrix elements column-wise (i.e.,
+                                      // all values in one column, then all
+                                      // in the next, etc.), whereas the
+                                      // FullMatrix stores them row-wise.
+                                      // We ignore that difference, and
+                                      // give our row-wise data to LAPACK,
+                                      // let LAPACK build the product of
+                                      // transpose matrices, and read the
+                                      // result as if it were row-wise
+                                      // again. In other words, we calculate 
+                                      // (B^T A^T)^T, which is AB.
+
+      const int m = src.n();
+      const int n = this->n();
+      const int k = this->m();
+      const char *trans = "t";
+      const char *notrans = "n";
+
+      const number alpha = 1.;
+      const number beta = (adding == true) ? 1. : 0.;
+
+                                      // Use the LAPACK function getrf for 
+                                      // calculating the LU factorization.
+      internal::gemm_switch::gemm(notrans, trans, &m, &n, &k, &alpha, &src(0,0), &m,
+                                 this->val, &n, &beta, &dst(0,0), &m);
+
+      return;
+    }
+
+#endif
+
+
   if (!adding)
     for (unsigned int i=0; i<n(); i++)
       for (unsigned int j=0; j<src.n(); j++)
@@ -1392,14 +1534,17 @@ FullMatrix<number>::gauss_jordan ()
   Assert (!this->empty(), ExcEmptyMatrix());  
   Assert (this->n_cols() == this->n_rows(), ExcNotQuadratic());
 
-                                  // see if we can use Lapack
-                                  // algorithms for this and if the
-                                  // type for 'number' works for us:
+                                  // see if we can use Lapack algorithms
+                                  // for this and if the type for 'number'
+                                  // works for us (it is usually not
+                                  // efficient to use Lapack for very small
+                                  // matrices):
 #if defined(HAVE_DGETRF_) && defined (HAVE_SGETRF_) && \
     defined(HAVE_DGETRI_) && defined (HAVE_SGETRI_)
   if (types_are_equal<number,double>::value
       ||
       types_are_equal<number,float>::value)
+    if (this->n_cols() > 15)
     {
                                       // In case we have the LAPACK functions 
                                       // getrf and getri detected at configure, 

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.