*/
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;
};
/*@}*/
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
*/
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
/**
* Return a reference to the MPI
* communicator object in use with
- * this vector.
+ * this matrix.
*/
const MPI_Comm & get_mpi_communicator () const;
* communicator object in use with
* this matrix.
*/
- const MPI_Comm & get_mpi_communicator () const;
+ virtual const MPI_Comm & get_mpi_communicator () const;
/** @addtogroup Exceptions
* @{ */
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:
/**
// $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
AssertThrow (ierr == 0, ExcPETScError(ierr));
}
+
+
+ const MPI_Comm &
+ FullMatrix::get_mpi_communicator () const
+ {
+ static const MPI_Comm communicator = MPI_COMM_SELF;
+ return communicator;
+ }
}
}
+
+ 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
// $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
+ 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,
--- /dev/null
+//---------------------------- 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;
+ };
+}
--- /dev/null
+//---------------------------- 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;
+ };
+}
--- /dev/null
+
+DEAL::7.15267e+09 7.15267e+09
+DEAL::32 32
+DEAL::6.84337e+09 6.84337e+09
+DEAL::30 32
+DEAL::OK
--- /dev/null
+
+DEAL::7.15267e+09 7.15267e+09
+DEAL::32 32
+DEAL::6.71272e+09 6.71272e+09
+DEAL::29 32
+DEAL::OK
projects / dealii-svn.git / commitdiff