Actual source code: test44.c
slepc-3.22.1 2024-10-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 using shell matrices.\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: /*
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);
40: /*
41: User context for shell matrices
42: */
43: typedef struct {
44: PetscScalar a,b,c,d;
45: } CTX_SHELL;
47: int main(int argc,char **argv)
48: {
49: Mat H,R,C; /* problem matrices */
50: EPS eps; /* eigenproblem solver context */
51: PetscReal lev;
52: PetscInt n=24,i,nconv;
53: PetscMPIInt size;
54: PetscBool terse,checkorthog;
55: Vec t,*x,*y;
56: CTX_SHELL *ctx;
58: PetscFunctionBeginUser;
59: PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
60: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size));
61: PetscCheck(size==1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only");
63: PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
64: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nShell Bethe-Salpeter eigenproblem, n=%" PetscInt_FMT "\n\n",n));
66: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
67: Generate the shell problem matrices R and C
68: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
70: 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: ctx->a = -0.1;
77: ctx->b = 1.0;
78: ctx->d = 2.0;
79: #endif
80: ctx->c = 4.5;
82: PetscCall(MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,(void*)ctx,&R));
83: PetscCall(MatShellSetOperation(R,MATOP_MULT,(void(*)(void))MatMult_R));
84: PetscCall(MatShellSetOperation(R,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_R));
85: PetscCall(MatSetOption(R,MAT_HERMITIAN,PETSC_TRUE));
86: PetscCall(MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,(void*)ctx,&C));
87: PetscCall(MatShellSetOperation(C,MATOP_MULT,(void(*)(void))MatMult_C));
88: PetscCall(MatShellSetOperation(C,MATOP_MULT_HERMITIAN_TRANSPOSE,(void(*)(void))MatMultHermitianTranspose_C));
89: PetscCall(MatSetOption(C,MAT_SYMMETRIC,PETSC_TRUE));
91: PetscCall(MatCreateBSE(R,C,&H));
92: PetscCall(MatDestroy(&R));
93: PetscCall(MatDestroy(&C));
95: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
96: Create the eigensolver and set various options
97: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
99: PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
100: PetscCall(EPSSetOperators(eps,H,NULL));
101: PetscCall(EPSSetProblemType(eps,EPS_BSE));
102: PetscCall(EPSSetWhichEigenpairs(eps,EPS_SMALLEST_MAGNITUDE));
103: PetscCall(EPSSetFromOptions(eps));
105: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106: Solve the eigensystem and display solution
107: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109: PetscCall(EPSSolve(eps));
111: /* show detailed info unless -terse option is given by user */
112: PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
113: if (terse) PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
114: else {
115: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
116: PetscCall(EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD));
117: PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
118: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
119: }
121: /* check bi-orthogonality */
122: PetscCall(PetscOptionsHasName(NULL,NULL,"-checkorthog",&checkorthog));
123: PetscCall(EPSGetConverged(eps,&nconv));
124: if (checkorthog && nconv>0) {
125: PetscCall(MatCreateVecs(H,&t,NULL));
126: PetscCall(VecDuplicateVecs(t,nconv,&x));
127: PetscCall(VecDuplicateVecs(t,nconv,&y));
128: for (i=0;i<nconv;i++) {
129: PetscCall(EPSGetEigenvector(eps,i,x[i],NULL));
130: PetscCall(EPSGetLeftEigenvector(eps,i,y[i],NULL));
131: }
132: PetscCall(VecCheckOrthogonality(x,nconv,y,nconv,NULL,NULL,&lev));
133: if (lev<100*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD," Level of bi-orthogonality of eigenvectors < 100*eps\n\n"));
134: else PetscCall(PetscPrintf(PETSC_COMM_WORLD," Level of bi-orthogonality of eigenvectors: %g\n\n",(double)lev));
135: PetscCall(VecDestroy(&t));
136: PetscCall(VecDestroyVecs(nconv,&x));
137: PetscCall(VecDestroyVecs(nconv,&y));
138: }
140: PetscCall(EPSDestroy(&eps));
141: PetscCall(MatDestroy(&H));
142: PetscCall(PetscFree(ctx));
143: PetscCall(SlepcFinalize());
144: return 0;
145: }
147: /*
148: Matrix-vector y = R*x.
150: R = pentadiag{a,b,c,conj(b),conj(a)}
151: */
152: PetscErrorCode MatMult_R(Mat R,Vec x,Vec y)
153: {
154: CTX_SHELL *ctx;
155: PetscInt n,i;
156: const PetscScalar *px;
157: PetscScalar *py;
159: PetscFunctionBeginUser;
160: PetscCall(MatShellGetContext(R,&ctx));
161: PetscCall(MatGetSize(R,NULL,&n));
162: PetscCall(VecGetArrayRead(x,&px));
163: PetscCall(VecGetArray(y,&py));
164: for (i=0;i<n;i++) {
165: py[i] = ctx->c*px[i];
166: if (i>1) py[i] += ctx->a*px[i-2];
167: if (i>0) py[i] += ctx->b*px[i-1];
168: if (i<n-1) py[i] += PetscConj(ctx->b)*px[i+1];
169: if (i<n-2) py[i] += PetscConj(ctx->a)*px[i+2];
170: }
171: PetscCall(VecRestoreArrayRead(x,&px));
172: PetscCall(VecRestoreArray(y,&py));
173: PetscFunctionReturn(PETSC_SUCCESS);
174: }
176: /*
177: Matrix-vector y = R^T*x.
179: Only needed to compute the residuals.
180: */
181: PetscErrorCode MatMultTranspose_R(Mat R,Vec x,Vec y)
182: {
183: Vec w;
185: PetscFunctionBeginUser;
186: PetscCall(VecDuplicate(x,&w));
187: PetscCall(VecCopy(x,w));
188: PetscCall(VecConjugate(w));
189: PetscCall(MatMult_R(R,w,y));
190: PetscCall(VecConjugate(y));
191: PetscCall(VecDestroy(&w));
192: PetscFunctionReturn(PETSC_SUCCESS);
193: }
195: /*
196: Matrix-vector y = C*x.
198: C = tridiag{b,d,b}
199: */
200: PetscErrorCode MatMult_C(Mat C,Vec x,Vec y)
201: {
202: CTX_SHELL *ctx;
203: PetscInt n,i;
204: const PetscScalar *px;
205: PetscScalar *py;
207: PetscFunctionBeginUser;
208: PetscCall(MatShellGetContext(C,&ctx));
209: PetscCall(MatGetSize(C,NULL,&n));
210: PetscCall(VecGetArrayRead(x,&px));
211: PetscCall(VecGetArray(y,&py));
212: for (i=0;i<n;i++) {
213: py[i] = ctx->d*px[i];
214: if (i>0) py[i] += ctx->b*px[i-1];
215: if (i<n-1) py[i] += ctx->b*px[i+1];
216: }
217: PetscCall(VecRestoreArrayRead(x,&px));
218: PetscCall(VecRestoreArray(y,&py));
219: PetscFunctionReturn(PETSC_SUCCESS);
220: }
222: /*
223: Matrix-vector y = C^H*x.
225: Only needed to compute the residuals.
226: */
227: PetscErrorCode MatMultHermitianTranspose_C(Mat C,Vec x,Vec y)
228: {
229: Vec w;
231: PetscFunctionBeginUser;
232: PetscCall(VecDuplicate(x,&w));
233: PetscCall(VecCopy(x,w));
234: PetscCall(VecConjugate(w));
235: PetscCall(MatMult_C(C,w,y));
236: PetscCall(VecConjugate(y));
237: PetscCall(VecDestroy(&w));
238: PetscFunctionReturn(PETSC_SUCCESS);
239: }
241: /*TEST
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
253: TEST*/