Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
6 : This file is part of SLEPc.
7 : SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9 : */
10 :
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";
14 :
15 : #include <slepceps.h>
16 :
17 : /*
18 : This example computes eigenvalues of a matrix
19 :
20 : H = [ R C
21 : -C^H -R^T ],
22 :
23 : where R is Hermitian and C is complex symmetric. In particular, R and C have the
24 : following Toeplitz structure:
25 :
26 : R = pentadiag{a,b,c,conj(b),conj(a)}
27 : C = tridiag{b,d,b}
28 :
29 : where a,b,d are complex scalars, and c is real.
30 : */
31 :
32 18 : int main(int argc,char **argv)
33 : {
34 18 : Mat H,R,C; /* problem matrices */
35 18 : EPS eps; /* eigenproblem solver context */
36 18 : PetscScalar a,b,c,d;
37 18 : PetscReal lev;
38 18 : PetscInt n=24,Istart,Iend,i,nconv;
39 18 : PetscBool terse,checkorthog;
40 18 : Vec t,*x,*y;
41 :
42 18 : PetscFunctionBeginUser;
43 18 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
44 :
45 18 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
46 18 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nBethe-Salpeter eigenproblem, n=%" PetscInt_FMT "\n\n",n));
47 :
48 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
49 : Compute the problem matrices R and C
50 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
51 :
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 18 : a = -0.1;
58 18 : b = 1.0;
59 18 : d = 2.0;
60 : #endif
61 18 : c = 4.5;
62 :
63 18 : PetscCall(MatCreate(PETSC_COMM_WORLD,&R));
64 18 : PetscCall(MatSetSizes(R,PETSC_DECIDE,PETSC_DECIDE,n,n));
65 18 : PetscCall(MatSetFromOptions(R));
66 :
67 18 : PetscCall(MatCreate(PETSC_COMM_WORLD,&C));
68 18 : PetscCall(MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,n,n));
69 18 : PetscCall(MatSetFromOptions(C));
70 :
71 18 : PetscCall(MatGetOwnershipRange(R,&Istart,&Iend));
72 774 : for (i=Istart;i<Iend;i++) {
73 756 : if (i>1) PetscCall(MatSetValue(R,i,i-2,a,INSERT_VALUES));
74 756 : if (i>0) PetscCall(MatSetValue(R,i,i-1,b,INSERT_VALUES));
75 756 : PetscCall(MatSetValue(R,i,i,c,INSERT_VALUES));
76 756 : if (i<n-1) PetscCall(MatSetValue(R,i,i+1,PetscConj(b),INSERT_VALUES));
77 756 : if (i<n-2) PetscCall(MatSetValue(R,i,i+2,PetscConj(a),INSERT_VALUES));
78 : }
79 :
80 18 : PetscCall(MatGetOwnershipRange(C,&Istart,&Iend));
81 774 : for (i=Istart;i<Iend;i++) {
82 756 : if (i>0) PetscCall(MatSetValue(C,i,i-1,b,INSERT_VALUES));
83 756 : PetscCall(MatSetValue(C,i,i,d,INSERT_VALUES));
84 756 : if (i<n-1) PetscCall(MatSetValue(C,i,i+1,b,INSERT_VALUES));
85 : }
86 :
87 18 : PetscCall(MatAssemblyBegin(R,MAT_FINAL_ASSEMBLY));
88 18 : PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY));
89 18 : PetscCall(MatAssemblyEnd(R,MAT_FINAL_ASSEMBLY));
90 18 : PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY));
91 :
92 18 : PetscCall(MatCreateBSE(R,C,&H));
93 :
94 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
95 : Create the eigensolver and set various options
96 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
97 :
98 18 : PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
99 18 : PetscCall(EPSSetOperators(eps,H,NULL));
100 18 : PetscCall(EPSSetProblemType(eps,EPS_BSE));
101 18 : PetscCall(EPSSetFromOptions(eps));
102 :
103 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104 : Solve the eigensystem and display solution
105 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
106 :
107 18 : PetscCall(EPSSolve(eps));
108 :
109 : /* show detailed info unless -terse option is given by user */
110 18 : PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
111 18 : if (terse) PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
112 : else {
113 0 : PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
114 0 : PetscCall(EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD));
115 0 : PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
116 0 : PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
117 : }
118 :
119 : /* check bi-orthogonality */
120 18 : PetscCall(PetscOptionsHasName(NULL,NULL,"-checkorthog",&checkorthog));
121 18 : PetscCall(EPSGetConverged(eps,&nconv));
122 18 : if (checkorthog && nconv>0) {
123 15 : PetscCall(MatCreateVecs(H,&t,NULL));
124 15 : PetscCall(VecDuplicateVecs(t,nconv,&x));
125 15 : PetscCall(VecDuplicateVecs(t,nconv,&y));
126 289 : for (i=0;i<nconv;i++) {
127 274 : PetscCall(EPSGetEigenvector(eps,i,x[i],NULL));
128 274 : PetscCall(EPSGetLeftEigenvector(eps,i,y[i],NULL));
129 : }
130 15 : PetscCall(VecCheckOrthogonality(x,nconv,y,nconv,NULL,NULL,&lev));
131 15 : if (lev<100*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD," Level of bi-orthogonality of eigenvectors < 100*eps\n\n"));
132 0 : else PetscCall(PetscPrintf(PETSC_COMM_WORLD," Level of bi-orthogonality of eigenvectors: %g\n\n",(double)lev));
133 15 : PetscCall(VecDestroy(&t));
134 15 : PetscCall(VecDestroyVecs(nconv,&x));
135 15 : PetscCall(VecDestroyVecs(nconv,&y));
136 : }
137 :
138 18 : PetscCall(EPSDestroy(&eps));
139 18 : PetscCall(MatDestroy(&R));
140 18 : PetscCall(MatDestroy(&C));
141 18 : PetscCall(MatDestroy(&H));
142 18 : PetscCall(SlepcFinalize());
143 : return 0;
144 : }
145 :
146 : /*TEST
147 :
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
158 :
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
180 :
181 : testset:
182 : 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
183 : test:
184 : suffix: 2
185 : requires: double complex
186 : test:
187 : suffix: 2_real
188 : requires: double !complex
189 :
190 : TEST*/
|
