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[] = "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";
14 :
15 : #include <slepceps.h>
16 :
17 : /*
18 : User-defined routines
19 : */
20 : PetscErrorCode ComputeResidualNorm(Mat,Mat,PetscBool,PetscScalar,PetscScalar,Vec,Vec,Vec*,PetscReal*);
21 :
22 1 : int main(int argc,char **argv)
23 : {
24 1 : Mat A,B;
25 1 : EPS eps;
26 1 : EPSType type;
27 1 : PetscInt i,nconv;
28 1 : PetscBool twosided,flg;
29 1 : PetscReal nrmr,nrml=0.0,re,im,lev;
30 1 : PetscScalar *kr,*ki;
31 1 : Vec t,*xr,*xi,*yr,*yi,*z;
32 1 : char filename[PETSC_MAX_PATH_LEN];
33 1 : PetscViewer viewer;
34 :
35 1 : PetscFunctionBeginUser;
36 1 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
37 :
38 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
39 : Load the matrices that define the eigensystem, Ax=kBx
40 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
41 :
42 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nGeneralized eigenproblem stored in file.\n\n"));
43 1 : PetscCall(PetscOptionsGetString(NULL,NULL,"-f1",filename,sizeof(filename),&flg));
44 1 : PetscCheck(flg,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Must indicate a file name for matrix A with the -f1 option");
45 :
46 : #if defined(PETSC_USE_COMPLEX)
47 1 : 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 1 : PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
52 1 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
53 1 : PetscCall(MatSetFromOptions(A));
54 1 : PetscCall(MatLoad(A,viewer));
55 1 : PetscCall(PetscViewerDestroy(&viewer));
56 :
57 1 : PetscCall(PetscOptionsGetString(NULL,NULL,"-f2",filename,sizeof(filename),&flg));
58 1 : if (flg) {
59 1 : PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer));
60 1 : PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
61 1 : PetscCall(MatSetFromOptions(B));
62 1 : PetscCall(MatLoad(B,viewer));
63 1 : PetscCall(PetscViewerDestroy(&viewer));
64 : } else {
65 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Matrix B was not provided, setting B=I\n\n"));
66 0 : B = NULL;
67 : }
68 1 : PetscCall(MatCreateVecs(A,NULL,&t));
69 :
70 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
71 : Create the eigensolver and set various options
72 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
73 :
74 1 : PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
75 1 : PetscCall(EPSSetOperators(eps,A,B));
76 :
77 : /* use a two-sided algorithm to compute left eigenvectors as well */
78 1 : PetscCall(EPSSetTwoSided(eps,PETSC_TRUE));
79 :
80 : /* allow user to change settings at run time */
81 1 : PetscCall(EPSSetFromOptions(eps));
82 1 : PetscCall(EPSGetTwoSided(eps,&twosided));
83 :
84 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
85 : Solve the eigensystem
86 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
87 :
88 1 : PetscCall(EPSSolve(eps));
89 :
90 : /*
91 : Optional: Get some information from the solver and display it
92 : */
93 1 : PetscCall(EPSGetType(eps,&type));
94 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type));
95 :
96 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97 : Display solution and clean up
98 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
99 :
100 : /*
101 : Get number of converged approximate eigenpairs
102 : */
103 1 : PetscCall(EPSGetConverged(eps,&nconv));
104 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of converged eigenpairs: %" PetscInt_FMT "\n\n",nconv));
105 1 : PetscCall(PetscMalloc2(nconv,&kr,nconv,&ki));
106 1 : PetscCall(VecDuplicateVecs(t,3,&z));
107 1 : PetscCall(VecDuplicateVecs(t,nconv,&xr));
108 1 : PetscCall(VecDuplicateVecs(t,nconv,&xi));
109 1 : if (twosided) {
110 1 : PetscCall(VecDuplicateVecs(t,nconv,&yr));
111 1 : PetscCall(VecDuplicateVecs(t,nconv,&yi));
112 : }
113 :
114 1 : if (nconv>0) {
115 : /*
116 : Display eigenvalues and relative errors
117 : */
118 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,
119 : " k ||Ax-kBx|| ||y'A-y'Bk||\n"
120 : " ---------------- ------------------ ------------------\n"));
121 :
122 7 : 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 6 : PetscCall(EPSGetEigenpair(eps,i,&kr[i],&ki[i],xr[i],xi[i]));
128 6 : if (twosided) PetscCall(EPSGetLeftEigenvector(eps,i,yr[i],yi[i]));
129 : /*
130 : Compute the residual norms associated to each eigenpair
131 : */
132 6 : PetscCall(ComputeResidualNorm(A,B,PETSC_FALSE,kr[i],ki[i],xr[i],xi[i],z,&nrmr));
133 6 : if (twosided) PetscCall(ComputeResidualNorm(A,B,PETSC_TRUE,kr[i],ki[i],yr[i],yi[i],z,&nrml));
134 :
135 : #if defined(PETSC_USE_COMPLEX)
136 6 : re = PetscRealPart(kr[i]);
137 6 : im = PetscImaginaryPart(kr[i]);
138 : #else
139 : re = kr[i];
140 : im = ki[i];
141 : #endif
142 6 : if (im!=0.0) PetscCall(PetscPrintf(PETSC_COMM_WORLD," %8f%+8fi %12g %12g\n",(double)re,(double)im,(double)nrmr,(double)nrml));
143 6 : else PetscCall(PetscPrintf(PETSC_COMM_WORLD," %12f %12g %12g\n",(double)re,(double)nrmr,(double)nrml));
144 : }
145 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n"));
146 : /*
147 : Check bi-orthogonality of eigenvectors
148 : */
149 1 : if (twosided) {
150 1 : PetscCall(VecCheckOrthogonality(xr,nconv,yr,nconv,B,NULL,&lev));
151 1 : if (lev<100*PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD," Level of bi-orthogonality of eigenvectors < 100*eps\n\n"));
152 0 : else PetscCall(PetscPrintf(PETSC_COMM_WORLD," Level of bi-orthogonality of eigenvectors: %g\n\n",(double)lev));
153 : }
154 : }
155 :
156 1 : PetscCall(EPSDestroy(&eps));
157 1 : PetscCall(MatDestroy(&A));
158 1 : PetscCall(MatDestroy(&B));
159 1 : PetscCall(VecDestroy(&t));
160 1 : PetscCall(PetscFree2(kr,ki));
161 1 : PetscCall(VecDestroyVecs(3,&z));
162 1 : PetscCall(VecDestroyVecs(nconv,&xr));
163 1 : PetscCall(VecDestroyVecs(nconv,&xi));
164 1 : if (twosided) {
165 1 : PetscCall(VecDestroyVecs(nconv,&yr));
166 1 : PetscCall(VecDestroyVecs(nconv,&yi));
167 : }
168 1 : PetscCall(SlepcFinalize());
169 : return 0;
170 : }
171 :
172 : /*
173 : ComputeResidualNorm - Computes the norm of the residual vector
174 : associated with an eigenpair.
175 :
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 12 : PetscErrorCode ComputeResidualNorm(Mat A,Mat B,PetscBool trans,PetscScalar kr,PetscScalar ki,Vec xr,Vec xi,Vec *z,PetscReal *norm)
183 : {
184 12 : Vec u,w=NULL;
185 12 : PetscScalar alpha;
186 : #if !defined(PETSC_USE_COMPLEX)
187 : Vec v;
188 : PetscReal ni,nr;
189 : #endif
190 12 : PetscErrorCode (*matmult)(Mat,Vec,Vec) = trans? MatMultHermitianTranspose: MatMult;
191 :
192 12 : PetscFunctionBegin;
193 12 : u = z[0];
194 12 : if (B) w = z[2];
195 :
196 : #if !defined(PETSC_USE_COMPLEX)
197 : v = z[1];
198 : if (ki == 0 || PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
199 : #endif
200 12 : PetscCall((*matmult)(A,xr,u)); /* u=A*x */
201 12 : if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
202 12 : if (B) PetscCall((*matmult)(B,xr,w)); /* w=B*x */
203 : else w = xr;
204 12 : alpha = trans? -PetscConj(kr): -kr;
205 12 : PetscCall(VecAXPY(u,alpha,w)); /* u=A*x-k*B*x */
206 : }
207 12 : 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 12 : PetscFunctionReturn(PETSC_SUCCESS);
230 : }
231 :
232 : /*TEST
233 :
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
246 :
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
253 :
254 : TEST*/
|
