Hi all, I have a very large eigenvalue problem of the form T(\lambda).x = 0. The eigenvalues appear in a complicated way, and I must use a matrix-free approach to compute the products T.x and T’.x.
I am trying to implement in SLEPc/NEP. To get started, I have defined a much
smaller and simpler system of the form
A.x - \lambda x = 0 where A is a 10x10 matrix. This is of course a simple
standard eigenvalue problem, but I am using it as a surrogate to understand how
to use NEP.
I have set the problem up using shell matrices (as that is my ultimate goal).
The full code is attached, but here is a smaller snippet of code:
!.... Create matrix-free operators for A and B
PetscCall(MatCreateShell(PETSC_COMM_SELF,n,n,PETSC_DETERMINE,PETSC_DETERMINE,
PETSC_NULL_INTEGER, A, ierr))
PetscCall(MatCreateShell(PETSC_COMM_SELF,n,n,PETSC_DETERMINE,PETSC_DETERMINE,
PETSC_NULL_INTEGER, B, ierr))
PetscCall(MatShellSetOperation(A, MATOP_MULT, MatMult_A, ierr))
PetscCall(MatShellSetOperation(B, MATOP_MULT, MatMult_B, ierr))
!.... Create nonlinear eigensolver
PetscCall(NEPCreate(PETSC_COMM_SELF, nep, ierr))
!.... Set the problem type
PetscCall(NEPSetProblemType(nep, NEP_GENERAL, ierr))
!
!.... set the solver type
PetscCall(NEPSetType(nep, NEPNLEIGS, ierr))
!
!.... Set functions and Jacobians for NEP
PetscCall(NEPSetFunction(nep, A, A, MyNEPFunction, PETSC_NULL_INTEGER,
ierr))
PetscCall(NEPSetJacobian(nep, B, MyNEPJacobian, PETSC_NULL_INTEGER,
ierr))
The code runs, calls MyNEPFunction and MatMult_A multiple times, sweeping over
the prescribed RG range, but crashes before it ever calls MyNEPJacobian or
MatMult_B. The NEP viewer and error messages are attached.
Any help on getting this problem properly set up would be greatly appreciated.
Kenneth Hall
ATTACHMENTS:
test_nep.f90
code_output
code_output
Description: code_output
test_nep.f90
Description: test_nep.f90
