Actual source code: ex55.c

slepc-main 2024-11-15
Report Typos and Errors
  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;
 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:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 95:                 Create the eigensolver and set various options
 96:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 98:   PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
 99:   PetscCall(EPSSetOperators(eps,H,NULL));
100:   PetscCall(EPSSetProblemType(eps,EPS_BSE));
101:   PetscCall(EPSSetFromOptions(eps));

103:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104:                  Solve the eigensystem and display solution
105:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

107:   PetscCall(EPSSolve(eps));

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

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

138:   PetscCall(EPSDestroy(&eps));
139:   PetscCall(MatDestroy(&R));
140:   PetscCall(MatDestroy(&C));
141:   PetscCall(MatDestroy(&H));
142:   PetscCall(SlepcFinalize());
143:   return 0;
144: }

146: /*TEST

148:    testset:
149:       args: -eps_nev 4 -eps_ncv 16 -eps_krylovschur_bse_type {{shao gruning projectedbse}} -terse -checkorthog
150:       filter: sed -e "s/17496/17495/g" | sed -e "s/38566/38567/g"
151:       nsize: {{1 2}}
152:       test:
153:          suffix: 1
154:          requires: complex
155:       test:
156:          suffix: 1_real
157:          requires: !complex

159:    testset:
160:       args: -eps_nev 4 -eps_ncv 16 -eps_krylovschur_bse_type {{shao gruning projectedbse}} -st_type sinvert -terse
161:       filter: sed -e "s/17496/17495/g" | sed -e "s/38566/38567/g"
162:       test:
163:          suffix: 1_sinvert
164:          requires: complex
165:       test:
166:          nsize: 4
167:          args: -mat_type scalapack
168:          suffix: 1_sinvert_scalapack
169:          requires: complex scalapack
170:          output_file: output/ex55_1_sinvert.out
171:       test:
172:          suffix: 1_real_sinvert
173:          requires: !complex
174:       test:
175:          nsize: 4
176:          args: -mat_type scalapack
177:          suffix: 1_real_sinvert_scalapack
178:          requires: !complex scalapack
179:          output_file: output/ex55_1_real_sinvert.out

181: TEST*/