[petsc-users] SLEPc/NEP for shell matrice T(lambda) and T'(lambda)

[Jose E. Roman]

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 <kenneth.c.hall at duke.edu> 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 <jroman at dsic.upv.es>
> Date: Friday, October 6, 2023 at 7:01 AM
> To: Kenneth C Hall <kenneth.c.hall at duke.edu>
> Cc: petsc-users at mcs.anl.gov <petsc-users at mcs.anl.gov>
> 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 <kenneth.c.hall at duke.edu> 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>
>