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 using shell matrices.\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 : /*
33 : User-defined routines
34 : */
35 : PetscErrorCode MatMult_R(Mat R,Vec x,Vec y);
36 : PetscErrorCode MatMultTranspose_R(Mat R,Vec x,Vec y);
37 : PetscErrorCode MatMult_C(Mat C,Vec x,Vec y);
38 : PetscErrorCode MatMultHermitianTranspose_C(Mat C,Vec x,Vec y);
39 :
40 : /*
41 : User context for shell matrices
42 : */
43 : typedef struct {
44 : PetscScalar a,b,c,d;
45 : } CTX_SHELL;
46 :
47 2 : int main(int argc,char **argv)
48 : {
49 2 : Mat H,R,C; /* problem matrices */
50 2 : EPS eps; /* eigenproblem solver context */
51 2 : PetscReal lev;
52 2 : PetscInt n=24,i,nconv;
53 2 : PetscMPIInt size;
54 2 : PetscBool terse,checkorthog;
55 2 : Vec t,*x,*y;
56 2 : CTX_SHELL *ctx;
57 :
58 2 : PetscFunctionBeginUser;
59 2 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
60 2 : PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size));
61 2 : PetscCheck(size==1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only");
62 :
63 2 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
64 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nShell Bethe-Salpeter eigenproblem, n=%" PetscInt_FMT "\n\n",n));
65 :
66 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
67 : Generate the shell problem matrices R and C
68 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
69 :
70 2 : PetscCall(PetscNew(&ctx));
71 : #if defined(PETSC_USE_COMPLEX)
72 : ctx->a = PetscCMPLX(-0.1,0.2);
73 : ctx->b = PetscCMPLX(1.0,0.5);
74 : ctx->d = PetscCMPLX(2.0,0.2);
75 : #else
76 2 : ctx->a = -0.1;
77 2 : ctx->b = 1.0;
78 2 : ctx->d = 2.0;
79 : #endif
80 2 : ctx->c = 4.5;
81 :
82 2 : PetscCall(MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,(void*)ctx,&R));
83 2 : PetscCall(MatShellSetOperation(R,MATOP_MULT,(void(*)(void))MatMult_R));
84 2 : PetscCall(MatShellSetOperation(R,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_R));
85 2 : PetscCall(MatSetOption(R,MAT_HERMITIAN,PETSC_TRUE));
86 2 : PetscCall(MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,(void*)ctx,&C));
87 2 : PetscCall(MatShellSetOperation(C,MATOP_MULT,(void(*)(void))MatMult_C));
88 2 : PetscCall(MatShellSetOperation(C,MATOP_MULT_HERMITIAN_TRANSPOSE,(void(*)(void))MatMultHermitianTranspose_C));
89 2 : PetscCall(MatSetOption(C,MAT_SYMMETRIC,PETSC_TRUE));
90 :
91 2 : PetscCall(MatCreateBSE(R,C,&H));
92 2 : PetscCall(MatDestroy(&R));
93 2 : PetscCall(MatDestroy(&C));
94 :
95 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
96 : Create the eigensolver and set various options
97 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
98 :
99 2 : PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
100 2 : PetscCall(EPSSetOperators(eps,H,NULL));
101 2 : PetscCall(EPSSetProblemType(eps,EPS_BSE));
102 2 : PetscCall(EPSSetWhichEigenpairs(eps,EPS_SMALLEST_MAGNITUDE));
103 2 : PetscCall(EPSSetFromOptions(eps));
104 :
105 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106 : Solve the eigensystem and display solution
107 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
108 :
109 2 : PetscCall(EPSSolve(eps));
110 :
111 : /* show detailed info unless -terse option is given by user */
112 2 : PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
113 2 : if (terse) PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
114 : else {
115 0 : PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
116 0 : PetscCall(EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD));
117 0 : PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
118 0 : PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
119 : }
120 :
121 : /* check bi-orthogonality */
122 2 : PetscCall(PetscOptionsHasName(NULL,NULL,"-checkorthog",&checkorthog));
123 2 : PetscCall(EPSGetConverged(eps,&nconv));
124 2 : if (checkorthog && nconv>0) {
125 2 : PetscCall(MatCreateVecs(H,&t,NULL));
126 2 : PetscCall(VecDuplicateVecs(t,nconv,&x));
127 2 : PetscCall(VecDuplicateVecs(t,nconv,&y));
128 30 : for (i=0;i<nconv;i++) {
129 28 : PetscCall(EPSGetEigenvector(eps,i,x[i],NULL));
130 28 : PetscCall(EPSGetLeftEigenvector(eps,i,y[i],NULL));
131 : }
132 2 : PetscCall(VecCheckOrthogonality(x,nconv,y,nconv,NULL,NULL,&lev));
133 2 : if (lev<100*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD," Level of bi-orthogonality of eigenvectors < 100*eps\n\n"));
134 0 : else PetscCall(PetscPrintf(PETSC_COMM_WORLD," Level of bi-orthogonality of eigenvectors: %g\n\n",(double)lev));
135 2 : PetscCall(VecDestroy(&t));
136 2 : PetscCall(VecDestroyVecs(nconv,&x));
137 2 : PetscCall(VecDestroyVecs(nconv,&y));
138 : }
139 :
140 2 : PetscCall(EPSDestroy(&eps));
141 2 : PetscCall(MatDestroy(&H));
142 2 : PetscCall(PetscFree(ctx));
143 2 : PetscCall(SlepcFinalize());
144 : return 0;
145 : }
146 :
147 : /*
148 : Matrix-vector y = R*x.
149 :
150 : R = pentadiag{a,b,c,conj(b),conj(a)}
151 : */
152 148 : PetscErrorCode MatMult_R(Mat R,Vec x,Vec y)
153 : {
154 148 : CTX_SHELL *ctx;
155 148 : PetscInt n,i;
156 148 : const PetscScalar *px;
157 148 : PetscScalar *py;
158 :
159 148 : PetscFunctionBeginUser;
160 148 : PetscCall(MatShellGetContext(R,&ctx));
161 148 : PetscCall(MatGetSize(R,NULL,&n));
162 148 : PetscCall(VecGetArrayRead(x,&px));
163 148 : PetscCall(VecGetArray(y,&py));
164 3700 : for (i=0;i<n;i++) {
165 3552 : py[i] = ctx->c*px[i];
166 3552 : if (i>1) py[i] += ctx->a*px[i-2];
167 3552 : if (i>0) py[i] += ctx->b*px[i-1];
168 3552 : if (i<n-1) py[i] += PetscConj(ctx->b)*px[i+1];
169 3552 : if (i<n-2) py[i] += PetscConj(ctx->a)*px[i+2];
170 : }
171 148 : PetscCall(VecRestoreArrayRead(x,&px));
172 148 : PetscCall(VecRestoreArray(y,&py));
173 148 : PetscFunctionReturn(PETSC_SUCCESS);
174 : }
175 :
176 : /*
177 : Matrix-vector y = R^T*x.
178 :
179 : Only needed to compute the residuals.
180 : */
181 8 : PetscErrorCode MatMultTranspose_R(Mat R,Vec x,Vec y)
182 : {
183 8 : Vec w;
184 :
185 8 : PetscFunctionBeginUser;
186 8 : PetscCall(VecDuplicate(x,&w));
187 8 : PetscCall(VecCopy(x,w));
188 8 : PetscCall(VecConjugate(w));
189 8 : PetscCall(MatMult_R(R,w,y));
190 8 : PetscCall(VecConjugate(y));
191 8 : PetscCall(VecDestroy(&w));
192 8 : PetscFunctionReturn(PETSC_SUCCESS);
193 : }
194 :
195 : /*
196 : Matrix-vector y = C*x.
197 :
198 : C = tridiag{b,d,b}
199 : */
200 148 : PetscErrorCode MatMult_C(Mat C,Vec x,Vec y)
201 : {
202 148 : CTX_SHELL *ctx;
203 148 : PetscInt n,i;
204 148 : const PetscScalar *px;
205 148 : PetscScalar *py;
206 :
207 148 : PetscFunctionBeginUser;
208 148 : PetscCall(MatShellGetContext(C,&ctx));
209 148 : PetscCall(MatGetSize(C,NULL,&n));
210 148 : PetscCall(VecGetArrayRead(x,&px));
211 148 : PetscCall(VecGetArray(y,&py));
212 3700 : for (i=0;i<n;i++) {
213 3552 : py[i] = ctx->d*px[i];
214 3552 : if (i>0) py[i] += ctx->b*px[i-1];
215 3552 : if (i<n-1) py[i] += ctx->b*px[i+1];
216 : }
217 148 : PetscCall(VecRestoreArrayRead(x,&px));
218 148 : PetscCall(VecRestoreArray(y,&py));
219 148 : PetscFunctionReturn(PETSC_SUCCESS);
220 : }
221 :
222 : /*
223 : Matrix-vector y = C^H*x.
224 :
225 : Only needed to compute the residuals.
226 : */
227 0 : PetscErrorCode MatMultHermitianTranspose_C(Mat C,Vec x,Vec y)
228 : {
229 0 : Vec w;
230 :
231 0 : PetscFunctionBeginUser;
232 0 : PetscCall(VecDuplicate(x,&w));
233 0 : PetscCall(VecCopy(x,w));
234 0 : PetscCall(VecConjugate(w));
235 0 : PetscCall(MatMult_C(C,w,y));
236 0 : PetscCall(VecConjugate(y));
237 0 : PetscCall(VecDestroy(&w));
238 0 : PetscFunctionReturn(PETSC_SUCCESS);
239 : }
240 :
241 : /*TEST
242 :
243 : testset:
244 : args: -eps_nev 4 -eps_ncv 16 -eps_krylovschur_bse_type {{shao gruning}} -terse -checkorthog
245 : filter: sed -e "s/17496/17495/g" | sed -e "s/32172/32173/g" | sed -e "s/38566/38567/g"
246 : test:
247 : suffix: 1
248 : requires: complex
249 : test:
250 : suffix: 1_real
251 : requires: !complex
252 :
253 : TEST*/
|
