From: wolf Date: Fri, 5 May 2000 17:10:43 +0000 (+0000) Subject: Change matrix_norm to matrix_norm_square and implement matrix_scalar_product. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9312a309349301ca9746b8c0385166cc9b4fde86;p=dealii-svn.git Change matrix_norm to matrix_norm_square and implement matrix_scalar_product. git-svn-id: https://svn.dealii.org/trunk@2815 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/sparse_matrix.2.templates b/deal.II/lac/include/lac/sparse_matrix.2.templates index 0755a0dcce..3e1d6fab19 100644 --- a/deal.II/lac/include/lac/sparse_matrix.2.templates +++ b/deal.II/lac/include/lac/sparse_matrix.2.templates @@ -33,7 +33,11 @@ template void SparseMatrix::Tvmult_add (Vector &, const Vector &) const; template TYPE2 -SparseMatrix::matrix_norm (const Vector &) const; +SparseMatrix::matrix_norm_square (const Vector &) const; + +template TYPE2 +SparseMatrix::matrix_scalar_product (const Vector &, + const Vector &) const; template TYPE2 SparseMatrix::residual (Vector &, const Vector &, diff --git a/deal.II/lac/include/lac/sparse_matrix.h b/deal.II/lac/include/lac/sparse_matrix.h index 1757d3d8c1..7a880eab2f 100644 --- a/deal.II/lac/include/lac/sparse_matrix.h +++ b/deal.II/lac/include/lac/sparse_matrix.h @@ -369,9 +369,10 @@ class SparseMatrix : public Subscriptor void Tvmult_add (Vector& dst, const Vector& 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 - somenumber matrix_norm (const Vector &v) const; + somenumber matrix_norm_square (const Vector &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 + void threaded_matrix_norm_square (const Vector &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 - void threaded_matrix_norm (const Vector &v, - const unsigned int begin_row, - const unsigned int end_row, - somenumber *partial_sum) const; + void threaded_matrix_scalar_product (const Vector &u, + const Vector &v, + const unsigned int begin_row, + const unsigned int end_row, + somenumber *partial_sum) const; /** * Version of #residual# which diff --git a/deal.II/lac/include/lac/sparse_matrix.templates.h b/deal.II/lac/include/lac/sparse_matrix.templates.h index 5492ec15f0..e11a1193e2 100644 --- a/deal.II/lac/include/lac/sparse_matrix.templates.h +++ b/deal.II/lac/include/lac/sparse_matrix.templates.h @@ -337,7 +337,7 @@ SparseMatrix::Tvmult_add (Vector& dst, const Vector template somenumber -SparseMatrix::matrix_norm (const Vector& v) const +SparseMatrix::matrix_norm_square (const Vector& v) const { Assert (cols != 0, ExcMatrixNotInitialized()); Assert (val != 0, ExcMatrixNotInitialized()); @@ -361,7 +361,7 @@ SparseMatrix::matrix_norm (const Vector& v) const for (unsigned int i=0; i:: - template threaded_matrix_norm) + template threaded_matrix_norm_square) .collect_args (this, v, n_rows * i / n_threads, n_rows * (i+1) / n_threads, @@ -394,13 +394,14 @@ SparseMatrix::matrix_norm (const Vector& v) const }; + template template void -SparseMatrix::threaded_matrix_norm (const Vector &v, - const unsigned int begin_row, - const unsigned int end_row, - somenumber *partial_sum) const +SparseMatrix::threaded_matrix_norm_square (const Vector &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::threaded_matrix_norm (const Vector &v, }; + +template +template +somenumber +SparseMatrix::matrix_scalar_product (const Vector& u, + const Vector& 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 partial_sums (n_threads, 0); + ACE_Thread_Manager thread_manager; + // spawn some jobs... + for (unsigned int i=0; i:: + template threaded_matrix_scalar_product) + .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; rowrowstart[row+1]]; + while (val_ptr != val_end_of_row) + s += *val_ptr++ * v(*colnum_ptr++); + + sum += s* u(row); + }; + + return sum; +}; + + + +template +template +void +SparseMatrix::threaded_matrix_scalar_product (const Vector &u, + const Vector &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; rowrowstart[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 number SparseMatrix::l1_norm () const { @@ -559,6 +657,7 @@ SparseMatrix::threaded_residual (Vector &dst, }; + template template void @@ -584,6 +683,7 @@ SparseMatrix::precondition_Jacobi (Vector& dst, }; + template template void