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*/