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\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: int main(int argc,char **argv)
     33: {
     34:   Mat            H,R,C;      /* problem matrices */
     35:   EPS            eps;        /* eigenproblem solver context */
     36:   PetscScalar    a,b,c,d;
     37:   PetscReal      lev;
     38:   PetscInt       n=24,Istart,Iend,i,nconv;
     39:   PetscBool      terse,checkorthog,nest=PETSC_FALSE;
     40:   Vec            t,*x,*y;

     42:   PetscFunctionBeginUser;
     43:   PetscCall(SlepcInitialize(&argc,&argv,NULL,help));

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

     48:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     49:                Compute the problem matrices R and C
     50:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     52: #if defined(PETSC_USE_COMPLEX)
     53:   a = PetscCMPLX(-0.1,0.2);
     54:   b = PetscCMPLX(1.0,0.5);
     55:   d = PetscCMPLX(2.0,0.2);
     56: #else
     57:   a = -0.1;
     58:   b = 1.0;
     59:   d = 2.0;
     60: #endif
     61:   c = 4.5;

     63:   PetscCall(MatCreate(PETSC_COMM_WORLD,&R));
     64:   PetscCall(MatSetSizes(R,PETSC_DECIDE,PETSC_DECIDE,n,n));
     65:   PetscCall(MatSetFromOptions(R));

     67:   PetscCall(MatCreate(PETSC_COMM_WORLD,&C));
     68:   PetscCall(MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,n,n));
     69:   PetscCall(MatSetFromOptions(C));

     71:   PetscCall(MatGetOwnershipRange(R,&Istart,&Iend));
     72:   for (i=Istart;i<Iend;i++) {
     73:     if (i>1) PetscCall(MatSetValue(R,i,i-2,a,INSERT_VALUES));
     74:     if (i>0) PetscCall(MatSetValue(R,i,i-1,b,INSERT_VALUES));
     75:     PetscCall(MatSetValue(R,i,i,c,INSERT_VALUES));
     76:     if (i<n-1) PetscCall(MatSetValue(R,i,i+1,PetscConj(b),INSERT_VALUES));
     77:     if (i<n-2) PetscCall(MatSetValue(R,i,i+2,PetscConj(a),INSERT_VALUES));
     78:   }

     80:   PetscCall(MatGetOwnershipRange(C,&Istart,&Iend));
     81:   for (i=Istart;i<Iend;i++) {
     82:     if (i>0) PetscCall(MatSetValue(C,i,i-1,b,INSERT_VALUES));
     83:     PetscCall(MatSetValue(C,i,i,d,INSERT_VALUES));
     84:     if (i<n-1) PetscCall(MatSetValue(C,i,i+1,b,INSERT_VALUES));
     85:   }

     87:   PetscCall(MatAssemblyBegin(R,MAT_FINAL_ASSEMBLY));
     88:   PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY));
     89:   PetscCall(MatAssemblyEnd(R,MAT_FINAL_ASSEMBLY));
     90:   PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY));

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

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

     98:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     99:                 Create the eigensolver and set various options
    100:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    102:   PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
    103:   PetscCall(EPSSetOperators(eps,H,NULL));
    104:   PetscCall(EPSSetProblemType(eps,EPS_BSE));
    105:   PetscCall(EPSSetFromOptions(eps));

    107:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    108:                  Solve the eigensystem and display solution
    109:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    111:   PetscCall(EPSSolve(eps));

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

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

    142:   PetscCall(EPSDestroy(&eps));
    143:   PetscCall(MatDestroy(&R));
    144:   PetscCall(MatDestroy(&C));
    145:   PetscCall(MatDestroy(&H));
    146:   PetscCall(SlepcFinalize());
    147:   return 0;
    148: }

    150: /*TEST

    152:    testset:
    153:       args: -eps_nev 4 -eps_ncv 16 -eps_krylovschur_bse_type {{shao gruning projectedbse}} -terse -checkorthog -nest {{0 1}}
    154:       filter: sed -e "s/17496/17495/g" | sed -e "s/38566/38567/g" | sed -e "s/32172/32173/g"
    155:       nsize: {{1 2}}
    156:       test:
    157:          suffix: 1
    158:          requires: complex
    159:       test:
    160:          suffix: 1_real
    161:          requires: !complex
    162:       test:
    163:          suffix: 1_dense
    164:          args: -mat_type dense
    165:          requires: complex
    166:          output_file: output/ex55_1.out
    167:       test:
    168:          suffix: 1_cuda
    169:          args: -mat_type aijcusparse
    170:          requires: cuda complex
    171:          output_file: output/ex55_1.out
    172:       test:
    173:          suffix: 1_real_cuda
    174:          args: -mat_type aijcusparse
    175:          requires: cuda !complex
    176:          output_file: output/ex55_1_real.out
    177:       test:
    178:          suffix: 1_hip
    179:          args: -mat_type aijhipsparse
    180:          requires: hip complex
    181:          output_file: output/ex55_1.out
    182:       test:
    183:          suffix: 1_real_hip
    184:          args: -mat_type aijhipsparse
    185:          requires: hip !complex
    186:          output_file: output/ex55_1_real.out

    188:    testset:
    189:       args: -eps_nev 4 -eps_ncv 16 -eps_krylovschur_bse_type {{shao gruning projectedbse}} -st_type sinvert -terse
    190:       filter: sed -e "s/17496/17495/g" | sed -e "s/38566/38567/g"
    191:       test:
    192:          suffix: 1_sinvert
    193:          requires: complex
    194:       test:
    195:          nsize: 4
    196:          args: -mat_type scalapack
    197:          suffix: 1_sinvert_scalapack
    198:          requires: complex scalapack
    199:          output_file: output/ex55_1_sinvert.out
    200:       test:
    201:          suffix: 1_real_sinvert
    202:          requires: !complex
    203:       test:
    204:          nsize: 4
    205:          args: -mat_type scalapack
    206:          suffix: 1_real_sinvert_scalapack
    207:          requires: !complex scalapack
    208:          output_file: output/ex55_1_real_sinvert.out

    210:    testset:
    211:       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
    212:       test:
    213:          suffix: 2
    214:          requires: double complex
    215:       test:
    216:          suffix: 2_real
    217:          requires: double !complex

    219:    testset:
    220:       args: -eps_nev 28 -eps_ncv 18 -terse
    221:       test:
    222:          suffix: 3
    223:          requires: !complex !single

    225: TEST*/