Dear Roy and John and all,
Please find attached my next version.
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.
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).
Thank you for the hint that a function pointer may not be casted into
a void*. I did not know this before. Anyway, this is no longer
necessary since I am now working with a pointer to a class.
Again, everything compiles well now but is not yet tested. Comments
are welcome. I'll keep you informed about my progress.
By the way, why is the .depend file contained in the repository?
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 3088)
+++ include/numerics/petsc_matrix.h (Arbeitskopie)
@@ -292,7 +292,13 @@
*/
Mat mat () { libmesh_assert (_mat != NULL); return _mat; }
-
+ /**
+ * Multiplies the matrix with \p arg and stores the result in \p
+ * dest.
+ */
+ virtual void vector_mult (NumericVector<T>& dest,
+ const NumericVector<T>& arg) const;
+
protected:
/**
Index: include/numerics/laspack_matrix.h
===================================================================
--- include/numerics/laspack_matrix.h (Revision 3088)
+++ include/numerics/laspack_matrix.h (Arbeitskopie)
@@ -272,7 +272,13 @@
*/
void print_personal(std::ostream& os=std::cout) const { this->print(os); }
-
+ /**
+ * Multiplies the matrix with \p arg and stores the result in \p
+ * dest.
+ */
+ virtual void vector_mult (NumericVector<T>& dest,
+ const NumericVector<T>& arg) const;
+
private:
/**
Index: include/numerics/linear_solver.h
===================================================================
--- include/numerics/linear_solver.h (Revision 3088)
+++ 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/petsc_linear_solver.h
===================================================================
--- include/numerics/petsc_linear_solver.h (Revision 3088)
+++ 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/laspack_linear_solver.h
===================================================================
--- include/numerics/laspack_linear_solver.h (Revision 3088)
+++ 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/sparse_matrix.h
===================================================================
--- include/numerics/sparse_matrix.h (Revision 3088)
+++ 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,13 @@
true); // true means REUSE submatrix
}
+ /**
+ * Multiplies the matrix with \p arg and stores the result in \p
+ * dest.
+ */
+ virtual void vector_mult (NumericVector<T>& dest,
+ const NumericVector<T>& arg) const = 0;
+
protected:
/**
Index: include/solvers/linear_implicit_system.h
===================================================================
--- include/solvers/linear_implicit_system.h (Revision 3088)
+++ 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,17 @@
*/
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);
+
protected:
/**
@@ -137,6 +149,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 3088)
+++ 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 3088)
+++ 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/laspack_matrix.C
===================================================================
--- src/numerics/laspack_matrix.C (Revision 3088)
+++ src/numerics/laspack_matrix.C (Arbeitskopie)
@@ -194,6 +194,16 @@
}
+
+template <typename T>
+void LaspackMatrix<T>::vector_mult (NumericVector<T>& dest,
+ const NumericVector<T>& arg) const
+{
+ libmesh_not_implemented();
+}
+
+
+
//------------------------------------------------------------------
// Explicit instantiations
template class LaspackMatrix<Number>;
Index: src/numerics/petsc_matrix.C
===================================================================
--- src/numerics/petsc_matrix.C (Revision 3088)
+++ 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,43 @@
+template <typename T>
+void PetscMatrix<T>::vector_mult (NumericVector<T>& dest,
+ const NumericVector<T>& arg) const
+{
+#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();
+#else
+ /* Check for matching dimensions. */
+ libmesh_assert(dest.size()==this->m());
+ libmesh_assert(arg.size()==this->n());
+
+ /* Cast the vectors to PetscVector and check whether that was
+ successful. */
+ PetscVector<T>* dest_petsc = dynamic_cast<PetscVector<T>*>(&dest);
+ const PetscVector<T>* arg_petsc = dynamic_cast<const PetscVector<T>*>(&arg);
+ libmesh_assert(dest_petsc != NULL);
+ libmesh_assert(arg_petsc != NULL);
+
+ /* Call PETSc function. Since PETSc doesn't eat const arguments,
+ this requires const_cast. */
+ int ierr = 0;
+ ierr = MatMult(const_cast<PetscMatrix<T>*>(this)->mat(),
+ const_cast<PetscVector<T>*>(arg_petsc)->vec(),
+ dest_petsc->vec());
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+#endif
+
+}
+
+
+
//------------------------------------------------------------------
// Explicit instantiations
template class PetscMatrix<Number>;
Index: src/solvers/linear_implicit_system.C
===================================================================
--- src/solvers/linear_implicit_system.C (Revision 3088)
+++ 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 3088)
+++ .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
