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
* 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
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
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());
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,
};
+
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.;
};
+
+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
{
};
+
template <typename number>
template <typename somenumber>
void
};
+
template <typename number>
template <typename somenumber>
void