--- /dev/null
+//---------------------------- 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 ---------------------------*/
--- /dev/null
+//---------------------------------------------------------------------------
+// $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
projects / dealii-svn.git / commitdiff