]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Merge in the changes made for multithreaded linear algebra into the main branch....
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 10 Sep 1999 09:15:22 +0000 (09:15 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 10 Sep 1999 09:15:22 +0000 (09:15 +0000)
git-svn-id: https://svn.dealii.org/trunk@1736 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/Todo
deal.II/lac/include/lac/sparse_matrix.h
deal.II/lac/include/lac/sparse_matrix.templates.h

index fa889abe22465dd2abeda2bcee64cca99d8f7767..77a2b39e12e2c83ade2f9156438fd02f38562a47 100644 (file)
@@ -37,4 +37,9 @@ Use the commented-out version in PreconditionBlock::invert_diagblocks
 Comment on the modification of the right hand side in the Vanka
   smoother (which is done to yield a homogeneous boundary value
   problem on each patch)
-  
\ No newline at end of file
+  
+
+Someone should take a look at the handling of the sparsity variable
+  in the sparse matrix: it is subscribed to twice (explicitely and 
+  through the subscriptor); the explicit subscription is redundant,
+  but I don't have the time to look into that now.
index 6e47f26412e512bd2e859da956823b7e3d77db1d..2a295495a711ddf76365595d4113dea1d8aa37d0 100644 (file)
@@ -15,6 +15,8 @@
 #include <base/subscriptor.h>
 #include <base/smartpointer.h>
 
+#include <utility>
+
 
 //forward declarations
 template <typename number> class Vector;
@@ -885,7 +887,12 @@ class SparseMatrix : public Subscriptor
                                      * function equals the matrix norm with
                                      * respect to the mass matrix of the vector
                                      * representing the nodal values of the
-                                     * finite element function.
+                                     * finite element function. Note that
+                                     * even though the function's name might
+                                     * suggest something different, for historic
+                                     * reasons not the norm but its square
+                                     * is returned, as defined above by
+                                     * the scalar product.
                                      *
                                      * Note the order in which the matrix
                                      * appears. For non-symmetric matrices
@@ -1137,7 +1144,48 @@ class SparseMatrix : public Subscriptor
                                      * function.
                                      */
     unsigned int max_len;
-    
+
+                                    /**
+                                     * Version of #vmult# which only performs
+                                     * its actions on the region defined by
+                                     * #[begin_row,end_row)#. This function
+                                     * is called by #vmult# in the case
+                                     * of enabled multithreading.
+                                     */
+    template <typename somenumber>
+    void * threaded_vmult (Vector<somenumber>       &dst,
+                          const Vector<somenumber> &src,
+                          const unsigned int        begin_row,
+                          const unsigned int        end_row) const;
+
+                                    /**
+                                     * Version of #matrix_norm# which only
+                                     * performs its actions on the region
+                                     * defined by #[begin_row,end_row)#. This
+                                     * function is called by #matrix_norm# in
+                                     * the case of enabled multithreading.
+                                     */
+    template <typename somenumber>
+    void * threaded_matrix_norm (const Vector<somenumber> &v,
+                                const unsigned int        begin_row,
+                                const unsigned int        end_row,
+                                somenumber               *partial_sum) const;
+
+                                    /**
+                                     * Version of #residual# which only
+                                     * performs its actions on the region
+                                     * defined by #[begin_row,end_row)#. This
+                                     * function is called by #residual# in
+                                     * the case of enabled multithreading.
+                                     */
+    template <typename somenumber>
+    void * threaded_residual (Vector<somenumber>       &dst,
+                             const Vector<somenumber> &u,
+                             const Vector<somenumber> &b,
+                             const unsigned int        begin_row,
+                             const unsigned int        end_row,
+                             somenumber               *partial_norm) const;
+
 
                                     // make all other sparse matrices
                                     // friends
index a8c60cc9de0630d4f59d3efea5294e202d80581d..ec663eee8694b4c685c4135e9b254fd7738f7b27 100644 (file)
 #include <iomanip>
 #include <algorithm>
 
+#ifdef DEAL_II_USE_MT
+#  include <vector>
+#  include <numeric>
+
+#  include <base/thread_manager.h>
+
+#  define NTHREADS 4
+#endif
+
+
 
 
 
@@ -69,6 +79,7 @@ SparseMatrix<number>::~SparseMatrix ()
 {
   if (cols != 0)
     cols->unsubscribe();
+  cols = 0;
   
   if (val != 0)
     delete[] val;
@@ -123,6 +134,8 @@ template <typename number>
 void
 SparseMatrix<number>::clear ()
 {
+  if (cols != 0)
+    cols->unsubscribe ();
   cols = 0;
   if (val) delete[] val;
   val = 0;
@@ -199,8 +212,53 @@ SparseMatrix<number>::vmult (Vector<somenumber>& dst, const Vector<somenumber>&
   Assert (val != 0, ExcMatrixNotInitialized());
   Assert(m() == dst.size(), ExcDimensionsDontMatch(m(),dst.size()));
   Assert(n() == src.size(), ExcDimensionsDontMatch(n(),src.size()));
-
+  
   const unsigned int n_rows = m();
+
+#ifdef DEAL_II_USE_MT
+                                  // in MT mode: start new threads only
+                                  // if the matrix is sufficiently large.
+                                  // the limit is mostly artificial
+  if (n_rows/NTHREADS > 2000)
+    {
+      const unsigned int n_threads = NTHREADS;
+
+      ThreadManager thread_manager;
+      
+      const ThreadManager::Mem_Fun_Data4<const SparseMatrix<number>,
+                                        Vector<somenumber> &,
+                                        const Vector<somenumber> &,
+                                        unsigned int,
+                                        unsigned int> 
+       mem_fun_data_all (this, dst, src, 0, 0,
+                         &SparseMatrix<number>::template threaded_vmult<somenumber>  );
+      vector<ThreadManager::Mem_Fun_Data4<const SparseMatrix<number>,
+                                          Vector<somenumber> &,
+                                          const Vector<somenumber> &,
+                                          unsigned int,
+                                          unsigned int> >
+       mem_fun_data(n_threads, mem_fun_data_all);
+  
+                                      // spawn some jobs...
+      for (unsigned int i=0; i<n_threads; ++i)
+       {
+                                          // compute the range of rows
+                                          // they are to serve
+         mem_fun_data[i].arg3 = n_rows * i / n_threads;
+         mem_fun_data[i].arg4 = n_rows * (i+1) / n_threads;
+         
+         thread_manager.spawn (&mem_fun_data[i]);
+       };
+      
+                                      // ... and wait until they're finished
+      thread_manager.wait ();
+
+      return;
+    };
+#endif
+
+                                  // if not in MT mode or size<2000
+                                  // do it in an oldfashioned way
   const number *val_ptr = &val[cols->rowstart[0]];
   const int *colnum_ptr = &cols->colnums[cols->rowstart[0]];
   somenumber   *dst_ptr = &dst(0);
@@ -215,6 +273,38 @@ SparseMatrix<number>::vmult (Vector<somenumber>& dst, const Vector<somenumber>&
 };
 
 
+
+template <typename number>
+template <typename somenumber>
+void *
+SparseMatrix<number>::threaded_vmult (Vector<somenumber>       &dst,
+                                     const Vector<somenumber> &src,
+                                     const unsigned int        begin_row,
+                                     const unsigned int        end_row) const
+{
+#ifdef DEAL_II_USE_MT
+  const number *val_ptr = &val[cols->rowstart[begin_row]];
+  const int *colnum_ptr = &cols->colnums[cols->rowstart[begin_row]];
+  somenumber   *dst_ptr = &dst(begin_row);
+  for (unsigned int row=begin_row; row<end_row; ++row)
+    {
+      somenumber s = 0.;
+      const number *const val_end_of_row = &val[cols->rowstart[row+1]];
+      while (val_ptr != val_end_of_row)
+       s += *val_ptr++ * src(*colnum_ptr++);
+      *dst_ptr++ = s;
+    };
+#else
+                                  // this function should not be called
+                                  // when not in parallel mode.
+  Assert (false, ExcInternalError());
+#endif
+
+  return 0;
+};
+
+
+
 template <typename number>
 template <typename somenumber>
 void
@@ -296,8 +386,58 @@ SparseMatrix<number>::matrix_norm (const Vector<somenumber>& v) const
   Assert(m() == v.size(), ExcDimensionsDontMatch(m(),v.size()));
   Assert(n() == v.size(), ExcDimensionsDontMatch(n(),v.size()));
 
-  somenumber sum = 0.;
   const unsigned int n_rows = m();
+#ifdef DEAL_II_USE_MT
+                                  // if in MT mode and size sufficiently
+                                  // large: do it in parallel; the limit
+                                  // is mostly artificial
+  if (n_rows/NTHREADS > 2000)
+    {
+      const unsigned int n_threads = NTHREADS;
+
+      ThreadManager thread_manager;
+      
+      const ThreadManager::Mem_Fun_Data4<const SparseMatrix<number>,
+                                        const Vector<somenumber> &,
+                                        unsigned int,
+                                        unsigned int,
+                                        somenumber *> 
+       mem_fun_data_all (this, v, 0, 0, 0,
+                         &SparseMatrix<number>::template threaded_matrix_norm<somenumber>  );
+      vector<ThreadManager::Mem_Fun_Data4<const SparseMatrix<number>,
+                                          const Vector<somenumber> &,
+                                          unsigned int,
+                                          unsigned int,
+                                         somenumber *> >
+       mem_fun_data(n_threads, mem_fun_data_all);
+
+                                      // space for the norms of
+                                      // the different parts
+      vector<somenumber> partial_sums (n_threads, 0);
+      
+                                      // spawn some jobs...
+      for (unsigned int i=0; i<n_threads; ++i)
+       {
+                                          // compute the range of rows
+                                          // they are to serve
+         mem_fun_data[i].arg2 = n_rows * i / n_threads;
+         mem_fun_data[i].arg3 = n_rows * (i+1) / n_threads;
+         mem_fun_data[i].arg4 = &partial_sums[i];
+                 
+         thread_manager.spawn (&mem_fun_data[i]);
+       };
+      
+                                      // ... and wait until they're finished
+      thread_manager.wait ();
+                                      // accumulate the partial results
+      return accumulate (partial_sums.begin(),
+                        partial_sums.end(),
+                        0.);
+    };
+#endif
+                                  // if not in MT mode or the matrix is
+                                  // too small: do it one-by-one
+  somenumber sum = 0.;
   const number *val_ptr = &val[cols->rowstart[0]];
   const int *colnum_ptr = &cols->colnums[cols->rowstart[0]];
   for (unsigned int row=0; row<n_rows; ++row)
@@ -314,6 +454,38 @@ SparseMatrix<number>::matrix_norm (const Vector<somenumber>& v) const
 };
 
 
+
+template <typename number>
+template <typename somenumber>
+void *
+SparseMatrix<number>::threaded_matrix_norm (const Vector<somenumber> &v,
+                                           const unsigned int        begin_row,
+                                           const unsigned int        end_row,
+                                           somenumber               *partial_sum) const
+{
+#ifdef DEAL_II_USE_MT
+  somenumber sum = 0.;
+  const number *val_ptr = &val[cols->rowstart[begin_row]];
+  const int *colnum_ptr = &cols->colnums[cols->rowstart[begin_row]];
+  for (unsigned int row=begin_row; row<end_row; ++row)
+    {
+      somenumber s = 0.;
+      const number *val_end_of_row = &val[cols->rowstart[row+1]];
+      while (val_ptr != val_end_of_row)
+       s += *val_ptr++ * v(*colnum_ptr++);
+
+      sum += s* v(row);
+    };
+  *partial_sum = sum;
+  
+#else
+                                  // function should not have been called
+  Assert (false, ExcInternalError());
+#endif
+  return 0;
+};
+
+
 template <typename number>
 number SparseMatrix<number>::l1_norm () const
 {
@@ -356,9 +528,9 @@ number SparseMatrix<number>::linfty_norm () const
 template <typename number>
 template <typename somenumber>
 somenumber
-SparseMatrix<number>::residual (Vector<somenumber>dst,
-                               const Vector<somenumber>u,
-                               const Vector<somenumber>b) const
+SparseMatrix<number>::residual (Vector<somenumber>       &dst,
+                               const Vector<somenumber> &u,
+                               const Vector<somenumber> &b) const
 {
   Assert (cols != 0, ExcMatrixNotInitialized());
   Assert (val != 0, ExcMatrixNotInitialized());
@@ -366,11 +538,65 @@ SparseMatrix<number>::residual (Vector<somenumber>& dst,
   Assert(m() == b.size(), ExcDimensionsDontMatch(m(),b.size()));
   Assert(n() == u.size(), ExcDimensionsDontMatch(n(),u.size()));
 
-  somenumber s,norm=0.;   
+  const unsigned int n_rows = m();
+#ifdef DEAL_II_USE_MT
+                                  // if in MT mode and size sufficiently
+                                  // large: do it in parallel; the limit
+                                  // is mostly artificial
+  if (n_rows/NTHREADS > 2000)
+    {
+      const unsigned int n_threads = NTHREADS;
+
+      ThreadManager thread_manager;
+      
+      const ThreadManager::Mem_Fun_Data6<const SparseMatrix<number>,
+                                        Vector<somenumber> &,        // dst
+                                        const Vector<somenumber> &,  // u
+                                        const Vector<somenumber> &,  // b
+                                        unsigned int,        // begin_row
+                                        unsigned int,        // end_row
+                                        somenumber *>            // partial norm
+       mem_fun_data_all (this, dst, u, b, 0, 0, 0,
+                         &SparseMatrix<number>::template threaded_residual<somenumber>  );
+      vector<ThreadManager::Mem_Fun_Data6<const SparseMatrix<number>,
+                                          Vector<somenumber> &,
+                                          const Vector<somenumber> &,
+                                          const Vector<somenumber> &,
+                                          unsigned int,
+                                          unsigned int,
+                                         somenumber *> >
+       mem_fun_data(n_threads, mem_fun_data_all);
+
+                                      // space for the square norms of
+                                      // the different parts
+      vector<somenumber> partial_norms (n_threads, 0);
+      
+                                      // spawn some jobs...
+      for (unsigned int i=0; i<n_threads; ++i)
+       {
+                                          // compute the range of rows
+                                          // they are to serve
+         mem_fun_data[i].arg4 = n_rows * i / n_threads;
+         mem_fun_data[i].arg5 = n_rows * (i+1) / n_threads;
+         mem_fun_data[i].arg6 = &partial_norms[i];
+                 
+         thread_manager.spawn (&mem_fun_data[i]);
+       };
+      
+                                      // ... and wait until they're finished
+      thread_manager.wait ();
+                                      // accumulate the partial results
+      return sqrt(accumulate (partial_norms.begin(),
+                             partial_norms.end(),
+                             0.));
+    };
+#endif
   
-  for (unsigned int i=0;i<m();i++)
+  somenumber norm=0.;   
+  
+  for (unsigned int i=0; i<n_rows; ++i)
     {
-      s = b(i);
+      somenumber s = b(i);
       for (unsigned int j=cols->rowstart[i]; j<cols->rowstart[i+1] ;j++)
        {
          int p = cols->colnums[j];
@@ -383,6 +609,40 @@ SparseMatrix<number>::residual (Vector<somenumber>& dst,
 }
 
 
+template <typename number>
+template <typename somenumber>
+void *
+SparseMatrix<number>::threaded_residual (Vector<somenumber>       &dst,
+                                        const Vector<somenumber> &u,
+                                        const Vector<somenumber> &b,
+                                        const unsigned int        begin_row,
+                                        const unsigned int        end_row,
+                                        somenumber               *partial_norm) const
+{
+#ifdef DEAL_II_USE_MT
+  somenumber norm=0.;   
+  
+  for (unsigned int i=begin_row; i<end_row; ++i)
+    {
+      somenumber s = b(i);
+      for (unsigned int j=cols->rowstart[i]; j<cols->rowstart[i+1] ;j++)
+       {
+         int p = cols->colnums[j];
+         s -= val[j] * u(p);
+       }
+      dst(i) = s;
+      norm += dst(i)*dst(i);
+    };
+
+  *partial_norm = norm;
+#else
+  Assert (false, ExcInternalError());
+#endif
+
+  return 0;
+};
+
+
 
 template <typename number>
 template <typename somenumber>

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.