]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Give the matrix base class a way to access the MPI communicator. Implement clear_row...
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 9 Aug 2005 18:18:12 +0000 (18:18 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 9 Aug 2005 18:18:12 +0000 (18:18 +0000)
git-svn-id: https://svn.dealii.org/trunk@11267 0785d39b-7218-0410-832d-ea1e28bc413d

12 files changed:
deal.II/lac/include/lac/petsc_full_matrix.h
deal.II/lac/include/lac/petsc_matrix_base.h
deal.II/lac/include/lac/petsc_parallel_block_sparse_matrix.h
deal.II/lac/include/lac/petsc_parallel_sparse_matrix.h
deal.II/lac/include/lac/petsc_sparse_matrix.h
deal.II/lac/source/petsc_full_matrix.cc
deal.II/lac/source/petsc_matrix_base.cc
deal.II/lac/source/petsc_sparse_matrix.cc
tests/bits/petsc_66.cc [new file with mode: 0644]
tests/bits/petsc_67.cc [new file with mode: 0644]
tests/results/i686-pc-linux-gnu+gcc3.2/bits/petsc_66.output [new file with mode: 0644]
tests/results/i686-pc-linux-gnu+gcc3.2/bits/petsc_67.output [new file with mode: 0644]

index 8d6006a3cd8415d9440531116692b0f66b3e3902..112bd9f99de4dbd4d9ffc480674cf6892de4d211 100644 (file)
@@ -48,6 +48,15 @@ namespace PETScWrappers
                                         */
       FullMatrix (const unsigned int m,
                   const unsigned int n);
+
+                                       /**
+                                        * Return a reference to the MPI
+                                        * communicator object in use with this
+                                        * matrix. Since this is a sequential
+                                        * matrix, it returns the MPI_COMM_SELF
+                                        * communicator.
+                                        */
+      virtual const MPI_Comm & get_mpi_communicator () const;
   };
   
 /*@}*/
index 643b7642bba166e0cc320e1894147e2218490533..5731f48d98ecd6545ef11f7e7a737a71ec2814f4 100644 (file)
@@ -362,6 +362,36 @@ namespace PETScWrappers
                 const unsigned int j,
                 const PetscScalar value);
 
+                                       /**
+                                        * Remove all elements from this #row
+                                        * by setting them to zero. The
+                                        * function does not modify the number
+                                        * of allocated nonzero entries, it
+                                        * only sets some entries to zero. It
+                                        * may drop them from the sparsity
+                                        * pattern, though (but retains the
+                                        * allocated memory in case new entries
+                                        * are again added later).
+                                        *
+                                        * This operation is used in
+                                        * eliminating constraints (e.g. due to
+                                        * hanging nodes) and makes sure that
+                                        * we can write this modification to
+                                        * the matrix without having to read
+                                        * entries (such as the locations of
+                                        * non-zero elements) from it --
+                                        * without this operation, removing
+                                        * constraints on parallel matrices is
+                                        * a rather complicated procedure.
+                                        */
+      void clear_row (const unsigned int row);
+
+                                       /**
+                                        * Same as clear_row(), except that it
+                                        * works on a number of rows at once.
+                                        */
+      void clear_rows (const std::vector<unsigned int> &rows);
+      
                                        /**
                                         * PETSc matrices store their own
                                         * sparsity patterns. So, in analogy to
@@ -480,6 +510,14 @@ namespace PETScWrappers
                                        */
       bool in_local_range (const unsigned int index) const;
 
+                                       /**
+                                        * Return a reference to the MPI
+                                        * communicator object in use with this
+                                        * matrix. This function has to be
+                                        * implemented in derived classes.
+                                        */
+      virtual const MPI_Comm & get_mpi_communicator () const = 0;
+
                                        /**
                                         * Return the number of nonzero
                                         * elements of this
index 2d088fd343118db553a9928c49008753c27bb454..f9f4029987594d02d42bea782be8c7894258355f 100644 (file)
@@ -264,7 +264,7 @@ namespace PETScWrappers
                                          /**
                                           * Return a reference to the MPI
                                           * communicator object in use with
-                                          * this vector.
+                                          * this matrix.
                                           */
         const MPI_Comm & get_mpi_communicator () const;
 
index 345eb7ab037544545ae5f9f5af63de1f4b557802..b92229e9ddd216ebdb0c3cb9c847fb50134835cb 100644 (file)
@@ -365,7 +365,7 @@ namespace PETScWrappers
                                           * communicator object in use with
                                           * this matrix.
                                           */
-        const MPI_Comm & get_mpi_communicator () const;
+        virtual const MPI_Comm & get_mpi_communicator () const;
         
                                     /** @addtogroup Exceptions
                                      * @{ */
index 4c8063e44c681fd5614e4fba850866f841807837..4708f792afd23def340b51b47190cfc99c06908c 100644 (file)
@@ -277,6 +277,15 @@ namespace PETScWrappers
       void reinit (const CompressedSparsityPattern &sparsity_pattern,
                    const bool                       preset_nonzero_locations = false);
 
+                                       /**
+                                        * Return a reference to the MPI
+                                        * communicator object in use with this
+                                        * matrix. Since this is a sequential
+                                        * matrix, it returns the MPI_COMM_SELF
+                                        * communicator.
+                                        */
+      virtual const MPI_Comm & get_mpi_communicator () const;
+
     private:
 
                                        /**
index 46db710ca62832766bb7ad8da763b02502eedbe4..8eb15456430c7bbde5ce8db4cb3dc8c35c487f4c 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2004 by the deal.II authors
+//    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
@@ -29,6 +29,14 @@ namespace PETScWrappers
 
     AssertThrow (ierr == 0, ExcPETScError(ierr));
   }
+
+
+  const MPI_Comm &
+  FullMatrix::get_mpi_communicator () const
+  {
+    static const MPI_Comm communicator = MPI_COMM_SELF;
+    return communicator;
+  }
 }
 
 
index 30e80aa1d35148f7961e4511f653e1b41fbedb58..9243c06410c01ac9d8469760a3821d685edee2b5 100644 (file)
@@ -196,6 +196,56 @@ namespace PETScWrappers
   }
 
 
+
+  void
+  MatrixBase::clear_row (const unsigned int row)
+  {
+    compress ();
+
+                                     // now set all the entries of this row to
+                                     // zero
+    IS index_set;
+    const PetscInt petsc_row = row;
+    ISCreateGeneral (get_mpi_communicator(), 1, &petsc_row, &index_set);
+    
+    static const PetscScalar zero = 0;
+    const int ierr
+      = MatZeroRows(matrix, index_set, &zero);
+    
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+    compress ();
+  }
+
+
+
+  void
+  MatrixBase::clear_rows (const std::vector<unsigned int> &rows)
+  {
+    compress ();
+
+                                     // now set all the entries of these rows
+                                     // to zero
+    IS index_set;
+    const std::vector<PetscInt> petsc_rows (rows.begin(), rows.end());
+
+                                     // call the functions. note that we have
+                                     // to call them even if #rows is empty,
+                                     // since this is a collective operation
+    ISCreateGeneral (get_mpi_communicator(), rows.size(),
+                     &petsc_rows[0], &index_set);
+    
+    static const PetscScalar zero = 0;
+    const int ierr
+      = MatZeroRows(matrix, index_set, &zero);
+    
+    AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+    compress ();
+  }
+  
+  
+
   PetscScalar
   MatrixBase::el (const unsigned int i,
                   const unsigned int j) const
index 701b7bdb8457f5161dbeb2be8bb8266a5a656bd7..0dbe4e6b5d589ed03830d3f603cf44037aa89ad2 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2004 by the deal.II authors
+//    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
@@ -115,6 +115,15 @@ namespace PETScWrappers
 
 
   
+  const MPI_Comm &
+  SparseMatrix::get_mpi_communicator () const
+  {
+    static const MPI_Comm communicator = MPI_COMM_SELF;
+    return communicator;
+  }
+
+
+
   void
   SparseMatrix::do_reinit (const unsigned int m,
                            const unsigned int n,
diff --git a/tests/bits/petsc_66.cc b/tests/bits/petsc_66.cc
new file mode 100644 (file)
index 0000000..f139763
--- /dev/null
@@ -0,0 +1,135 @@
+//----------------------------  petsc_66.cc  ---------------------------
+//    $Id$
+//    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_66.cc  ---------------------------
+
+
+// check PETScWrappers::MatrixBase::clear_row ()
+
+#include "../tests.h"
+#include <lac/petsc_sparse_matrix.h>
+#include <lac/vector.h>
+
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+
+void test (PETScWrappers::MatrixBase &m)
+{
+  Assert (m.m() != 0, ExcInternalError());
+  Assert (m.n() != 0, ExcInternalError());
+
+                                   // build a tri-diagonal pattern
+  double norm_sqr = 0;
+  unsigned int nnz = 0;
+  const unsigned int N = m.m();
+  for (unsigned int i=0; i<N; ++i)
+    {
+      if (i>=5)
+        {
+          const double s = rand();
+          m.add (i,i-5, s);
+          norm_sqr += s*s;
+          ++nnz;
+        }
+      
+      if (i<N-5)
+        {
+          const double s = rand();
+          m.add (i,i+5, s);
+          norm_sqr += s*s;
+          ++nnz;
+        }
+      
+      const double s = rand();
+      m.add (i,i,s);
+      norm_sqr += s*s;
+      ++nnz;
+    }
+  m.compress ();
+  
+  deallog << m.frobenius_norm() << ' ' << std::sqrt (norm_sqr)
+          << std::endl;
+  deallog << m.n_nonzero_elements() << ' ' << nnz << std::endl;
+
+  Assert (std::fabs (m.frobenius_norm() - std::sqrt(norm_sqr))
+          < std::fabs (std::sqrt (norm_sqr)),
+          ExcInternalError());
+  Assert (m.n_nonzero_elements()-nnz == 0, ExcInternalError());
+
+                                   // now remove the entries of row N/2
+  for (unsigned int i=0; i<N; ++i)
+    {
+      const double s = m.el(N/2,i);
+      norm_sqr -= s*s;
+    }
+  m.clear_row (N/2);
+  
+  deallog << m.frobenius_norm() << ' ' << std::sqrt (norm_sqr)
+          << std::endl;
+  deallog << m.n_nonzero_elements() << ' ' << nnz << std::endl;
+
+  Assert (std::fabs (m.frobenius_norm() - std::sqrt(norm_sqr))
+          < std::fabs (std::sqrt (norm_sqr)),
+          ExcInternalError());
+
+                                   // make sure that zeroing out rows does at
+                                   // least not add new nonzero entries (it
+                                   // may remove some, though)
+  Assert (m.n_nonzero_elements() <= nnz, ExcInternalError());
+  
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main (int argc,char **argv) 
+{
+  std::ofstream logfile("petsc_66.output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  try
+    {
+      PetscInitialize(&argc,&argv,0,0);
+      {
+        PETScWrappers::SparseMatrix v (14,14,3);
+        test (v);
+      }
+      PetscFinalize();
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+               << exc.what() << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      
+      return 1;
+    }
+  catch (...) 
+    {
+      std::cerr << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      return 1;
+    };
+}
diff --git a/tests/bits/petsc_67.cc b/tests/bits/petsc_67.cc
new file mode 100644 (file)
index 0000000..002745f
--- /dev/null
@@ -0,0 +1,141 @@
+//----------------------------  petsc_67.cc  ---------------------------
+//    $Id$
+//    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_67.cc  ---------------------------
+
+
+// check PETScWrappers::MatrixBase::clear_rows ()
+
+#include "../tests.h"
+#include <lac/petsc_sparse_matrix.h>
+#include <lac/vector.h>
+
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+
+void test (PETScWrappers::MatrixBase &m)
+{
+  Assert (m.m() != 0, ExcInternalError());
+  Assert (m.n() != 0, ExcInternalError());
+
+                                   // build a tri-diagonal pattern
+  double norm_sqr = 0;
+  unsigned int nnz = 0;
+  const unsigned int N = m.m();
+  for (unsigned int i=0; i<N; ++i)
+    {
+      if (i>=5)
+        {
+          const double s = rand();
+          m.add (i,i-5, s);
+          norm_sqr += s*s;
+          ++nnz;
+        }
+      
+      if (i<N-5)
+        {
+          const double s = rand();
+          m.add (i,i+5, s);
+          norm_sqr += s*s;
+          ++nnz;
+        }
+      
+      const double s = rand();
+      m.add (i,i,s);
+      norm_sqr += s*s;
+      ++nnz;
+    }
+  m.compress ();
+  
+  deallog << m.frobenius_norm() << ' ' << std::sqrt (norm_sqr)
+          << std::endl;
+  deallog << m.n_nonzero_elements() << ' ' << nnz << std::endl;
+
+  Assert (std::fabs (m.frobenius_norm() - std::sqrt(norm_sqr))
+          < std::fabs (std::sqrt (norm_sqr)),
+          ExcInternalError());
+  Assert (m.n_nonzero_elements()-nnz == 0, ExcInternalError());
+
+                                   // now remove the entries of rows N/2 and N/3
+  for (unsigned int i=0; i<N; ++i)
+    {
+      const double s = m.el(N/2,i);
+      norm_sqr -= s*s;
+    }
+  for (unsigned int i=0; i<N; ++i)
+    {
+      const double s = m.el(N/3,i);
+      norm_sqr -= s*s;
+    }
+  const unsigned int rows[2] = { N/3, N/2 };
+  m.clear_rows (std::vector<unsigned int>(&rows[0], &rows[2]));
+  
+  deallog << m.frobenius_norm() << ' ' << std::sqrt (norm_sqr)
+          << std::endl;
+  deallog << m.n_nonzero_elements() << ' ' << nnz << std::endl;
+
+  Assert (std::fabs (m.frobenius_norm() - std::sqrt(norm_sqr))
+          < std::fabs (std::sqrt (norm_sqr)),
+          ExcInternalError());
+
+                                   // make sure that zeroing out rows does at
+                                   // least not add new nonzero entries (it
+                                   // may remove some, though)
+  Assert (m.n_nonzero_elements() <= nnz, ExcInternalError());
+  
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main (int argc,char **argv) 
+{
+  std::ofstream logfile("petsc_67.output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  try
+    {
+      PetscInitialize(&argc,&argv,0,0);
+      {
+        PETScWrappers::SparseMatrix v (14,14,3);
+        test (v);
+      }
+      PetscFinalize();
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+               << exc.what() << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      
+      return 1;
+    }
+  catch (...) 
+    {
+      std::cerr << std::endl << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+               << "Aborting!" << std::endl
+               << "----------------------------------------------------"
+               << std::endl;
+      return 1;
+    };
+}
diff --git a/tests/results/i686-pc-linux-gnu+gcc3.2/bits/petsc_66.output b/tests/results/i686-pc-linux-gnu+gcc3.2/bits/petsc_66.output
new file mode 100644 (file)
index 0000000..3a786b7
--- /dev/null
@@ -0,0 +1,6 @@
+
+DEAL::7.15267e+09 7.15267e+09
+DEAL::32 32
+DEAL::6.84337e+09 6.84337e+09
+DEAL::30 32
+DEAL::OK
diff --git a/tests/results/i686-pc-linux-gnu+gcc3.2/bits/petsc_67.output b/tests/results/i686-pc-linux-gnu+gcc3.2/bits/petsc_67.output
new file mode 100644 (file)
index 0000000..8a0f539
--- /dev/null
@@ -0,0 +1,6 @@
+
+DEAL::7.15267e+09 7.15267e+09
+DEAL::32 32
+DEAL::6.71272e+09 6.71272e+09
+DEAL::29 32
+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.