Actual source code: ex55.c
slepc-3.22.0 2024-09-28
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*/