]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Change matrix_norm to matrix_norm_square and implement matrix_scalar_product.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 5 May 2000 17:10:43 +0000 (17:10 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Fri, 5 May 2000 17:10:43 +0000 (17:10 +0000)
git-svn-id: https://svn.dealii.org/trunk@2815 0785d39b-7218-0410-832d-ea1e28bc413d

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

index 0755a0dcce48f71d1561795dc0c1714b3c3d28d0..3e1d6fab198759650e968d34366204bac0b45e8e 100644 (file)
@@ -33,7 +33,11 @@ template void SparseMatrix<TYPEMAT>::Tvmult_add (Vector<TYPE2> &,
                                                const Vector<TYPE2> &) const;
 
 template TYPE2
-SparseMatrix<TYPEMAT>::matrix_norm (const Vector<TYPE2> &) const;
+SparseMatrix<TYPEMAT>::matrix_norm_square (const Vector<TYPE2> &) const;
+
+template TYPE2
+SparseMatrix<TYPEMAT>::matrix_scalar_product (const Vector<TYPE2> &,
+                                             const Vector<TYPE2> &) const;
 
 template TYPE2 SparseMatrix<TYPEMAT>::residual (Vector<TYPE2> &,
                                               const Vector<TYPE2> &,
index 1757d3d8c1db61e3a915f1cbee946532c98ce8bf..7a880eab2f36657431394abc4957332d6b0a678c 100644 (file)
@@ -369,9 +369,10 @@ class SparseMatrix : public Subscriptor
     void Tvmult_add (Vector<somenumber>& dst, const Vector<somenumber>& src) const;
   
                                     /**
-                                     * Return the norm of the vector
-                                     * $v$ with respect to the norm
-                                     * induced by this matrix,
+                                     * Return the square of the norm
+                                     * of the vector $v$ with respect
+                                     * to the norm induced by this
+                                     * matrix,
                                      * i.e. $\left(v,Mv\right)$. This
                                      * is useful, e.g. in the finite
                                      * element context, where the
@@ -380,19 +381,13 @@ class SparseMatrix : public Subscriptor
                                      * respect to the mass matrix of
                                      * the vector representing the
                                      * nodal values of the 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.
+                                     * element function.
                                      *
                                      * Obviously, the matrix needs to
                                      * be square for this operation.
                                      */
     template <typename somenumber>
-    somenumber matrix_norm (const Vector<somenumber> &v) const;
+    somenumber matrix_norm_square (const Vector<somenumber> &v) const;
 
                                     /**
                                      * Compute the matrix scalar
@@ -693,19 +688,39 @@ class SparseMatrix : public Subscriptor
                         const unsigned int        end_row) const;
 
                                     /**
-                                     * Version of #matrix_norm# which
+                                     * Version of
+                                     * #matrix_norm_square# 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.
+                                     * #matrix_norm_square# in the
+                                     * case of enabled
+                                     * multithreading.
+                                     */
+    template <typename somenumber>
+    void threaded_matrix_norm_square (const Vector<somenumber> &v,
+                                     const unsigned int        begin_row,
+                                     const unsigned int        end_row,
+                                     somenumber               *partial_sum) const;
+
+                                    /**
+                                     * Version of
+                                     * #matrix_scalar_product# which
+                                     * only performs its actions on
+                                     * the region defined by
+                                     * #[begin_row,end_row)#. This
+                                     * function is called by
+                                     * #matrix_scalar_product# 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;
+    void threaded_matrix_scalar_product (const Vector<somenumber> &u,
+                                        const Vector<somenumber> &v,
+                                        const unsigned int        begin_row,
+                                        const unsigned int        end_row,
+                                        somenumber               *partial_sum) const;
 
                                     /**
                                      * Version of #residual# which
index 5492ec15f00054339110f196f5aaa6bda7d203a2..e11a1193e24b881b8b630ea82f1ae6b907b36140 100644 (file)
@@ -337,7 +337,7 @@ SparseMatrix<number>::Tvmult_add (Vector<somenumber>& dst, const Vector<somenumb
 template <typename number>
 template <typename somenumber>
 somenumber
-SparseMatrix<number>::matrix_norm (const Vector<somenumber>& v) const
+SparseMatrix<number>::matrix_norm_square (const Vector<somenumber>& v) const
 {
   Assert (cols != 0, ExcMatrixNotInitialized());
   Assert (val != 0, ExcMatrixNotInitialized());
@@ -361,7 +361,7 @@ SparseMatrix<number>::matrix_norm (const Vector<somenumber>& v) const
       for (unsigned int i=0; i<n_threads; ++i)
        Threads::spawn (thread_manager,
                        Threads::encapsulate (&SparseMatrix<number>::
-                                             template threaded_matrix_norm<somenumber>)
+                                             template threaded_matrix_norm_square<somenumber>)
                        .collect_args (this, v,
                                       n_rows * i / n_threads,
                                       n_rows * (i+1) / n_threads,
@@ -394,13 +394,14 @@ 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
+SparseMatrix<number>::threaded_matrix_norm_square (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.;
@@ -424,6 +425,103 @@ SparseMatrix<number>::threaded_matrix_norm (const Vector<somenumber> &v,
 };
 
 
+
+template <typename number>
+template <typename somenumber>
+somenumber
+SparseMatrix<number>::matrix_scalar_product (const Vector<somenumber>& u,
+                                            const Vector<somenumber>& v) const
+{
+  Assert (cols != 0, ExcMatrixNotInitialized());
+  Assert (val != 0, ExcMatrixNotInitialized());
+  Assert(m() == u.size(), ExcDimensionsDontMatch(m(),u.size()));
+  Assert(n() == v.size(), ExcDimensionsDontMatch(n(),v.size()));
+
+  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/multithread_info.n_default_threads > 2000)
+    {
+      const unsigned int n_threads = multithread_info.n_default_threads;
+
+                                      // space for the norms of
+                                      // the different parts
+      vector<somenumber> partial_sums (n_threads, 0);
+      ACE_Thread_Manager thread_manager;
+                                      // spawn some jobs...
+      for (unsigned int i=0; i<n_threads; ++i)
+       Threads::spawn (thread_manager,
+                       Threads::encapsulate (&SparseMatrix<number>::
+                                             template threaded_matrix_scalar_product<somenumber>)
+                       .collect_args (this, u, v,
+                                      n_rows * i / n_threads,
+                                      n_rows * (i+1) / n_threads,
+                                      &partial_sums[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 unsigned int *colnum_ptr = &cols->colnums[cols->rowstart[0]];
+  for (unsigned int row=0; row<n_rows; ++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* u(row);
+    };
+
+  return sum;
+};
+
+
+
+template <typename number>
+template <typename somenumber>
+void
+SparseMatrix<number>::threaded_matrix_scalar_product (const Vector<somenumber> &u,
+                                                     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 unsigned 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* u(row);
+    };
+  *partial_sum = sum;
+  
+#else
+                                  // function should not have been called
+  Assert (false, ExcInternalError());
+#endif
+};
+
+
+
+
+
 template <typename number>
 number SparseMatrix<number>::l1_norm () const
 {
@@ -559,6 +657,7 @@ SparseMatrix<number>::threaded_residual (Vector<somenumber>       &dst,
 };
 
 
+
 template <typename number>
 template <typename somenumber>
 void
@@ -584,6 +683,7 @@ SparseMatrix<number>::precondition_Jacobi (Vector<somenumber>& dst,
 };
 
 
+
 template <typename number>
 template <typename somenumber>
 void

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.