Actual source code: test29.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  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*/