From 6fbfff772afc4a54873ef28449306b0a35d54c99 Mon Sep 17 00:00:00 2001 From: guido Date: Fri, 11 Oct 2002 12:13:12 +0000 Subject: [PATCH] improvements in copy_from, conjugate_add git-svn-id: https://svn.dealii.org/trunk@6631 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/sparse_matrix_ez.h | 92 ++++++++++++++-------- 1 file changed, 60 insertions(+), 32 deletions(-) diff --git a/deal.II/lac/include/lac/sparse_matrix_ez.h b/deal.II/lac/include/lac/sparse_matrix_ez.h index c892790bd8..340b79ccb6 100644 --- a/deal.II/lac/include/lac/sparse_matrix_ez.h +++ b/deal.II/lac/include/lac/sparse_matrix_ez.h @@ -644,16 +644,16 @@ class SparseMatrixEZ : public Subscriptor * Obviously, the matrix needs to * be square for this operation. */ - template - somenumber matrix_norm_square (const Vector &v) const; +// template +// somenumber matrix_norm_square (const Vector &v) const; /** * Compute the matrix scalar * product $\left(u,Mv\right)$. */ - template - somenumber matrix_scalar_product (const Vector &u, - const Vector &v) const; +// template +// somenumber matrix_scalar_product (const Vector &u, +// const Vector &v) const; /** * Frobenius-norm of the matrix. @@ -731,6 +731,15 @@ class SparseMatrixEZ : public Subscriptor const Vector &src, const number om = 1.) const; + /** + * Add the matrix @p{A} + * conjugated by @p{B}, that is, + * $B A B^T$ to this object. + */ + template + void conjugate_add (const MATRIXA& A, + const MATRIXB& B); + /** * STL-like iterator with the * first entry. @@ -774,7 +783,7 @@ class SparseMatrixEZ : public Subscriptor * nonzero entry of the matrix * per line. */ - void print (std::ostream &out) const; +// void print (std::ostream &out) const; /** * Print the matrix in the usual @@ -815,12 +824,12 @@ class SparseMatrixEZ : public Subscriptor * may produce @em{large} amounts of * output if applied to a large matrix! */ - void print_formatted (std::ostream &out, - const unsigned int precision = 3, - const bool scientific = true, - const unsigned int width = 0, - const char *zero_string = " ", - const double denominator = 1.) const; +// void print_formatted (std::ostream &out, +// const unsigned int precision = 3, +// const bool scientific = true, +// const unsigned int width = 0, +// const char *zero_string = " ", +// const double denominator = 1.) const; /** * Determine an estimate for the @@ -914,24 +923,6 @@ class SparseMatrixEZ : public Subscriptor const unsigned int end_row, somenumber *partial_sum) const; - /** - * Version of @p{residual} which - * only performs its actions on - * the region defined by - * @p{[begin_row,end_row)} (these - * numbers are the components of - * @p{interval}). This function is - * called by @p{residual} in the - * case of enabled - * multithreading. - */ - template - void threaded_residual (Vector &dst, - const Vector &u, - const Vector &b, - const std::pair interval, - somenumber *partial_norm) const; - /** * Number of columns. This is * used to check vector @@ -1346,7 +1337,8 @@ SparseMatrixEZ::copy_from (const MATRIX& M) while (start != final) { - set(start->row(), start->column(), start->value()); + if (start->value() != 0.) + set(start->row(), start->column(), start->value()); ++start; } return *this; @@ -1361,16 +1353,52 @@ SparseMatrixEZ::add_scaled (number factor, { Assert (M.m() == m(), ExcDimensionMismatch(M.m(), m())); Assert (M.n() == n(), ExcDimensionMismatch(M.n(), n())); + + if (factor == 0.) + return; typename MATRIX::const_iterator start = M.begin(); const typename MATRIX::const_iterator final = M.end(); while (start != final) { - add(start->row(), start->column(), factor * start->value()); + if (start->value() != 0.) + add(start->row(), start->column(), factor * start->value()); ++start; } } + +template +template +void +SparseMatrixEZ::conjugate_add (const MATRIXA& A, + const MATRIXB& B) +{ +// Compute the result +// r_ij = \sum_kl b_ik b_jl a_kl + + typename MATRIXB::const_iterator b1 = B.begin(); + typename MATRIXB::const_iterator b2 = B.begin(); + const typename MATRIXB::const_iterator b_final = B.end(); + while (b1 != b_final) + { + const unsigned int i = b->row(); + const unsigned int k = b->column(); + while (b2 != b_final) + { + const unsigned int j = b->row(); + const unsigned int l = b->column(); + + const typename MATRIXA::value_type a = A.el(k,l); + + if (a != 0.) + add (i, j, a * b1->value() * b2->value()); + ++b2; + } + ++b1; + } +} + #endif /*---------------------------- sparse_matrix.h ---------------------------*/ -- 2.39.5