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