Actual source code: test29.c
slepc-3.22.2 2024-12-02
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*);
22: int main(int argc,char **argv)
23: {
24: Mat A,B;
25: EPS eps;
26: EPSType type;
27: PetscInt i,nconv;
28: PetscBool twosided,flg;
29: PetscReal nrmr,nrml=0.0,re,im,lev;
30: PetscScalar *kr,*ki;
31: Vec t,*xr,*xi,*yr,*yi,*z;
32: char filename[PETSC_MAX_PATH_LEN];
33: PetscViewer viewer;
35: PetscFunctionBeginUser;
36: PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
38: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
39: Load the matrices that define the eigensystem, Ax=kBx
40: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
42: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nGeneralized eigenproblem stored in file.\n\n"));
43: PetscCall(PetscOptionsGetString(NULL,NULL,"-f1",filename,sizeof(filename),&flg));
44: PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must indicate a file name for matrix A with the -f1 option");
46: #if defined(PETSC_USE_COMPLEX)
47: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Reading COMPLEX matrices from binary files...\n"));
48: #else
49: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Reading REAL matrices from binary files...\n"));
50: #endif
51: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
52: PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
53: PetscCall(MatSetFromOptions(A));
54: PetscCall(MatLoad(A,viewer));
55: PetscCall(PetscViewerDestroy(&viewer));
57: PetscCall(PetscOptionsGetString(NULL,NULL,"-f2",filename,sizeof(filename),&flg));
58: if (flg) {
59: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
60: PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
61: PetscCall(MatSetFromOptions(B));
62: PetscCall(MatLoad(B,viewer));
63: PetscCall(PetscViewerDestroy(&viewer));
64: } else {
65: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Matrix B was not provided, setting B=I\n\n"));
66: B = NULL;
67: }
68: PetscCall(MatCreateVecs(A,NULL,&t));
70: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
71: Create the eigensolver and set various options
72: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
74: PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
75: PetscCall(EPSSetOperators(eps,A,B));
77: /* use a two-sided algorithm to compute left eigenvectors as well */
78: PetscCall(EPSSetTwoSided(eps,PETSC_TRUE));
80: /* allow user to change settings at run time */
81: PetscCall(EPSSetFromOptions(eps));
82: PetscCall(EPSGetTwoSided(eps,&twosided));
84: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
85: Solve the eigensystem
86: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
88: PetscCall(EPSSolve(eps));
90: /*
91: Optional: Get some information from the solver and display it
92: */
93: PetscCall(EPSGetType(eps,&type));
94: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type));
96: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97: Display solution and clean up
98: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
100: /*
101: Get number of converged approximate eigenpairs
102: */
103: PetscCall(EPSGetConverged(eps,&nconv));
104: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of converged eigenpairs: %" PetscInt_FMT "\n\n",nconv));
105: PetscCall(PetscMalloc2(nconv,&kr,nconv,&ki));
106: PetscCall(VecDuplicateVecs(t,3,&z));
107: PetscCall(VecDuplicateVecs(t,nconv,&xr));
108: PetscCall(VecDuplicateVecs(t,nconv,&xi));
109: if (twosided) {
110: PetscCall(VecDuplicateVecs(t,nconv,&yr));
111: PetscCall(VecDuplicateVecs(t,nconv,&yi));
112: }
114: if (nconv>0) {
115: /*
116: Display eigenvalues and relative errors
117: */
118: PetscCall(PetscPrintf(PETSC_COMM_WORLD,
119: " k ||Ax-kBx|| ||y'A-y'Bk||\n"
120: " ---------------- ------------------ ------------------\n"));
122: for (i=0;i<nconv;i++) {
123: /*
124: Get converged eigenpairs: i-th eigenvalue is stored in kr (real part) and
125: ki (imaginary part)
126: */
127: PetscCall(EPSGetEigenpair(eps,i,&kr[i],&ki[i],xr[i],xi[i]));
128: if (twosided) PetscCall(EPSGetLeftEigenvector(eps,i,yr[i],yi[i]));
129: /*
130: Compute the residual norms associated to each eigenpair
131: */
132: PetscCall(ComputeResidualNorm(A,B,PETSC_FALSE,kr[i],ki[i],xr[i],xi[i],z,&nrmr));
133: if (twosided) PetscCall(ComputeResidualNorm(A,B,PETSC_TRUE,kr[i],ki[i],yr[i],yi[i],z,&nrml));
135: #if defined(PETSC_USE_COMPLEX)
136: re = PetscRealPart(kr[i]);
137: im = PetscImaginaryPart(kr[i]);
138: #else
139: re = kr[i];
140: im = ki[i];
141: #endif
142: if (im!=0.0) PetscCall(PetscPrintf(PETSC_COMM_WORLD," %8f%+8fi %12g %12g\n",(double)re,(double)im,(double)nrmr,(double)nrml));
143: else PetscCall(PetscPrintf(PETSC_COMM_WORLD," %12f %12g %12g\n",(double)re,(double)nrmr,(double)nrml));
144: }
145: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n"));
146: /*
147: Check bi-orthogonality of eigenvectors
148: */
149: if (twosided) {
150: PetscCall(VecCheckOrthogonality(xr,nconv,yr,nconv,B,NULL,&lev));
151: if (lev<100*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD," Level of bi-orthogonality of eigenvectors < 100*eps\n\n"));
152: else PetscCall(PetscPrintf(PETSC_COMM_WORLD," Level of bi-orthogonality of eigenvectors: %g\n\n",(double)lev));
153: }
154: }
156: PetscCall(EPSDestroy(&eps));
157: PetscCall(MatDestroy(&A));
158: PetscCall(MatDestroy(&B));
159: PetscCall(VecDestroy(&t));
160: PetscCall(PetscFree2(kr,ki));
161: PetscCall(VecDestroyVecs(3,&z));
162: PetscCall(VecDestroyVecs(nconv,&xr));
163: PetscCall(VecDestroyVecs(nconv,&xi));
164: if (twosided) {
165: PetscCall(VecDestroyVecs(nconv,&yr));
166: PetscCall(VecDestroyVecs(nconv,&yi));
167: }
168: PetscCall(SlepcFinalize());
169: return 0;
170: }
172: /*
173: ComputeResidualNorm - Computes the norm of the residual vector
174: associated with an eigenpair.
176: Input Parameters:
177: trans - whether A' must be used instead of A
178: kr,ki - eigenvalue
179: xr,xi - eigenvector
180: z - three work vectors (the second one not referenced in complex scalars)
181: */
182: PetscErrorCode ComputeResidualNorm(Mat A,Mat B,PetscBool trans,PetscScalar kr,PetscScalar ki,Vec xr,Vec xi,Vec *z,PetscReal *norm)
183: {
184: Vec u,w=NULL;
185: PetscScalar alpha;
186: #if !defined(PETSC_USE_COMPLEX)
187: Vec v;
188: PetscReal ni,nr;
189: #endif
190: PetscErrorCode (*matmult)(Mat,Vec,Vec) = trans? MatMultHermitianTranspose: MatMult;
192: PetscFunctionBegin;
193: u = z[0];
194: if (B) w = z[2];
196: #if !defined(PETSC_USE_COMPLEX)
197: v = z[1];
198: if (ki == 0 || PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
199: #endif
200: PetscCall((*matmult)(A,xr,u)); /* u=A*x */
201: if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
202: if (B) PetscCall((*matmult)(B,xr,w)); /* w=B*x */
203: else w = xr;
204: alpha = trans? -PetscConj(kr): -kr;
205: PetscCall(VecAXPY(u,alpha,w)); /* u=A*x-k*B*x */
206: }
207: PetscCall(VecNorm(u,NORM_2,norm));
208: #if !defined(PETSC_USE_COMPLEX)
209: } else {
210: PetscCall((*matmult)(A,xr,u)); /* u=A*xr */
211: if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
212: if (B) PetscCall((*matmult)(B,xr,v)); /* v=B*xr */
213: else PetscCall(VecCopy(xr,v));
214: PetscCall(VecAXPY(u,-kr,v)); /* u=A*xr-kr*B*xr */
215: if (B) PetscCall((*matmult)(B,xi,w)); /* w=B*xi */
216: else w = xi;
217: PetscCall(VecAXPY(u,trans?-ki:ki,w)); /* u=A*xr-kr*B*xr+ki*B*xi */
218: }
219: PetscCall(VecNorm(u,NORM_2,&nr));
220: PetscCall((*matmult)(A,xi,u)); /* u=A*xi */
221: if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
222: PetscCall(VecAXPY(u,-kr,w)); /* u=A*xi-kr*B*xi */
223: PetscCall(VecAXPY(u,trans?ki:-ki,v)); /* u=A*xi-kr*B*xi-ki*B*xr */
224: }
225: PetscCall(VecNorm(u,NORM_2,&ni));
226: *norm = SlepcAbsEigenvalue(nr,ni);
227: }
228: #endif
229: PetscFunctionReturn(PETSC_SUCCESS);
230: }
232: /*TEST
234: testset:
235: 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
236: filter: grep -v "method" | sed -e "s/[+-]0\.0*i//g" | sed -e "s/[0-9]\.[0-9]*e[+-]\([0-9]*\)/removed/g"
237: requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
238: test:
239: suffix: 1
240: test:
241: suffix: 1_rqi
242: args: -eps_type power -eps_power_shift_type rayleigh -eps_nev 2 -eps_target -2000
243: test:
244: suffix: 1_rqi_singular
245: args: -eps_type power -eps_power_shift_type rayleigh -eps_nev 1 -eps_target -195500
247: test:
248: suffix: 2
249: args: -f1 ${DATAFILESPATH}/matrices/complex/mhd1280a.petsc -f2 ${DATAFILESPATH}/matrices/complex/mhd1280b.petsc -eps_nev 6 -eps_tol 1e-11
250: 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"
251: requires: double complex datafilespath !defined(PETSC_USE_64BIT_INDICES)
252: timeoutfactor: 2
254: TEST*/