This works fine for standard eigenvalue problems -- a modifed version of ex16 using a ShellMatrix works fine. But I haven't yet got it working for generalized eigenproblems. I've attached a patch (slepc-patch) for the changes I've made to EigenSolver and SlepcEigenSolver, as well as a modified version of ex17 that uses shell matrices (in this case, both A and B are shell matrices -- I've also tried it when only matrix_A is a shell matrix, but I get the same error). The gdb stack trace when the modified ex17 crashes is pasted below.
If anyone has some insights into what is going wrong here (and why the generalized eigenvalue problem behaves differently to the standard eigenvalue problem), that'd be a big help!
Best wishes,
Dave
Program received signal SIGABRT, Aborted.
[Switching to Thread 0x21aa770 (LWP 17661)]
0x00110416 in __kernel_vsyscall ()
(gdb) bt
#0 0x00110416 in __kernel_vsyscall ()
#1 0x040ed660 in raise () from /lib/libc.so.6
#2 0x040ef028 in abort () from /lib/libc.so.6
#3 0x01f11771 in PetscAbortErrorHandler (line=283,
fun=0x485104 "EPSSetOperators", file=0x484f51 "setup.c",
dir=0x484f59 "src/eps/interface/", n=64, p=1,
mess=0xbfd6dd80 "Invalid Pointer to Object: Parameter # 1", ctx=0x0)
at errabort.c:62
#4 0x01f12f0f in PetscError (line=283, func=0x485104 "EPSSetOperators",
file=0x484f51 "setup.c", dir=0x484f59 "src/eps/interface/", n=64, p=1,
mess=0x484f88 "Invalid Pointer to Object: Parameter # %d") at err.c:482
#5 0x003f6352 in EPSSetOperators (eps=0x65702d65, A=0x8a8ecf8, B=0x8a8f868)
at setup.c:283
warning: (Internal error: pc 0x1920dbc in read in psymtab, but not in
symtab.)
#6 0x01920dbd in SlepcEigenSolver<double>::_solve_generalized_helper (this=warning: (Internal error: pc 0x1920dbc in read in psymtab, but not in symtab.)
warning: (Internal error: pc 0x1920c14 in read in psymtab, but not in symtab.)
warning: (Internal error: pc 0x1920dbc in read in psymtab, but not in symtab.)
0x8a411d0, mat_A=warning: (Internal error: pc 0x1920dbc in read in psymtab, but not in symtab.)
0x8a8ecf8, mat_B=warning: (Internal error: pc 0x1920dbc in read in psymtab, but not in symtab.)
0x8a8f868, nev=warning: (Internal error: pc 0x1920dbc in read in psymtab, but not in symtab.)
5, ncv=warning: (Internal error: pc 0x1920dbc in read in psymtab, but not in symtab.)
15, tol=warning: (Internal error: pc 0x1920dbc in read in psymtab, but not in symtab.)
9.9999999999999887e-11, m_its=warning: (Internal error: pc 0x1920dbc in read in psymtab, but not in symtab.)
1000) at src/numerics/slepc_eigen_solver.C:367warning: (Internal error: pc 0x1922558 in read in psymtab, but not in symtab.)
warning: (Internal error: pc 0x1922557 in read in psymtab, but not in symtab.)
warning: (Internal error: pc 0x1922557 in read in psymtab, but not in symtab.)
#7 0x01922558 in SlepcEigenSolver<double>::solve_generalized (this=warning: (Internal error: pc 0x1922557 in read in psymtab, but not in symtab.)
warning: (Internal error: pc 0x1922274 in read in psymtab, but not in symtab.)
warning: (Internal error: pc 0x1922557 in read in psymtab, but not in symtab.)
0x8a411d0,shell_matrix_A=warning: (Internal error: pc 0x1922557 in read in psymtab, but not in symtab.)
@0xbfd6ebc0, shell_matrix_B=warning: (Internal error: pc 0x1922557 in read in psymtab, but not in symtab.)
@0xbfd6ebb8, nev=warning: (Internal error: pc 0x1922557 in read in psymtab, but not in symtab.)
5, ncv=warning: (Internal error: pc 0x1922557 in read in psymtab, but not in symtab.)
15, tol=warning: (Internal error: pc 0x1922557 in read in psymtab, but not in symtab.)
9.9999999999999887e-11, m_its=warning: (Internal error: pc 0x1922557 in read in psymtab, but not in symtab.)
1000) at src/numerics/slepc_eigen_solver.C:336
#8 0x0809499e in ShellEigenSystem::shell_solve (this=0x8a41070)
at shell_eigen_system.C:93
#9 0x0806e422 in main (argc=4, argv=0xbfd6f0d4) at ex17.C:168
Index: include/numerics/eigen_solver.h
===================================================================
--- include/numerics/eigen_solver.h (revision 3351)
+++ include/numerics/eigen_solver.h (working copy)
@@ -38,6 +38,7 @@
// forward declarations
template <typename T> class AutoPtr;
template <typename T> class SparseMatrix;
+template <typename T> class ShellMatrix;
template <typename T> class NumericVector;
@@ -62,7 +63,7 @@
virtual ~EigenSolver ();
/**
- * Builds a \p EigenSolver using the linear solver package specified by
+ * Builds an \p EigenSolver using the linear solver package specified by
* \p solver_package
*/
static AutoPtr<EigenSolver<T> > build(const SolverPackage solver_package =
@@ -120,9 +121,9 @@
{_position_of_spectrum= pos;}
/**
- * Solves the standard eigen problem and returns the
- * number of converged eigenpairs and the number
- * of iterations.
+ * Solves the standard eigen problem when matrix_A is a
+ * \p SparseMatrix, and returns the number of converged
+ * eigenpairs and the number of iterations.
*/
virtual std::pair<unsigned int, unsigned int> solve_standard
(SparseMatrix<T> &matrix_A,
int nev,
@@ -130,10 +131,21 @@
const double
tol,
const unsigned
int m_its) = 0;
+ /**
+ * Solves the standard eigen problem when matrix_A is a
+ * \p ShellMatrix, and returns the number of converged
+ * eigenpairs and the number of iterations.
+ */
+ virtual std::pair<unsigned int, unsigned int> solve_standard (ShellMatrix<T>
&matrix_A,
+ int nev,
+ int ncv,
+ const double
tol,
+ const unsigned
int m_its) = 0;
/**
- * Solves the generalized eigen problem and returns the
+ * Solves the generalized eigen problem when both matrix_A
+ * and matrix_B are of type \p SparseMatrix and returns the
* number of converged eigenpairs and the number
* of iterations.
*/
@@ -144,6 +156,16 @@
const double
tol,
const
unsigned int m_its) = 0;
+ /**
+ * Solves the generalized eigen problem when both matrix_A
+ * and matrix_B are of type ShellMatrix.
+ */
+ virtual std::pair<unsigned int, unsigned int> solve_generalized
(ShellMatrix<T> &matrix_A,
+
ShellMatrix<T> &matrix_B,
+ int nev,
+ int ncv,
+ const double
tol,
+ const
unsigned int m_its) = 0;
/**
Index: include/numerics/slepc_eigen_solver.h
===================================================================
--- include/numerics/slepc_eigen_solver.h (revision 3351)
+++ include/numerics/slepc_eigen_solver.h (working copy)
@@ -96,7 +96,7 @@
/**
* This function calls the SLEPc solver to compute
- * the eigenpairs of matrix matrix_A. \p nev is
+ * the eigenpairs of the SparseMatrix matrix_A. \p nev is
* the number of eigenpairs to be computed and
* \p ncv is the number of basis vectors to be
* used in the solution procedure. Return values
@@ -109,19 +109,30 @@
int ncv,
const double tol,
const unsigned int
m_its);
+
+ /**
+ * Same as above except that matrix_A is a ShellMatrix
+ * in this case.
+ */
+ std::pair<unsigned int, unsigned int> solve_standard (ShellMatrix<T>
&shell_matrix,
+ int nev,
+ int ncv,
+ const double tol,
+ const unsigned int
m_its);
-/**
+ /**
* This function calls the SLEPc solver to compute
- * the eigenpairs of matrix matrix_A provided the matrix B.
- * in case of a generalized eigenproblem \p nev is
- * the number of eigenpairs to be computed and
- * \p ncv is the number of basis vectors to be
+ * the eigenpairs for the generalized eigenproblem
+ * defined by the matrix_A and matrix_B,
+ * which are of type SparseMatrix. The argument
+ * \p nev is the number of eigenpairs to be computed
+ * and \p ncv is the number of basis vectors to be
* used in the solution procedure. Return values
* are the number of converged eigen values and the
* number of the iterations carried out by the eigen
* solver.
- */
+ */
std::pair<unsigned int, unsigned int> solve_generalized(SparseMatrix<T>
&matrix_A,
SparseMatrix<T>
&matrix_B,
int nev,
@@ -129,6 +140,18 @@
const double tol,
const unsigned int
m_its);
+
+ /**
+ * Same as above except that matrix_A and matrix_B are of
+ * type ShellMatrix.
+ */
+ std::pair<unsigned int, unsigned int> solve_generalized(ShellMatrix<T>
&matrix_A,
+ ShellMatrix<T>
&matrix_B,
+ int nev,
+ int ncv,
+ const double tol,
+ const unsigned int
m_its);
+
/**
@@ -150,6 +173,25 @@
private:
/**
+ * Helper function that actually performs the standard eigensolve.
+ */
+ std::pair<unsigned int, unsigned int> _solve_standard_helper (Mat mat,
+ int nev,
+ int ncv,
+ const double tol,
+ const unsigned int
m_its);
+
+ /**
+ * Helper function that actually performs the generalized eigensolve.
+ */
+ std::pair<unsigned int, unsigned int> _solve_generalized_helper (Mat mat_A,
+ Mat mat_B,
+ int nev,
+ int ncv,
+ const double
tol,
+ const unsigned
int m_its);
+
+ /**
* Tells Slepc to use the user-specified solver stored in
* \p _eigen_solver_type
*/
@@ -166,13 +208,20 @@
* stored in \p _position_of_spectrum
*/
void set_slepc_position_of_spectrum();
-
-
+
/**
+ * Internal function if shell matrix mode is used, this just
+ * calls the shell matrix's matrix multiplication function.
+ * See PetscLinearSolver for a similar implementation.
+ */
+ static PetscErrorCode _petsc_shell_matrix_mult(Mat mat, Vec arg, Vec dest);
+
+ /**
* Eigenproblem solver context
*/
EPS _eps;
+
};
Index: src/numerics/slepc_eigen_solver.C
===================================================================
--- src/numerics/slepc_eigen_solver.C (revision 3351)
+++ src/numerics/slepc_eigen_solver.C (working copy)
@@ -29,6 +29,7 @@
// Local Includes
#include "libmesh_logging.h"
#include "slepc_eigen_solver.h"
+#include "shell_matrix.h"
/*----------------------- functions ----------------------------------*/
@@ -50,7 +51,7 @@
#else
// Krylov-Schur showed up as of Slepc 2.3.2
this->_eigen_solver_type = KRYLOVSCHUR;
-#endif
+#endif
}
}
@@ -87,15 +88,76 @@
const double tol, // solver
tolerance
const unsigned int m_its) // maximum
number of iterations
{
- START_LOG("solve_standard()", "SlepcEigenSolver");
+// START_LOG("solve_standard()", "SlepcEigenSolver");
this->init ();
// Make sure the SparseMatrix passed in is really a PetscMatrix
PetscMatrix<T>* matrix_A = libmesh_cast_ptr<PetscMatrix<T>*>(&matrix_A_in);
+ // Close the matrix and vectors in case this wasn't already done.
+ matrix_A->close ();
+
+ // just for debugging, remove this
+// char mat_file[] = "matA.petsc";
+// PetscViewer petsc_viewer;
+// ierr = PetscViewerBinaryOpen(libMesh::COMM_WORLD, mat_file,
PETSC_FILE_CREATE, &petsc_viewer);
+// CHKERRABORT(libMesh::COMM_WORLD,ierr);
+// ierr = MatView(matrix_A->mat(),petsc_viewer);
+// CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ return _solve_standard_helper(matrix_A->mat(), nev, ncv, tol, m_its);
+}
+
+
+template <typename T>
+std::pair<unsigned int, unsigned int>
+SlepcEigenSolver<T>::solve_standard (ShellMatrix<T> &shell_matrix,
+ int nev, // number of
requested eigenpairs
+ int ncv, // number of
basis vectors
+ const double tol, // solver
tolerance
+ const unsigned int m_its) // maximum
number of iterations
+{
+ this->init ();
+
int ierr=0;
+ // Prepare the matrix.
+ Mat mat;
+ ierr = MatCreateShell(libMesh::COMM_WORLD,
+ shell_matrix.m(), // Specify the number of local rows
+ shell_matrix.n(), // Specify the number of local columns
+ PETSC_DETERMINE, // Can't specify global and not local??
+ PETSC_DETERMINE,
+ 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));
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+
+ return _solve_standard_helper(mat, nev, ncv, tol, m_its);
+}
+
+template <typename T>
+std::pair<unsigned int, unsigned int>
+SlepcEigenSolver<T>::_solve_standard_helper
+ (Mat mat,
+ int nev, // number of
requested eigenpairs
+ int ncv, // number of
basis vectors
+ const double tol, // solver
tolerance
+ const unsigned int m_its) // maximum
number of iterations
+{
+ START_LOG("solve_standard()", "SlepcEigenSolver");
+
+ int ierr=0;
+
// converged eigen pairs and number of iterations
int nconv=0;
int its=0;
@@ -108,19 +170,8 @@
PetscScalar kr, ki;
#endif
- // Close the matrix and vectors in case this wasn't already done.
- matrix_A->close ();
-
- // just for debugging, remove this
-// char mat_file[] = "matA.petsc";
-// PetscViewer petsc_viewer;
-// ierr = PetscViewerBinaryOpen(libMesh::COMM_WORLD, mat_file,
PETSC_FILE_CREATE, &petsc_viewer);
-// CHKERRABORT(libMesh::COMM_WORLD,ierr);
-// ierr = MatView(matrix_A->mat(),petsc_viewer);
-// CHKERRABORT(libMesh::COMM_WORLD,ierr);
-
// Set operators.
- ierr = EPSSetOperators (_eps, matrix_A->mat(), PETSC_NULL);
+ ierr = EPSSetOperators (_eps, mat, PETSC_NULL);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
//set the problem type and the position of the spectrum
@@ -202,6 +253,7 @@
CHKERRABORT(libMesh::COMM_WORLD,ierr);
#endif // DEBUG
+
STOP_LOG("solve_standard()", "SlepcEigenSolver");
// return the number of converged eigenpairs
@@ -223,21 +275,87 @@
const double tol, // solver
tolerance
const unsigned int m_its) // maximum
number of iterations
{
- START_LOG("solve_generalized()", "SlepcEigenSolver");
+ this->init ();
- this->init ();
-
// Make sure the data passed in are really of Petsc types
PetscMatrix<T>* matrix_A = libmesh_cast_ptr<PetscMatrix<T>*>(&matrix_A_in);
PetscMatrix<T>* matrix_B = libmesh_cast_ptr<PetscMatrix<T>*>(&matrix_B_in);
+ // Close the matrix and vectors in case this wasn't already done.
+ matrix_A->close ();
+ matrix_B->close ();
+
+
+ return _solve_generalized_helper (matrix_A->mat(), matrix_B->mat(), nev,
ncv, tol, m_its);
+}
+
+
+template <typename T>
+std::pair<unsigned int, unsigned int>
+SlepcEigenSolver<T>::solve_generalized (ShellMatrix<T> &shell_matrix_A,
+ ShellMatrix<T> &shell_matrix_B,
+ int nev, // number of
requested eigenpairs
+ int ncv, // number of
basis vectors
+ const double tol, // solver
tolerance
+ const unsigned int m_its) // maximum
number of iterations
+{
int ierr=0;
+ // Prepare the matrix.
+ Mat mat_A;
+ ierr = MatCreateShell(libMesh::COMM_WORLD,
+ shell_matrix_A.m(), // Specify the number of local rows
+ shell_matrix_A.n(), // Specify the number of local columns
+ PETSC_DETERMINE, // Can't specify global and not local??
+ PETSC_DETERMINE,
+ const_cast<void*>(static_cast<const void*>(&shell_matrix_A)),
+ &mat_A);
+
+ Mat mat_B;
+ ierr = MatCreateShell(libMesh::COMM_WORLD,
+ shell_matrix_B.m(), // Specify the number of local rows
+ shell_matrix_B.n(), // Specify the number of local columns
+ PETSC_DETERMINE, // Can't specify global and not local??
+ PETSC_DETERMINE,
+ const_cast<void*>(static_cast<const void*>(&shell_matrix_B)),
+ &mat_B);
+
+ /* 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_A,MATOP_MULT,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult));
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ ierr =
MatShellSetOperation(mat_B,MATOP_MULT,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult));
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ return _solve_generalized_helper (mat_A, mat_B, nev, ncv, tol, m_its);
+}
+
+
+
+template <typename T>
+std::pair<unsigned int, unsigned int>
+SlepcEigenSolver<T>::_solve_generalized_helper (Mat mat_A,
+ Mat mat_B,
+ int nev, //
number of requested eigenpairs
+ int ncv, //
number of basis vectors
+ const double tol, //
solver tolerance
+ const unsigned int m_its) //
maximum number of iterations
+{
+ START_LOG("solve_generalized()", "SlepcEigenSolver");
+
+ int ierr=0;
+
// converged eigen pairs and number of iterations
int nconv=0;
int its=0;
-#ifdef DEBUG
+#ifdef DEBUG
// The relative error.
PetscReal error, re, im;
@@ -245,20 +363,8 @@
PetscScalar kr, ki;
#endif
- // Close the matrix and vectors in case this wasn't already done.
- matrix_A->close ();
- matrix_B->close ();
-
- // just for debugging, remove this
-// char mat_file[] = "matA.petsc";
-// PetscViewer petsc_viewer;
-// ierr = PetscViewerBinaryOpen(libMesh::COMM_WORLD, mat_file,
PETSC_FILE_CREATE, &petsc_viewer);
-// CHKERRABORT(libMesh::COMM_WORLD,ierr);
-// ierr = MatView(matrix_A->mat(),petsc_viewer);
-// CHKERRABORT(libMesh::COMM_WORLD,ierr);
-
// Set operators.
- ierr = EPSSetOperators (_eps, matrix_A->mat(),matrix_B->mat());
+ ierr = EPSSetOperators (_eps, mat_A, mat_B);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
//set the problem type and the position of the spectrum
@@ -495,9 +601,29 @@
return error;
}
+template <typename T>
+PetscErrorCode SlepcEigenSolver<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;
+}
+
+
//------------------------------------------------------------------
// Explicit instantiations
template class SlepcEigenSolver<Number>;
ex17-shell-matrix.tar.gz
Description: GNU Zip compressed data
------------------------------------------------------------------------------ This SF.net email is sponsored by: High Quality Requirements in a Collaborative Environment. Download a free trial of Rational Requirements Composer Now! http://p.sf.net/sfu/www-ibm-com
_______________________________________________ Libmesh-devel mailing list [email protected] https://lists.sourceforge.net/lists/listinfo/libmesh-devel
