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

[Kenneth C Hall]

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>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20231006/03f6c029/attachment.html>