Actual source code: test7.c

      1: /*
      2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
      4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

      6:    This file is part of SLEPc.
      7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
      8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      9: */

     11: static char help[] = "Test the NLEIGS solver with shell matrices.\n\n"
     12:   "This is based on ex27.\n"
     13:   "The command line options are:\n"
     14:   "  -n <n>, where <n> = matrix dimension.\n"
     15:   "  -split <0/1>, to select the split form in the problem definition (enabled by default).\n";

     17: /*
     18:    Solve T(lambda)x=0 using NLEIGS solver
     19:       with T(lambda) = -D+sqrt(lambda)*I
     20:       where D is the Laplacian operator in 1 dimension
     21:       and with the interpolation interval [.01,16]
     22: */

     24: #include <slepcnep.h>

     26: /* User-defined routines */
     27: PetscErrorCode FormFunction(NEP,PetscScalar,Mat,Mat,void*);
     28: PetscErrorCode ComputeSingularities(NEP,PetscInt*,PetscScalar*,void*);
     29: PetscErrorCode MatMult_A0(Mat,Vec,Vec);
     30: PetscErrorCode MatGetDiagonal_A0(Mat,Vec);
     31: PetscErrorCode MatDuplicate_A0(Mat,MatDuplicateOption,Mat*);
     32: PetscErrorCode MatMult_A1(Mat,Vec,Vec);
     33: PetscErrorCode MatGetDiagonal_A1(Mat,Vec);
     34: PetscErrorCode MatDuplicate_A1(Mat,MatDuplicateOption,Mat*);
     35: PetscErrorCode MatMult_F(Mat,Vec,Vec);
     36: PetscErrorCode MatGetDiagonal_F(Mat,Vec);
     37: PetscErrorCode MatDuplicate_F(Mat,MatDuplicateOption,Mat*);
     38: PetscErrorCode MatDestroy_F(Mat);

     40: typedef struct {
     41:   PetscScalar t;  /* square root of lambda */
     42: } MatCtx;

     44: int main(int argc,char **argv)
     45: {
     46:   NEP            nep;
     47:   KSP            *ksp;
     48:   PC             pc;
     49:   Mat            F,A[2];
     50:   NEPType        type;
     51:   PetscInt       i,n=100,nev,its,nsolve;
     52:   PetscReal      keep,tol=PETSC_SQRT_MACHINE_EPSILON/10;
     53:   RG             rg;
     54:   FN             f[2];
     55:   PetscBool      terse,flg,lock,split=PETSC_TRUE;
     56:   PetscScalar    coeffs;
     57:   MatCtx         *ctx;

     59:   PetscFunctionBeginUser;
     60:   PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
     61:   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
     62:   PetscCall(PetscOptionsGetBool(NULL,NULL,"-split",&split,NULL));
     63:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nSquare root eigenproblem, n=%" PetscInt_FMT "%s\n\n",n,split?" (in split form)":""));

     65:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     66:      Create NEP context, configure NLEIGS with appropriate options
     67:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     69:   PetscCall(NEPCreate(PETSC_COMM_WORLD,&nep));
     70:   PetscCall(NEPSetType(nep,NEPNLEIGS));
     71:   PetscCall(NEPNLEIGSSetSingularitiesFunction(nep,ComputeSingularities,NULL));
     72:   PetscCall(NEPGetRG(nep,&rg));
     73:   PetscCall(RGSetType(rg,RGINTERVAL));
     74: #if defined(PETSC_USE_COMPLEX)
     75:   PetscCall(RGIntervalSetEndpoints(rg,0.01,16.0,-0.001,0.001));
     76: #else
     77:   PetscCall(RGIntervalSetEndpoints(rg,0.01,16.0,0,0));
     78: #endif
     79:   PetscCall(NEPSetTarget(nep,1.1));
     80:   PetscCall(NEPNLEIGSGetKSPs(nep,&nsolve,&ksp));
     81:   for (i=0;i<nsolve;i++) {
     82:    PetscCall(KSPSetType(ksp[i],KSPBICG));
     83:    PetscCall(KSPGetPC(ksp[i],&pc));
     84:    PetscCall(PCSetType(pc,PCJACOBI));
     85:    PetscCall(KSPSetTolerances(ksp[i],tol,PETSC_CURRENT,PETSC_CURRENT,PETSC_CURRENT));
     86:   }

     88:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     89:      Define the nonlinear problem
     90:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     92:   if (split) {
     93:     /* Create matrix A0 (tridiagonal) */
     94:     PetscCall(MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,NULL,&A[0]));
     95:     PetscCall(MatShellSetOperation(A[0],MATOP_MULT,(PetscErrorCodeFn*)MatMult_A0));
     96:     PetscCall(MatShellSetOperation(A[0],MATOP_MULT_TRANSPOSE,(PetscErrorCodeFn*)MatMult_A0));
     97:     PetscCall(MatShellSetOperation(A[0],MATOP_GET_DIAGONAL,(PetscErrorCodeFn*)MatGetDiagonal_A0));
     98:     PetscCall(MatShellSetOperation(A[0],MATOP_DUPLICATE,(PetscErrorCodeFn*)MatDuplicate_A0));

    100:     /* Create matrix A0 (identity) */
    101:     PetscCall(MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,NULL,&A[1]));
    102:     PetscCall(MatShellSetOperation(A[1],MATOP_MULT,(PetscErrorCodeFn*)MatMult_A1));
    103:     PetscCall(MatShellSetOperation(A[1],MATOP_MULT_TRANSPOSE,(PetscErrorCodeFn*)MatMult_A1));
    104:     PetscCall(MatShellSetOperation(A[1],MATOP_GET_DIAGONAL,(PetscErrorCodeFn*)MatGetDiagonal_A1));
    105:     PetscCall(MatShellSetOperation(A[1],MATOP_DUPLICATE,(PetscErrorCodeFn*)MatDuplicate_A1));

    107:     /* Define functions for the split form */
    108:     PetscCall(FNCreate(PETSC_COMM_WORLD,&f[0]));
    109:     PetscCall(FNSetType(f[0],FNRATIONAL));
    110:     coeffs = 1.0;
    111:     PetscCall(FNRationalSetNumerator(f[0],1,&coeffs));
    112:     PetscCall(FNCreate(PETSC_COMM_WORLD,&f[1]));
    113:     PetscCall(FNSetType(f[1],FNSQRT));
    114:     PetscCall(NEPSetSplitOperator(nep,2,A,f,SUBSET_NONZERO_PATTERN));
    115:   } else {
    116:     /* Callback form: create shell matrix for F=A0+sqrt(lambda)*A1  */
    117:     PetscCall(PetscNew(&ctx));
    118:     PetscCall(MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,(void*)ctx,&F));
    119:     PetscCall(MatShellSetOperation(F,MATOP_MULT,(PetscErrorCodeFn*)MatMult_F));
    120:     PetscCall(MatShellSetOperation(F,MATOP_MULT_TRANSPOSE,(PetscErrorCodeFn*)MatMult_F));
    121:     PetscCall(MatShellSetOperation(F,MATOP_GET_DIAGONAL,(PetscErrorCodeFn*)MatGetDiagonal_F));
    122:     PetscCall(MatShellSetOperation(F,MATOP_DUPLICATE,(PetscErrorCodeFn*)MatDuplicate_F));
    123:     PetscCall(MatShellSetOperation(F,MATOP_DESTROY,(PetscErrorCodeFn*)MatDestroy_F));
    124:     /* Set Function evaluation routine */
    125:     PetscCall(NEPSetFunction(nep,F,F,FormFunction,NULL));
    126:   }

    128:   /* Set solver parameters at runtime */
    129:   PetscCall(NEPSetFromOptions(nep));

    131:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    132:                       Solve the eigensystem
    133:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    134:   PetscCall(NEPSolve(nep));
    135:   PetscCall(NEPGetType(nep,&type));
    136:   PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n",type));
    137:   PetscCall(NEPGetDimensions(nep,&nev,NULL,NULL));
    138:   PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
    139:   PetscCall(PetscObjectTypeCompare((PetscObject)nep,NEPNLEIGS,&flg));
    140:   if (flg) {
    141:     PetscCall(NEPNLEIGSGetRestart(nep,&keep));
    142:     PetscCall(NEPNLEIGSGetLocking(nep,&lock));
    143:     PetscCall(NEPNLEIGSGetInterpolation(nep,&tol,&its));
    144:     PetscCall(PetscPrintf(PETSC_COMM_WORLD," Restart factor is %3.2f",(double)keep));
    145:     if (lock) PetscCall(PetscPrintf(PETSC_COMM_WORLD," (locking activated)"));
    146:     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n Divided diferences with tol=%6.2g maxit=%" PetscInt_FMT "\n",(double)tol,its));
    147:     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n"));
    148:   }

    150:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    151:                     Display solution and clean up
    152:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    154:   /* show detailed info unless -terse option is given by user */
    155:   PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
    156:   if (terse) PetscCall(NEPErrorView(nep,NEP_ERROR_RELATIVE,NULL));
    157:   else {
    158:     PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
    159:     PetscCall(NEPConvergedReasonView(nep,PETSC_VIEWER_STDOUT_WORLD));
    160:     PetscCall(NEPErrorView(nep,NEP_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
    161:     PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
    162:   }
    163:   PetscCall(NEPDestroy(&nep));
    164:   if (split) {
    165:     PetscCall(MatDestroy(&A[0]));
    166:     PetscCall(MatDestroy(&A[1]));
    167:     PetscCall(FNDestroy(&f[0]));
    168:     PetscCall(FNDestroy(&f[1]));
    169:   } else PetscCall(MatDestroy(&F));
    170:   PetscCall(SlepcFinalize());
    171:   return 0;
    172: }

    174: /*
    175:    FormFunction - Computes Function matrix  T(lambda)
    176: */
    177: PetscErrorCode FormFunction(NEP nep,PetscScalar lambda,Mat fun,Mat B,void *ctx)
    178: {
    179:   MatCtx         *ctxF;

    181:   PetscFunctionBeginUser;
    182:   PetscCall(MatShellGetContext(fun,&ctxF));
    183:   ctxF->t = PetscSqrtScalar(lambda);
    184:   PetscFunctionReturn(PETSC_SUCCESS);
    185: }

    187: /*
    188:    ComputeSingularities - Computes maxnp points (at most) in the complex plane where
    189:    the function T(.) is not analytic.

    191:    In this case, we discretize the singularity region (-inf,0)~(-1e+6,-1e-6)
    192: */
    193: PetscErrorCode ComputeSingularities(NEP nep,PetscInt *maxnp,PetscScalar *xi,void *pt)
    194: {
    195:   PetscReal h;
    196:   PetscInt  i;

    198:   PetscFunctionBeginUser;
    199:   h = 12.0/(*maxnp-1);
    200:   xi[0] = -1e-6; xi[*maxnp-1] = -1e+6;
    201:   for (i=1;i<*maxnp-1;i++) xi[i] = -PetscPowReal(10,-6+h*i);
    202:   PetscFunctionReturn(PETSC_SUCCESS);
    203: }

    205: /* -------------------------------- A0 ----------------------------------- */

    207: PetscErrorCode MatMult_A0(Mat A,Vec x,Vec y)
    208: {
    209:   PetscInt          i,n;
    210:   PetscMPIInt       rank,size,next,prev;
    211:   const PetscScalar *px;
    212:   PetscScalar       *py,upper=0.0,lower=0.0;
    213:   MPI_Comm          comm;

    215:   PetscFunctionBeginUser;
    216:   PetscCall(PetscObjectGetComm((PetscObject)A,&comm));
    217:   PetscCallMPI(MPI_Comm_size(comm,&size));
    218:   PetscCallMPI(MPI_Comm_rank(comm,&rank));
    219:   next = rank==size-1? MPI_PROC_NULL: rank+1;
    220:   prev = rank==0? MPI_PROC_NULL: rank-1;

    222:   PetscCall(VecGetArrayRead(x,&px));
    223:   PetscCall(VecGetArray(y,&py));
    224:   PetscCall(VecGetLocalSize(x,&n));

    226:   PetscCallMPI(MPI_Sendrecv(px,1,MPIU_SCALAR,prev,0,&lower,1,MPIU_SCALAR,next,0,comm,MPI_STATUS_IGNORE));
    227:   PetscCallMPI(MPI_Sendrecv(px+n-1,1,MPIU_SCALAR,next,0,&upper,1,MPIU_SCALAR,prev,0,comm,MPI_STATUS_IGNORE));

    229:   py[0] = upper-2.0*px[0]+px[1];
    230:   for (i=1;i<n-1;i++) py[i] = px[i-1]-2.0*px[i]+px[i+1];
    231:   py[n-1] = px[n-2]-2.0*px[n-1]+lower;
    232:   PetscCall(VecRestoreArrayRead(x,&px));
    233:   PetscCall(VecRestoreArray(y,&py));
    234:   PetscFunctionReturn(PETSC_SUCCESS);
    235: }

    237: PetscErrorCode MatGetDiagonal_A0(Mat A,Vec diag)
    238: {
    239:   PetscFunctionBeginUser;
    240:   PetscCall(VecSet(diag,-2.0));
    241:   PetscFunctionReturn(PETSC_SUCCESS);
    242: }

    244: PetscErrorCode MatDuplicate_A0(Mat A,MatDuplicateOption op,Mat *B)
    245: {
    246:   PetscInt       m,n,M,N;
    247:   MPI_Comm       comm;

    249:   PetscFunctionBegin;
    250:   PetscCall(MatGetSize(A,&M,&N));
    251:   PetscCall(MatGetLocalSize(A,&m,&n));
    252:   PetscCall(PetscObjectGetComm((PetscObject)A,&comm));
    253:   PetscCall(MatCreateShell(comm,m,n,M,N,NULL,B));
    254:   PetscCall(MatShellSetOperation(*B,MATOP_MULT,(PetscErrorCodeFn*)MatMult_A0));
    255:   PetscCall(MatShellSetOperation(*B,MATOP_MULT_TRANSPOSE,(PetscErrorCodeFn*)MatMult_A0));
    256:   PetscCall(MatShellSetOperation(*B,MATOP_GET_DIAGONAL,(PetscErrorCodeFn*)MatGetDiagonal_A0));
    257:   PetscCall(MatShellSetOperation(*B,MATOP_DUPLICATE,(PetscErrorCodeFn*)MatDuplicate_A0));
    258:   PetscFunctionReturn(PETSC_SUCCESS);
    259: }

    261: /* -------------------------------- A1 ----------------------------------- */

    263: PetscErrorCode MatMult_A1(Mat A,Vec x,Vec y)
    264: {
    265:   PetscFunctionBeginUser;
    266:   PetscCall(VecCopy(x,y));
    267:   PetscFunctionReturn(PETSC_SUCCESS);
    268: }

    270: PetscErrorCode MatGetDiagonal_A1(Mat A,Vec diag)
    271: {
    272:   PetscFunctionBeginUser;
    273:   PetscCall(VecSet(diag,1.0));
    274:   PetscFunctionReturn(PETSC_SUCCESS);
    275: }

    277: PetscErrorCode MatDuplicate_A1(Mat A,MatDuplicateOption op,Mat *B)
    278: {
    279:   PetscInt       m,n,M,N;
    280:   MPI_Comm       comm;

    282:   PetscFunctionBegin;
    283:   PetscCall(MatGetSize(A,&M,&N));
    284:   PetscCall(MatGetLocalSize(A,&m,&n));
    285:   PetscCall(PetscObjectGetComm((PetscObject)A,&comm));
    286:   PetscCall(MatCreateShell(comm,m,n,M,N,NULL,B));
    287:   PetscCall(MatShellSetOperation(*B,MATOP_MULT,(PetscErrorCodeFn*)MatMult_A1));
    288:   PetscCall(MatShellSetOperation(*B,MATOP_MULT_TRANSPOSE,(PetscErrorCodeFn*)MatMult_A1));
    289:   PetscCall(MatShellSetOperation(*B,MATOP_GET_DIAGONAL,(PetscErrorCodeFn*)MatGetDiagonal_A1));
    290:   PetscCall(MatShellSetOperation(*B,MATOP_DUPLICATE,(PetscErrorCodeFn*)MatDuplicate_A1));
    291:   PetscFunctionReturn(PETSC_SUCCESS);
    292: }

    294: /* -------------------------------- F ----------------------------------- */

    296: PetscErrorCode MatMult_F(Mat A,Vec x,Vec y)
    297: {
    298:   PetscInt          i,n;
    299:   PetscMPIInt       rank,size,next,prev;
    300:   const PetscScalar *px;
    301:   PetscScalar       *py,d,upper=0.0,lower=0.0;
    302:   MatCtx            *ctx;
    303:   MPI_Comm          comm;

    305:   PetscFunctionBeginUser;
    306:   PetscCall(PetscObjectGetComm((PetscObject)A,&comm));
    307:   PetscCallMPI(MPI_Comm_size(comm,&size));
    308:   PetscCallMPI(MPI_Comm_rank(comm,&rank));
    309:   next = rank==size-1? MPI_PROC_NULL: rank+1;
    310:   prev = rank==0? MPI_PROC_NULL: rank-1;

    312:   PetscCall(MatShellGetContext(A,&ctx));
    313:   PetscCall(VecGetArrayRead(x,&px));
    314:   PetscCall(VecGetArray(y,&py));
    315:   PetscCall(VecGetLocalSize(x,&n));

    317:   PetscCallMPI(MPI_Sendrecv(px,1,MPIU_SCALAR,prev,0,&lower,1,MPIU_SCALAR,next,0,comm,MPI_STATUS_IGNORE));
    318:   PetscCallMPI(MPI_Sendrecv(px+n-1,1,MPIU_SCALAR,next,0,&upper,1,MPIU_SCALAR,prev,0,comm,MPI_STATUS_IGNORE));

    320:   d = -2.0+ctx->t;
    321:   py[0] = upper+d*px[0]+px[1];
    322:   for (i=1;i<n-1;i++) py[i] = px[i-1]+d*px[i]+px[i+1];
    323:   py[n-1] = px[n-2]+d*px[n-1]+lower;
    324:   PetscCall(VecRestoreArrayRead(x,&px));
    325:   PetscCall(VecRestoreArray(y,&py));
    326:   PetscFunctionReturn(PETSC_SUCCESS);
    327: }

    329: PetscErrorCode MatGetDiagonal_F(Mat A,Vec diag)
    330: {
    331:   MatCtx         *ctx;

    333:   PetscFunctionBeginUser;
    334:   PetscCall(MatShellGetContext(A,&ctx));
    335:   PetscCall(VecSet(diag,-2.0+ctx->t));
    336:   PetscFunctionReturn(PETSC_SUCCESS);
    337: }

    339: PetscErrorCode MatDuplicate_F(Mat A,MatDuplicateOption op,Mat *B)
    340: {
    341:   MatCtx         *actx,*bctx;
    342:   PetscInt       m,n,M,N;
    343:   MPI_Comm       comm;

    345:   PetscFunctionBegin;
    346:   PetscCall(MatShellGetContext(A,&actx));
    347:   PetscCall(MatGetSize(A,&M,&N));
    348:   PetscCall(MatGetLocalSize(A,&m,&n));
    349:   PetscCall(PetscNew(&bctx));
    350:   bctx->t = actx->t;
    351:   PetscCall(PetscObjectGetComm((PetscObject)A,&comm));
    352:   PetscCall(MatCreateShell(comm,m,n,M,N,(void*)bctx,B));
    353:   PetscCall(MatShellSetOperation(*B,MATOP_MULT,(PetscErrorCodeFn*)MatMult_F));
    354:   PetscCall(MatShellSetOperation(*B,MATOP_MULT_TRANSPOSE,(PetscErrorCodeFn*)MatMult_F));
    355:   PetscCall(MatShellSetOperation(*B,MATOP_GET_DIAGONAL,(PetscErrorCodeFn*)MatGetDiagonal_F));
    356:   PetscCall(MatShellSetOperation(*B,MATOP_DUPLICATE,(PetscErrorCodeFn*)MatDuplicate_F));
    357:   PetscCall(MatShellSetOperation(*B,MATOP_DESTROY,(PetscErrorCodeFn*)MatDestroy_F));
    358:   PetscFunctionReturn(PETSC_SUCCESS);
    359: }

    361: PetscErrorCode MatDestroy_F(Mat A)
    362: {
    363:   MatCtx         *ctx;

    365:   PetscFunctionBegin;
    366:   PetscCall(MatShellGetContext(A,&ctx));
    367:   PetscCall(PetscFree(ctx));
    368:   PetscFunctionReturn(PETSC_SUCCESS);
    369: }

    371: /*TEST

    373:    testset:
    374:       nsize: {{1 2}}
    375:       args: -nep_nev 3 -nep_tol 1e-8 -terse
    376:       filter: sed -e "s/[+-]0\.0*i//g"
    377:       requires: !single
    378:       test:
    379:          suffix: 1
    380:          args: -nep_nleigs_locking 0 -nep_nleigs_interpolation_degree 90 -nep_nleigs_interpolation_tol 1e-8 -nep_nleigs_restart 0.4
    381:       test:
    382:          suffix: 2
    383:          args: -split 0

    385: TEST*/