From: wolf Date: Tue, 9 Aug 2005 18:18:12 +0000 (+0000) Subject: Give the matrix base class a way to access the MPI communicator. Implement clear_row... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b9e5b29e038ed6ab4e54cd12ddbc0a313c443399;p=dealii-svn.git Give the matrix base class a way to access the MPI communicator. Implement clear_row and clear_rows for the matrix classes. git-svn-id: https://svn.dealii.org/trunk@11267 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/petsc_full_matrix.h b/deal.II/lac/include/lac/petsc_full_matrix.h index 8d6006a3cd..112bd9f99d 100644 --- a/deal.II/lac/include/lac/petsc_full_matrix.h +++ b/deal.II/lac/include/lac/petsc_full_matrix.h @@ -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; }; /*@}*/ diff --git a/deal.II/lac/include/lac/petsc_matrix_base.h b/deal.II/lac/include/lac/petsc_matrix_base.h index 643b7642bb..5731f48d98 100644 --- a/deal.II/lac/include/lac/petsc_matrix_base.h +++ b/deal.II/lac/include/lac/petsc_matrix_base.h @@ -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 &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 diff --git a/deal.II/lac/include/lac/petsc_parallel_block_sparse_matrix.h b/deal.II/lac/include/lac/petsc_parallel_block_sparse_matrix.h index 2d088fd343..f9f4029987 100644 --- a/deal.II/lac/include/lac/petsc_parallel_block_sparse_matrix.h +++ b/deal.II/lac/include/lac/petsc_parallel_block_sparse_matrix.h @@ -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; diff --git a/deal.II/lac/include/lac/petsc_parallel_sparse_matrix.h b/deal.II/lac/include/lac/petsc_parallel_sparse_matrix.h index 345eb7ab03..b92229e9dd 100644 --- a/deal.II/lac/include/lac/petsc_parallel_sparse_matrix.h +++ b/deal.II/lac/include/lac/petsc_parallel_sparse_matrix.h @@ -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 * @{ */ diff --git a/deal.II/lac/include/lac/petsc_sparse_matrix.h b/deal.II/lac/include/lac/petsc_sparse_matrix.h index 4c8063e44c..4708f792af 100644 --- a/deal.II/lac/include/lac/petsc_sparse_matrix.h +++ b/deal.II/lac/include/lac/petsc_sparse_matrix.h @@ -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: /** diff --git a/deal.II/lac/source/petsc_full_matrix.cc b/deal.II/lac/source/petsc_full_matrix.cc index 46db710ca6..8eb1545643 100644 --- a/deal.II/lac/source/petsc_full_matrix.cc +++ b/deal.II/lac/source/petsc_full_matrix.cc @@ -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; + } } diff --git a/deal.II/lac/source/petsc_matrix_base.cc b/deal.II/lac/source/petsc_matrix_base.cc index 30e80aa1d3..9243c06410 100644 --- a/deal.II/lac/source/petsc_matrix_base.cc +++ b/deal.II/lac/source/petsc_matrix_base.cc @@ -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 &rows) + { + compress (); + + // now set all the entries of these rows + // to zero + IS index_set; + const std::vector 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 diff --git a/deal.II/lac/source/petsc_sparse_matrix.cc b/deal.II/lac/source/petsc_sparse_matrix.cc index 701b7bdb84..0dbe4e6b5d 100644 --- a/deal.II/lac/source/petsc_sparse_matrix.cc +++ b/deal.II/lac/source/petsc_sparse_matrix.cc @@ -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 index 0000000000..f139763f1f --- /dev/null +++ b/tests/bits/petsc_66.cc @@ -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 +#include + +#include +#include +#include + + +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=5) + { + const double s = rand(); + m.add (i,i-5, s); + norm_sqr += s*s; + ++nnz; + } + + if (i +#include + +#include +#include +#include + + +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=5) + { + const double s = rand(); + m.add (i,i-5, s); + norm_sqr += s*s; + ++nnz; + } + + if (i(&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 index 0000000000..3a786b7f71 --- /dev/null +++ b/tests/results/i686-pc-linux-gnu+gcc3.2/bits/petsc_66.output @@ -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 index 0000000000..8a0f539840 --- /dev/null +++ b/tests/results/i686-pc-linux-gnu+gcc3.2/bits/petsc_67.output @@ -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