Actual source code: epsdefault.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: */
 10: /*
 11:    This file contains some simple default routines for common operations
 12: */

 14: #include <slepc/private/epsimpl.h>   /*I "slepceps.h" I*/
 15: #include <slepcvec.h>

 17: PetscErrorCode EPSBackTransform_Default(EPS eps)
 18: {

 22:   STBackTransform(eps->st,eps->nconv,eps->eigr,eps->eigi);
 23:   return(0);
 24: }

 26: /*
 27:   EPSComputeVectors_Hermitian - Copies the Lanczos vectors as eigenvectors
 28:   using purification for generalized eigenproblems.
 29:  */
 30: PetscErrorCode EPSComputeVectors_Hermitian(EPS eps)
 31: {
 33:   PetscInt       i;
 34:   PetscScalar    dot;
 35:   PetscReal      norm;
 36:   PetscBool      iscayley;
 37:   Mat            B;
 38:   Vec            w,x;

 41:   if (eps->purify) {
 42:     EPS_Purify(eps,eps->nconv);
 43:     for (i=0;i<eps->nconv;i++) {
 44:       BVNormColumn(eps->V,i,NORM_2,&norm);
 45:       BVScaleColumn(eps->V,i,1.0/norm);
 46:     }
 47:   } else {
 48:     /* In the case of Cayley transform, eigenvectors need to be B-normalized */
 49:     PetscObjectTypeCompare((PetscObject)eps->st,STCAYLEY,&iscayley);
 50:     if (iscayley && eps->isgeneralized) {
 51:       STGetMatrix(eps->st,1,&B);
 52:       MatCreateVecs(B,NULL,&w);
 53:       for (i=0;i<eps->nconv;i++) {
 54:         BVGetColumn(eps->V,i,&x);
 55:         MatMult(B,x,w);
 56:         VecDot(w,x,&dot);
 57:         VecScale(x,1.0/PetscSqrtScalar(dot));
 58:         BVRestoreColumn(eps->V,i,&x);
 59:       }
 60:       VecDestroy(&w);
 61:     }
 62:   }
 63:   return(0);
 64: }

 66: /*
 67:   EPSComputeVectors_Indefinite - similar to the Schur version but
 68:   for indefinite problems
 69:  */
 70: PetscErrorCode EPSComputeVectors_Indefinite(EPS eps)
 71: {
 73:   PetscInt       n,i;
 74:   Mat            X;
 75:   Vec            v;
 76: #if !defined(PETSC_USE_COMPLEX)
 77:   Vec            v1;
 78: #endif

 81:   DSGetDimensions(eps->ds,&n,NULL,NULL,NULL,NULL);
 82:   DSVectors(eps->ds,DS_MAT_X,NULL,NULL);
 83:   DSGetMat(eps->ds,DS_MAT_X,&X);
 84:   BVSetActiveColumns(eps->V,0,n);
 85:   BVMultInPlace(eps->V,X,0,n);
 86:   MatDestroy(&X);

 88:   /* purification */
 89:   if (eps->purify) { EPS_Purify(eps,eps->nconv); }

 91:   /* normalization */
 92:   for (i=0;i<n;i++) {
 93: #if !defined(PETSC_USE_COMPLEX)
 94:     if (eps->eigi[i] != 0.0) {
 95:       BVGetColumn(eps->V,i,&v);
 96:       BVGetColumn(eps->V,i+1,&v1);
 97:       VecNormalizeComplex(v,v1,PETSC_TRUE,NULL);
 98:       BVRestoreColumn(eps->V,i,&v);
 99:       BVRestoreColumn(eps->V,i+1,&v1);
100:       i++;
101:     } else
102: #endif
103:     {
104:       BVGetColumn(eps->V,i,&v);
105:       VecNormalize(v,NULL);
106:       BVRestoreColumn(eps->V,i,&v);
107:     }
108:   }
109:   return(0);
110: }

112: /*
113:   EPSComputeVectors_Twosided - Adjust left eigenvectors in generalized problems: y = B^-* y.
114:  */
115: PetscErrorCode EPSComputeVectors_Twosided(EPS eps)
116: {
118:   PetscInt       i;
119:   Vec            w,y;

122:   if (!eps->twosided || !eps->isgeneralized) return(0);
123:   EPSSetWorkVecs(eps,1);
124:   w = eps->work[0];
125:   for (i=0;i<eps->nconv;i++) {
126:     BVCopyVec(eps->W,i,w);
127:     VecConjugate(w);
128:     BVGetColumn(eps->W,i,&y);
129:     STMatSolveTranspose(eps->st,w,y);
130:     VecConjugate(y);
131:     BVRestoreColumn(eps->W,i,&y);
132:   }
133:   return(0);
134: }

136: /*
137:   EPSComputeVectors_Schur - Compute eigenvectors from the vectors
138:   provided by the eigensolver. This version is intended for solvers
139:   that provide Schur vectors. Given the partial Schur decomposition
140:   OP*V=V*T, the following steps are performed:
141:       1) compute eigenvectors of T: T*Z=Z*D
142:       2) compute eigenvectors of OP: X=V*Z
143:  */
144: PetscErrorCode EPSComputeVectors_Schur(EPS eps)
145: {
147:   PetscInt       n,i;
148:   Mat            Z;
149:   Vec            z,v;
150: #if !defined(PETSC_USE_COMPLEX)
151:   Vec            v1;
152: #endif

155:   if (eps->ishermitian) {
156:     if (eps->isgeneralized && !eps->ispositive) {
157:        EPSComputeVectors_Indefinite(eps);
158:     } else {
159:       EPSComputeVectors_Hermitian(eps);
160:     }
161:     return(0);
162:   }
163:   DSGetDimensions(eps->ds,&n,NULL,NULL,NULL,NULL);

165:   /* right eigenvectors */
166:   DSVectors(eps->ds,DS_MAT_X,NULL,NULL);

168:   /* V = V * Z */
169:   DSGetMat(eps->ds,DS_MAT_X,&Z);
170:   BVSetActiveColumns(eps->V,0,n);
171:   BVMultInPlace(eps->V,Z,0,n);
172:   MatDestroy(&Z);

174:   /* Purify eigenvectors */
175:   if (eps->purify) { EPS_Purify(eps,eps->nconv); }

177:   /* Fix eigenvectors if balancing was used */
178:   if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
179:     for (i=0;i<n;i++) {
180:       BVGetColumn(eps->V,i,&z);
181:       VecPointwiseDivide(z,z,eps->D);
182:       BVRestoreColumn(eps->V,i,&z);
183:     }
184:   }

186:   /* normalize eigenvectors (when using purification or balancing) */
187:   if (eps->purify || (eps->balance!=EPS_BALANCE_NONE && eps->D)) {
188:     for (i=0;i<n;i++) {
189: #if !defined(PETSC_USE_COMPLEX)
190:       if (eps->eigi[i] != 0.0) {
191:         BVGetColumn(eps->V,i,&v);
192:         BVGetColumn(eps->V,i+1,&v1);
193:         VecNormalizeComplex(v,v1,PETSC_TRUE,NULL);
194:         BVRestoreColumn(eps->V,i,&v);
195:         BVRestoreColumn(eps->V,i+1,&v1);
196:         i++;
197:       } else
198: #endif
199:       {
200:         BVGetColumn(eps->V,i,&v);
201:         VecNormalize(v,NULL);
202:         BVRestoreColumn(eps->V,i,&v);
203:       }
204:     }
205:   }

207:   /* left eigenvectors from eps->dsts */
208:   if (eps->twosided) {
209:     DSVectors(eps->dsts,DS_MAT_X,NULL,NULL);
210:     /* W = W * Z */
211:     DSGetMat(eps->dsts,DS_MAT_X,&Z);
212:     BVSetActiveColumns(eps->W,0,n);
213:     BVMultInPlace(eps->W,Z,0,n);
214:     MatDestroy(&Z);
215:     /* Fix left eigenvectors if balancing was used */
216:     if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
217:       for (i=0;i<n;i++) {
218:         BVGetColumn(eps->W,i,&z);
219:         VecPointwiseMult(z,z,eps->D);
220:         BVRestoreColumn(eps->W,i,&z);
221:       }
222:     }
223:     EPSComputeVectors_Twosided(eps);
224:     /* normalize */
225:     for (i=0;i<eps->nconv;i++) {
226: #if !defined(PETSC_USE_COMPLEX)
227:       if (eps->eigi[i]!=0.0) {   /* first eigenvalue of a complex conjugate pair */
228:         BVGetColumn(eps->W,i,&v);
229:         BVGetColumn(eps->W,i+1,&v1);
230:         VecNormalizeComplex(v,v1,PETSC_TRUE,NULL);
231:         BVRestoreColumn(eps->W,i,&v);
232:         BVRestoreColumn(eps->W,i+1,&v1);
233:         i++;
234:       } else   /* real eigenvalue */
235: #endif
236:       {
237:         BVGetColumn(eps->W,i,&v);
238:         VecNormalize(v,NULL);
239:         BVRestoreColumn(eps->W,i,&v);
240:       }
241:     }
242: #if !defined(PETSC_USE_COMPLEX)
243:     for (i=0;i<eps->nconv-1;i++) {
244:       if (eps->eigi[i] != 0.0) {
245:         if (eps->eigi[i] > 0.0) { BVScaleColumn(eps->W,i+1,-1.0); }
246:         i++;
247:       }
248:     }
249: #endif
250:   }
251:   return(0);
252: }

254: /*@
255:    EPSSetWorkVecs - Sets a number of work vectors into an EPS object.

257:    Collective on EPS

259:    Input Parameters:
260: +  eps - eigensolver context
261: -  nw  - number of work vectors to allocate

263:    Developers Note:
264:    This is SLEPC_EXTERN because it may be required by user plugin EPS
265:    implementations.

267:    Level: developer
268: @*/
269: PetscErrorCode EPSSetWorkVecs(EPS eps,PetscInt nw)
270: {
272:   Vec            t;

277:   if (nw <= 0) SETERRQ1(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"nw must be > 0: nw = %D",nw);
278:   if (eps->nwork < nw) {
279:     VecDestroyVecs(eps->nwork,&eps->work);
280:     eps->nwork = nw;
281:     BVGetColumn(eps->V,0,&t);
282:     VecDuplicateVecs(t,nw,&eps->work);
283:     BVRestoreColumn(eps->V,0,&t);
284:     PetscLogObjectParents(eps,nw,eps->work);
285:   }
286:   return(0);
287: }

289: /*
290:   EPSSetWhichEigenpairs_Default - Sets the default value for which,
291:   depending on the ST.
292:  */
293: PetscErrorCode EPSSetWhichEigenpairs_Default(EPS eps)
294: {
295:   PetscBool      target;

299:   PetscObjectTypeCompareAny((PetscObject)eps->st,&target,STSINVERT,STCAYLEY,"");
300:   if (target) eps->which = EPS_TARGET_MAGNITUDE;
301:   else eps->which = EPS_LARGEST_MAGNITUDE;
302:   return(0);
303: }

305: /*
306:   EPSConvergedRelative - Checks convergence relative to the eigenvalue.
307: */
308: PetscErrorCode EPSConvergedRelative(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
309: {
310:   PetscReal w;

313:   w = SlepcAbsEigenvalue(eigr,eigi);
314:   *errest = res/w;
315:   return(0);
316: }

318: /*
319:   EPSConvergedAbsolute - Checks convergence absolutely.
320: */
321: PetscErrorCode EPSConvergedAbsolute(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
322: {
324:   *errest = res;
325:   return(0);
326: }

328: /*
329:   EPSConvergedNorm - Checks convergence relative to the eigenvalue and
330:   the matrix norms.
331: */
332: PetscErrorCode EPSConvergedNorm(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
333: {
334:   PetscReal w;

337:   w = SlepcAbsEigenvalue(eigr,eigi);
338:   *errest = res / (eps->nrma + w*eps->nrmb);
339:   return(0);
340: }

342: /*@C
343:    EPSStoppingBasic - Default routine to determine whether the outer eigensolver
344:    iteration must be stopped.

346:    Collective on EPS

348:    Input Parameters:
349: +  eps    - eigensolver context obtained from EPSCreate()
350: .  its    - current number of iterations
351: .  max_it - maximum number of iterations
352: .  nconv  - number of currently converged eigenpairs
353: .  nev    - number of requested eigenpairs
354: -  ctx    - context (not used here)

356:    Output Parameter:
357: .  reason - result of the stopping test

359:    Notes:
360:    A positive value of reason indicates that the iteration has finished successfully
361:    (converged), and a negative value indicates an error condition (diverged). If
362:    the iteration needs to be continued, reason must be set to EPS_CONVERGED_ITERATING
363:    (zero).

365:    EPSStoppingBasic() will stop if all requested eigenvalues are converged, or if
366:    the maximum number of iterations has been reached.

368:    Use EPSSetStoppingTest() to provide your own test instead of using this one.

370:    Level: advanced

372: .seealso: EPSSetStoppingTest(), EPSConvergedReason, EPSGetConvergedReason()
373: @*/
374: PetscErrorCode EPSStoppingBasic(EPS eps,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nev,EPSConvergedReason *reason,void *ctx)
375: {

379:   *reason = EPS_CONVERGED_ITERATING;
380:   if (nconv >= nev) {
381:     PetscInfo2(eps,"Linear eigensolver finished successfully: %D eigenpairs converged at iteration %D\n",nconv,its);
382:     *reason = EPS_CONVERGED_TOL;
383:   } else if (its >= max_it) {
384:     *reason = EPS_DIVERGED_ITS;
385:     PetscInfo1(eps,"Linear eigensolver iteration reached maximum number of iterations (%D)\n",its);
386:   }
387:   return(0);
388: }

390: /*
391:   EPSComputeRitzVector - Computes the current Ritz vector.

393:   Simple case (complex scalars or real scalars with Zi=NULL):
394:     x = V*Zr  (V is a basis of nv vectors, Zr has length nv)

396:   Split case:
397:     x = V*Zr  y = V*Zi  (Zr and Zi have length nv)
398: */
399: PetscErrorCode EPSComputeRitzVector(EPS eps,PetscScalar *Zr,PetscScalar *Zi,BV V,Vec x,Vec y)
400: {
402:   PetscInt       l,k;
403:   PetscReal      norm;
404: #if !defined(PETSC_USE_COMPLEX)
405:   Vec            z;
406: #endif

409:   /* compute eigenvector */
410:   BVGetActiveColumns(V,&l,&k);
411:   BVSetActiveColumns(V,0,k);
412:   BVMultVec(V,1.0,0.0,x,Zr);

414:   /* purify eigenvector if necessary */
415:   if (eps->purify) {
416:     STApply(eps->st,x,y);
417:     if (eps->ishermitian) {
418:       BVNormVec(eps->V,y,NORM_2,&norm);
419:     } else {
420:       VecNorm(y,NORM_2,&norm);
421:     }
422:     VecScale(y,1.0/norm);
423:     VecCopy(y,x);
424:   }
425:   /* fix eigenvector if balancing is used */
426:   if (!eps->ishermitian && eps->balance!=EPS_BALANCE_NONE && eps->D) {
427:     VecPointwiseDivide(x,x,eps->D);
428:   }
429: #if !defined(PETSC_USE_COMPLEX)
430:   /* compute imaginary part of eigenvector */
431:   if (Zi) {
432:     BVMultVec(V,1.0,0.0,y,Zi);
433:     if (eps->ispositive) {
434:       BVCreateVec(V,&z);
435:       STApply(eps->st,y,z);
436:       VecNorm(z,NORM_2,&norm);
437:       VecScale(z,1.0/norm);
438:       VecCopy(z,y);
439:       VecDestroy(&z);
440:     }
441:     if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
442:       VecPointwiseDivide(y,y,eps->D);
443:     }
444:   } else
445: #endif
446:   { VecSet(y,0.0); }

448:   /* normalize eigenvectors (when using balancing) */
449:   if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
450: #if !defined(PETSC_USE_COMPLEX)
451:     if (Zi) {
452:       VecNormalizeComplex(x,y,PETSC_TRUE,NULL);
453:     } else
454: #endif
455:     {
456:       VecNormalize(x,NULL);
457:     }
458:   }
459:   BVSetActiveColumns(V,l,k);
460:   return(0);
461: }

463: /*
464:   EPSBuildBalance_Krylov - uses a Krylov subspace method to compute the
465:   diagonal matrix to be applied for balancing in non-Hermitian problems.
466: */
467: PetscErrorCode EPSBuildBalance_Krylov(EPS eps)
468: {
469:   Vec               z,p,r;
470:   PetscInt          i,j;
471:   PetscReal         norma;
472:   PetscScalar       *pz,*pD;
473:   const PetscScalar *pr,*pp;
474:   PetscRandom       rand;
475:   PetscErrorCode    ierr;

478:   EPSSetWorkVecs(eps,3);
479:   BVGetRandomContext(eps->V,&rand);
480:   r = eps->work[0];
481:   p = eps->work[1];
482:   z = eps->work[2];
483:   VecSet(eps->D,1.0);

485:   for (j=0;j<eps->balance_its;j++) {

487:     /* Build a random vector of +-1's */
488:     VecSetRandom(z,rand);
489:     VecGetArray(z,&pz);
490:     for (i=0;i<eps->nloc;i++) {
491:       if (PetscRealPart(pz[i])<0.5) pz[i]=-1.0;
492:       else pz[i]=1.0;
493:     }
494:     VecRestoreArray(z,&pz);

496:     /* Compute p=DA(D\z) */
497:     VecPointwiseDivide(r,z,eps->D);
498:     STApply(eps->st,r,p);
499:     VecPointwiseMult(p,p,eps->D);
500:     if (eps->balance == EPS_BALANCE_TWOSIDE) {
501:       if (j==0) {
502:         /* Estimate the matrix inf-norm */
503:         VecAbs(p);
504:         VecMax(p,NULL,&norma);
505:       }
506:       /* Compute r=D\(A'Dz) */
507:       VecPointwiseMult(z,z,eps->D);
508:       STApplyTranspose(eps->st,z,r);
509:       VecPointwiseDivide(r,r,eps->D);
510:     }

512:     /* Adjust values of D */
513:     VecGetArrayRead(r,&pr);
514:     VecGetArrayRead(p,&pp);
515:     VecGetArray(eps->D,&pD);
516:     for (i=0;i<eps->nloc;i++) {
517:       if (eps->balance == EPS_BALANCE_TWOSIDE) {
518:         if (PetscAbsScalar(pp[i])>eps->balance_cutoff*norma && pr[i]!=0.0)
519:           pD[i] *= PetscSqrtReal(PetscAbsScalar(pr[i]/pp[i]));
520:       } else {
521:         if (pp[i]!=0.0) pD[i] /= PetscAbsScalar(pp[i]);
522:       }
523:     }
524:     VecRestoreArrayRead(r,&pr);
525:     VecRestoreArrayRead(p,&pp);
526:     VecRestoreArray(eps->D,&pD);
527:   }
528:   return(0);
529: }