Actual source code: test44.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[] = "Eigenvalue problem with Bethe-Salpeter structure using shell matrices.\n\n"
     12:   "The command line options are:\n"
     13:   "  -n <n>, where <n> = dimension of the blocks.\n\n";

     15: #include <slepceps.h>

     17: /*
     18:    This example computes eigenvalues of a matrix

     20:         H = [  R    C
     21:               -C^H -R^T ],

     23:    where R is Hermitian and C is complex symmetric. In particular, R and C have the
     24:    following Toeplitz structure:

     26:         R = pentadiag{a,b,c,conj(b),conj(a)}
     27:         C = tridiag{b,d,b}

     29:    where a,b,d are complex scalars, and c is real.
     30: */

     32: /*
     33:    User-defined routines
     34: */
     35: PetscErrorCode MatMult_R(Mat R,Vec x,Vec y);
     36: PetscErrorCode MatMultTranspose_R(Mat R,Vec x,Vec y);
     37: PetscErrorCode MatMult_C(Mat C,Vec x,Vec y);
     38: PetscErrorCode MatMultHermitianTranspose_C(Mat C,Vec x,Vec y);

     40: /*
     41:    User context for shell matrices
     42: */
     43: typedef struct {
     44:   PetscScalar a,b,c,d;
     45: } CTX_SHELL;

     47: int main(int argc,char **argv)
     48: {
     49:   Mat            H,R,C;      /* problem matrices */
     50:   EPS            eps;        /* eigenproblem solver context */
     51:   PetscReal      lev;
     52:   PetscInt       n=24,i,nconv;
     53:   PetscMPIInt    size;
     54:   PetscBool      terse,checkorthog;
     55:   Vec            t,*x,*y;
     56:   CTX_SHELL      *ctx;

     58:   PetscFunctionBeginUser;
     59:   PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
     60:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size));
     61:   PetscCheck(size==1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only");

     63:   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
     64:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nShell Bethe-Salpeter eigenproblem, n=%" PetscInt_FMT "\n\n",n));

     66:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     67:                Generate the shell problem matrices R and C
     68:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     70:   PetscCall(PetscNew(&ctx));
     71: #if defined(PETSC_USE_COMPLEX)
     72:   ctx->a = PetscCMPLX(-0.1,0.2);
     73:   ctx->b = PetscCMPLX(1.0,0.5);
     74:   ctx->d = PetscCMPLX(2.0,0.2);
     75: #else
     76:   ctx->a = -0.1;
     77:   ctx->b = 1.0;
     78:   ctx->d = 2.0;
     79: #endif
     80:   ctx->c = 4.5;

     82:   PetscCall(MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,(void*)ctx,&R));
     83:   PetscCall(MatShellSetOperation(R,MATOP_MULT,(PetscErrorCodeFn*)MatMult_R));
     84:   PetscCall(MatShellSetOperation(R,MATOP_MULT_TRANSPOSE,(PetscErrorCodeFn*)MatMultTranspose_R));
     85:   PetscCall(MatSetOption(R,MAT_HERMITIAN,PETSC_TRUE));
     86:   PetscCall(MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,(void*)ctx,&C));
     87:   PetscCall(MatShellSetOperation(C,MATOP_MULT,(PetscErrorCodeFn*)MatMult_C));
     88:   PetscCall(MatShellSetOperation(C,MATOP_MULT_HERMITIAN_TRANSPOSE,(PetscErrorCodeFn*)MatMultHermitianTranspose_C));
     89:   PetscCall(MatSetOption(C,MAT_SYMMETRIC,PETSC_TRUE));

     91:   PetscCall(MatCreateBSE(R,C,&H));
     92:   PetscCall(MatDestroy(&R));
     93:   PetscCall(MatDestroy(&C));

     95:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     96:                 Create the eigensolver and set various options
     97:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     99:   PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
    100:   PetscCall(EPSSetOperators(eps,H,NULL));
    101:   PetscCall(EPSSetProblemType(eps,EPS_BSE));
    102:   PetscCall(EPSSetWhichEigenpairs(eps,EPS_SMALLEST_MAGNITUDE));
    103:   PetscCall(EPSSetFromOptions(eps));

    105:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    106:                  Solve the eigensystem and display solution
    107:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    109:   PetscCall(EPSSolve(eps));

    111:   /* show detailed info unless -terse option is given by user */
    112:   PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
    113:   if (terse) PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
    114:   else {
    115:     PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
    116:     PetscCall(EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD));
    117:     PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
    118:     PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
    119:   }

    121:   /* check bi-orthogonality */
    122:   PetscCall(PetscOptionsHasName(NULL,NULL,"-checkorthog",&checkorthog));
    123:   PetscCall(EPSGetConverged(eps,&nconv));
    124:   if (checkorthog && nconv>0) {
    125:     PetscCall(MatCreateVecs(H,&t,NULL));
    126:     PetscCall(VecDuplicateVecs(t,nconv,&x));
    127:     PetscCall(VecDuplicateVecs(t,nconv,&y));
    128:     for (i=0;i<nconv;i++) {
    129:       PetscCall(EPSGetEigenvector(eps,i,x[i],NULL));
    130:       PetscCall(EPSGetLeftEigenvector(eps,i,y[i],NULL));
    131:     }
    132:     PetscCall(VecCheckOrthogonality(x,nconv,y,nconv,NULL,NULL,&lev));
    133:     if (lev<100*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD," Level of bi-orthogonality of eigenvectors < 100*eps\n\n"));
    134:     else PetscCall(PetscPrintf(PETSC_COMM_WORLD," Level of bi-orthogonality of eigenvectors: %g\n\n",(double)lev));
    135:     PetscCall(VecDestroy(&t));
    136:     PetscCall(VecDestroyVecs(nconv,&x));
    137:     PetscCall(VecDestroyVecs(nconv,&y));
    138:   }

    140:   PetscCall(EPSDestroy(&eps));
    141:   PetscCall(MatDestroy(&H));
    142:   PetscCall(PetscFree(ctx));
    143:   PetscCall(SlepcFinalize());
    144:   return 0;
    145: }

    147: /*
    148:     Matrix-vector y = R*x.

    150:     R = pentadiag{a,b,c,conj(b),conj(a)}
    151:  */
    152: PetscErrorCode MatMult_R(Mat R,Vec x,Vec y)
    153: {
    154:   CTX_SHELL         *ctx;
    155:   PetscInt          n,i;
    156:   const PetscScalar *px;
    157:   PetscScalar       *py;

    159:   PetscFunctionBeginUser;
    160:   PetscCall(MatShellGetContext(R,&ctx));
    161:   PetscCall(MatGetSize(R,NULL,&n));
    162:   PetscCall(VecGetArrayRead(x,&px));
    163:   PetscCall(VecGetArray(y,&py));
    164:   for (i=0;i<n;i++) {
    165:     py[i] = ctx->c*px[i];
    166:     if (i>1) py[i] += ctx->a*px[i-2];
    167:     if (i>0) py[i] += ctx->b*px[i-1];
    168:     if (i<n-1) py[i] += PetscConj(ctx->b)*px[i+1];
    169:     if (i<n-2) py[i] += PetscConj(ctx->a)*px[i+2];
    170:   }
    171:   PetscCall(VecRestoreArrayRead(x,&px));
    172:   PetscCall(VecRestoreArray(y,&py));
    173:   PetscFunctionReturn(PETSC_SUCCESS);
    174: }

    176: /*
    177:     Matrix-vector y = R^T*x.

    179:     Only needed to compute the residuals.
    180:  */
    181: PetscErrorCode MatMultTranspose_R(Mat R,Vec x,Vec y)
    182: {
    183:   Vec w;

    185:   PetscFunctionBeginUser;
    186:   PetscCall(VecDuplicate(x,&w));
    187:   PetscCall(VecCopy(x,w));
    188:   PetscCall(VecConjugate(w));
    189:   PetscCall(MatMult_R(R,w,y));
    190:   PetscCall(VecConjugate(y));
    191:   PetscCall(VecDestroy(&w));
    192:   PetscFunctionReturn(PETSC_SUCCESS);
    193: }

    195: /*
    196:     Matrix-vector y = C*x.

    198:     C = tridiag{b,d,b}
    199:  */
    200: PetscErrorCode MatMult_C(Mat C,Vec x,Vec y)
    201: {
    202:   CTX_SHELL         *ctx;
    203:   PetscInt          n,i;
    204:   const PetscScalar *px;
    205:   PetscScalar       *py;

    207:   PetscFunctionBeginUser;
    208:   PetscCall(MatShellGetContext(C,&ctx));
    209:   PetscCall(MatGetSize(C,NULL,&n));
    210:   PetscCall(VecGetArrayRead(x,&px));
    211:   PetscCall(VecGetArray(y,&py));
    212:   for (i=0;i<n;i++) {
    213:     py[i] = ctx->d*px[i];
    214:     if (i>0) py[i] += ctx->b*px[i-1];
    215:     if (i<n-1) py[i] += ctx->b*px[i+1];
    216:   }
    217:   PetscCall(VecRestoreArrayRead(x,&px));
    218:   PetscCall(VecRestoreArray(y,&py));
    219:   PetscFunctionReturn(PETSC_SUCCESS);
    220: }

    222: /*
    223:     Matrix-vector y = C^H*x.

    225:     Only needed to compute the residuals.
    226:  */
    227: PetscErrorCode MatMultHermitianTranspose_C(Mat C,Vec x,Vec y)
    228: {
    229:   Vec w;

    231:   PetscFunctionBeginUser;
    232:   PetscCall(VecDuplicate(x,&w));
    233:   PetscCall(VecCopy(x,w));
    234:   PetscCall(VecConjugate(w));
    235:   PetscCall(MatMult_C(C,w,y));
    236:   PetscCall(VecConjugate(y));
    237:   PetscCall(VecDestroy(&w));
    238:   PetscFunctionReturn(PETSC_SUCCESS);
    239: }

    241: /*TEST

    243:    testset:
    244:       args: -eps_nev 4 -eps_ncv 16 -eps_krylovschur_bse_type {{shao gruning}} -terse -checkorthog
    245:       filter: sed -e "s/17496/17495/g" | sed -e "s/32172/32173/g" | sed -e "s/38566/38567/g"
    246:       test:
    247:          suffix: 1
    248:          requires: complex
    249:       test:
    250:          suffix: 1_real
    251:          requires: !complex

    253: TEST*/