]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Fix SparseMatrix::{T,}mmult in two different regards.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 30 Oct 2011 14:32:26 +0000 (14:32 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 30 Oct 2011 14:32:26 +0000 (14:32 +0000)
git-svn-id: https://svn.dealii.org/trunk@24702 0785d39b-7218-0410-832d-ea1e28bc413d

12 files changed:
deal.II/doc/news/changes.h
deal.II/include/deal.II/lac/sparse_matrix.h
deal.II/include/deal.II/lac/sparse_matrix.templates.h
deal.II/source/lac/sparse_matrix.inst.in
tests/lac/sparse_matrix_Tmmult_01.cc [new file with mode: 0644]
tests/lac/sparse_matrix_Tmmult_01/cmp/generic [new file with mode: 0644]
tests/lac/sparse_matrix_Tmmult_02.cc [new file with mode: 0644]
tests/lac/sparse_matrix_Tmmult_02/cmp/generic [new file with mode: 0644]
tests/lac/sparse_matrix_mmult_01.cc [new file with mode: 0644]
tests/lac/sparse_matrix_mmult_01/cmp/generic [new file with mode: 0644]
tests/lac/sparse_matrix_mmult_02.cc [new file with mode: 0644]
tests/lac/sparse_matrix_mmult_02/cmp/generic [new file with mode: 0644]

index ff36636034dff29dbb8af4b390eb0fc3f6664989..d66ce3821b0bfdbac0b1f4b9a4d22fc50f9525ba 100644 (file)
@@ -48,10 +48,20 @@ enabled due to a missing include file in file
 <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)
 
index 13ef275822fd3341c7eb239dd9e4c1fe3f72c69f..46274071b98a5a376124d51a2a571301ecfecba8 100644 (file)
@@ -1,7 +1,7 @@
 //---------------------------------------------------------------------------
 //    $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
@@ -1504,12 +1504,19 @@ class SparseMatrix : public virtual Subscriptor
                                      * 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;
 
                                     /**
@@ -1545,12 +1552,19 @@ class SparseMatrix : public virtual Subscriptor
                                      * 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;
 
 //@}
index 1c899618cc159f5d520270a0bfc5257adfc38efc..0550a859ffaad6e87f8d6de50d398316087a75db 100644 (file)
@@ -897,11 +897,11 @@ SparseMatrix<number>::matrix_scalar_product (const Vector<somenumber>& u,
 
 
 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;
@@ -916,6 +916,17 @@ SparseMatrix<number>::mmult (SparseMatrix<numberC>       &C,
                                   // 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 =
@@ -1034,11 +1045,11 @@ SparseMatrix<number>::mmult (SparseMatrix<numberC>       &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;
@@ -1053,8 +1064,19 @@ SparseMatrix<number>::Tmmult (SparseMatrix<numberC>       &C,
                                   // 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();
index f5b3b0eabdd13b9fac8aabfa9c8389ce5b556b26..2ce3a691ab0f2439620c1745b394bdd44439960b 100644 (file)
@@ -2,7 +2,7 @@
 //    $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
@@ -26,22 +26,22 @@ for (S1, S2 : REAL_SCALARS)
     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);
   }
@@ -117,7 +117,7 @@ for (S1, S2 : REAL_SCALARS)
                     const S1) const;
     template void SparseMatrix<S1>::
       SSOR_step<S2> (Vector<S2> &,
-                    const Vector<S2> &, 
+                    const Vector<S2> &,
                     const S1) const;
   }
 
@@ -134,14 +134,14 @@ for (S1, S2, S3 : REAL_SCALARS;
       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;
   }
 
 
@@ -160,7 +160,7 @@ for (S1, S2, S3, S4: REAL_SCALARS)
 //     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,
@@ -233,7 +233,7 @@ for (S1, S2, S3, S4: REAL_SCALARS)
 //                  const S1) const;
 //     template void SparseMatrix<S1>::
 //       SSOR_step<S2> (Vector<S2> &,
-//                  const Vector<S2> &, 
+//                  const Vector<S2> &,
 //                  const S1) const;
 //   }
 
diff --git a/tests/lac/sparse_matrix_Tmmult_01.cc b/tests/lac/sparse_matrix_Tmmult_01.cc
new file mode 100644 (file)
index 0000000..5cfb459
--- /dev/null
@@ -0,0 +1,92 @@
+//----------------------------  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);
+}
diff --git a/tests/lac/sparse_matrix_Tmmult_01/cmp/generic b/tests/lac/sparse_matrix_Tmmult_01/cmp/generic
new file mode 100644 (file)
index 0000000..8b3b075
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::OK
+DEAL::OK
diff --git a/tests/lac/sparse_matrix_Tmmult_02.cc b/tests/lac/sparse_matrix_Tmmult_02.cc
new file mode 100644 (file)
index 0000000..3a62974
--- /dev/null
@@ -0,0 +1,97 @@
+//----------------------------  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);
+}
diff --git a/tests/lac/sparse_matrix_Tmmult_02/cmp/generic b/tests/lac/sparse_matrix_Tmmult_02/cmp/generic
new file mode 100644 (file)
index 0000000..8b3b075
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::OK
+DEAL::OK
diff --git a/tests/lac/sparse_matrix_mmult_01.cc b/tests/lac/sparse_matrix_mmult_01.cc
new file mode 100644 (file)
index 0000000..706f1e6
--- /dev/null
@@ -0,0 +1,92 @@
+//----------------------------  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);
+}
diff --git a/tests/lac/sparse_matrix_mmult_01/cmp/generic b/tests/lac/sparse_matrix_mmult_01/cmp/generic
new file mode 100644 (file)
index 0000000..8b3b075
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::OK
+DEAL::OK
diff --git a/tests/lac/sparse_matrix_mmult_02.cc b/tests/lac/sparse_matrix_mmult_02.cc
new file mode 100644 (file)
index 0000000..6aaf2f9
--- /dev/null
@@ -0,0 +1,97 @@
+//----------------------------  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);
+}
diff --git a/tests/lac/sparse_matrix_mmult_02/cmp/generic b/tests/lac/sparse_matrix_mmult_02/cmp/generic
new file mode 100644 (file)
index 0000000..8b3b075
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::OK
+DEAL::OK

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.