Actual source code: epsdefault.c

slepc-3.16.0 2021-09-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2021, 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>
 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:   PetscBool      iscayley,indef;
 34:   Mat            B,C;

 37:   if (eps->purify) {
 38:     EPS_Purify(eps,eps->nconv);
 39:     BVNormalize(eps->V,NULL);
 40:   } else {
 41:     /* In the case of Cayley transform, eigenvectors need to be B-normalized */
 42:     PetscObjectTypeCompare((PetscObject)eps->st,STCAYLEY,&iscayley);
 43:     if (iscayley && eps->isgeneralized) {
 44:       STGetMatrix(eps->st,1,&B);
 45:       BVGetMatrix(eps->V,&C,&indef);
 46:       if (indef) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"The inner product should not be indefinite");
 47:       BVSetMatrix(eps->V,B,PETSC_FALSE);
 48:       BVNormalize(eps->V,NULL);
 49:       BVSetMatrix(eps->V,C,PETSC_FALSE);  /* restore original matrix */
 50:     }
 51:   }
 52:   return(0);
 53: }

 55: /*
 56:   EPSComputeVectors_Indefinite - similar to the Schur version but
 57:   for indefinite problems
 58:  */
 59: PetscErrorCode EPSComputeVectors_Indefinite(EPS eps)
 60: {
 62:   PetscInt       n;
 63:   Mat            X;

 66:   DSGetDimensions(eps->ds,&n,NULL,NULL,NULL);
 67:   DSVectors(eps->ds,DS_MAT_X,NULL,NULL);
 68:   DSGetMat(eps->ds,DS_MAT_X,&X);
 69:   BVMultInPlace(eps->V,X,0,n);
 70:   MatDestroy(&X);

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

 75:   /* normalization */
 76:   BVNormalize(eps->V,eps->eigi);
 77:   return(0);
 78: }

 80: /*
 81:   EPSComputeVectors_Twosided - Adjust left eigenvectors in generalized problems: y = B^-* y.
 82:  */
 83: PetscErrorCode EPSComputeVectors_Twosided(EPS eps)
 84: {
 86:   PetscInt       i;
 87:   Vec            w,y;

 90:   if (!eps->twosided || !eps->isgeneralized) return(0);
 91:   EPSSetWorkVecs(eps,1);
 92:   w = eps->work[0];
 93:   for (i=0;i<eps->nconv;i++) {
 94:     BVCopyVec(eps->W,i,w);
 95:     VecConjugate(w);
 96:     BVGetColumn(eps->W,i,&y);
 97:     STMatSolveTranspose(eps->st,w,y);
 98:     VecConjugate(y);
 99:     BVRestoreColumn(eps->W,i,&y);
100:   }
101:   return(0);
102: }

104: /*
105:   EPSComputeVectors_Schur - Compute eigenvectors from the vectors
106:   provided by the eigensolver. This version is intended for solvers
107:   that provide Schur vectors. Given the partial Schur decomposition
108:   OP*V=V*T, the following steps are performed:
109:       1) compute eigenvectors of T: T*Z=Z*D
110:       2) compute eigenvectors of OP: X=V*Z
111:  */
112: PetscErrorCode EPSComputeVectors_Schur(EPS eps)
113: {
115:   PetscInt       i;
116:   Mat            Z;
117:   Vec            z;

120:   if (eps->ishermitian) {
121:     if (eps->isgeneralized && !eps->ispositive) {
122:       EPSComputeVectors_Indefinite(eps);
123:     } else {
124:       EPSComputeVectors_Hermitian(eps);
125:     }
126:     return(0);
127:   }

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

132:   /* V = V * Z */
133:   DSGetMat(eps->ds,DS_MAT_X,&Z);
134:   BVMultInPlace(eps->V,Z,0,eps->nconv);
135:   MatDestroy(&Z);

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

140:   /* Fix eigenvectors if balancing was used */
141:   if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
142:     for (i=0;i<eps->nconv;i++) {
143:       BVGetColumn(eps->V,i,&z);
144:       VecPointwiseDivide(z,z,eps->D);
145:       BVRestoreColumn(eps->V,i,&z);
146:     }
147:   }

149:   /* normalize eigenvectors (when using purification or balancing) */
150:   if (eps->purify || (eps->balance!=EPS_BALANCE_NONE && eps->D)) {
151:     BVNormalize(eps->V,eps->eigi);
152:   }

154:   /* left eigenvectors */
155:   if (eps->twosided) {
156:     DSVectors(eps->ds,DS_MAT_Y,NULL,NULL);
157:     /* W = W * Z */
158:     DSGetMat(eps->ds,DS_MAT_Y,&Z);
159:     BVMultInPlace(eps->W,Z,0,eps->nconv);
160:     MatDestroy(&Z);
161:     /* Fix left eigenvectors if balancing was used */
162:     if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
163:       for (i=0;i<eps->nconv;i++) {
164:         BVGetColumn(eps->W,i,&z);
165:         VecPointwiseMult(z,z,eps->D);
166:         BVRestoreColumn(eps->W,i,&z);
167:       }
168:     }
169:     EPSComputeVectors_Twosided(eps);
170:     /* normalize */
171:     BVNormalize(eps->W,eps->eigi);
172: #if !defined(PETSC_USE_COMPLEX)
173:     for (i=0;i<eps->nconv-1;i++) {
174:       if (eps->eigi[i] != 0.0) {
175:         if (eps->eigi[i] > 0.0) { BVScaleColumn(eps->W,i+1,-1.0); }
176:         i++;
177:       }
178:     }
179: #endif
180:   }
181:   return(0);
182: }

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

187:    Collective on eps

189:    Input Parameters:
190: +  eps - eigensolver context
191: -  nw  - number of work vectors to allocate

193:    Developers Note:
194:    This is SLEPC_EXTERN because it may be required by user plugin EPS
195:    implementations.

197:    Level: developer
198: @*/
199: PetscErrorCode EPSSetWorkVecs(EPS eps,PetscInt nw)
200: {
202:   Vec            t;

207:   if (nw <= 0) SETERRQ1(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"nw must be > 0: nw = %D",nw);
208:   if (eps->nwork < nw) {
209:     VecDestroyVecs(eps->nwork,&eps->work);
210:     eps->nwork = nw;
211:     BVGetColumn(eps->V,0,&t);
212:     VecDuplicateVecs(t,nw,&eps->work);
213:     BVRestoreColumn(eps->V,0,&t);
214:     PetscLogObjectParents(eps,nw,eps->work);
215:   }
216:   return(0);
217: }

219: /*
220:   EPSSetWhichEigenpairs_Default - Sets the default value for which,
221:   depending on the ST.
222:  */
223: PetscErrorCode EPSSetWhichEigenpairs_Default(EPS eps)
224: {
225:   PetscBool      target;

229:   PetscObjectTypeCompareAny((PetscObject)eps->st,&target,STSINVERT,STCAYLEY,"");
230:   if (target) eps->which = EPS_TARGET_MAGNITUDE;
231:   else eps->which = EPS_LARGEST_MAGNITUDE;
232:   return(0);
233: }

235: /*
236:   EPSConvergedRelative - Checks convergence relative to the eigenvalue.
237: */
238: PetscErrorCode EPSConvergedRelative(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
239: {
240:   PetscReal w;

243:   w = SlepcAbsEigenvalue(eigr,eigi);
244:   *errest = res/w;
245:   return(0);
246: }

248: /*
249:   EPSConvergedAbsolute - Checks convergence absolutely.
250: */
251: PetscErrorCode EPSConvergedAbsolute(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
252: {
254:   *errest = res;
255:   return(0);
256: }

258: /*
259:   EPSConvergedNorm - Checks convergence relative to the eigenvalue and
260:   the matrix norms.
261: */
262: PetscErrorCode EPSConvergedNorm(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
263: {
264:   PetscReal w;

267:   w = SlepcAbsEigenvalue(eigr,eigi);
268:   *errest = res / (eps->nrma + w*eps->nrmb);
269:   return(0);
270: }

272: /*@C
273:    EPSStoppingBasic - Default routine to determine whether the outer eigensolver
274:    iteration must be stopped.

276:    Collective on eps

278:    Input Parameters:
279: +  eps    - eigensolver context obtained from EPSCreate()
280: .  its    - current number of iterations
281: .  max_it - maximum number of iterations
282: .  nconv  - number of currently converged eigenpairs
283: .  nev    - number of requested eigenpairs
284: -  ctx    - context (not used here)

286:    Output Parameter:
287: .  reason - result of the stopping test

289:    Notes:
290:    A positive value of reason indicates that the iteration has finished successfully
291:    (converged), and a negative value indicates an error condition (diverged). If
292:    the iteration needs to be continued, reason must be set to EPS_CONVERGED_ITERATING
293:    (zero).

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

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

300:    Level: advanced

302: .seealso: EPSSetStoppingTest(), EPSConvergedReason, EPSGetConvergedReason()
303: @*/
304: PetscErrorCode EPSStoppingBasic(EPS eps,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nev,EPSConvergedReason *reason,void *ctx)
305: {

309:   *reason = EPS_CONVERGED_ITERATING;
310:   if (nconv >= nev) {
311:     PetscInfo2(eps,"Linear eigensolver finished successfully: %D eigenpairs converged at iteration %D\n",nconv,its);
312:     *reason = EPS_CONVERGED_TOL;
313:   } else if (its >= max_it) {
314:     *reason = EPS_DIVERGED_ITS;
315:     PetscInfo1(eps,"Linear eigensolver iteration reached maximum number of iterations (%D)\n",its);
316:   }
317:   return(0);
318: }

320: /*
321:   EPSComputeRitzVector - Computes the current Ritz vector.

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

326:   Split case:
327:     x = V*Zr  y = V*Zi  (Zr and Zi have length nv)
328: */
329: PetscErrorCode EPSComputeRitzVector(EPS eps,PetscScalar *Zr,PetscScalar *Zi,BV V,Vec x,Vec y)
330: {
332:   PetscInt       l,k;
333:   PetscReal      norm;
334: #if !defined(PETSC_USE_COMPLEX)
335:   Vec            z;
336: #endif

339:   /* compute eigenvector */
340:   BVGetActiveColumns(V,&l,&k);
341:   BVSetActiveColumns(V,0,k);
342:   BVMultVec(V,1.0,0.0,x,Zr);

344:   /* purify eigenvector if necessary */
345:   if (eps->purify) {
346:     STApply(eps->st,x,y);
347:     if (eps->ishermitian) {
348:       BVNormVec(eps->V,y,NORM_2,&norm);
349:     } else {
350:       VecNorm(y,NORM_2,&norm);
351:     }
352:     VecScale(y,1.0/norm);
353:     VecCopy(y,x);
354:   }
355:   /* fix eigenvector if balancing is used */
356:   if (!eps->ishermitian && eps->balance!=EPS_BALANCE_NONE && eps->D) {
357:     VecPointwiseDivide(x,x,eps->D);
358:   }
359: #if !defined(PETSC_USE_COMPLEX)
360:   /* compute imaginary part of eigenvector */
361:   if (Zi) {
362:     BVMultVec(V,1.0,0.0,y,Zi);
363:     if (eps->ispositive) {
364:       BVCreateVec(V,&z);
365:       STApply(eps->st,y,z);
366:       VecNorm(z,NORM_2,&norm);
367:       VecScale(z,1.0/norm);
368:       VecCopy(z,y);
369:       VecDestroy(&z);
370:     }
371:     if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
372:       VecPointwiseDivide(y,y,eps->D);
373:     }
374:   } else
375: #endif
376:   { VecSet(y,0.0); }

378:   /* normalize eigenvectors (when using balancing) */
379:   if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
380: #if !defined(PETSC_USE_COMPLEX)
381:     if (Zi) {
382:       VecNormalizeComplex(x,y,PETSC_TRUE,NULL);
383:     } else
384: #endif
385:     {
386:       VecNormalize(x,NULL);
387:     }
388:   }
389:   BVSetActiveColumns(V,l,k);
390:   return(0);
391: }

393: /*
394:   EPSBuildBalance_Krylov - uses a Krylov subspace method to compute the
395:   diagonal matrix to be applied for balancing in non-Hermitian problems.
396: */
397: PetscErrorCode EPSBuildBalance_Krylov(EPS eps)
398: {
399:   Vec               z,p,r;
400:   PetscInt          i,j;
401:   PetscReal         norma;
402:   PetscScalar       *pz,*pD;
403:   const PetscScalar *pr,*pp;
404:   PetscRandom       rand;
405:   PetscErrorCode    ierr;

408:   EPSSetWorkVecs(eps,3);
409:   BVGetRandomContext(eps->V,&rand);
410:   r = eps->work[0];
411:   p = eps->work[1];
412:   z = eps->work[2];
413:   VecSet(eps->D,1.0);

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

417:     /* Build a random vector of +-1's */
418:     VecSetRandom(z,rand);
419:     VecGetArray(z,&pz);
420:     for (i=0;i<eps->nloc;i++) {
421:       if (PetscRealPart(pz[i])<0.5) pz[i]=-1.0;
422:       else pz[i]=1.0;
423:     }
424:     VecRestoreArray(z,&pz);

426:     /* Compute p=DA(D\z) */
427:     VecPointwiseDivide(r,z,eps->D);
428:     STApply(eps->st,r,p);
429:     VecPointwiseMult(p,p,eps->D);
430:     if (eps->balance == EPS_BALANCE_TWOSIDE) {
431:       if (j==0) {
432:         /* Estimate the matrix inf-norm */
433:         VecAbs(p);
434:         VecMax(p,NULL,&norma);
435:       }
436:       /* Compute r=D\(A'Dz) */
437:       VecPointwiseMult(z,z,eps->D);
438:       STApplyHermitianTranspose(eps->st,z,r);
439:       VecPointwiseDivide(r,r,eps->D);
440:     }

442:     /* Adjust values of D */
443:     VecGetArrayRead(r,&pr);
444:     VecGetArrayRead(p,&pp);
445:     VecGetArray(eps->D,&pD);
446:     for (i=0;i<eps->nloc;i++) {
447:       if (eps->balance == EPS_BALANCE_TWOSIDE) {
448:         if (PetscAbsScalar(pp[i])>eps->balance_cutoff*norma && pr[i]!=0.0)
449:           pD[i] *= PetscSqrtReal(PetscAbsScalar(pr[i]/pp[i]));
450:       } else {
451:         if (pp[i]!=0.0) pD[i] /= PetscAbsScalar(pp[i]);
452:       }
453:     }
454:     VecRestoreArrayRead(r,&pr);
455:     VecRestoreArrayRead(p,&pp);
456:     VecRestoreArray(eps->D,&pD);
457:   }
458:   return(0);
459: }