Jose, Thanks for this. I will try this and report back.
Kenneth From: Jose E. Roman <[email protected]> Date: Friday, October 6, 2023 at 7:01 AM To: Kenneth C Hall <[email protected]> Cc: [email protected] <[email protected]> Subject: Re: [petsc-users] SLEPc/NEP for shell matrice T(lambda) and T'(lambda) I am getting an error in a different place than you. I started to debug, but don't have much time at the moment. Can you try something? Comparing to ex21.c, I see that a difference that may be relevant is the MATOP_DUPLICATE operation. Can you try defining it for your A matrix? Note: If you plan to use the NLEIGS solver, there is no need to define the derivative T' so you can skip the call to NEPSetJacobian(). Jose > El 6 oct 2023, a las 0:37, Kenneth C Hall <[email protected]> escribió: > > 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><test_nep.f90>
