add PETScWrappers::*Matrix::add(other, factor) with test. Patch by Jose Javier Munoz...

git-svn-id: https://svn.dealii.org/trunk@28478 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h
index a5dd2d7f2fb7b9b4380437dab90ba2e0c0a5a3d8..4f791bf6afa78e07eaec25377e6551a0fd5a8c77 100644 (file)
--- a/deal.II/doc/news/changes.h
+++ b/deal.II/doc/news/changes.h
@@ -208,6 +208,11 @@ DoFHandler, in particular removal of specializations.
 <h3>Specific improvements</h3>
 
 <ol>
+<li> New: PETScWrappers::*Matrix::add(other, factor) to
+add a scaled other matrix to the current matrix.
+<br>
+(Jose Javier Munoz Criollo, 2013/02/19)
+
 <li> New: GridGenerator::extrude_triangulation() allows
 you to extrude a 2d mesh to turn it into a 3d mesh.
 <br>

diff --git a/deal.II/include/deal.II/lac/petsc_matrix_base.h b/deal.II/include/deal.II/lac/petsc_matrix_base.h
index a8f137b8460dc5be4f8dbdf8c7f16a59b1f7c7cf..bc81805f3f55fe87a884328cf2ba20559599f046 100644 (file)
--- a/deal.II/include/deal.II/lac/petsc_matrix_base.h
+++ b/deal.II/include/deal.II/lac/petsc_matrix_base.h
@@ -926,6 +926,13 @@ namespace PETScWrappers
      */
     MatrixBase &operator /= (const PetscScalar factor);
 
+    /**
+     * Add the matrix @p other scaled by the factor @p factor to the current
+     * matrix.
+     */
+    MatrixBase & add (const MatrixBase &other,
+                     const PetscScalar factor);
+
     /**
      * Matrix-vector multiplication:
      * let <i>dst = M*src</i> with
@@ -1700,7 +1707,7 @@ namespace PETScWrappers
 #ifndef PETSC_USE_64BIT_INDICES
     if (elide_zero_values == false)
       {
-        col_index_ptr = (int *)col_indices;
+        col_index_ptr = (int*)col_indices;
         col_value_ptr = values;
         n_columns = n_cols;
       }

diff --git a/deal.II/source/lac/petsc_matrix_base.cc b/deal.II/source/lac/petsc_matrix_base.cc
index c8a3d456c1a2a2d0687ab7521a8a9b302ec2cb68..f41fbd9f21bb49cf1893028799b2b137c1ef62e6 100644 (file)
--- a/deal.II/source/lac/petsc_matrix_base.cc
+++ b/deal.II/source/lac/petsc_matrix_base.cc
@@ -478,6 +478,19 @@ namespace PETScWrappers
   }
 
 
+  MatrixBase &
+  MatrixBase::add (const MatrixBase &other,
+                  const PetscScalar factor)
+  {
+    const int ierr = MatAXPY (matrix, factor,
+                             other, DIFFERENT_NONZERO_PATTERN);
+
+    Assert (ierr == 0, ExcPETScError(ierr));
+
+    return *this;
+  }
+  
+
   void
   MatrixBase::vmult (VectorBase       &dst,
                      const VectorBase &src) const

diff --git a/tests/petsc/sparse_matrix_01.cc b/tests/petsc/sparse_matrix_01.cc
new file mode 100644 (file)
index 0000000..d91d62e
--- /dev/null
+++ b/tests/petsc/sparse_matrix_01.cc
@@ -0,0 +1,79 @@
+//----------------------------  petsc_sparse_matrix_vector_07.cc  ---------------------------
+//    $Id: sparse_matrix_vector_07.cc 24924 2012-01-25 12:35:17Z kormann $
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 2005 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.
+//
+//----------------------------  petsc_sparse_matrix_vector_07.cc  ---------------------------
+
+
+// check SparseMatrix::add(other, factor)
+
+#include "../tests.h"
+#include <deal.II/lac/petsc_vector.h>
+#include <deal.II/lac/petsc_sparse_matrix.h>
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+
+void test ()
+{
+  const unsigned int s = 10;
+  PETScWrappers::SparseMatrix m(s,s,s);
+  for (unsigned int i=0; i<m.m(); ++i)
+    for (unsigned int j=0; j<=i; ++j)
+        m.set (i,j, i+2*j);
+      
+  m.compress ();
+
+  PETScWrappers::SparseMatrix m2(s,s,s);
+  m2.set(0,1,5.0);
+  for (unsigned int i=0; i<m2.m(); ++i)
+    m2.set(i,i,1.0+ i);
+  m2.compress();
+
+  // we now print the matrix m,
+  // print after adding m2, and then subtract m2 again
+  // to get the original matrix back.
+  
+  deallog << "before: " << m(0,1) << std::endl;
+  for (unsigned int i=0;i<s;++i)
+    deallog << m(i,i) << " ";
+  deallog << std::endl;
+
+  m.add(m2,1.0);
+
+  deallog << "after: " << m(0,1) << std::endl;
+  for (unsigned int i=0;i<s;++i)
+    deallog << m(i,i) << " ";
+  deallog << std::endl;
+
+  m.add(m2,-1.0);
+
+  deallog << "back to original: " << m(0,1) << std::endl;
+  for (unsigned int i=0;i<s;++i)
+    deallog << m(i,i) << " ";
+  deallog << std::endl;
+  
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main (int argc, char **argv)
+{
+  std::ofstream logfile("sparse_matrix_01/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  PetscInitialize(&argc,&argv,0,0);
+  test ();
+  PetscFinalize ();
+}

diff --git a/tests/petsc/sparse_matrix_01/cmp/generic b/tests/petsc/sparse_matrix_01/cmp/generic
new file mode 100644 (file)
index 0000000..0250cad
--- /dev/null
+++ b/tests/petsc/sparse_matrix_01/cmp/generic
@@ -0,0 +1,8 @@
+
+DEAL::before: 0
+DEAL::0 3.00000 6.00000 9.00000 12.0000 15.0000 18.0000 21.0000 24.0000 27.0000 
+DEAL::after: 5.00000
+DEAL::1.00000 5.00000 9.00000 13.0000 17.0000 21.0000 25.0000 29.0000 33.0000 37.0000 
+DEAL::back to original: 0
+DEAL::0 3.00000 6.00000 9.00000 12.0000 15.0000 18.0000 21.0000 24.0000 27.0000 
+DEAL::OK