Actual source code: epsdefault.c

slepc-3.10.1 2018-10-23
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2018, 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:   PetscScalar    tmp;
 79:   PetscReal      norm,normi;
 80: #endif

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

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

 93:   /* normalization */
 94:   for (i=0;i<n;i++) {
 95: #if !defined(PETSC_USE_COMPLEX)
 96:     if (eps->eigi[i] != 0.0) {
 97:       BVGetColumn(eps->V,i,&v);
 98:       BVGetColumn(eps->V,i+1,&v1);
 99:       VecNorm(v,NORM_2,&norm);
100:       VecNorm(v1,NORM_2,&normi);
101:       tmp = 1.0 / SlepcAbsEigenvalue(norm,normi);
102:       VecScale(v,tmp);
103:       VecScale(v1,tmp);
104:       BVRestoreColumn(eps->V,i,&v);
105:       BVRestoreColumn(eps->V,i+1,&v1);
106:       i++;
107:     } else
108: #endif
109:     {
110:       BVGetColumn(eps->V,i,&v);
111:       VecNormalize(v,NULL);
112:       BVRestoreColumn(eps->V,i,&v);
113:     }
114:   }
115:   return(0);
116: }

118: /*
119:   EPSComputeVectors_Schur - Compute eigenvectors from the vectors
120:   provided by the eigensolver. This version is intended for solvers
121:   that provide Schur vectors. Given the partial Schur decomposition
122:   OP*V=V*T, the following steps are performed:
123:       1) compute eigenvectors of T: T*Z=Z*D
124:       2) compute eigenvectors of OP: X=V*Z
125:  */
126: PetscErrorCode EPSComputeVectors_Schur(EPS eps)
127: {
129:   PetscInt       n,i;
130:   Mat            Z;
131:   Vec            z,v;
132: #if !defined(PETSC_USE_COMPLEX)
133:   Vec            v1;
134:   PetscScalar    tmp;
135:   PetscReal      norm,normi;
136: #endif

139:   if (eps->ishermitian) {
140:     if (eps->isgeneralized && !eps->ispositive) {
141:        EPSComputeVectors_Indefinite(eps);
142:     } else {
143:       EPSComputeVectors_Hermitian(eps);
144:     }
145:     return(0);
146:   }
147:   DSGetDimensions(eps->ds,&n,NULL,NULL,NULL,NULL);

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

152:   /* V = V * Z */
153:   DSGetMat(eps->ds,DS_MAT_X,&Z);
154:   BVSetActiveColumns(eps->V,0,n);
155:   BVMultInPlace(eps->V,Z,0,n);
156:   MatDestroy(&Z);

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

161:   /* Fix eigenvectors if balancing was used */
162:   if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
163:     for (i=0;i<n;i++) {
164:       BVGetColumn(eps->V,i,&z);
165:       VecPointwiseDivide(z,z,eps->D);
166:       BVRestoreColumn(eps->V,i,&z);
167:     }
168:   }

170:   /* normalize eigenvectors (when using purification or balancing) */
171:   if (eps->purify || (eps->balance!=EPS_BALANCE_NONE && eps->D)) {
172:     for (i=0;i<n;i++) {
173: #if !defined(PETSC_USE_COMPLEX)
174:       if (eps->eigi[i] != 0.0) {
175:         BVGetColumn(eps->V,i,&v);
176:         BVGetColumn(eps->V,i+1,&v1);
177:         VecNorm(v,NORM_2,&norm);
178:         VecNorm(v1,NORM_2,&normi);
179:         tmp = 1.0 / SlepcAbsEigenvalue(norm,normi);
180:         VecScale(v,tmp);
181:         VecScale(v1,tmp);
182:         BVRestoreColumn(eps->V,i,&v);
183:         BVRestoreColumn(eps->V,i+1,&v1);
184:         i++;
185:       } else
186: #endif
187:       {
188:         BVGetColumn(eps->V,i,&v);
189:         VecNormalize(v,NULL);
190:         BVRestoreColumn(eps->V,i,&v);
191:       }
192:     }
193:   }
194:   return(0);
195: }

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

200:    Collective on EPS

202:    Input Parameters:
203: +  eps - eigensolver context
204: -  nw  - number of work vectors to allocate

206:    Developers Note:
207:    This is PETSC_EXTERN because it may be required by user plugin EPS
208:    implementations.

210:    Level: developer
211: @*/
212: PetscErrorCode EPSSetWorkVecs(EPS eps,PetscInt nw)
213: {
215:   Vec            t;

220:   if (nw <= 0) SETERRQ1(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"nw must be > 0: nw = %D",nw);
221:   if (eps->nwork < nw) {
222:     VecDestroyVecs(eps->nwork,&eps->work);
223:     eps->nwork = nw;
224:     BVGetColumn(eps->V,0,&t);
225:     VecDuplicateVecs(t,nw,&eps->work);
226:     BVRestoreColumn(eps->V,0,&t);
227:     PetscLogObjectParents(eps,nw,eps->work);
228:   }
229:   return(0);
230: }

232: /*
233:   EPSSetWhichEigenpairs_Default - Sets the default value for which,
234:   depending on the ST.
235:  */
236: PetscErrorCode EPSSetWhichEigenpairs_Default(EPS eps)
237: {
238:   PetscBool      target;

242:   PetscObjectTypeCompareAny((PetscObject)eps->st,&target,STSINVERT,STCAYLEY,"");
243:   if (target) eps->which = EPS_TARGET_MAGNITUDE;
244:   else eps->which = EPS_LARGEST_MAGNITUDE;
245:   return(0);
246: }

248: /*
249:   EPSConvergedRelative - Checks convergence relative to the eigenvalue.
250: */
251: PetscErrorCode EPSConvergedRelative(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
252: {
253:   PetscReal w;

256:   w = SlepcAbsEigenvalue(eigr,eigi);
257:   *errest = res/w;
258:   return(0);
259: }

261: /*
262:   EPSConvergedAbsolute - Checks convergence absolutely.
263: */
264: PetscErrorCode EPSConvergedAbsolute(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
265: {
267:   *errest = res;
268:   return(0);
269: }

271: /*
272:   EPSConvergedNorm - Checks convergence relative to the eigenvalue and
273:   the matrix norms.
274: */
275: PetscErrorCode EPSConvergedNorm(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
276: {
277:   PetscReal w;

280:   w = SlepcAbsEigenvalue(eigr,eigi);
281:   *errest = res / (eps->nrma + w*eps->nrmb);
282:   return(0);
283: }

285: /*@C
286:    EPSStoppingBasic - Default routine to determine whether the outer eigensolver
287:    iteration must be stopped.

289:    Collective on EPS

291:    Input Parameters:
292: +  eps    - eigensolver context obtained from EPSCreate()
293: .  its    - current number of iterations
294: .  max_it - maximum number of iterations
295: .  nconv  - number of currently converged eigenpairs
296: .  nev    - number of requested eigenpairs
297: -  ctx    - context (not used here)

299:    Output Parameter:
300: .  reason - result of the stopping test

302:    Notes:
303:    A positive value of reason indicates that the iteration has finished successfully
304:    (converged), and a negative value indicates an error condition (diverged). If
305:    the iteration needs to be continued, reason must be set to EPS_CONVERGED_ITERATING
306:    (zero).

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

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

313:    Level: advanced

315: .seealso: EPSSetStoppingTest(), EPSConvergedReason, EPSGetConvergedReason()
316: @*/
317: PetscErrorCode EPSStoppingBasic(EPS eps,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nev,EPSConvergedReason *reason,void *ctx)
318: {

322:   *reason = EPS_CONVERGED_ITERATING;
323:   if (nconv >= nev) {
324:     PetscInfo2(eps,"Linear eigensolver finished successfully: %D eigenpairs converged at iteration %D\n",nconv,its);
325:     *reason = EPS_CONVERGED_TOL;
326:   } else if (its >= max_it) {
327:     *reason = EPS_DIVERGED_ITS;
328:     PetscInfo1(eps,"Linear eigensolver iteration reached maximum number of iterations (%D)\n",its);
329:   }
330:   return(0);
331: }

333: /*
334:   EPSComputeRitzVector - Computes the current Ritz vector.

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

339:   Split case:
340:     x = V*Zr  y = V*Zi  (Zr and Zi have length nv)
341: */
342: PetscErrorCode EPSComputeRitzVector(EPS eps,PetscScalar *Zr,PetscScalar *Zi,BV V,Vec x,Vec y)
343: {
345:   PetscInt       l,k;
346:   PetscReal      norm;
347: #if !defined(PETSC_USE_COMPLEX)
348:   Vec            z;
349:   PetscReal      normi;
350:   PetscScalar    tmp;
351: #endif

354:   /* compute eigenvector */
355:   BVGetActiveColumns(V,&l,&k);
356:   BVSetActiveColumns(V,0,k);
357:   BVMultVec(V,1.0,0.0,x,Zr);

359:   /* purify eigenvector if necessary */
360:   if (eps->purify) {
361:     STApply(eps->st,x,y);
362:     if (eps->ishermitian) {
363:       BVNormVec(eps->V,y,NORM_2,&norm);
364:     } else {
365:       VecNorm(y,NORM_2,&norm);
366:     }
367:     VecScale(y,1.0/norm);
368:     VecCopy(y,x);
369:   }
370:   /* fix eigenvector if balancing is used */
371:   if (!eps->ishermitian && eps->balance!=EPS_BALANCE_NONE && eps->D) {
372:     VecPointwiseDivide(x,x,eps->D);
373:   }
374: #if !defined(PETSC_USE_COMPLEX)
375:   /* compute imaginary part of eigenvector */
376:   if (Zi) {
377:     BVMultVec(V,1.0,0.0,y,Zi);
378:     if (eps->ispositive) {
379:       BVCreateVec(V,&z);
380:       STApply(eps->st,y,z);
381:       VecNorm(z,NORM_2,&norm);
382:       VecScale(z,1.0/norm);
383:       VecCopy(z,y);
384:       VecDestroy(&z);
385:     }
386:     if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
387:       VecPointwiseDivide(y,y,eps->D);
388:     }
389:   } else
390: #endif
391:   { VecSet(y,0.0); }

393:   /* normalize eigenvectors (when using balancing) */
394:   if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
395: #if !defined(PETSC_USE_COMPLEX)
396:     if (Zi) {
397:       VecNorm(x,NORM_2,&norm);
398:       VecNorm(y,NORM_2,&normi);
399:       tmp = 1.0 / SlepcAbsEigenvalue(norm,normi);
400:       VecScale(x,tmp);
401:       VecScale(y,tmp);
402:     } else
403: #endif
404:     {
405:       VecNormalize(x,NULL);
406:     }
407:   }
408:   BVSetActiveColumns(V,l,k);
409:   return(0);
410: }

412: /*
413:   EPSBuildBalance_Krylov - uses a Krylov subspace method to compute the
414:   diagonal matrix to be applied for balancing in non-Hermitian problems.
415: */
416: PetscErrorCode EPSBuildBalance_Krylov(EPS eps)
417: {
418:   Vec               z,p,r;
419:   PetscInt          i,j;
420:   PetscReal         norma;
421:   PetscScalar       *pz,*pD;
422:   const PetscScalar *pr,*pp;
423:   PetscRandom       rand;
424:   PetscErrorCode    ierr;

427:   EPSSetWorkVecs(eps,3);
428:   BVGetRandomContext(eps->V,&rand);
429:   r = eps->work[0];
430:   p = eps->work[1];
431:   z = eps->work[2];
432:   VecSet(eps->D,1.0);

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

436:     /* Build a random vector of +-1's */
437:     VecSetRandom(z,rand);
438:     VecGetArray(z,&pz);
439:     for (i=0;i<eps->nloc;i++) {
440:       if (PetscRealPart(pz[i])<0.5) pz[i]=-1.0;
441:       else pz[i]=1.0;
442:     }
443:     VecRestoreArray(z,&pz);

445:     /* Compute p=DA(D\z) */
446:     VecPointwiseDivide(r,z,eps->D);
447:     STApply(eps->st,r,p);
448:     VecPointwiseMult(p,p,eps->D);
449:     if (eps->balance == EPS_BALANCE_TWOSIDE) {
450:       if (j==0) {
451:         /* Estimate the matrix inf-norm */
452:         VecAbs(p);
453:         VecMax(p,NULL,&norma);
454:       }
455:       /* Compute r=D\(A'Dz) */
456:       VecPointwiseMult(z,z,eps->D);
457:       STApplyTranspose(eps->st,z,r);
458:       VecPointwiseDivide(r,r,eps->D);
459:     }

461:     /* Adjust values of D */
462:     VecGetArrayRead(r,&pr);
463:     VecGetArrayRead(p,&pp);
464:     VecGetArray(eps->D,&pD);
465:     for (i=0;i<eps->nloc;i++) {
466:       if (eps->balance == EPS_BALANCE_TWOSIDE) {
467:         if (PetscAbsScalar(pp[i])>eps->balance_cutoff*norma && pr[i]!=0.0)
468:           pD[i] *= PetscSqrtReal(PetscAbsScalar(pr[i]/pp[i]));
469:       } else {
470:         if (pp[i]!=0.0) pD[i] /= PetscAbsScalar(pp[i]);
471:       }
472:     }
473:     VecRestoreArrayRead(r,&pr);
474:     VecRestoreArrayRead(p,&pp);
475:     VecRestoreArray(eps->D,&pD);
476:   }
477:   return(0);
478: }