Actual source code: ex10.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[] = "Illustrates the use of shell spectral transformations.\n\n"
     12:   "The problem to be solved is the same as ex1.c and"
     13:   "corresponds to the Laplacian operator in 1 dimension.\n\n"
     14:   "The command line options are:\n"
     15:   "  -n <n>, where <n> = number of grid subdivisions = matrix dimension.\n\n";

     17: #include <slepceps.h>

     19: /* Define context for user-provided spectral transformation */
     20: typedef struct {
     21:   KSP        ksp;
     22: } SampleShellST;

     24: /* Declare routines for user-provided spectral transformation */
     25: PetscErrorCode STCreate_User(SampleShellST**);
     26: PetscErrorCode STSetUp_User(SampleShellST*,ST);
     27: PetscErrorCode STApply_User(ST,Vec,Vec);
     28: PetscErrorCode STApplyTranspose_User(ST,Vec,Vec);
     29: #if defined(PETSC_USE_COMPLEX)
     30: PetscErrorCode STApplyHermitianTranspose_User(ST,Vec,Vec);
     31: #endif
     32: PetscErrorCode STBackTransform_User(ST,PetscInt,PetscScalar*,PetscScalar*);
     33: PetscErrorCode STDestroy_User(ST);

     35: int main (int argc,char **argv)
     36: {
     37:   Mat            A;               /* operator matrix */
     38:   EPS            eps;             /* eigenproblem solver context */
     39:   ST             st;              /* spectral transformation context */
     40:   SampleShellST  *shell;          /* user-defined spectral transform context */
     41:   EPSType        type;
     42:   PetscInt       n=30,i,Istart,Iend,nev;
     43:   PetscBool      isShell,terse;
     44:   PetscBool      set_ht=PETSC_FALSE;

     46:   PetscFunctionBeginUser;
     47:   PetscCall(SlepcInitialize(&argc,&argv,NULL,help));

     49:   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
     50:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n1-D Laplacian Eigenproblem (shell-enabled), n=%" PetscInt_FMT "\n\n",n));
     51:   PetscCall(PetscOptionsGetBool(NULL,NULL,"-set_ht",&set_ht,NULL));

     53:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     54:      Compute the operator matrix that defines the eigensystem, Ax=kx
     55:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     57:   PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
     58:   PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n,n));
     59:   PetscCall(MatSetFromOptions(A));

     61:   PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
     62:   for (i=Istart;i<Iend;i++) {
     63:     if (i>0) PetscCall(MatSetValue(A,i,i-1,-1.0,INSERT_VALUES));
     64:     if (i<n-1) PetscCall(MatSetValue(A,i,i+1,-1.0,INSERT_VALUES));
     65:     PetscCall(MatSetValue(A,i,i,2.0,INSERT_VALUES));
     66:   }
     67:   PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
     68:   PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));

     70:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     71:                 Create the eigensolver and set various options
     72:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     74:   /*
     75:      Create eigensolver context
     76:   */
     77:   PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));

     79:   /*
     80:      Set operators. In this case, it is a standard eigenvalue problem
     81:   */
     82:   PetscCall(EPSSetOperators(eps,A,NULL));
     83:   PetscCall(EPSSetProblemType(eps,EPS_HEP));

     85:   /*
     86:      Set solver parameters at runtime
     87:   */
     88:   PetscCall(EPSSetFromOptions(eps));

     90:   /*
     91:      Initialize shell spectral transformation if selected by user
     92:   */
     93:   PetscCall(EPSGetST(eps,&st));
     94:   PetscCall(PetscObjectTypeCompare((PetscObject)st,STSHELL,&isShell));
     95:   if (isShell) {
     96:     /* Change sorting criterion since this ST example computes values
     97:        closest to 0 */
     98:     PetscCall(EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL));

    100:     /* (Required) Create a context for the user-defined spectral transform;
    101:        this context can be defined to contain any application-specific data. */
    102:     PetscCall(STCreate_User(&shell));
    103:     PetscCall(STShellSetContext(st,shell));
    104:     PetscCall(STShellSetDestroy(st,STDestroy_User));

    106:     /* (Required) Set the user-defined routine for applying the operator */
    107:     PetscCall(STShellSetApply(st,STApply_User));

    109:     /* (Optional) Set the user-defined routine for applying the transposed operator */
    110:     PetscCall(STShellSetApplyTranspose(st,STApplyTranspose_User));

    112:     /* (Optional) Set the user-defined routine for applying the conjugate-transposed operator */
    113: #if defined(PETSC_USE_COMPLEX)
    114:     if (set_ht) PetscCall(STShellSetApplyHermitianTranspose(st,STApplyHermitianTranspose_User));
    115: #endif

    117:     /* (Optional) Set the user-defined routine for back-transformation */
    118:     PetscCall(STShellSetBackTransform(st,STBackTransform_User));

    120:     /* (Optional) Set a name for the transformation, used for STView() */
    121:     PetscCall(PetscObjectSetName((PetscObject)st,"MyTransformation"));

    123:     /* (Optional) Do any setup required for the new transformation */
    124:     PetscCall(STSetUp_User(shell,st));
    125:   }

    127:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    128:                       Solve the eigensystem
    129:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    131:   PetscCall(EPSSolve(eps));

    133:   /*
    134:      Optional: Get some information from the solver and display it
    135:   */
    136:   PetscCall(EPSGetType(eps,&type));
    137:   PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type));
    138:   PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
    139:   PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));

    141:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    142:                     Display solution and clean up
    143:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    145:   /* show detailed info unless -terse option is given by user */
    146:   PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
    147:   if (terse) PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
    148:   else {
    149:     PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
    150:     PetscCall(EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD));
    151:     PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
    152:     PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
    153:   }
    154:   PetscCall(EPSDestroy(&eps));
    155:   PetscCall(MatDestroy(&A));
    156:   PetscCall(SlepcFinalize());
    157:   return 0;
    158: }

    160: /***********************************************************************/
    161: /*     Routines for a user-defined shell spectral transformation       */
    162: /***********************************************************************/

    164: /*
    165:    STCreate_User - This routine creates a user-defined
    166:    spectral transformation context.

    168:    Output Parameter:
    169: .  shell - user-defined spectral transformation context
    170: */
    171: PetscErrorCode STCreate_User(SampleShellST **shell)
    172: {
    173:   SampleShellST  *newctx;

    175:   PetscFunctionBeginUser;
    176:   PetscCall(PetscNew(&newctx));
    177:   PetscCall(KSPCreate(PETSC_COMM_WORLD,&newctx->ksp));
    178:   PetscCall(KSPAppendOptionsPrefix(newctx->ksp,"st_"));
    179:   *shell = newctx;
    180:   PetscFunctionReturn(PETSC_SUCCESS);
    181: }
    182: /* ------------------------------------------------------------------- */
    183: /*
    184:    STSetUp_User - This routine sets up a user-defined
    185:    spectral transformation context.

    187:    Input Parameters:
    188: +  shell - user-defined spectral transformation context
    189: -  st    - spectral transformation context containing the operator matrices

    191:    Output Parameter:
    192: .  shell - fully set up user-defined transformation context

    194:    Notes:
    195:    In this example, the user-defined transformation is simply OP=A^-1.
    196:    Therefore, the eigenpairs converge in reversed order. The KSP object
    197:    used for the solution of linear systems with A is handled via the
    198:    user-defined context SampleShellST.
    199: */
    200: PetscErrorCode STSetUp_User(SampleShellST *shell,ST st)
    201: {
    202:   Mat            A;

    204:   PetscFunctionBeginUser;
    205:   PetscCall(STGetMatrix(st,0,&A));
    206:   PetscCall(KSPSetOperators(shell->ksp,A,A));
    207:   PetscCall(KSPSetFromOptions(shell->ksp));
    208:   PetscFunctionReturn(PETSC_SUCCESS);
    209: }
    210: /* ------------------------------------------------------------------- */
    211: /*
    212:    STApply_User - This routine demonstrates the use of a
    213:    user-provided spectral transformation.

    215:    Input Parameters:
    216: +  st - spectral transformation context
    217: -  x - input vector

    219:    Output Parameter:
    220: .  y - output vector

    222:    Notes:
    223:    The transformation implemented in this code is just OP=A^-1 and
    224:    therefore it is of little use, merely as an example of working with
    225:    a STSHELL.
    226: */
    227: PetscErrorCode STApply_User(ST st,Vec x,Vec y)
    228: {
    229:   SampleShellST  *shell;

    231:   PetscFunctionBeginUser;
    232:   PetscCall(STShellGetContext(st,&shell));
    233:   PetscCall(KSPSolve(shell->ksp,x,y));
    234:   PetscFunctionReturn(PETSC_SUCCESS);
    235: }
    236: /* ------------------------------------------------------------------- */
    237: /*
    238:    STApplyTranspose_User - This is not required unless using a two-sided
    239:    eigensolver.

    241:    Input Parameters:
    242: +  st - spectral transformation context
    243: -  x - input vector

    245:    Output Parameter:
    246: .  y - output vector
    247: */
    248: PetscErrorCode STApplyTranspose_User(ST st,Vec x,Vec y)
    249: {
    250:   SampleShellST  *shell;

    252:   PetscFunctionBeginUser;
    253:   PetscCall(STShellGetContext(st,&shell));
    254:   PetscCall(KSPSolveTranspose(shell->ksp,x,y));
    255:   PetscFunctionReturn(PETSC_SUCCESS);
    256: }
    257: #if defined(PETSC_USE_COMPLEX)
    258: /* ------------------------------------------------------------------- */
    259: /*
    260:    STApplyHermitianTranspose_User - This is not required unless using a two-sided
    261:    eigensolver in complex scalars.

    263:    Input Parameters:
    264: +  st - spectral transformation context
    265: -  x - input vector

    267:    Output Parameter:
    268: .  y - output vector
    269: */
    270: PetscErrorCode STApplyHermitianTranspose_User(ST st,Vec x,Vec y)
    271: {
    272:   SampleShellST  *shell;
    273:   Vec            w;

    275:   PetscFunctionBeginUser;
    276:   PetscCall(STShellGetContext(st,&shell));
    277:   PetscCall(VecDuplicate(x,&w));
    278:   PetscCall(VecCopy(x,w));
    279:   PetscCall(VecConjugate(w));
    280:   PetscCall(KSPSolveTranspose(shell->ksp,w,y));
    281:   PetscCall(VecConjugate(y));
    282:   PetscCall(VecDestroy(&w));
    283:   PetscFunctionReturn(PETSC_SUCCESS);
    284: }
    285: #endif
    286: /* ------------------------------------------------------------------- */
    287: /*
    288:    STBackTransform_User - This routine demonstrates the use of a
    289:    user-provided spectral transformation.

    291:    Input Parameters:
    292: +  st - spectral transformation context
    293: -  n  - number of eigenvalues to transform

    295:    Input/Output Parameters:
    296: +  eigr - pointer to real part of eigenvalues
    297: -  eigi - pointer to imaginary part of eigenvalues

    299:    Notes:
    300:    This code implements the back transformation of eigenvalues in
    301:    order to retrieve the eigenvalues of the original problem. In this
    302:    example, simply set k_i = 1/k_i.
    303: */
    304: PetscErrorCode STBackTransform_User(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
    305: {
    306:   PetscInt j;

    308:   PetscFunctionBeginUser;
    309:   for (j=0;j<n;j++) {
    310:     eigr[j] = 1.0 / eigr[j];
    311:   }
    312:   PetscFunctionReturn(PETSC_SUCCESS);
    313: }
    314: /* ------------------------------------------------------------------- */
    315: /*
    316:    STDestroy_User - This routine destroys a user-defined
    317:    spectral transformation context.

    319:    Input Parameter:
    320: .  st - spectral transformation context
    321: */
    322: PetscErrorCode STDestroy_User(ST st)
    323: {
    324:   SampleShellST *shell;

    326:   PetscFunctionBeginUser;
    327:   PetscCall(STShellGetContext(st,&shell));
    328:   PetscCall(KSPDestroy(&shell->ksp));
    329:   PetscCall(PetscFree(shell));
    330:   PetscFunctionReturn(PETSC_SUCCESS);
    331: }

    333: /*TEST

    335:    testset:
    336:       args: -eps_nev 5 -eps_non_hermitian -terse
    337:       output_file: output/ex10_1.out
    338:       test:
    339:          suffix: 1_sinvert
    340:          args: -st_type sinvert
    341:          requires: !single
    342:       test:
    343:          suffix: 1_sinvert_twoside
    344:          args: -st_type sinvert -eps_balance twoside -set_ht {{0 1}}
    345:          requires: !single
    346:       test:
    347:          suffix: 1_shell
    348:          args: -st_type shell
    349:          requires: !single
    350:       test:
    351:          suffix: 1_shell_twoside
    352:          args: -st_type shell -eps_balance twoside -set_ht {{0 1}}

    354: TEST*/