Actual source code: test4.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 RII solver with a user-provided KSP.\n\n"
     12:   "This is a simplified version of ex20.\n"
     13:   "The command line options are:\n"
     14:   "  -n <n>, where <n> = number of grid subdivisions.\n";

     16: /*
     17:    Solve 1-D PDE
     18:             -u'' = lambda*u
     19:    on [0,1] subject to
     20:             u(0)=0, u'(1)=u(1)*lambda*kappa/(kappa-lambda)
     21: */

     23: #include <slepcnep.h>

     25: /*
     26:    User-defined routines
     27: */
     28: PetscErrorCode FormFunction(NEP,PetscScalar,Mat,Mat,void*);
     29: PetscErrorCode FormJacobian(NEP,PetscScalar,Mat,void*);

     31: /*
     32:    User-defined application context
     33: */
     34: typedef struct {
     35:   PetscScalar kappa;   /* ratio between stiffness of spring and attached mass */
     36:   PetscReal   h;       /* mesh spacing */
     37: } ApplicationCtx;

     39: int main(int argc,char **argv)
     40: {
     41:   NEP            nep;
     42:   KSP            ksp;
     43:   PC             pc;
     44:   Mat            F,J;
     45:   ApplicationCtx ctx;
     46:   PetscInt       n=128,lag,its;
     47:   PetscBool      terse,flg,cct,herm;
     48:   PetscReal      thres;

     50:   PetscFunctionBeginUser;
     51:   PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
     52:   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
     53:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n1-D Nonlinear Eigenproblem, n=%" PetscInt_FMT "\n\n",n));
     54:   ctx.h = 1.0/(PetscReal)n;
     55:   ctx.kappa = 1.0;

     57:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     58:         Create a standalone KSP with appropriate settings
     59:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     61:   PetscCall(KSPCreate(PETSC_COMM_WORLD,&ksp));
     62:   PetscCall(KSPSetType(ksp,KSPBCGS));
     63:   PetscCall(KSPGetPC(ksp,&pc));
     64:   PetscCall(PCSetType(pc,PCBJACOBI));
     65:   PetscCall(KSPSetFromOptions(ksp));

     67:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     68:                Prepare nonlinear eigensolver context
     69:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     71:   PetscCall(NEPCreate(PETSC_COMM_WORLD,&nep));

     73:   /* Create Function and Jacobian matrices; set evaluation routines */
     74:   PetscCall(MatCreate(PETSC_COMM_WORLD,&F));
     75:   PetscCall(MatSetSizes(F,PETSC_DECIDE,PETSC_DECIDE,n,n));
     76:   PetscCall(MatSetFromOptions(F));
     77:   PetscCall(MatSeqAIJSetPreallocation(F,3,NULL));
     78:   PetscCall(MatMPIAIJSetPreallocation(F,3,NULL,1,NULL));
     79:   PetscCall(NEPSetFunction(nep,F,F,FormFunction,&ctx));

     81:   PetscCall(MatCreate(PETSC_COMM_WORLD,&J));
     82:   PetscCall(MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n));
     83:   PetscCall(MatSetFromOptions(J));
     84:   PetscCall(MatSeqAIJSetPreallocation(J,3,NULL));
     85:   PetscCall(MatMPIAIJSetPreallocation(F,3,NULL,1,NULL));
     86:   PetscCall(NEPSetJacobian(nep,J,FormJacobian,&ctx));

     88:   PetscCall(NEPSetType(nep,NEPRII));
     89:   PetscCall(NEPRIISetKSP(nep,ksp));
     90:   PetscCall(NEPRIISetMaximumIterations(nep,6));
     91:   PetscCall(NEPRIISetConstCorrectionTol(nep,PETSC_TRUE));
     92:   PetscCall(NEPRIISetHermitian(nep,PETSC_TRUE));
     93:   PetscCall(NEPSetFromOptions(nep));

     95:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     96:                       Solve the eigensystem
     97:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     99:   PetscCall(NEPSolve(nep));
    100:   PetscCall(PetscObjectTypeCompare((PetscObject)nep,NEPRII,&flg));
    101:   if (flg) {
    102:     PetscCall(NEPRIIGetMaximumIterations(nep,&its));
    103:     PetscCall(NEPRIIGetLagPreconditioner(nep,&lag));
    104:     PetscCall(NEPRIIGetDeflationThreshold(nep,&thres));
    105:     PetscCall(NEPRIIGetConstCorrectionTol(nep,&cct));
    106:     PetscCall(NEPRIIGetHermitian(nep,&herm));
    107:     PetscCall(PetscPrintf(PETSC_COMM_WORLD," Maximum inner iterations of RII is %" PetscInt_FMT "\n",its));
    108:     PetscCall(PetscPrintf(PETSC_COMM_WORLD," Preconditioner rebuilt every %" PetscInt_FMT " iterations\n",lag));
    109:     if (thres>0.0) PetscCall(PetscPrintf(PETSC_COMM_WORLD," Using deflation threshold=%g\n",(double)thres));
    110:     if (cct) PetscCall(PetscPrintf(PETSC_COMM_WORLD," Using a constant correction tolerance\n"));
    111:     if (herm) PetscCall(PetscPrintf(PETSC_COMM_WORLD," Hermitian version of scalar equation\n"));
    112:     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n"));
    113:   }

    115:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    116:                     Display solution and clean up
    117:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    119:   /* show detailed info unless -terse option is given by user */
    120:   PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
    121:   if (terse) PetscCall(NEPErrorView(nep,NEP_ERROR_RELATIVE,NULL));
    122:   else {
    123:     PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
    124:     PetscCall(NEPConvergedReasonView(nep,PETSC_VIEWER_STDOUT_WORLD));
    125:     PetscCall(NEPErrorView(nep,NEP_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
    126:     PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
    127:   }

    129:   PetscCall(NEPDestroy(&nep));
    130:   PetscCall(KSPDestroy(&ksp));
    131:   PetscCall(MatDestroy(&F));
    132:   PetscCall(MatDestroy(&J));
    133:   PetscCall(SlepcFinalize());
    134:   return 0;
    135: }

    137: /* ------------------------------------------------------------------- */
    138: /*
    139:    FormFunction - Computes Function matrix  T(lambda)

    141:    Input Parameters:
    142: .  nep    - the NEP context
    143: .  lambda - the scalar argument
    144: .  ctx    - optional user-defined context, as set by NEPSetFunction()

    146:    Output Parameters:
    147: .  fun - Function matrix
    148: .  B   - optionally different preconditioning matrix
    149: */
    150: PetscErrorCode FormFunction(NEP nep,PetscScalar lambda,Mat fun,Mat B,void *ctx)
    151: {
    152:   ApplicationCtx *user = (ApplicationCtx*)ctx;
    153:   PetscScalar    A[3],c,d;
    154:   PetscReal      h;
    155:   PetscInt       i,n,j[3],Istart,Iend;
    156:   PetscBool      FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;

    158:   PetscFunctionBeginUser;
    159:   /*
    160:      Compute Function entries and insert into matrix
    161:   */
    162:   PetscCall(MatGetSize(fun,&n,NULL));
    163:   PetscCall(MatGetOwnershipRange(fun,&Istart,&Iend));
    164:   if (Istart==0) FirstBlock=PETSC_TRUE;
    165:   if (Iend==n) LastBlock=PETSC_TRUE;
    166:   h = user->h;
    167:   c = user->kappa/(lambda-user->kappa);
    168:   d = n;

    170:   /*
    171:      Interior grid points
    172:   */
    173:   for (i=(FirstBlock? Istart+1: Istart);i<(LastBlock? Iend-1: Iend);i++) {
    174:     j[0] = i-1; j[1] = i; j[2] = i+1;
    175:     A[0] = A[2] = -d-lambda*h/6.0; A[1] = 2.0*(d-lambda*h/3.0);
    176:     PetscCall(MatSetValues(fun,1,&i,3,j,A,INSERT_VALUES));
    177:   }

    179:   /*
    180:      Boundary points
    181:   */
    182:   if (FirstBlock) {
    183:     i = 0;
    184:     j[0] = 0; j[1] = 1;
    185:     A[0] = 2.0*(d-lambda*h/3.0); A[1] = -d-lambda*h/6.0;
    186:     PetscCall(MatSetValues(fun,1,&i,2,j,A,INSERT_VALUES));
    187:   }

    189:   if (LastBlock) {
    190:     i = n-1;
    191:     j[0] = n-2; j[1] = n-1;
    192:     A[0] = -d-lambda*h/6.0; A[1] = d-lambda*h/3.0+lambda*c;
    193:     PetscCall(MatSetValues(fun,1,&i,2,j,A,INSERT_VALUES));
    194:   }

    196:   /*
    197:      Assemble matrix
    198:   */
    199:   PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
    200:   PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
    201:   if (fun != B) {
    202:     PetscCall(MatAssemblyBegin(fun,MAT_FINAL_ASSEMBLY));
    203:     PetscCall(MatAssemblyEnd(fun,MAT_FINAL_ASSEMBLY));
    204:   }
    205:   PetscFunctionReturn(PETSC_SUCCESS);
    206: }

    208: /* ------------------------------------------------------------------- */
    209: /*
    210:    FormJacobian - Computes Jacobian matrix  T'(lambda)

    212:    Input Parameters:
    213: .  nep    - the NEP context
    214: .  lambda - the scalar argument
    215: .  ctx    - optional user-defined context, as set by NEPSetJacobian()

    217:    Output Parameters:
    218: .  jac - Jacobian matrix
    219: .  B   - optionally different preconditioning matrix
    220: */
    221: PetscErrorCode FormJacobian(NEP nep,PetscScalar lambda,Mat jac,void *ctx)
    222: {
    223:   ApplicationCtx *user = (ApplicationCtx*)ctx;
    224:   PetscScalar    A[3],c;
    225:   PetscReal      h;
    226:   PetscInt       i,n,j[3],Istart,Iend;
    227:   PetscBool      FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;

    229:   PetscFunctionBeginUser;
    230:   /*
    231:      Compute Jacobian entries and insert into matrix
    232:   */
    233:   PetscCall(MatGetSize(jac,&n,NULL));
    234:   PetscCall(MatGetOwnershipRange(jac,&Istart,&Iend));
    235:   if (Istart==0) FirstBlock=PETSC_TRUE;
    236:   if (Iend==n) LastBlock=PETSC_TRUE;
    237:   h = user->h;
    238:   c = user->kappa/(lambda-user->kappa);

    240:   /*
    241:      Interior grid points
    242:   */
    243:   for (i=(FirstBlock? Istart+1: Istart);i<(LastBlock? Iend-1: Iend);i++) {
    244:     j[0] = i-1; j[1] = i; j[2] = i+1;
    245:     A[0] = A[2] = -h/6.0; A[1] = -2.0*h/3.0;
    246:     PetscCall(MatSetValues(jac,1,&i,3,j,A,INSERT_VALUES));
    247:   }

    249:   /*
    250:      Boundary points
    251:   */
    252:   if (FirstBlock) {
    253:     i = 0;
    254:     j[0] = 0; j[1] = 1;
    255:     A[0] = -2.0*h/3.0; A[1] = -h/6.0;
    256:     PetscCall(MatSetValues(jac,1,&i,2,j,A,INSERT_VALUES));
    257:   }

    259:   if (LastBlock) {
    260:     i = n-1;
    261:     j[0] = n-2; j[1] = n-1;
    262:     A[0] = -h/6.0; A[1] = -h/3.0-c*c;
    263:     PetscCall(MatSetValues(jac,1,&i,2,j,A,INSERT_VALUES));
    264:   }

    266:   /*
    267:      Assemble matrix
    268:   */
    269:   PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY));
    270:   PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY));
    271:   PetscFunctionReturn(PETSC_SUCCESS);
    272: }

    274: /*TEST

    276:    test:
    277:       suffix: 1
    278:       args: -nep_target 21 -nep_rii_lag_preconditioner 2 -terse
    279:       requires: !single

    281: TEST*/