Kenneth, The MatDuplicate issue should be fixed in the following MR https://gitlab.com/petsc/petsc/-/merge_requests/6912
Note that the NLEIGS solver internally uses MatDuplicate for creating multiple copies of the shell matrix, each one with its own value of lambda. Hence your implementation of the shell matrix is not appropriate, since you have a single global lambda within the module. I have attempted to write a Fortran example that duplicates the lambda correctly (see the MR), but does not work yet. Jose > El 6 oct 2023, a las 22:28, Kenneth C Hall <[email protected]> escribió: > > Jose, > > Unfortunately, I was unable to implement the MATOP_DUPLICATE operation in > fortran (and I do not know enough c to work in c). Here is the error message > I get: > > [0]PETSC ERROR: #1 MatShellSetOperation_Fortran() at > /Users/hall/Documents/Fortran_Codes/Packages/petsc/src/mat/impls/shell/ftn-custom/zshellf.c:283 > [0]PETSC ERROR: #2 src/test_nep.f90:62 > > When I look at zshellf.c, MATOP_DUPLICATE is not one of the supported > operations. See below. > > Kenneth > > > /** > > > > * Subset of MatOperation that is supported by the Fortran wrappers. > > > > */ > enum FortranMatOperation { > FORTRAN_MATOP_MULT = 0, > FORTRAN_MATOP_MULT_ADD = 1, > FORTRAN_MATOP_MULT_TRANSPOSE = 2, > FORTRAN_MATOP_MULT_TRANSPOSE_ADD = 3, > FORTRAN_MATOP_SOR = 4, > FORTRAN_MATOP_TRANSPOSE = 5, > FORTRAN_MATOP_GET_DIAGONAL = 6, > FORTRAN_MATOP_DIAGONAL_SCALE = 7, > FORTRAN_MATOP_ZERO_ENTRIES = 8, > FORTRAN_MATOP_AXPY = 9, > FORTRAN_MATOP_SHIFT = 10, > FORTRAN_MATOP_DIAGONAL_SET = 11, > FORTRAN_MATOP_DESTROY = 12, > FORTRAN_MATOP_VIEW = 13, > FORTRAN_MATOP_CREATE_VECS = 14, > FORTRAN_MATOP_GET_DIAGONAL_BLOCK = 15, > FORTRAN_MATOP_COPY = 16, > FORTRAN_MATOP_SCALE = 17, > FORTRAN_MATOP_SET_RANDOM = 18, > FORTRAN_MATOP_ASSEMBLY_BEGIN = 19, > FORTRAN_MATOP_ASSEMBLY_END = 20, > FORTRAN_MATOP_SIZE = 21 > }; > > > 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> >
