Dear John/Roy/all,
As Roy predicted, preconditioning is tricky. I extended the
ShellMatrix class (and all the derived classes) by a get_diagonal()
method, thus allowing at least JACOBI preconditioning, which for the
moment seems to be sufficient for my example. This required to add an
analogous function to the SparseMatrix class and a pointwise_mult()
function to NumericVector.
A patch with the current version is attached.
It works for my example on one processor, where "works" means "does
not crash and gives results that might or might not be correct".
I am currently unable to test it in parallel since there is some
trouble with the cluster I've been using. I hope this will be
resolved soon.
I'll keep you informed.
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/laspack_vector.h
===================================================================
--- include/numerics/laspack_vector.h (Revision 3096)
+++ include/numerics/laspack_vector.h (Arbeitskopie)
@@ -371,6 +371,13 @@
void localize_to_one (std::vector<T>& v_local,
const unsigned int proc_id=0) const;
+ /**
+ * Computes the pointwise (i.e. component-wise) product of \p vec1
+ * and \p vec2 and stores the result in \p *this.
+ */
+ virtual void pointwise_mult (const NumericVector<T>& vec1,
+ const NumericVector<T>& vec2);
+
private:
/**
Index: include/numerics/petsc_vector.h
===================================================================
--- include/numerics/petsc_vector.h (Revision 3096)
+++ include/numerics/petsc_vector.h (Arbeitskopie)
@@ -391,6 +391,13 @@
const unsigned int proc_id=0) const;
/**
+ * Computes the pointwise (i.e. component-wise) product of \p vec1
+ * and \p vec2 and stores the result in \p *this.
+ */
+ virtual void pointwise_mult (const NumericVector<T>& vec1,
+ const NumericVector<T>& vec2);
+
+ /**
* Print the contents of the vector in Matlab
* format. Optionally prints the
* matrix to the file named \p name. If \p name
Index: include/numerics/petsc_matrix.h
===================================================================
--- include/numerics/petsc_matrix.h (Revision 3096)
+++ include/numerics/petsc_matrix.h (Arbeitskopie)
@@ -281,6 +281,11 @@
void print_matlab(const std::string name="NULL") const;
/**
+ * Copies the diagonal part of the matrix into \p dest.
+ */
+ virtual void get_diagonal (NumericVector<T>& dest) const;
+
+ /**
* Swaps the raw PETSc matrix context pointers.
*/
void swap (PetscMatrix<T> &);
@@ -292,7 +297,6 @@
*/
Mat mat () { libmesh_assert (_mat != NULL); return _mat; }
-
protected:
/**
Index: include/numerics/laspack_matrix.h
===================================================================
--- include/numerics/laspack_matrix.h (Revision 3096)
+++ include/numerics/laspack_matrix.h (Arbeitskopie)
@@ -272,7 +272,11 @@
*/
void print_personal(std::ostream& os=std::cout) const { this->print(os); }
-
+ /**
+ * Copies the diagonal part of the matrix into \p dest.
+ */
+ virtual void get_diagonal (NumericVector<T>& dest) const;
+
private:
/**
Index: include/numerics/linear_solver.h
===================================================================
--- include/numerics/linear_solver.h (Revision 3096)
+++ 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/distributed_vector.h
===================================================================
--- include/numerics/distributed_vector.h (Revision 3096)
+++ include/numerics/distributed_vector.h (Arbeitskopie)
@@ -365,6 +365,13 @@
void localize_to_one (std::vector<T>& v_local,
const unsigned int proc_id=0) const;
+ /**
+ * Computes the pointwise (i.e. component-wise) product of \p vec1
+ * and \p vec2 and stores the result in \p *this.
+ */
+ virtual void pointwise_mult (const NumericVector<T>& vec1,
+ const NumericVector<T>& vec2);
+
private:
/**
Index: include/numerics/petsc_linear_solver.h
===================================================================
--- include/numerics/petsc_linear_solver.h (Revision 3096)
+++ 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,16 @@
*/
void set_petsc_preconditioner_type ();
+ /**
+ * Internal function if shell matrix mode is used.
+ */
+ static PetscErrorCode _petsc_shell_matrix_mult(Mat mat, Vec arg, Vec dest);
+
+ /**
+ * Internal function if shell matrix mode is used.
+ */
+ static PetscErrorCode _petsc_shell_matrix_get_diagonal(Mat mat, 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 3096)
+++ 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 3096)
+++ 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;
/**
@@ -303,7 +304,7 @@
* sparse matrix format. Optionally prints the
* matrix to the file named \p name. If \p name
* is not specified it is dumped to the screen.
- */
+x */
virtual void print_matlab(const std::string name="NULL") const
{
std::cerr << "ERROR: Not Implemented in base class yet!" << std::endl;
@@ -342,6 +343,24 @@
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;
+
+ /**
+ * Copies the diagonal part of the matrix into \p dest.
+ */
+ virtual void get_diagonal (NumericVector<T>& dest) const = 0;
+
protected:
/**
Index: include/numerics/numeric_vector.h
===================================================================
--- include/numerics/numeric_vector.h (Revision 3096)
+++ 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
@@ -402,6 +410,13 @@
const Real threshold = TOLERANCE) const;
/**
+ * Computes the pointwise (i.e. component-wise) product of \p vec1
+ * and \p vec2 and stores the result in \p *this.
+ */
+ virtual void pointwise_mult (const NumericVector<T>& vec1,
+ const NumericVector<T>& vec2) = 0;
+
+ /**
* Prints the local contents of the vector to the screen.
*/
virtual void print(std::ostream& os=std::cout) const;
Index: include/solvers/linear_implicit_system.h
===================================================================
--- include/solvers/linear_implicit_system.h (Revision 3096)
+++ 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 3096)
+++ 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,101 @@
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_mult));
+ ierr =
MatShellSetOperation(mat,MATOP_GET_DIAGONAL,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_get_diagonal));
+ 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 +744,54 @@
}
+
+template <typename T>
+PetscErrorCode PetscLinearSolver<T>::_petsc_shell_matrix_mult(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 around 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;
+}
+
+
+
+template <typename T>
+PetscErrorCode PetscLinearSolver<T>::_petsc_shell_matrix_get_diagonal(Mat mat,
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 around the vector. */
+ PetscVector<T> dest_global(dest);
+
+ /* Call the user function. */
+ shell_matrix.get_diagonal(dest_global);
+
+ return ierr;
+}
+
+
+
//------------------------------------------------------------------
// Explicit instantiations
template class PetscLinearSolver<Number>;
Index: src/numerics/laspack_linear_solver.C
===================================================================
--- src/numerics/laspack_linear_solver.C (Revision 3096)
+++ 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 3096)
+++ 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 3096)
+++ 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_vector.C
===================================================================
--- src/numerics/petsc_vector.C (Revision 3096)
+++ src/numerics/petsc_vector.C (Arbeitskopie)
@@ -879,6 +879,38 @@
template <typename T>
+void PetscVector<T>::pointwise_mult (const NumericVector<T>& vec1,
+ const NumericVector<T>& vec2)
+{
+ int ierr = 0;
+
+ // Convert arguments to PetscVector*.
+ const PetscVector<T>* vec1_petsc = dynamic_cast<const
PetscVector<T>*>(&vec1);
+ const PetscVector<T>* vec2_petsc = dynamic_cast<const
PetscVector<T>*>(&vec2);
+ libmesh_assert (vec1_petsc != NULL);
+ libmesh_assert (vec2_petsc != NULL);
+
+ // Call PETSc function.
+
+#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
+
+ ierr = VecPointwiseMult(this->vec(),
+ const_cast<PetscVector<T>*>(vec1_petsc)->vec(),
+ const_cast<PetscVector<T>*>(vec2_petsc)->vec());
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+#endif
+
+}
+
+
+
+template <typename T>
void PetscVector<T>::print_matlab (const std::string name) const
{
libmesh_assert (this->initialized());
Index: src/numerics/laspack_vector.C
===================================================================
--- src/numerics/laspack_vector.C (Revision 3096)
+++ src/numerics/laspack_vector.C (Arbeitskopie)
@@ -387,6 +387,15 @@
+template <typename T>
+void LaspackVector<T>::pointwise_mult (const NumericVector<T>& vec1,
+ const NumericVector<T>& vec2)
+{
+ libmesh_not_implemented();
+}
+
+
+
template <typename T>
Real LaspackVector<T>::max() const
{
Index: src/numerics/laspack_matrix.C
===================================================================
--- src/numerics/laspack_matrix.C (Revision 3096)
+++ src/numerics/laspack_matrix.C (Arbeitskopie)
@@ -194,6 +194,15 @@
}
+
+template <typename T>
+void LaspackMatrix<T>::get_diagonal (NumericVector<T>& dest) const
+{
+ libmesh_not_implemented();
+}
+
+
+
//------------------------------------------------------------------
// Explicit instantiations
template class LaspackMatrix<Number>;
Index: src/numerics/petsc_matrix.C
===================================================================
--- src/numerics/petsc_matrix.C (Revision 3096)
+++ 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,33 @@
+template <typename T>
+void PetscMatrix<T>::get_diagonal (NumericVector<T>& dest) const
+{
+ int ierr=0;
+ // Convert vector to PetscVector.
+ PetscVector<T>* petsc_dest = dynamic_cast<PetscVector<T>*>(&dest);
+ libmesh_assert(petsc_vector != NULL);
+ // Call PETSc function.
+#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
+
+ // Needs a const_cast since PETSc does not work with const.
+ ierr =
MatGetDiagonal(const_cast<PetscMatrix<T>*>(this)->mat(),petsc_dest->vec());
CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+#endif
+
+}
+
+
+
//------------------------------------------------------------------
// Explicit instantiations
template class PetscMatrix<Number>;
Index: src/numerics/distributed_vector.C
===================================================================
--- src/numerics/distributed_vector.C (Revision 3096)
+++ src/numerics/distributed_vector.C (Arbeitskopie)
@@ -534,6 +534,16 @@
}
+
+template <typename T>
+void DistributedVector<T>::pointwise_mult (const NumericVector<T>& vec1,
+ const NumericVector<T>& vec2)
+{
+ libmesh_not_implemented();
+}
+
+
+
//--------------------------------------------------------------
// Explicit instantiations
template class DistributedVector<Number>;
Index: src/solvers/linear_implicit_system.C
===================================================================
--- src/solvers/linear_implicit_system.C (Revision 3096)
+++ 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);
+ }
+ else 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 3096)
+++ .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
