Dear all,
thank you very much for all your suggestions.
I started implementing the thing now. Find attached a patch that
contains my first attempt. It compiles well at the moment, but it
most probably will not work as expected, because I did not have the
time to test it yet, and I am not too much experienced in PETSc
programming.
However, I would appreciate your feedback concerning my API
suggestion. Let me explain this to you:
Coarse scale: LinearImplicitSystem gets a method attach_shell_matrix()
that the user should call if he wants to do "matrix free" linear
solving. What then happens is that LinearImplicitSystem::solve()
notices that there is a shell matrix and then calls a different solver
function.
Fine scale: LinearSolver gets another solve() method that eats a shell
matrix function rather than a SparseMatrix object. This calls the
respective PETSc functions to do a matrix free solve.
It is a little bit tricky not to get confused with the contexts that
these shell functions need: The user supplied shell function has a
reference to EquationSystems; that should be enough as a context. The
shell function that LinearImplicitSystem passes to LinearSolver cannot
do this because (a) LinearSolver does not know about EquationSystems,
and (b) it wouldn't know which system is the right one. Hence, I use
a pointer to *this here, casted to void*. The shell function that is
passed to PETSc needs both the user shell function and user context
(which is the pointer to the LinearImplicitSystem of course), hence I
construct a std::vector<void*>* which I again cast into void*. I
agree that this is somehow confusing, but I think it should work.
Any comments are welcome. I will test the thing in the next days (by
implementing the application that I need it for) and let you know when
it works.
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/linear_solver.h
===================================================================
--- include/numerics/linear_solver.h (Revision 3087)
+++ include/numerics/linear_solver.h (Arbeitskopie)
@@ -140,7 +140,23 @@
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 (void
shell_matrix(NumericVector<T>& dest,
+ const
NumericVector<T>& arg,
+ void* ctx),
+ void* ctx,
+ 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 3087)
+++ include/numerics/petsc_linear_solver.h (Arbeitskopie)
@@ -122,6 +122,19 @@
const unsigned int m_its);
/**
+ * This function solves a system whose matrix is a shell matrix.
+ */
+ std::pair<unsigned int, Real>
+ solve (void shell_matrix(NumericVector<T>& dest,
+ const NumericVector<T>& arg,
+ void* ctx),
+ void* ctx,
+ 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 +182,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 3087)
+++ include/numerics/laspack_linear_solver.h (Arbeitskopie)
@@ -95,6 +95,19 @@
const unsigned int m_its);
/**
+ * This function solves a system whose matrix is a shell matrix.
+ */
+ std::pair<unsigned int, Real>
+ solve (void shell_matrix(NumericVector<T>& dest,
+ const NumericVector<T>& arg,
+ void* ctx),
+ void* ctx,
+ 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/solvers/linear_implicit_system.h
===================================================================
--- include/solvers/linear_implicit_system.h (Revision 3087)
+++ include/solvers/linear_implicit_system.h (Arbeitskopie)
@@ -125,6 +125,19 @@
*/
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 (void fptr(NumericVector<Number>& dest,
+ const NumericVector<Number>& arg,
+ EquationSystems& es));
+
protected:
/**
@@ -137,6 +150,22 @@
* The final residual for the linear system Ax=b.
*/
Real _final_linear_residual;
+
+ /**
+ * User supplies shell matrix function or \p NULL if no such
+ * function is used.
+ */
+ void(*_shell_matrix)(NumericVector<Number>& dest,
+ const NumericVector<Number>& arg,
+ EquationSystems& es);
+
+private:
+ /**
+ * Internal function for being passed to \p LinearSolver.
+ */
+ static void _internal_shell_matrix(NumericVector<Number>& dest,
+ const NumericVector<Number>& arg,
+ void* ctx);
};
Index: src/numerics/petsc_linear_solver.C
===================================================================
--- src/numerics/petsc_linear_solver.C (Revision 3087)
+++ src/numerics/petsc_linear_solver.C (Arbeitskopie)
@@ -413,6 +413,92 @@
template <typename T>
+std::pair<unsigned int, Real>
+PetscLinearSolver<T>::solve (void shell_matrix(NumericVector<T>& dest,
+ const NumericVector<T>& arg,
+ void* ctx),
+ void* ctx,
+ NumericVector<T>& solution_in,
+ NumericVector<T>& rhs_in,
+ const double tol,
+ const unsigned int m_its)
+{
+ 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 ();
+
+ /* Make a context vector for PETSc. It consits of void pointers to
+ (0) the user context, (1) the shell matrix. */
+ std::vector<void*> inner_ctx;
+ inner_ctx.push_back(ctx);
+ inner_ctx.push_back(reinterpret_cast<void*>(shell_matrix));
+
+ // 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(),
+ &inner_ctx,
+ &mat);
+ 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);
+}
+
+
+
+template <typename T>
void PetscLinearSolver<T>::get_residual_history(std::vector<double>& hist)
{
int ierr = 0;
@@ -648,6 +734,36 @@
}
+template <typename T>
+PetscErrorCode PetscLinearSolver<T>::_petsc_shell_matrix(Mat mat, Vec arg, Vec
dest)
+{
+ /* Get the matrix context. */
+ int ierr=0;
+ void* inner_ctx;
+ ierr = MatShellGetContext(mat,&inner_ctx);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ /* Get user context and used shell function. */
+ void* ctx = (*reinterpret_cast<std::vector<void*>*>(inner_ctx))[0];
+ void(*shell_matrix)(NumericVector<T>& dest,
+ const NumericVector<T>& arg,
+ void* ctx)=
+ reinterpret_cast<void (*)(NumericVector<double>&, const
NumericVector<double>&, void*)>
+ (
+ (*reinterpret_cast<std::vector<void*>*>(inner_ctx))[1]
+ );
+
+ /* Make \p NumericVector instances arount the vectors. */
+ PetscVector<T> arg_global(arg);
+ PetscVector<T> dest_global(dest);
+
+ /* Call the user function. */
+ shell_matrix(dest_global,arg_global,ctx);
+
+ return ierr;
+}
+
+
//------------------------------------------------------------------
// Explicit instantiations
template class PetscLinearSolver<Number>;
Index: src/numerics/laspack_linear_solver.C
===================================================================
--- src/numerics/laspack_linear_solver.C (Revision 3087)
+++ src/numerics/laspack_linear_solver.C (Arbeitskopie)
@@ -284,6 +284,23 @@
template <typename T>
+std::pair<unsigned int, Real>
+LaspackLinearSolver<T>::solve (void shell_matrix(NumericVector<T>& dest,
+ const NumericVector<T>& arg,
+ void* ctx),
+ void* ctx,
+ 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/solvers/linear_implicit_system.C
===================================================================
--- src/solvers/linear_implicit_system.C (Revision 3087)
+++ 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(&_internal_shell_matrix, this, *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,30 @@
// Update the system after the solve
this->update();
}
+
+
+
+void LinearImplicitSystem::attach_shell_matrix (void
fptr(NumericVector<Number>& dest,
+ const
NumericVector<Number>& arg,
+ EquationSystems& es))
+{
+ _shell_matrix = fptr;
+}
+
+
+
+void LinearImplicitSystem::_internal_shell_matrix(NumericVector<Number>& dest,
+ const NumericVector<Number>&
arg,
+ void* ctx)
+{
+ /* Get the "*this" instance. */
+ LinearImplicitSystem& s = *reinterpret_cast<LinearImplicitSystem*>(ctx);
+
+ /* Make sure there is really a shell function in that system. */
+ libmesh_assert(s._shell_matrix!=NULL);
+
+ /* Call it. */
+ s._shell_matrix(dest, arg, s.get_equation_systems());
+}
+
+
-------------------------------------------------------------------------
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
