From df1d1e35c7574f8e3354c9a12d45b8ecea4a66ac Mon Sep 17 00:00:00 2001 From: wolf Date: Wed, 24 Mar 2004 00:30:09 +0000 Subject: [PATCH] Start of parallel matrices. git-svn-id: https://svn.dealii.org/trunk@8857 0785d39b-7218-0410-832d-ea1e28bc413d --- .../lac/petsc_parallel_sparse_matrix.h | 197 ++++++++++++++++++ .../source/petsc_parallel_sparse_matrix.cc | 191 +++++++++++++++++ 2 files changed, 388 insertions(+) create mode 100644 deal.II/lac/include/lac/petsc_parallel_sparse_matrix.h create mode 100644 deal.II/lac/source/petsc_parallel_sparse_matrix.cc 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 index 0000000000..41cb1051de --- /dev/null +++ b/deal.II/lac/include/lac/petsc_parallel_sparse_matrix.h @@ -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 +#include + +#ifdef DEAL_II_USE_PETSC + +#include + +#include + + +//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 &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 &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 &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 index 0000000000..a445c2ab24 --- /dev/null +++ b/deal.II/lac/source/petsc_parallel_sparse_matrix.cc @@ -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 +#include + +#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 &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 &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 &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 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 -- 2.39.5