<h3>Specific improvements</h3>
<ol>
+<li> Fixed: SparseMatrix::mmult and SpareMatrix::Tmmult had a number of
+issues that are now fixed: (i) rebuilding the sparsity pattern was allowed
+even if several of the matrices involved in these operations shared a
+sparsity pattern; (ii) the functions had a vector argument that had a default
+value but the default value could not be used because it wasn't used in a
+template context deducible by the compiler.
+<br>
+(Wolfgang Bangerth, 2011/10/30)
-<li> New: parallel::distributed::Triangulation<dim>::mesh_reconstruction_after_repartitioning
-setting which is necessary for save()/load() to be deterministic. Otherwise the matrix assembly
-is done in a different order depending on the order of old refinements.
+<li> New:
+parallel::distributed::Triangulation<dim>::mesh_reconstruction_after_repartitioning
+setting which is necessary for save()/load() to be deterministic. Otherwise
+the matrix assembly is done in a different order depending on the order of old
+refinements.
<br>
(Timo Heister, 2011/10/26)
//---------------------------------------------------------------------------
// $Id$
//
-// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009, 2010 by the deal.II authors
+// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009, 2010, 2011 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
* and instead uses the sparsity
* pattern stored in <tt>C</tt>. In
* that case, make sure that it really
- * fits.
+ * fits. The default is to rebuild the
+ * sparsity pattern.
+ *
+ * @note Rebuilding the sparsity pattern
+ * requires changing it. This means that
+ * all other matrices that are associated
+ * with this sparsity pattern will
+ * then have invalid entries.
*/
- template <typename numberB, typename numberC, typename numberV>
+ template <typename numberB, typename numberC>
void mmult (SparseMatrix<numberC> &C,
const SparseMatrix<numberB> &B,
- const Vector<numberV> &V = Vector<numberV>(),
+ const Vector<number> &V = Vector<number>(),
const bool rebuild_sparsity_pattern = true) const;
/**
* and instead uses the sparsity
* pattern stored in <tt>C</tt>. In
* that case, make sure that it really
- * fits.
+ * fits. The default is to rebuild the
+ * sparsity pattern.
+ *
+ * @note Rebuilding the sparsity pattern
+ * requires changing it. This means that
+ * all other matrices that are associated
+ * with this sparsity pattern will
+ * then have invalid entries.
*/
- template <typename numberB, typename numberC, typename numberV>
+ template <typename numberB, typename numberC>
void Tmmult (SparseMatrix<numberC> &C,
const SparseMatrix<numberB> &B,
- const Vector<numberV> &V = Vector<numberV>(),
+ const Vector<number> &V = Vector<number>(),
const bool rebuild_sparsity_pattern = true) const;
//@}
template <typename number>
-template <typename numberB, typename numberC, typename numberV>
+template <typename numberB, typename numberC>
void
SparseMatrix<number>::mmult (SparseMatrix<numberC> &C,
const SparseMatrix<numberB> &B,
- const Vector<numberV> &V,
+ const Vector<number> &V,
const bool rebuild_sparsity_C) const
{
const bool use_vector = V.size() == n() ? true : false;
// clear previous content of C
if (rebuild_sparsity_C == true)
{
+ // we are about to change the sparsity
+ // pattern of C. this can not work if
+ // either A or B use the same sparsity
+ // pattern
+ Assert (&C.get_sparsity_pattern() != &this->get_sparsity_pattern(),
+ ExcMessage ("Can't use the same sparsity pattern for "
+ "different matrices if it is to be rebuilt."));
+ Assert (&C.get_sparsity_pattern() != &B.get_sparsity_pattern(),
+ ExcMessage ("Can't use the same sparsity pattern for "
+ "different matrices if it is to be rebuilt."));
+
// need to change the sparsity pattern of
// C, so cast away const-ness.
SparsityPattern & sp_C =
template <typename number>
-template <typename numberB, typename numberC, typename numberV>
+template <typename numberB, typename numberC>
void
SparseMatrix<number>::Tmmult (SparseMatrix<numberC> &C,
const SparseMatrix<numberB> &B,
- const Vector<numberV> &V,
+ const Vector<number> &V,
const bool rebuild_sparsity_C) const
{
const bool use_vector = V.size() == m() ? true : false;
// clear previous content of C
if (rebuild_sparsity_C == true)
{
- // need to change the sparsity pattern of
- // C, so cast away const-ness.
+ // we are about to change the sparsity
+ // pattern of C. this can not work if
+ // either A or B use the same sparsity
+ // pattern
+ Assert (&C.get_sparsity_pattern() != &this->get_sparsity_pattern(),
+ ExcMessage ("Can't use the same sparsity pattern for "
+ "different matrices if it is to be rebuilt."));
+ Assert (&C.get_sparsity_pattern() != &B.get_sparsity_pattern(),
+ ExcMessage ("Can't use the same sparsity pattern for "
+ "different matrices if it is to be rebuilt."));
+
+ // need to change the sparsity pattern of
+ // C, so cast away const-ness.
SparsityPattern & sp_C =
*(const_cast<SparsityPattern *>(&C.get_sparsity_pattern()));
C.clear();
// $Id$
// Version: $Name$
//
-// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2006, 2007, 2009 by the deal.II authors
+// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2006, 2007, 2009, 2011 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
template SparseMatrix<S1> &
SparseMatrix<S1>::copy_from<S2> (const SparseMatrix<S2> &);
- template
+ template
void SparseMatrix<S1>::copy_from<S2> (const FullMatrix<S2> &);
template void SparseMatrix<S1>::add<S2> (const S1,
const SparseMatrix<S2> &);
template void SparseMatrix<S1>::add<S2> (const unsigned int,
- const unsigned int,
- const unsigned int *,
+ const unsigned int,
+ const unsigned int *,
const S2 *,
- const bool,
+ const bool,
const bool);
template void SparseMatrix<S1>::set<S2> (const unsigned int,
- const unsigned int,
- const unsigned int *,
+ const unsigned int,
+ const unsigned int *,
const S2 *,
const bool);
}
const S1) const;
template void SparseMatrix<S1>::
SSOR_step<S2> (Vector<S2> &,
- const Vector<S2> &,
+ const Vector<S2> &,
const S1) const;
}
Tvmult_add (V1<S2> &, const V2<S3> &) const;
}
-for (S1, S2, S3, S4: REAL_SCALARS)
+for (S1, S2, S3: REAL_SCALARS)
{
template void SparseMatrix<S1>::
- mmult (SparseMatrix<S2> &, const SparseMatrix<S3> &, const Vector<S4>&,
+ mmult (SparseMatrix<S2> &, const SparseMatrix<S3> &, const Vector<S1>&,
const bool) const;
template void SparseMatrix<S1>::
- Tmmult (SparseMatrix<S2> &, const SparseMatrix<S3> &, const Vector<S4>&,
- const bool) const;
+ Tmmult (SparseMatrix<S2> &, const SparseMatrix<S3> &, const Vector<S1>&,
+ const bool) const;
}
// template SparseMatrix<S1> &
// SparseMatrix<S1>::copy_from<S2> (const SparseMatrix<S2> &);
-// template
+// template
// void SparseMatrix<S1>::copy_from<S2> (const FullMatrix<S2> &);
// template void SparseMatrix<S1>::add<S2> (const S1,
// const S1) const;
// template void SparseMatrix<S1>::
// SSOR_step<S2> (Vector<S2> &,
-// const Vector<S2> &,
+// const Vector<S2> &,
// const S1) const;
// }
--- /dev/null
+//---------------------------- sparse_matrix_Tmmult_01.cc,v ---------------------------
+// $Id$
+//
+// Copyright (C) 2011 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- sparse_matrix_Tmmult_01.cc,v ---------------------------
+
+// check SparseMatrix::Tmmult. this function has a default argument
+// that could previously not be instantiated because it was in a
+// non-deduced context but that should not be possible to omit
+//
+// this test checks SparseMatrix::Tmmult without additional arguments
+
+#include "../tests.h"
+#include <cmath>
+#include <fstream>
+#include <iomanip>
+#include <iomanip>
+#include <cstdlib>
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/sparsity_pattern.h>
+#include <deal.II/lac/sparse_matrix.h>
+
+
+void test (const unsigned int n)
+{
+ // Create some random full matrices in the
+ // data structures of a sparse matrix
+ SparsityPattern sp (n,n);
+ SparsityPattern C_sp (n,n);
+ for (unsigned int i=0;i<n;++i)
+ for (unsigned int j=0;j<n;++j)
+ {
+ sp.add (i,j);
+ C_sp.add (i,j);
+ }
+ sp.compress ();
+ C_sp.compress ();
+
+ SparseMatrix<double> A(sp);
+ SparseMatrix<double> B(sp);
+ SparseMatrix<double> C(C_sp);
+
+ for (unsigned int i=0;i<n;++i)
+ {
+ for (unsigned int j=0;j<n;++j)
+ A.set(i,j,std::rand());
+ for (unsigned int j=0;j<n;++j)
+ B.set(i,j,std::rand());
+ }
+
+ // now form the matrix-matrix product and
+ // initialize a test rhs
+ A.Tmmult (C, B);
+
+ Vector<double> x(n), y(n), z(n), tmp(n);
+ for (unsigned int j=0;j<n;++j)
+ x(j) = std::rand();
+
+ // then test for correctness
+ C.vmult (y, x);
+
+ B.vmult (tmp, x);
+ A.Tvmult (z, tmp);
+
+ y -= z;
+ Assert (y.l2_norm() <= 1e-12 * z.l2_norm(),
+ ExcInternalError());
+
+ deallog << "OK" << std::endl;
+}
+
+
+int
+main ()
+{
+ const std::string logname = JobIdentifier::base_name(__FILE__) + std::string("/output");
+ std::ofstream logfile(logname.c_str());
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ std::srand(3391466);
+
+ test(3);
+ test(7);
+}
--- /dev/null
+
+DEAL::OK
+DEAL::OK
--- /dev/null
+//---------------------------- sparse_matrix_Tmmult_02.cc,v ---------------------------
+// $Id$
+//
+// Copyright (C) 2011 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- sparse_matrix_Tmmult_02.cc,v ---------------------------
+
+// check SparseMatrix::Tmmult. this function has a default argument
+// that could previously not be instantiated because it was in a
+// non-deduced context but that should not be possible to omit
+//
+// this test checks SparseMatrix::Tmmult with additional arguments
+
+#include "../tests.h"
+#include <cmath>
+#include <fstream>
+#include <iomanip>
+#include <iomanip>
+#include <cstdlib>
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/sparsity_pattern.h>
+#include <deal.II/lac/sparse_matrix.h>
+
+
+void test (const unsigned int n)
+{
+ // Create some random full matrices in the
+ // data structures of a sparse matrix
+ SparsityPattern sp (n,n);
+ SparsityPattern C_sp (n,n);
+ for (unsigned int i=0;i<n;++i)
+ for (unsigned int j=0;j<n;++j)
+ {
+ sp.add (i,j);
+ C_sp.add (i,j);
+ }
+ sp.compress ();
+ C_sp.compress ();
+
+ SparseMatrix<double> A(sp);
+ SparseMatrix<double> B(sp);
+ SparseMatrix<double> C(C_sp);
+
+ for (unsigned int i=0;i<n;++i)
+ {
+ for (unsigned int j=0;j<n;++j)
+ A.set(i,j,std::rand());
+ for (unsigned int j=0;j<n;++j)
+ B.set(i,j,std::rand());
+ }
+
+ Vector<double> v(n);
+ for (unsigned int j=0;j<n;++j)
+ v(j) = std::rand();
+
+ // now form the matrix-matrix product and
+ // initialize a test rhs
+ A.Tmmult (C, B, v);
+
+ Vector<double> x(n), y(n), z(n), tmp(n);
+ for (unsigned int j=0;j<n;++j)
+ x(j) = std::rand();
+
+ // then test for correctness
+ C.vmult (y, x);
+
+ B.vmult (tmp, x);
+ tmp.scale (v);
+ A.Tvmult (z, tmp);
+
+ y -= z;
+ Assert (y.l2_norm() <= 1e-12 * z.l2_norm(),
+ ExcInternalError());
+
+ deallog << "OK" << std::endl;
+}
+
+
+int
+main ()
+{
+ const std::string logname = JobIdentifier::base_name(__FILE__) + std::string("/output");
+ std::ofstream logfile(logname.c_str());
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ std::srand(3391466);
+
+ test(3);
+ test(7);
+}
--- /dev/null
+
+DEAL::OK
+DEAL::OK
--- /dev/null
+//---------------------------- sparse_matrix_mmult_01.cc,v ---------------------------
+// $Id$
+//
+// Copyright (C) 2011 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- sparse_matrix_mmult_01.cc,v ---------------------------
+
+// check SparseMatrix::mmult. this function has a default argument
+// that could previously not be instantiated because it was in a
+// non-deduced context but that should not be possible to omit
+//
+// this test checks SparseMatrix::mmult without additional arguments
+
+#include "../tests.h"
+#include <cmath>
+#include <fstream>
+#include <iomanip>
+#include <iomanip>
+#include <cstdlib>
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/sparsity_pattern.h>
+#include <deal.II/lac/sparse_matrix.h>
+
+
+void test (const unsigned int n)
+{
+ // Create some random full matrices in the
+ // data structures of a sparse matrix
+ SparsityPattern sp (n,n);
+ SparsityPattern C_sp (n,n);
+ for (unsigned int i=0;i<n;++i)
+ for (unsigned int j=0;j<n;++j)
+ {
+ sp.add (i,j);
+ C_sp.add (i,j);
+ }
+ sp.compress ();
+ C_sp.compress ();
+
+ SparseMatrix<double> A(sp);
+ SparseMatrix<double> B(sp);
+ SparseMatrix<double> C(C_sp);
+
+ for (unsigned int i=0;i<n;++i)
+ {
+ for (unsigned int j=0;j<n;++j)
+ A.set(i,j,std::rand());
+ for (unsigned int j=0;j<n;++j)
+ B.set(i,j,std::rand());
+ }
+
+ // now form the matrix-matrix product and
+ // initialize a test rhs
+ A.mmult (C, B);
+
+ Vector<double> x(n), y(n), z(n), tmp(n);
+ for (unsigned int j=0;j<n;++j)
+ x(j) = std::rand();
+
+ // then test for correctness
+ C.vmult (y, x);
+
+ B.vmult (tmp, x);
+ A.vmult (z, tmp);
+
+ y -= z;
+ Assert (y.l2_norm() <= 1e-12 * z.l2_norm(),
+ ExcInternalError());
+
+ deallog << "OK" << std::endl;
+}
+
+
+int
+main ()
+{
+ const std::string logname = JobIdentifier::base_name(__FILE__) + std::string("/output");
+ std::ofstream logfile(logname.c_str());
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ std::srand(3391466);
+
+ test(3);
+ test(7);
+}
--- /dev/null
+
+DEAL::OK
+DEAL::OK
--- /dev/null
+//---------------------------- sparse_matrix_mmult_01.cc,v ---------------------------
+// $Id$
+//
+// Copyright (C) 2011 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- sparse_matrix_mmult_01.cc,v ---------------------------
+
+// check SparseMatrix::mmult. this function has a default argument
+// that could previously not be instantiated because it was in a
+// non-deduced context but that should not be possible to omit
+//
+// this test checks SparseMatrix::mmult with additional arguments
+
+#include "../tests.h"
+#include <cmath>
+#include <fstream>
+#include <iomanip>
+#include <iomanip>
+#include <cstdlib>
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/sparsity_pattern.h>
+#include <deal.II/lac/sparse_matrix.h>
+
+
+void test (const unsigned int n)
+{
+ // Create some random full matrices in the
+ // data structures of a sparse matrix
+ SparsityPattern sp (n,n);
+ SparsityPattern C_sp (n,n);
+ for (unsigned int i=0;i<n;++i)
+ for (unsigned int j=0;j<n;++j)
+ {
+ sp.add (i,j);
+ C_sp.add (i,j);
+ }
+ sp.compress ();
+ C_sp.compress ();
+
+ SparseMatrix<double> A(sp);
+ SparseMatrix<double> B(sp);
+ SparseMatrix<double> C(C_sp);
+
+ for (unsigned int i=0;i<n;++i)
+ {
+ for (unsigned int j=0;j<n;++j)
+ A.set(i,j,std::rand());
+ for (unsigned int j=0;j<n;++j)
+ B.set(i,j,std::rand());
+ }
+
+ Vector<double> v(n);
+ for (unsigned int j=0;j<n;++j)
+ v(j) = std::rand();
+
+ // now form the matrix-matrix product and
+ // initialize a test rhs
+ A.mmult (C, B, v);
+
+ Vector<double> x(n), y(n), z(n), tmp(n);
+ for (unsigned int j=0;j<n;++j)
+ x(j) = std::rand();
+
+ // then test for correctness
+ C.vmult (y, x);
+
+ B.vmult (tmp, x);
+ tmp.scale (v);
+ A.vmult (z, tmp);
+
+ y -= z;
+ Assert (y.l2_norm() <= 1e-12 * z.l2_norm(),
+ ExcInternalError());
+
+ deallog << "OK" << std::endl;
+}
+
+
+int
+main ()
+{
+ const std::string logname = JobIdentifier::base_name(__FILE__) + std::string("/output");
+ std::ofstream logfile(logname.c_str());
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ std::srand(3391466);
+
+ test(3);
+ test(7);
+}
--- /dev/null
+
+DEAL::OK
+DEAL::OK