Actual source code: ex55.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.\n\n"
     12:   "The command line options are:\n"
     13:   "  -n <n>, where <n> = dimension of the blocks.\n"
     14:   "  -R <filename>, R matrix.\n"
     15:   "  -C <filename>, C matrix.\n\n";

     17: #include <slepceps.h>

     19: /*
     20:    This example computes eigenvalues of a matrix

     22:         H = [  R    C
     23:               -C^H -R^T ],

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

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

     31:    where a,b,d are complex scalars, and c is real.

     33:    Alternatively, use -R and -C command-line options to load matrices from file.
     34: */

     36: int main(int argc,char **argv)
     37: {
     38:   Mat            H,R,C;      /* problem matrices */
     39:   EPS            eps;        /* eigenproblem solver context */
     40:   PetscScalar    a,b,c,d;
     41:   PetscReal      lev;
     42:   PetscInt       n=24,Istart,Iend,i,nconv;
     43:   PetscBool      flg,terse,checkorthog,nest=PETSC_FALSE;
     44:   Vec            t,*x,*y;
     45:   PetscViewer    viewer;
     46:   char           filename[PETSC_MAX_PATH_LEN];

     48:   PetscFunctionBeginUser;
     49:   PetscCall(SlepcInitialize(&argc,&argv,NULL,help));

     51:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     52:              Compute (or load) the problem matrices R and C
     53:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     55:   PetscCall(PetscOptionsHasName(NULL,NULL,"-R",&flg));

     57:   if (flg) {

     59:     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nBethe-Salpeter eigenproblem stored in file\n\n"));
     60:     PetscCall(PetscPrintf(PETSC_COMM_WORLD," Reading %s matrices from a binary file...\n",PetscDefined(USE_COMPLEX)?"COMPLEX":"REAL"));

     62:     PetscCall(PetscOptionsGetString(NULL,NULL,"-R",filename,sizeof(filename),&flg));
     63:     PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
     64:     PetscCall(MatCreate(PETSC_COMM_WORLD,&R));
     65:     PetscCall(MatSetFromOptions(R));
     66:     PetscCall(MatLoad(R,viewer));
     67:     PetscCall(PetscViewerDestroy(&viewer));

     69:     PetscCall(PetscOptionsGetString(NULL,NULL,"-C",filename,sizeof(filename),&flg));
     70:     PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must indicate file for both R and C matrices");
     71:     PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
     72:     PetscCall(MatCreate(PETSC_COMM_WORLD,&C));
     73:     PetscCall(MatSetFromOptions(C));
     74:     PetscCall(MatLoad(C,viewer));
     75:     PetscCall(PetscViewerDestroy(&viewer));

     77:   } else {

     79:     PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
     80:     PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nBethe-Salpeter eigenproblem, n=%" PetscInt_FMT "\n\n",n));

     82: #if defined(PETSC_USE_COMPLEX)
     83:     a = PetscCMPLX(-0.1,0.2);
     84:     b = PetscCMPLX(1.0,0.5);
     85:     d = PetscCMPLX(2.0,0.2);
     86: #else
     87:     a = -0.1;
     88:     b = 1.0;
     89:     d = 2.0;
     90: #endif
     91:     c = 4.5;

     93:     PetscCall(MatCreate(PETSC_COMM_WORLD,&R));
     94:     PetscCall(MatSetSizes(R,PETSC_DECIDE,PETSC_DECIDE,n,n));
     95:     PetscCall(MatSetFromOptions(R));

     97:     PetscCall(MatCreate(PETSC_COMM_WORLD,&C));
     98:     PetscCall(MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,n,n));
     99:     PetscCall(MatSetFromOptions(C));

    101:     PetscCall(MatGetOwnershipRange(R,&Istart,&Iend));
    102:     for (i=Istart;i<Iend;i++) {
    103:       if (i>1) PetscCall(MatSetValue(R,i,i-2,a,INSERT_VALUES));
    104:       if (i>0) PetscCall(MatSetValue(R,i,i-1,b,INSERT_VALUES));
    105:       PetscCall(MatSetValue(R,i,i,c,INSERT_VALUES));
    106:       if (i<n-1) PetscCall(MatSetValue(R,i,i+1,PetscConj(b),INSERT_VALUES));
    107:       if (i<n-2) PetscCall(MatSetValue(R,i,i+2,PetscConj(a),INSERT_VALUES));
    108:     }

    110:     PetscCall(MatGetOwnershipRange(C,&Istart,&Iend));
    111:     for (i=Istart;i<Iend;i++) {
    112:       if (i>0) PetscCall(MatSetValue(C,i,i-1,b,INSERT_VALUES));
    113:       PetscCall(MatSetValue(C,i,i,d,INSERT_VALUES));
    114:       if (i<n-1) PetscCall(MatSetValue(C,i,i+1,b,INSERT_VALUES));
    115:     }

    117:     PetscCall(MatAssemblyBegin(R,MAT_FINAL_ASSEMBLY));
    118:     PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY));
    119:     PetscCall(MatAssemblyEnd(R,MAT_FINAL_ASSEMBLY));
    120:     PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY));
    121:   }

    123:   PetscCall(MatCreateBSE(R,C,&H));

    125:   /* if you prefer, set the vector type so that MatCreateVecs() returns nested vectors */
    126:   PetscCall(PetscOptionsGetBool(NULL,NULL,"-nest",&nest,NULL));
    127:   if (nest) PetscCall(MatNestSetVecType(H,VECNEST));

    129:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    130:                 Create the eigensolver and set various options
    131:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    133:   PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
    134:   PetscCall(EPSSetOperators(eps,H,NULL));
    135:   PetscCall(EPSSetProblemType(eps,EPS_BSE));
    136:   PetscCall(EPSSetFromOptions(eps));

    138:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    139:                  Solve the eigensystem and display solution
    140:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    142:   PetscCall(EPSSolve(eps));

    144:   /* show detailed info unless -terse option is given by user */
    145:   PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
    146:   if (terse) PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
    147:   else {
    148:     PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
    149:     PetscCall(EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD));
    150:     PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
    151:     PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
    152:   }

    154:   /* check bi-orthogonality */
    155:   PetscCall(PetscOptionsHasName(NULL,NULL,"-checkorthog",&checkorthog));
    156:   PetscCall(EPSGetConverged(eps,&nconv));
    157:   if (checkorthog && nconv>0) {
    158:     PetscCall(MatCreateVecs(H,&t,NULL));
    159:     PetscCall(VecDuplicateVecs(t,nconv,&x));
    160:     PetscCall(VecDuplicateVecs(t,nconv,&y));
    161:     for (i=0;i<nconv;i++) {
    162:       PetscCall(EPSGetEigenvector(eps,i,x[i],NULL));
    163:       PetscCall(EPSGetLeftEigenvector(eps,i,y[i],NULL));
    164:     }
    165:     PetscCall(VecCheckOrthogonality(x,nconv,y,nconv,NULL,NULL,&lev));
    166:     if (lev<100*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD," Level of bi-orthogonality of eigenvectors < 100*eps\n\n"));
    167:     else PetscCall(PetscPrintf(PETSC_COMM_WORLD," Level of bi-orthogonality of eigenvectors: %g\n\n",(double)lev));
    168:     PetscCall(VecDestroy(&t));
    169:     PetscCall(VecDestroyVecs(nconv,&x));
    170:     PetscCall(VecDestroyVecs(nconv,&y));
    171:   }

    173:   PetscCall(EPSDestroy(&eps));
    174:   PetscCall(MatDestroy(&R));
    175:   PetscCall(MatDestroy(&C));
    176:   PetscCall(MatDestroy(&H));
    177:   PetscCall(SlepcFinalize());
    178:   return 0;
    179: }

    181: /*TEST

    183:    testset:
    184:       args: -eps_nev 4 -eps_ncv 16 -eps_krylovschur_bse_type {{shao gruning projectedbse}} -terse -checkorthog -nest {{0 1}}
    185:       filter: sed -e "s/17496/17495/g" | sed -e "s/38566/38567/g" | sed -e "s/32172/32173/g"
    186:       nsize: {{1 2}}
    187:       test:
    188:          suffix: 1
    189:          requires: complex
    190:       test:
    191:          suffix: 1_real
    192:          requires: !complex
    193:       test:
    194:          suffix: 1_dense
    195:          args: -mat_type dense
    196:          requires: complex
    197:          output_file: output/ex55_1.out
    198:       test:
    199:          suffix: 1_cuda
    200:          args: -mat_type aijcusparse
    201:          requires: cuda complex
    202:          output_file: output/ex55_1.out
    203:       test:
    204:          suffix: 1_real_cuda
    205:          args: -mat_type aijcusparse
    206:          requires: cuda !complex
    207:          output_file: output/ex55_1_real.out
    208:       test:
    209:          suffix: 1_hip
    210:          args: -mat_type aijhipsparse
    211:          requires: hip complex
    212:          output_file: output/ex55_1.out
    213:       test:
    214:          suffix: 1_real_hip
    215:          args: -mat_type aijhipsparse
    216:          requires: hip !complex
    217:          output_file: output/ex55_1_real.out

    219:    testset:
    220:       args: -eps_nev 4 -eps_ncv 16 -eps_krylovschur_bse_type {{shao gruning projectedbse}} -st_type sinvert -terse
    221:       filter: sed -e "s/17496/17495/g" | sed -e "s/38566/38567/g"
    222:       test:
    223:          suffix: 1_sinvert
    224:          requires: complex
    225:       test:
    226:          nsize: 4
    227:          args: -mat_type scalapack
    228:          suffix: 1_sinvert_scalapack
    229:          requires: complex scalapack !__float128
    230:          output_file: output/ex55_1_sinvert.out
    231:       test:
    232:          suffix: 1_real_sinvert
    233:          requires: !complex
    234:       test:
    235:          nsize: 4
    236:          args: -mat_type scalapack
    237:          suffix: 1_real_sinvert_scalapack
    238:          requires: !complex scalapack !__float128
    239:          output_file: output/ex55_1_real_sinvert.out

    241:    testset:
    242:       args: -n 90 -eps_threshold_absolute 2.5 -eps_ncv {{10 24}} -eps_krylovschur_bse_type {{shao gruning projectedbse}} -eps_tol 1e-14 -terse -checkorthog
    243:       test:
    244:          suffix: 2
    245:          requires: double complex
    246:       test:
    247:          suffix: 2_real
    248:          requires: double !complex

    250:    testset:
    251:       args: -eps_nev 28 -eps_ncv 18 -terse
    252:       test:
    253:          suffix: 3
    254:          requires: !complex !single

    256: TEST*/