Actual source code: test29.c
slepc-3.23.3 2025-09-08
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[] = "Illustrates the computation of left eigenvectors for generalized eigenproblems.\n\n"
12: "The command line options are:\n"
13: " -f1 <filename> -f2 <filename>, PETSc binary files containing A and B\n\n";
15: #include <slepceps.h>
17: /*
18: User-defined routines
19: */
20: PetscErrorCode ComputeResidualNorm(Mat,Mat,PetscBool,PetscScalar,PetscScalar,Vec,Vec,Vec*,PetscReal*);
21: PetscErrorCode CheckBiorthogonality(Vec*,Vec*,PetscScalar*,PetscInt,Mat,PetscReal*);
23: int main(int argc,char **argv)
24: {
25: Mat A,B;
26: EPS eps;
27: EPSType type;
28: PetscInt i,nconv;
29: PetscBool twosided,flg;
30: PetscReal nrmr,nrml=0.0,re,im,lev;
31: PetscScalar *kr,*ki;
32: Vec t,*xr,*xi,*yr,*yi,*z;
33: char filename[PETSC_MAX_PATH_LEN];
34: PetscViewer viewer;
36: PetscFunctionBeginUser;
37: PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
39: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
40: Load the matrices that define the eigensystem, Ax=kBx
41: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
43: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nGeneralized eigenproblem stored in file.\n\n"));
44: PetscCall(PetscOptionsGetString(NULL,NULL,"-f1",filename,sizeof(filename),&flg));
45: PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must indicate a file name for matrix A with the -f1 option");
47: #if defined(PETSC_USE_COMPLEX)
48: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Reading COMPLEX matrices from binary files...\n"));
49: #else
50: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Reading REAL matrices from binary files...\n"));
51: #endif
52: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
53: PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
54: PetscCall(MatSetFromOptions(A));
55: PetscCall(MatLoad(A,viewer));
56: PetscCall(PetscViewerDestroy(&viewer));
58: PetscCall(PetscOptionsGetString(NULL,NULL,"-f2",filename,sizeof(filename),&flg));
59: if (flg) {
60: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
61: PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
62: PetscCall(MatSetFromOptions(B));
63: PetscCall(MatLoad(B,viewer));
64: PetscCall(PetscViewerDestroy(&viewer));
65: } else {
66: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Matrix B was not provided, setting B=I\n\n"));
67: B = NULL;
68: }
69: PetscCall(MatCreateVecs(A,NULL,&t));
71: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
72: Create the eigensolver and set various options
73: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
75: PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
76: PetscCall(EPSSetOperators(eps,A,B));
78: /* use a two-sided algorithm to compute left eigenvectors as well */
79: PetscCall(EPSSetTwoSided(eps,PETSC_TRUE));
81: /* allow user to change settings at run time */
82: PetscCall(EPSSetFromOptions(eps));
83: PetscCall(EPSGetTwoSided(eps,&twosided));
85: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
86: Solve the eigensystem
87: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
89: PetscCall(EPSSolve(eps));
91: /*
92: Optional: Get some information from the solver and display it
93: */
94: PetscCall(EPSGetType(eps,&type));
95: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type));
97: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
98: Display solution and clean up
99: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
101: /*
102: Get number of converged approximate eigenpairs
103: */
104: PetscCall(EPSGetConverged(eps,&nconv));
105: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of converged eigenpairs: %" PetscInt_FMT "\n\n",nconv));
106: PetscCall(PetscMalloc2(nconv,&kr,nconv,&ki));
107: PetscCall(VecDuplicateVecs(t,3,&z));
108: PetscCall(VecDuplicateVecs(t,nconv,&xr));
109: PetscCall(VecDuplicateVecs(t,nconv,&xi));
110: if (twosided) {
111: PetscCall(VecDuplicateVecs(t,nconv,&yr));
112: PetscCall(VecDuplicateVecs(t,nconv,&yi));
113: }
115: if (nconv>0) {
116: /*
117: Display eigenvalues and relative errors
118: */
119: PetscCall(PetscPrintf(PETSC_COMM_WORLD,
120: " k ||Ax-kBx|| ||y'A-y'Bk||\n"
121: " ---------------- ------------------ ------------------\n"));
123: for (i=0;i<nconv;i++) {
124: /*
125: Get converged eigenpairs: i-th eigenvalue is stored in kr (real part) and
126: ki (imaginary part)
127: */
128: PetscCall(EPSGetEigenpair(eps,i,&kr[i],&ki[i],xr[i],xi[i]));
129: if (twosided) PetscCall(EPSGetLeftEigenvector(eps,i,yr[i],yi[i]));
130: /*
131: Compute the residual norms associated to each eigenpair
132: */
133: PetscCall(ComputeResidualNorm(A,B,PETSC_FALSE,kr[i],ki[i],xr[i],xi[i],z,&nrmr));
134: if (twosided) PetscCall(ComputeResidualNorm(A,B,PETSC_TRUE,kr[i],ki[i],yr[i],yi[i],z,&nrml));
136: #if defined(PETSC_USE_COMPLEX)
137: re = PetscRealPart(kr[i]);
138: im = PetscImaginaryPart(kr[i]);
139: #else
140: re = kr[i];
141: im = ki[i];
142: #endif
143: if (im!=0.0) PetscCall(PetscPrintf(PETSC_COMM_WORLD," %8f%+8fi %12g %12g\n",(double)re,(double)im,(double)nrmr,(double)nrml));
144: else PetscCall(PetscPrintf(PETSC_COMM_WORLD," %12f %12g %12g\n",(double)re,(double)nrmr,(double)nrml));
145: }
146: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n"));
147: /*
148: Check bi-orthogonality of eigenvectors
149: */
150: if (twosided) {
151: PetscCall(CheckBiorthogonality(xr,yr,ki,nconv,B,&lev));
152: if (lev<100*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD," Level of bi-orthogonality of eigenvectors < 100*eps\n\n"));
153: else PetscCall(PetscPrintf(PETSC_COMM_WORLD," Level of bi-orthogonality of eigenvectors: %g\n\n",(double)lev));
154: }
155: }
157: PetscCall(EPSDestroy(&eps));
158: PetscCall(MatDestroy(&A));
159: PetscCall(MatDestroy(&B));
160: PetscCall(VecDestroy(&t));
161: PetscCall(PetscFree2(kr,ki));
162: PetscCall(VecDestroyVecs(3,&z));
163: PetscCall(VecDestroyVecs(nconv,&xr));
164: PetscCall(VecDestroyVecs(nconv,&xi));
165: if (twosided) {
166: PetscCall(VecDestroyVecs(nconv,&yr));
167: PetscCall(VecDestroyVecs(nconv,&yi));
168: }
169: PetscCall(SlepcFinalize());
170: return 0;
171: }
173: /*
174: ComputeResidualNorm - Computes the norm of the residual vector
175: associated with an eigenpair.
177: Input Parameters:
178: trans - whether A' must be used instead of A
179: kr,ki - eigenvalue
180: xr,xi - eigenvector
181: z - three work vectors (the second one not referenced in complex scalars)
182: */
183: PetscErrorCode ComputeResidualNorm(Mat A,Mat B,PetscBool trans,PetscScalar kr,PetscScalar ki,Vec xr,Vec xi,Vec *z,PetscReal *norm)
184: {
185: Vec u,w=NULL;
186: PetscScalar alpha;
187: #if !defined(PETSC_USE_COMPLEX)
188: Vec v;
189: PetscReal ni,nr;
190: #endif
191: PetscErrorCode (*matmult)(Mat,Vec,Vec) = trans? MatMultHermitianTranspose: MatMult;
193: PetscFunctionBegin;
194: u = z[0];
195: if (B) w = z[2];
197: #if !defined(PETSC_USE_COMPLEX)
198: v = z[1];
199: if (ki == 0 || PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
200: #endif
201: PetscCall((*matmult)(A,xr,u)); /* u=A*x */
202: if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
203: if (B) PetscCall((*matmult)(B,xr,w)); /* w=B*x */
204: else w = xr;
205: alpha = trans? -PetscConj(kr): -kr;
206: PetscCall(VecAXPY(u,alpha,w)); /* u=A*x-k*B*x */
207: }
208: PetscCall(VecNorm(u,NORM_2,norm));
209: #if !defined(PETSC_USE_COMPLEX)
210: } else {
211: PetscCall((*matmult)(A,xr,u)); /* u=A*xr */
212: if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
213: if (B) PetscCall((*matmult)(B,xr,v)); /* v=B*xr */
214: else PetscCall(VecCopy(xr,v));
215: PetscCall(VecAXPY(u,-kr,v)); /* u=A*xr-kr*B*xr */
216: if (B) PetscCall((*matmult)(B,xi,w)); /* w=B*xi */
217: else w = xi;
218: PetscCall(VecAXPY(u,trans?-ki:ki,w)); /* u=A*xr-kr*B*xr+ki*B*xi */
219: }
220: PetscCall(VecNorm(u,NORM_2,&nr));
221: PetscCall((*matmult)(A,xi,u)); /* u=A*xi */
222: if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
223: PetscCall(VecAXPY(u,-kr,w)); /* u=A*xi-kr*B*xi */
224: PetscCall(VecAXPY(u,trans?ki:-ki,v)); /* u=A*xi-kr*B*xi-ki*B*xr */
225: }
226: PetscCall(VecNorm(u,NORM_2,&ni));
227: *norm = SlepcAbsEigenvalue(nr,ni);
228: }
229: #endif
230: PetscFunctionReturn(PETSC_SUCCESS);
231: }
233: /*
234: CheckOrthogonality - Similar to VecCheckOrthogonality but taking into account complex eigenvalues.
236: Input Parameters:
237: V - right eigenvectors
238: W - left eigenvectors
239: eigi - imaginary part of eigenvalues (only to check if the corresponding vector is complex)
240: n - number of V,W vectors
241: B - matrix defining the inner product
242: Output Parameter:
243: lev - level of bi-orthogonality
244: */
245: PetscErrorCode CheckBiorthogonality(Vec *V,Vec *W,PetscScalar *ki,PetscInt n,Mat B,PetscReal *lev)
246: {
247: PetscInt i,j;
248: PetscScalar val1;
249: PetscBool wcmplx;
250: Vec wr;
251: #if !defined(PETSC_USE_COMPLEX)
252: Vec wi;
253: PetscScalar val2,dotr,doti;
254: #endif
256: PetscFunctionBegin;
257: if (n<=0) PetscFunctionReturn(PETSC_SUCCESS);
258: PetscCall(VecDuplicate(V[0],&wr));
259: #if !defined(PETSC_USE_COMPLEX)
260: PetscCall(VecDuplicate(V[0],&wi));
261: #endif
262: *lev = 0.0;
263: for (i=0;i<n;i++) {
264: wcmplx = PETSC_FALSE;
265: PetscCall(MatMult(B,W[i],wr));
266: #if !defined(PETSC_USE_COMPLEX)
267: if (ki[i] != 0.0) {
268: PetscCall(MatMult(B,W[i+1],wi));
269: wcmplx = PETSC_TRUE;
270: }
271: #endif
272: for (j=0;j<n;j++) {
273: if (j==i || (wcmplx && j==i+1)) continue;
274: PetscCall(VecDot(wr,V[j],&val1));
275: #if !defined(PETSC_USE_COMPLEX)
276: dotr = val1;
277: if (ki[j] != 0.0) { /* V complex */
278: if (wcmplx) {
279: PetscCall(VecDot(wi,V[j+1],&val2));
280: dotr = val1+val2;
281: PetscCall(VecDot(wi,V[j],&val1));
282: PetscCall(VecDot(wr,V[j+1],&val2));
283: doti = val1-val2;
284: } else {
285: PetscCall(VecDot(wr,V[j+1],&val1));
286: doti = -val1;
287: }
288: } else if (wcmplx) {
289: PetscCall(VecDot(wi,V[j],&doti));
290: } else doti = 0.0;
291: *lev = PetscMax(*lev,SlepcAbsEigenvalue(dotr,doti));
292: if (ki[j] != 0.0) j++;
293: #else
294: *lev = PetscMax(*lev,PetscAbsScalar(val1));
295: #endif
296: }
297: if (wcmplx) i++;
298: }
299: PetscCall(VecDestroy(&wr));
300: #if !defined(PETSC_USE_COMPLEX)
301: PetscCall(VecDestroy(&wi));
302: #endif
303: PetscFunctionReturn(PETSC_SUCCESS);
304: }
306: /*TEST
308: testset:
309: args: -f1 ${SLEPC_DIR}/share/slepc/datafiles/matrices/bfw62a.petsc -f2 ${SLEPC_DIR}/share/slepc/datafiles/matrices/bfw62b.petsc -eps_nev 4 -st_type sinvert -eps_target -190000
310: filter: grep -v "method" | sed -e "s/[+-]0\.0*i//g" | sed -e "s/[0-9]\.[0-9]*e[+-]\([0-9]*\)/removed/g"
311: requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
312: test:
313: suffix: 1
314: test:
315: suffix: 1_cmplxvals
316: args: -eps_target -250000
317: test:
318: suffix: 1_rqi
319: args: -eps_type power -eps_power_shift_type rayleigh -eps_nev 2 -eps_target -2000
320: test:
321: suffix: 1_rqi_singular
322: args: -eps_type power -eps_power_shift_type rayleigh -eps_nev 1 -eps_target -195500
324: test:
325: suffix: 2
326: args: -f1 ${DATAFILESPATH}/matrices/complex/mhd1280a.petsc -f2 ${DATAFILESPATH}/matrices/complex/mhd1280b.petsc -eps_nev 6 -eps_tol 1e-11
327: filter: sed -e "s/-892/+892/" | sed -e "s/-759/+759/" | sed -e "s/-674/+674/" | sed -e "s/[0-9]\.[0-9]*e[+-]\([0-9]*\)/removed/g"
328: requires: double complex datafilespath !defined(PETSC_USE_64BIT_INDICES)
329: timeoutfactor: 2
331: test:
332: suffix: 3
333: args: -f1 ${SLEPC_DIR}/share/slepc/datafiles/matrices/bfw62a.petsc -f2 ${SLEPC_DIR}/share/slepc/datafiles/matrices/bfw62b.petsc -eps_nev 4 -st_type sinvert -eps_target -250000
334: filter: grep -v "method" | sed -e "s/[+-]0\.0*i//g" | sed -e "s/[0-9]\.[0-9]*e[+-]\([0-9]*\)/removed/g" | grep -v bi-orthog
335: requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
337: TEST*/