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