]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Start of parallel matrices.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 24 Mar 2004 00:30:09 +0000 (00:30 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 24 Mar 2004 00:30:09 +0000 (00:30 +0000)
git-svn-id: https://svn.dealii.org/trunk@8857 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/petsc_parallel_sparse_matrix.h [new file with mode: 0644]
deal.II/lac/source/petsc_parallel_sparse_matrix.cc [new file with mode: 0644]

diff --git a/deal.II/lac/include/lac/petsc_parallel_sparse_matrix.h b/deal.II/lac/include/lac/petsc_parallel_sparse_matrix.h
new file mode 100644 (file)
index 0000000..41cb105
--- /dev/null
@@ -0,0 +1,197 @@
+//----------------------------  petsc_parallel_sparse_matrix.h  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2004 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_parallel_sparse_matrix.h  ---------------------------
+#ifndef __deal2__petsc_parallel_sparse_matrix_h
+#define __deal2__petsc_parallel_sparse_matrix_h
+
+
+#include <base/config.h>
+#include <base/exceptions.h>
+
+#ifdef DEAL_II_USE_PETSC
+
+#include <lac/petsc_matrix_base.h>
+
+#include <vector>
+
+
+//TODO: How does a matrix or vector get its communicator if zero-initialized and later reinit'ed?
+
+namespace PETScWrappers
+{
+  namespace MPI
+  {
+    
+/**
+ * Implementation of a parallel sparse matrix class based on PETSC, with rows
+ * of the matrix distributed across an MPI network. All the functionality is
+ * actually in the base class, except for the calls to generate a parallel
+ * sparse matrix. This is possible since PETSc only works on an abstract
+ * matrix type and internally distributes to functions that do the actual work
+ * depending on the actual matrix type (much like using virtual
+ * functions). Only the functions creating a matrix of specific type differ,
+ * and are implemented in this particular class.
+ *
+ * @author Wolfgang Bangerth, 2004
+ */
+    class SparseMatrix : public MatrixBase
+    {
+      public:
+                                         /**
+                                          * Default constructor. Create an empty
+                                          * matrix.
+                                          */
+        SparseMatrix ();
+      
+                                         /**
+                                          * Create a sparse matrix of dimensions
+                                          * @arg m times @arg n, with an
+                                          * initial guess of
+                                          * @arg n_nonzero_per_row nonzero
+                                          * elements per row. PETSc is able to
+                                          * cope with the situation that more
+                                          * than this number of elements is
+                                          * later allocated for a row, but this
+                                          * involves copying data, and is thus
+                                          * expensive.
+                                          *
+                                          * The @p{is_symmetric} flag determines
+                                          * whether we should tell PETSc that
+                                          * the matrix is going to be symmetric
+                                          * (as indicated by the call
+                                          * @p{MatSetOption(mat,MAT_SYMMETRIC)}. Note
+                                          * that the PETSc documentation states
+                                          * that one cannot form an ILU
+                                          * decomposition of a matrix for which
+                                          * this flag has been set. The default
+                                          * value of this flag is @p{false}.
+                                          */
+        SparseMatrix (const MPI_Comm    &communicator,
+                      const unsigned int m,
+                      const unsigned int n,
+                      const unsigned int local_rows,
+                      const unsigned int n_nonzero_per_row,
+                      const bool         is_symmetric = false);
+
+                                         /**
+                                          * Initialize a rectangular matrix with
+                                          * @arg m rows and @arg n
+                                          * columns.  The maximal number of
+                                          * nonzero entries for each row
+                                          * separately is given by the
+                                          * @arg row_lengths array.
+                                          *
+                                          * Just as for the other constructors:
+                                          * PETSc is able to cope with the
+                                          * situation that more than this number
+                                          * of elements is later allocated for a
+                                          * row, but this involves copying data,
+                                          * and is thus expensive.
+                                          *
+                                          * The @p{is_symmetric} flag determines
+                                          * whether we should tell PETSc that
+                                          * the matrix is going to be symmetric
+                                          * (as indicated by the call
+                                          * @p{MatSetOption(mat,MAT_SYMMETRIC)}. Note
+                                          * that the PETSc documentation states
+                                          * that one cannot form an ILU
+                                          * decomposition of a matrix for which
+                                          * this flag has been set. The default
+                                          * value of this flag is @p{false}.
+                                          */
+        SparseMatrix (const MPI_Comm                  &communicator,
+                      const unsigned int               m,
+                      const unsigned int               n,
+                      const unsigned int local_rows,
+                      const std::vector<unsigned int> &row_lengths,
+                      const bool                       is_symmetric = false);
+
+                                         /**
+                                          * Set all matrix entries to zero, but
+                                          * retain the sparsity pattern. This
+                                          * function simply calls the respective
+                                          * function of the base class.
+                                          */
+        void reinit ();
+
+                                         /**
+                                          * Throw away the present matrix and
+                                          * generate one that has the same
+                                          * properties as if it were created by
+                                          * the constructor of this class with
+                                          * the same argument list as the
+                                          * present function.
+                                          */
+        void reinit (const MPI_Comm    &communicator,
+                     const unsigned int m,
+                     const unsigned int n,
+                     const unsigned int local_rows,
+                     const unsigned int n_nonzero_per_row,
+                     const bool         is_symmetric = false);
+
+                                         /**
+                                          * Throw away the present matrix and
+                                          * generate one that has the same
+                                          * properties as if it were created by
+                                          * the constructor of this class with
+                                          * the same argument list as the
+                                          * present function.
+                                          */
+        void reinit (const MPI_Comm                  &communicator,
+                     const unsigned int               m,
+                     const unsigned int               n,
+                     const unsigned int               local_rows,
+                     const std::vector<unsigned int> &row_lengths,
+                     const bool                       is_symmetric = false);
+
+      private:
+
+                                         /**
+                                          * Do the actual work for the
+                                          * respective reinit() function and the
+                                          * matching constructor, i.e. create a
+                                          * matrix. Getting rid of the previous
+                                          * matrix is left to the caller.
+                                          */
+        void do_reinit (const unsigned int m,
+                        const unsigned int n,
+                        const unsigned int local_rows,
+                        const unsigned int n_nonzero_per_row,
+                        const bool         is_symmetric = false);
+
+                                         /**
+                                          * Same as previous function.
+                                          */
+        void do_reinit (const unsigned int               m,
+                        const unsigned int               n,
+                        const unsigned int local_rows,
+                        const std::vector<unsigned int> &row_lengths,
+                        const bool                       is_symmetric = false);
+
+      private:
+                                         /**
+                                          * Copy of the communicator object to
+                                          * be used for this parallel vector.
+                                          */
+        MPI_Comm communicator;
+    };
+
+  }
+  
+}
+
+#endif // DEAL_II_USE_PETSC
+
+/*----------------------------   petsc_parallel_sparse_matrix.h     ---------------------------*/
+
+#endif
+/*----------------------------   petsc_parallel_sparse_matrix.h     ---------------------------*/
diff --git a/deal.II/lac/source/petsc_parallel_sparse_matrix.cc b/deal.II/lac/source/petsc_parallel_sparse_matrix.cc
new file mode 100644 (file)
index 0000000..a445c2a
--- /dev/null
@@ -0,0 +1,191 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2004 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.
+//
+//---------------------------------------------------------------------------
+
+
+#include <lac/petsc_parallel_sparse_matrix.h>
+#include <lac/petsc_vector.h>
+
+#ifdef DEAL_II_USE_PETSC
+
+
+namespace PETScWrappers
+{
+  namespace MPI
+  {
+    
+    SparseMatrix::SparseMatrix ()
+    {
+                                       // just like for vectors: since we
+                                       // create an empty matrix, we can as
+                                       // well make it sequential
+      const int m=0, n=0, n_nonzero_per_row=0;
+      const int ierr
+        = MatCreateSeqAIJ(PETSC_COMM_SELF, m, n, n_nonzero_per_row,
+                          0, &matrix);
+      AssertThrow (ierr == 0, ExcPETScError(ierr));
+    }
+
+
+
+    SparseMatrix::SparseMatrix (const MPI_Comm    &communicator,
+                                const unsigned int m,
+                                const unsigned int n,
+                                const unsigned int local_rows,
+                                const unsigned int n_nonzero_per_row,
+                                const bool         is_symmetric)
+                    :
+                    communicator (communicator)
+    {
+      do_reinit (m, n, local_rows, n_nonzero_per_row, is_symmetric);
+    }
+
+
+
+    SparseMatrix::SparseMatrix (const MPI_Comm    &communicator,
+                                const unsigned int m,
+                                const unsigned int n,
+                                const unsigned int local_rows,
+                                const std::vector<unsigned int> &row_lengths,
+                                const bool         is_symmetric)
+                    :
+                    communicator (communicator)
+    {
+      do_reinit (m, n, local_rows, row_lengths, is_symmetric);
+    }
+
+
+
+    void
+    SparseMatrix::reinit (const MPI_Comm    &communicator,
+                          const unsigned int m,
+                          const unsigned int n,
+                          const unsigned int local_rows,
+                          const unsigned int n_nonzero_per_row,
+                          const bool         is_symmetric)
+    {
+      this->communicator = communicator;
+
+                                       // get rid of old matrix and generate a
+                                       // new one
+      const int ierr = MatDestroy (matrix);
+      AssertThrow (ierr == 0, ExcPETScError(ierr));    
+    
+      do_reinit (m, n, local_rows, n_nonzero_per_row, is_symmetric);
+    }
+
+
+
+    void
+    SparseMatrix::reinit (const MPI_Comm    &communicator,
+                          const unsigned int m,
+                          const unsigned int n,
+                          const unsigned int local_rows,
+                          const std::vector<unsigned int> &row_lengths,
+                          const bool         is_symmetric)
+    {
+      this->communicator = communicator;
+      
+                                       // get rid of old matrix and generate a
+                                       // new one
+      const int ierr = MatDestroy (matrix);
+      AssertThrow (ierr == 0, ExcPETScError(ierr));    
+
+      do_reinit (m, n, local_rows, row_lengths, is_symmetric);
+    }  
+
+
+
+    void
+    SparseMatrix::do_reinit (const unsigned int m,
+                             const unsigned int n,
+                             const unsigned int local_rows,
+                             const unsigned int n_nonzero_per_row,
+                             const bool         is_symmetric)
+    {
+                                       // use the call sequence indicating only
+                                       // a maximal number of elements per row
+                                       // for all rows globally
+                                       //
+                                       // note that we partition the columns
+                                       // in the same way as we partition the
+                                       // rows, i.e. we assume that users
+                                       // specify the same number for
+                                       // local_rows as they do for local_size
+                                       // on parallel vectors
+      const int ierr
+        = MatCreateMPIAIJ(communicator, local_rows, local_rows, m, n,
+                          n_nonzero_per_row, 0, 0, 0,
+                          &matrix);
+      AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+                                       // set symmetric flag, if so requested
+      if (is_symmetric == true)
+        {
+          const int ierr
+            = MatSetOption (matrix, MAT_SYMMETRIC);
+          AssertThrow (ierr == 0, ExcPETScError(ierr));
+        }
+    }
+
+
+
+    void
+    SparseMatrix::do_reinit (const unsigned int m,
+                             const unsigned int n,
+                             const unsigned int local_rows,
+                             const std::vector<unsigned int> &row_lengths,
+                             const bool         is_symmetric)
+    {
+      Assert (row_lengths.size() == m,
+              ExcDimensionMismatch (row_lengths.size(), m));
+    
+                                       // use the call sequence indicating a
+                                       // maximal number of elements for each
+                                       // row individually. annoyingly, we
+                                       // always use unsigned ints for cases
+                                       // like this, while PETSc wants to see
+                                       // signed integers. so we have to
+                                       // convert, unless we want to play dirty
+                                       // tricks with conversions of pointers
+      const std::vector<signed int> int_row_lengths (row_lengths.begin(),
+                                                     row_lengths.end());
+                                       
+                                       // note that we partition the columns
+                                       // in the same way as we partition the
+                                       // rows, i.e. we assume that users
+                                       // specify the same number for
+                                       // local_rows as they do for local_size
+                                       // on parallel vectors
+      const int ierr
+        = MatCreateMPIAIJ(communicator, local_rows, local_rows, m, n,
+                          0, &int_row_lengths[0], 0, 0,
+                          &matrix);
+      AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+                                       // set symmetric flag, if so requested
+      if (is_symmetric == true)
+        {
+          const int ierr
+            = MatSetOption (matrix, MAT_SYMMETRIC);
+          AssertThrow (ierr == 0, ExcPETScError(ierr));
+        }    
+    }
+  }
+}
+
+
+#else
+// On gcc2.95 on Alpha OSF1, the native assembler does not like empty
+// files, so provide some dummy code
+namespace { void dummy () {} }
+#endif // DEAL_II_USE_PETSC

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.