Dear John,
On Thu, 9 Oct 2008, John Peterson wrote:
I enclosed all my PETSc calles with "#if PETSC_VERSION_LESS_THAN(2,3,3)" and
added error messages in the "#else" case since I am not familiar with
PETSc's API history. Whenever somebody might need the shell matrix
functionality together with older PETSc versions, he might want to implement
that himself.
OK, that may unnecessarily break your code when Petsc 2.3.4 (or
whatever their next version will be) comes out though. I'd do
instead:
#if PETSC_VERSION_LESS_THAN(2,3,2)
// suitable error or unimplemented message
#else
// your current code
#endif
Well, that's actually almost what I did, except that I compared to
(2,3,3) rather than (2,3,2), since I concluded from other usages of
that macro that it is a "<" comparison rather than a "<=". (I agree
that my mail could lead to the assumption that I had done it the wrong
way around.)
I switched the user API for shell matrix usage from a callback function to a
class, which in fact seems to be much more sensible. I made an abstract
base class "ShellMatrix" and three derived classes. One of the derived
classes just wraps a SparseMatrix. I noticed that SparseMatrix apparently
did not have a method to compute a matrix-vector product, so I added this
functionality to SparseMatrix and PetscMatrix (whereas calling
libmesh_not_implemented() in LaspackMatrix).
This functionality (matrix-vector product) is actually in the
NumericVector class, it's the add_vector routine which takes a
SparseMatrix. Enough people have been confused about this that we
should probably implement it in SparseMatrix as well, but I'd vote for
an implementation in terms of the existing NumericVector one (if
possible) rather than duplicating code.
Okay, done. Also, I did the other way round for ShellMatrix, i.e. I
added NumericVector::add_vector(ShellMatrix,NumericVector) which is
implemented in terms of ShellMatrix::vector_mult_add().
My application that tests the shell matrix is still under
construction.
Best Regards,
Tim
--
Dr. Tim Kroeger Phone +49-421-218-7710
[EMAIL PROTECTED], [EMAIL PROTECTED] Fax +49-421-218-4236
MeVis Research GmbH, Universitaetsallee 29, 28359 Bremen, Germany
Amtsgericht Bremen HRB 16222
Geschaeftsfuehrer: Prof. Dr. H.-O. Peitgen
Index: include/numerics/petsc_matrix.h
===================================================================
--- include/numerics/petsc_matrix.h (Revision 3090)
+++ include/numerics/petsc_matrix.h (Arbeitskopie)
@@ -292,7 +292,6 @@
*/
Mat mat () { libmesh_assert (_mat != NULL); return _mat; }
-
protected:
/**
Index: include/numerics/laspack_matrix.h
===================================================================
--- include/numerics/laspack_matrix.h (Revision 3090)
+++ include/numerics/laspack_matrix.h (Arbeitskopie)
@@ -272,7 +272,6 @@
*/
void print_personal(std::ostream& os=std::cout) const { this->print(os); }
-
private:
/**
Index: include/numerics/linear_solver.h
===================================================================
--- include/numerics/linear_solver.h (Revision 3090)
+++ include/numerics/linear_solver.h (Arbeitskopie)
@@ -37,6 +37,7 @@
template <typename T> class AutoPtr;
template <typename T> class SparseMatrix;
template <typename T> class NumericVector;
+template <typename T> class ShellMatrix;
@@ -140,7 +141,20 @@
const double, // Stopping
tolerance
const unsigned int) = 0; // N.
Iterations
+
+
/**
+ * This function solves a system whose matrix is a shell matrix.
+ */
+ virtual std::pair<unsigned int, Real> solve (const ShellMatrix<T>&
shell_matrix,
+ NumericVector<T>&, // Solution
vector
+ NumericVector<T>&, // RHS vector
+ const double, // Stopping
tolerance
+ const unsigned int) = 0; // N.
Iterations
+
+
+
+ /**
* Prints a useful message about why the latest linear solve
* con(di)verged.
*/
Index: include/numerics/laspack_linear_solver.h
===================================================================
--- include/numerics/laspack_linear_solver.h (Revision 3090)
+++ include/numerics/laspack_linear_solver.h (Arbeitskopie)
@@ -95,6 +95,16 @@
const unsigned int m_its);
/**
+ * This function solves a system whose matrix is a shell matrix.
+ */
+ std::pair<unsigned int, Real>
+ solve (const ShellMatrix<T>& shell_matrix,
+ NumericVector<T>& solution_in,
+ NumericVector<T>& rhs_in,
+ const double tol,
+ const unsigned int m_its);
+
+ /**
* Prints a useful message about why the latest linear solve
* con(di)verged.
*/
Index: include/numerics/petsc_linear_solver.h
===================================================================
--- include/numerics/petsc_linear_solver.h (Revision 3090)
+++ include/numerics/petsc_linear_solver.h (Arbeitskopie)
@@ -122,6 +122,16 @@
const unsigned int m_its);
/**
+ * This function solves a system whose matrix is a shell matrix.
+ */
+ std::pair<unsigned int, Real>
+ solve (const ShellMatrix<T>& shell_matrix,
+ NumericVector<T>& solution_in,
+ NumericVector<T>& rhs_in,
+ const double tol,
+ const unsigned int m_its);
+
+ /**
* Returns the raw PETSc preconditioner context pointer. This allows
* you to specify the PCShellSetApply() and PCShellSetSetUp() functions
* if you desire. Just don't do anything crazy like calling PCDestroy()!
@@ -169,6 +179,12 @@
*/
void set_petsc_preconditioner_type ();
+ /**
+ * Internal matrix multiplication function if shell matrix mode is
+ * used.
+ */
+ static PetscErrorCode _petsc_shell_matrix(Mat mat, Vec arg, Vec dest);
+
// SLES removed from >= PETSc 2.2.0
#if PETSC_VERSION_LESS_THAN(2,2,0)
Index: include/numerics/sparse_matrix.h
===================================================================
--- include/numerics/sparse_matrix.h (Revision 3090)
+++ include/numerics/sparse_matrix.h (Arbeitskopie)
@@ -39,6 +39,7 @@
template <typename T> inline std::ostream& operator << (std::ostream& os,
const SparseMatrix<T>& m);
class DofMap;
namespace SparsityPattern { class Graph; }
+template <typename T> class NumericVector;
/**
@@ -342,6 +343,19 @@
true); // true means REUSE submatrix
}
+ /**
+ * Multiplies the matrix with \p arg and stores the result in \p
+ * dest.
+ */
+ void vector_mult (NumericVector<T>& dest,
+ const NumericVector<T>& arg) const;
+
+ /**
+ * Multiplies the matrix with \p arg and adds the result to \p dest.
+ */
+ void vector_mult_add (NumericVector<T>& dest,
+ const NumericVector<T>& arg) const;
+
protected:
/**
Index: include/numerics/numeric_vector.h
===================================================================
--- include/numerics/numeric_vector.h (Revision 3090)
+++ include/numerics/numeric_vector.h (Arbeitskopie)
@@ -38,6 +38,7 @@
template <typename T> class NumericVector;
template <typename T> class DenseVector;
template <typename T> class SparseMatrix;
+template <typename T> class ShellMatrix;
/**
@@ -309,6 +310,13 @@
const SparseMatrix<T>&) = 0;
/**
+ * \f$U+=A*V\f$, add the product of a \p ShellMatrix \p A
+ * and a \p NumericVector \p V to \p this, where \p this=U.
+ */
+ void add_vector (const NumericVector<T>& v,
+ const ShellMatrix<T>& a);
+
+ /**
* \f$ U+=V \f$ where U and V are type
* DenseVector<T> and you
* want to specify WHERE to add
Index: include/solvers/linear_implicit_system.h
===================================================================
--- include/solvers/linear_implicit_system.h (Revision 3090)
+++ include/solvers/linear_implicit_system.h (Arbeitskopie)
@@ -30,6 +30,7 @@
// Forward Declarations
template <typename T> class LinearSolver;
+template <typename T> class ShellMatrix;
/**
@@ -125,6 +126,22 @@
*/
Real final_linear_residual() const { return _final_linear_residual; }
+ /**
+ * This function enables the user to provide a shell matrix, i.e. a
+ * matrix that is not stored element-wise, but as a function. When
+ * you register your shell matrix using this function, calling \p
+ * solve() will no longer use the \p matrix member but the
+ * registered shell matrix instead. You can reset this behaviour to
+ * its original state by supplying a \p NULL pointer to this
+ * function.
+ */
+ void attach_shell_matrix (ShellMatrix<Number>* shell_matrix);
+
+ /**
+ * Detaches a shell matrix. Same as \p attach_shell_matrix(NULL).
+ */
+ void detach_shell_matrix (void) { attach_shell_matrix(NULL); }
+
protected:
/**
@@ -137,6 +154,12 @@
* The final residual for the linear system Ax=b.
*/
Real _final_linear_residual;
+
+ /**
+ * User supplies shell matrix or \p NULL if no shell matrix is used.
+ */
+ ShellMatrix<Number>* _shell_matrix;
+
};
Index: src/numerics/petsc_linear_solver.C
===================================================================
--- src/numerics/petsc_linear_solver.C (Revision 3090)
+++ src/numerics/petsc_linear_solver.C (Arbeitskopie)
@@ -29,6 +29,7 @@
// Local Includes
#include "libmesh_logging.h"
#include "petsc_linear_solver.h"
+#include "shell_matrix.h"
@@ -413,6 +414,100 @@
template <typename T>
+std::pair<unsigned int, Real>
+PetscLinearSolver<T>::solve (const ShellMatrix<T>& shell_matrix,
+ NumericVector<T>& solution_in,
+ NumericVector<T>& rhs_in,
+ const double tol,
+ const unsigned int m_its)
+{
+
+#if PETSC_VERSION_LESS_THAN(2,3,3)
+
+ std::cout << "This method has been developed with PETSc 2.3.3. No one has
made it backwards compatible with older version of PETSc so far; however, it
might work without any change with some older version." << std::endl;
+ libmesh_error();
+ return std::make_pair(0,0.0);
+
+#else
+
+ START_LOG("solve()", "PetscLinearSolver");
+
+ libmesh_assert(shell_matrix!=NULL);
+
+ PetscVector<T>* solution = dynamic_cast<PetscVector<T>*>(&solution_in);
+ PetscVector<T>* rhs = dynamic_cast<PetscVector<T>*>(&rhs_in);
+
+ // We cast to pointers so we can be sure that they succeeded
+ // by comparing the result against NULL.
+ libmesh_assert(solution != NULL);
+ libmesh_assert(rhs != NULL);
+
+ this->init ();
+
+ int ierr=0;
+ int its=0, max_its = static_cast<int>(m_its);
+ PetscReal final_resid=0.;
+
+ // Close the matrices and vectors in case this wasn't already done.
+ solution->close ();
+ rhs->close ();
+
+ // Prepare the matrix.
+ Mat mat;
+ ierr = MatCreateShell(libMesh::COMM_WORLD,
+ rhs_in.local_size(),
+ solution_in.local_size(),
+ rhs_in.local_size(),
+ solution_in.local_size(),
+ const_cast<void*>(static_cast<const
void*>(&shell_matrix)),
+ &mat);
+ /* Note that the const_cast above is only necessary because PETSc
+ does not accept a const void*. Inside the member function
+ _petsc_shell_matrix() below, the pointer is casted back to a
+ const ShellMatrix<T>*. */
+
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ ierr =
MatShellSetOperation(mat,MATOP_MULT,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix));
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ // Set operators. The input matrix works as the preconditioning matrix
+ ierr = KSPSetOperators(_ksp, mat, mat,
+ SAME_NONZERO_PATTERN);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ // Set the tolerances for the iterative solver. Use the user-supplied
+ // tolerance for the relative residual & leave the others at default values.
+ ierr = KSPSetTolerances (_ksp, tol, PETSC_DEFAULT,
+ PETSC_DEFAULT, max_its);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ // Solve the linear system
+ ierr = KSPSolve (_ksp, rhs->vec(), solution->vec());
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ // Get the number of iterations required for convergence
+ ierr = KSPGetIterationNumber (_ksp, &its);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ // Get the norm of the final residual to return to the user.
+ ierr = KSPGetResidualNorm (_ksp, &final_resid);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ // Destroy the matrix.
+ ierr = MatDestroy(mat);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ STOP_LOG("solve()", "PetscLinearSolver");
+ // return the # of its. and the final residual norm.
+ return std::make_pair(its, final_resid);
+
+#endif
+
+}
+
+
+
+template <typename T>
void PetscLinearSolver<T>::get_residual_history(std::vector<double>& hist)
{
int ierr = 0;
@@ -648,6 +743,29 @@
}
+template <typename T>
+PetscErrorCode PetscLinearSolver<T>::_petsc_shell_matrix(Mat mat, Vec arg, Vec
dest)
+{
+ /* Get the matrix context. */
+ int ierr=0;
+ void* ctx;
+ ierr = MatShellGetContext(mat,&ctx);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ /* Get user shell matrix object. */
+ const ShellMatrix<T>& shell_matrix = *static_cast<const
ShellMatrix<T>*>(ctx);
+
+ /* Make \p NumericVector instances arount the vectors. */
+ PetscVector<T> arg_global(arg);
+ PetscVector<T> dest_global(dest);
+
+ /* Call the user function. */
+ shell_matrix.vector_mult(dest_global,arg_global);
+
+ return ierr;
+}
+
+
//------------------------------------------------------------------
// Explicit instantiations
template class PetscLinearSolver<Number>;
Index: src/numerics/laspack_linear_solver.C
===================================================================
--- src/numerics/laspack_linear_solver.C (Revision 3090)
+++ src/numerics/laspack_linear_solver.C (Arbeitskopie)
@@ -284,6 +284,20 @@
template <typename T>
+std::pair<unsigned int, Real>
+LaspackLinearSolver<T>::solve (const ShellMatrix<T>& shell_matrix,
+ NumericVector<T> &solution_in,
+ NumericVector<T> &rhs_in,
+ const double tol,
+ const unsigned int m_its)
+{
+ libmesh_not_implemented();
+ return std::make_pair(0,0.0);
+}
+
+
+
+template <typename T>
void LaspackLinearSolver<T>::set_laspack_preconditioner_type ()
{
switch (this->_preconditioner_type)
Index: src/numerics/sparse_matrix.C
===================================================================
--- src/numerics/sparse_matrix.C (Revision 3090)
+++ src/numerics/sparse_matrix.C (Arbeitskopie)
@@ -26,6 +26,7 @@
#include "laspack_matrix.h"
#include "petsc_matrix.h"
#include "trilinos_epetra_matrix.h"
+#include "numeric_vector.h"
//------------------------------------------------------------------
@@ -81,6 +82,27 @@
+template <typename T>
+void SparseMatrix<T>::vector_mult (NumericVector<T>& dest,
+ const NumericVector<T>& arg) const
+{
+ dest.zero();
+ this->vector_mult_add(dest,arg);
+}
+
+
+
+template <typename T>
+void SparseMatrix<T>::vector_mult_add (NumericVector<T>& dest,
+ const NumericVector<T>& arg) const
+{
+ /* This functionality is actually implemented in the \p
+ NumericVector class. */
+ dest.add_vector(arg,*this);
+}
+
+
+
//------------------------------------------------------------------
// Explicit instantiations
template class SparseMatrix<Number>;
Index: src/numerics/numeric_vector.C
===================================================================
--- src/numerics/numeric_vector.C (Revision 3090)
+++ src/numerics/numeric_vector.C (Arbeitskopie)
@@ -28,6 +28,7 @@
#include "laspack_vector.h"
#include "petsc_vector.h"
#include "trilinos_epetra_vector.h"
+#include "shell_matrix.h"
@@ -189,6 +190,15 @@
+template <typename T>
+void NumericVector<T>::add_vector (const NumericVector<T>& v,
+ const ShellMatrix<T>& a)
+{
+ a.vector_mult_add(*this,v);
+}
+
+
+
//------------------------------------------------------------------
// Explicit instantiations
template class NumericVector<Number>;
Index: src/numerics/petsc_matrix.C
===================================================================
--- src/numerics/petsc_matrix.C (Revision 3090)
+++ src/numerics/petsc_matrix.C (Arbeitskopie)
@@ -27,6 +27,7 @@
#include "petsc_matrix.h"
#include "dof_map.h"
#include "dense_matrix.h"
+#include "petsc_vector.h"
@@ -373,10 +374,6 @@
-
-
-
-
//------------------------------------------------------------------
// Explicit instantiations
template class PetscMatrix<Number>;
Index: src/numerics/laspack_matrix.C
===================================================================
--- src/numerics/laspack_matrix.C (Revision 3090)
+++ src/numerics/laspack_matrix.C (Arbeitskopie)
@@ -194,6 +194,7 @@
}
+
//------------------------------------------------------------------
// Explicit instantiations
template class LaspackMatrix<Number>;
Index: src/solvers/linear_implicit_system.C
===================================================================
--- src/solvers/linear_implicit_system.C (Revision 3090)
+++ src/solvers/linear_implicit_system.C (Arbeitskopie)
@@ -37,7 +37,8 @@
Parent (es, name, number),
linear_solver (LinearSolver<Number>::build()),
_n_linear_iterations (0),
- _final_linear_residual (1.e20)
+ _final_linear_residual (1.e20),
+ _shell_matrix(NULL)
{
}
@@ -95,13 +96,23 @@
const unsigned int maxits =
es.parameters.get<unsigned int>("linear solver maximum iterations");
- // Solve the linear system. Two cases:
- const std::pair<unsigned int, Real> rval =
- (this->have_matrix("Preconditioner")) ?
- // 1.) User-supplied preconditioner
- linear_solver->solve (*matrix, this->get_matrix("Preconditioner"),
*solution, *rhs, tol, maxits) :
- // 2.) Use system matrix for the preconditioner
- linear_solver->solve (*matrix, *solution, *rhs, tol, maxits);
+ // Solve the linear system. Three cases:
+ std::pair<unsigned int, Real> rval = std::make_pair(0,0.0);
+ if(_shell_matrix!=NULL)
+ {
+ // 1.) Shell matrix. Use *this as the context.
+ rval = linear_solver->solve(*_shell_matrix, *solution, *rhs, tol,
maxits);
+ }
+ if(this->have_matrix("Preconditioner"))
+ {
+ // 2.) No shell matrix, but with user-supplied preconditioner
+ rval = linear_solver->solve (*matrix,
this->get_matrix("Preconditioner"), *solution, *rhs, tol, maxits);
+ }
+ else
+ {
+ // 3.) No shell matrix, and use system matrix for the preconditioner
+ rval = linear_solver->solve (*matrix, *solution, *rhs, tol, maxits);
+ }
// Store the number of linear iterations required to
// solve and the final residual.
@@ -115,3 +126,13 @@
// Update the system after the solve
this->update();
}
+
+
+
+void LinearImplicitSystem::attach_shell_matrix (ShellMatrix<Number>*
shell_matrix)
+{
+ _shell_matrix = shell_matrix;
+}
+
+
+
Index: .depend
===================================================================
--- .depend (Revision 3090)
+++ .depend (Arbeitskopie)
@@ -6536,6 +6536,7 @@
include/numerics/petsc_macro.h\
include/numerics/petsc_matrix.h\
include/numerics/petsc_vector.h\
+ include/numerics/shell_matrix.h\
include/numerics/sparse_matrix.h\
include/parallel/threads.h\
include/utils/o_string_stream.h\
@@ -6570,8 +6571,10 @@
include/mesh/mesh_base.h\
include/numerics/dense_matrix.h\
include/numerics/dense_matrix_base.h\
+ include/numerics/numeric_vector.h\
include/numerics/petsc_macro.h\
include/numerics/petsc_matrix.h\
+ include/numerics/petsc_vector.h\
include/numerics/sparse_matrix.h\
include/numerics/type_vector.h\
include/numerics/vector_value.h\
@@ -6706,6 +6709,57 @@
include/utils/compare_types.h\
include/utils/o_string_stream.h\
include/utils/perf_log.h
+src/numerics/sparse_shell_matrix.$(obj-suffix):\
+ src/numerics/sparse_shell_matrix.C\
+ include/base/libmesh.h\
+ include/base/libmesh_base.h\
+ include/base/libmesh_common.h\
+ include/base/libmesh_config.h\
+ include/base/libmesh_exceptions.h\
+ include/base/libmesh_logging.h\
+ include/base/print_trace.h\
+ include/base/reference_counted_object.h\
+ include/base/reference_counter.h\
+ include/enums/enum_solver_package.h\
+ include/numerics/shell_matrix.h\
+ include/numerics/sparse_shell_matrix.h\
+ include/parallel/threads.h\
+ include/utils/o_string_stream.h\
+ include/utils/perf_log.h
+src/numerics/sum_shell_matrix.$(obj-suffix):\
+ src/numerics/sum_shell_matrix.C\
+ include/base/libmesh.h\
+ include/base/libmesh_base.h\
+ include/base/libmesh_common.h\
+ include/base/libmesh_config.h\
+ include/base/libmesh_exceptions.h\
+ include/base/libmesh_logging.h\
+ include/base/print_trace.h\
+ include/base/reference_counted_object.h\
+ include/base/reference_counter.h\
+ include/enums/enum_solver_package.h\
+ include/numerics/shell_matrix.h\
+ include/numerics/sum_shell_matrix.h\
+ include/parallel/threads.h\
+ include/utils/o_string_stream.h\
+ include/utils/perf_log.h
+src/numerics/tensor_shell_matrix.$(obj-suffix):\
+ src/numerics/tensor_shell_matrix.C\
+ include/base/libmesh.h\
+ include/base/libmesh_base.h\
+ include/base/libmesh_common.h\
+ include/base/libmesh_config.h\
+ include/base/libmesh_exceptions.h\
+ include/base/libmesh_logging.h\
+ include/base/print_trace.h\
+ include/base/reference_counted_object.h\
+ include/base/reference_counter.h\
+ include/enums/enum_solver_package.h\
+ include/numerics/shell_matrix.h\
+ include/numerics/tensor_shell_matrix.h\
+ include/parallel/threads.h\
+ include/utils/o_string_stream.h\
+ include/utils/perf_log.h
src/numerics/trilinos_aztec_linear_solver.$(obj-suffix):\
src/numerics/trilinos_aztec_linear_solver.C\
include/base/auto_ptr.h\
-------------------------------------------------------------------------
This SF.Net email is sponsored by the Moblin Your Move Developer's challenge
Build the coolest Linux based applications with Moblin SDK & win great prizes
Grand prize is a trip for two to an Open Source event anywhere in the world
http://moblin-contest.org/redirect.php?banner_id=100&url=/
_______________________________________________
Libmesh-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-users
