Actual source code: epsdefault.c

slepc-3.19.2 2023-09-05
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, 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: {
 19:   PetscFunctionBegin;
 20:   PetscCall(STBackTransform(eps->st,eps->nconv,eps->eigr,eps->eigi));
 21:   PetscFunctionReturn(PETSC_SUCCESS);
 22: }

 24: /*
 25:   EPSComputeVectors_Hermitian - Copies the Lanczos vectors as eigenvectors
 26:   using purification for generalized eigenproblems.
 27:  */
 28: PetscErrorCode EPSComputeVectors_Hermitian(EPS eps)
 29: {
 30:   PetscBool      iscayley,indef;
 31:   Mat            B,C;

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

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

 61:   PetscFunctionBegin;
 62:   PetscCall(DSGetDimensions(eps->ds,&n,NULL,NULL,NULL));
 63:   PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));
 64:   PetscCall(DSGetMat(eps->ds,DS_MAT_X,&X));
 65:   PetscCall(BVMultInPlace(eps->V,X,0,n));
 66:   PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&X));

 68:   /* purification */
 69:   if (eps->purify) PetscCall(EPS_Purify(eps,eps->nconv));

 71:   /* normalization */
 72:   PetscCall(BVNormalize(eps->V,eps->eigi));
 73:   PetscFunctionReturn(PETSC_SUCCESS);
 74: }

 76: /*
 77:   EPSComputeVectors_Twosided - Adjust left eigenvectors in generalized problems: y = B^-* y.
 78:  */
 79: PetscErrorCode EPSComputeVectors_Twosided(EPS eps)
 80: {
 81:   PetscInt       i;
 82:   Vec            w,y;

 84:   PetscFunctionBegin;
 85:   if (!eps->twosided || !eps->isgeneralized) PetscFunctionReturn(PETSC_SUCCESS);
 86:   PetscCall(EPSSetWorkVecs(eps,1));
 87:   w = eps->work[0];
 88:   for (i=0;i<eps->nconv;i++) {
 89:     PetscCall(BVCopyVec(eps->W,i,w));
 90:     PetscCall(VecConjugate(w));
 91:     PetscCall(BVGetColumn(eps->W,i,&y));
 92:     PetscCall(STMatSolveTranspose(eps->st,w,y));
 93:     PetscCall(VecConjugate(y));
 94:     PetscCall(BVRestoreColumn(eps->W,i,&y));
 95:   }
 96:   PetscFunctionReturn(PETSC_SUCCESS);
 97: }

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

113:   PetscFunctionBegin;
114:   if (eps->ishermitian) {
115:     if (eps->isgeneralized && !eps->ispositive) PetscCall(EPSComputeVectors_Indefinite(eps));
116:     else PetscCall(EPSComputeVectors_Hermitian(eps));
117:     PetscFunctionReturn(PETSC_SUCCESS);
118:   }

120:   /* right eigenvectors */
121:   PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));

123:   /* V = V * Z */
124:   PetscCall(DSGetMat(eps->ds,DS_MAT_X,&Z));
125:   PetscCall(BVMultInPlace(eps->V,Z,0,eps->nconv));
126:   PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&Z));

128:   /* Purify eigenvectors */
129:   if (eps->purify) PetscCall(EPS_Purify(eps,eps->nconv));

131:   /* Fix eigenvectors if balancing was used */
132:   if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
133:     for (i=0;i<eps->nconv;i++) {
134:       PetscCall(BVGetColumn(eps->V,i,&z));
135:       PetscCall(VecPointwiseDivide(z,z,eps->D));
136:       PetscCall(BVRestoreColumn(eps->V,i,&z));
137:     }
138:   }

140:   /* normalize eigenvectors (when using purification or balancing) */
141:   if (eps->purify || (eps->balance!=EPS_BALANCE_NONE && eps->D)) PetscCall(BVNormalize(eps->V,eps->eigi));

143:   /* left eigenvectors */
144:   if (eps->twosided) {
145:     PetscCall(DSVectors(eps->ds,DS_MAT_Y,NULL,NULL));
146:     /* W = W * Z */
147:     PetscCall(DSGetMat(eps->ds,DS_MAT_Y,&Z));
148:     PetscCall(BVMultInPlace(eps->W,Z,0,eps->nconv));
149:     PetscCall(DSRestoreMat(eps->ds,DS_MAT_Y,&Z));
150:     /* Fix left eigenvectors if balancing was used */
151:     if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
152:       for (i=0;i<eps->nconv;i++) {
153:         PetscCall(BVGetColumn(eps->W,i,&z));
154:         PetscCall(VecPointwiseMult(z,z,eps->D));
155:         PetscCall(BVRestoreColumn(eps->W,i,&z));
156:       }
157:     }
158:     PetscCall(EPSComputeVectors_Twosided(eps));
159:     /* normalize */
160:     PetscCall(BVNormalize(eps->W,eps->eigi));
161: #if !defined(PETSC_USE_COMPLEX)
162:     for (i=0;i<eps->nconv-1;i++) {
163:       if (eps->eigi[i] != 0.0) {
164:         if (eps->eigi[i] > 0.0) PetscCall(BVScaleColumn(eps->W,i+1,-1.0));
165:         i++;
166:       }
167:     }
168: #endif
169:   }
170:   PetscFunctionReturn(PETSC_SUCCESS);
171: }

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

176:    Collective

178:    Input Parameters:
179: +  eps - eigensolver context
180: -  nw  - number of work vectors to allocate

182:    Developer Notes:
183:    This is SLEPC_EXTERN because it may be required by user plugin EPS
184:    implementations.

186:    Level: developer

188: .seealso: EPSSetUp()
189: @*/
190: PetscErrorCode EPSSetWorkVecs(EPS eps,PetscInt nw)
191: {
192:   Vec            t;

194:   PetscFunctionBegin;
197:   PetscCheck(nw>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"nw must be > 0: nw = %" PetscInt_FMT,nw);
198:   if (eps->nwork < nw) {
199:     PetscCall(VecDestroyVecs(eps->nwork,&eps->work));
200:     eps->nwork = nw;
201:     PetscCall(BVGetColumn(eps->V,0,&t));
202:     PetscCall(VecDuplicateVecs(t,nw,&eps->work));
203:     PetscCall(BVRestoreColumn(eps->V,0,&t));
204:   }
205:   PetscFunctionReturn(PETSC_SUCCESS);
206: }

208: /*
209:   EPSSetWhichEigenpairs_Default - Sets the default value for which,
210:   depending on the ST.
211:  */
212: PetscErrorCode EPSSetWhichEigenpairs_Default(EPS eps)
213: {
214:   PetscBool      target;

216:   PetscFunctionBegin;
217:   PetscCall(PetscObjectTypeCompareAny((PetscObject)eps->st,&target,STSINVERT,STCAYLEY,""));
218:   if (target) eps->which = EPS_TARGET_MAGNITUDE;
219:   else eps->which = EPS_LARGEST_MAGNITUDE;
220:   PetscFunctionReturn(PETSC_SUCCESS);
221: }

223: /*
224:   EPSConvergedRelative - Checks convergence relative to the eigenvalue.
225: */
226: PetscErrorCode EPSConvergedRelative(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
227: {
228:   PetscReal w;

230:   PetscFunctionBegin;
231:   w = SlepcAbsEigenvalue(eigr,eigi);
232:   *errest = res/w;
233:   PetscFunctionReturn(PETSC_SUCCESS);
234: }

236: /*
237:   EPSConvergedAbsolute - Checks convergence absolutely.
238: */
239: PetscErrorCode EPSConvergedAbsolute(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
240: {
241:   PetscFunctionBegin;
242:   *errest = res;
243:   PetscFunctionReturn(PETSC_SUCCESS);
244: }

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

254:   PetscFunctionBegin;
255:   w = SlepcAbsEigenvalue(eigr,eigi);
256:   *errest = res / (eps->nrma + w*eps->nrmb);
257:   PetscFunctionReturn(PETSC_SUCCESS);
258: }

260: /*@C
261:    EPSStoppingBasic - Default routine to determine whether the outer eigensolver
262:    iteration must be stopped.

264:    Collective

266:    Input Parameters:
267: +  eps    - eigensolver context obtained from EPSCreate()
268: .  its    - current number of iterations
269: .  max_it - maximum number of iterations
270: .  nconv  - number of currently converged eigenpairs
271: .  nev    - number of requested eigenpairs
272: -  ctx    - context (not used here)

274:    Output Parameter:
275: .  reason - result of the stopping test

277:    Notes:
278:    A positive value of reason indicates that the iteration has finished successfully
279:    (converged), and a negative value indicates an error condition (diverged). If
280:    the iteration needs to be continued, reason must be set to EPS_CONVERGED_ITERATING
281:    (zero).

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

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

288:    Level: advanced

290: .seealso: EPSSetStoppingTest(), EPSConvergedReason, EPSGetConvergedReason()
291: @*/
292: PetscErrorCode EPSStoppingBasic(EPS eps,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nev,EPSConvergedReason *reason,void *ctx)
293: {
294:   PetscFunctionBegin;
295:   *reason = EPS_CONVERGED_ITERATING;
296:   if (nconv >= nev) {
297:     PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: %" PetscInt_FMT " eigenpairs converged at iteration %" PetscInt_FMT "\n",nconv,its));
298:     *reason = EPS_CONVERGED_TOL;
299:   } else if (its >= max_it) {
300:     *reason = EPS_DIVERGED_ITS;
301:     PetscCall(PetscInfo(eps,"Linear eigensolver iteration reached maximum number of iterations (%" PetscInt_FMT ")\n",its));
302:   }
303:   PetscFunctionReturn(PETSC_SUCCESS);
304: }

306: /*
307:   EPSComputeRitzVector - Computes the current Ritz vector.

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

312:   Split case:
313:     x = V*Zr  y = V*Zi  (Zr and Zi have length nv)
314: */
315: PetscErrorCode EPSComputeRitzVector(EPS eps,PetscScalar *Zr,PetscScalar *Zi,BV V,Vec x,Vec y)
316: {
317:   PetscInt       l,k;
318:   PetscReal      norm;
319: #if !defined(PETSC_USE_COMPLEX)
320:   Vec            z;
321: #endif

323:   PetscFunctionBegin;
324:   /* compute eigenvector */
325:   PetscCall(BVGetActiveColumns(V,&l,&k));
326:   PetscCall(BVSetActiveColumns(V,0,k));
327:   PetscCall(BVMultVec(V,1.0,0.0,x,Zr));

329:   /* purify eigenvector if necessary */
330:   if (eps->purify) {
331:     PetscCall(STApply(eps->st,x,y));
332:     if (eps->ishermitian) PetscCall(BVNormVec(eps->V,y,NORM_2,&norm));
333:     else PetscCall(VecNorm(y,NORM_2,&norm));
334:     PetscCall(VecScale(y,1.0/norm));
335:     PetscCall(VecCopy(y,x));
336:   }
337:   /* fix eigenvector if balancing is used */
338:   if (!eps->ishermitian && eps->balance!=EPS_BALANCE_NONE && eps->D) PetscCall(VecPointwiseDivide(x,x,eps->D));
339: #if !defined(PETSC_USE_COMPLEX)
340:   /* compute imaginary part of eigenvector */
341:   if (Zi) {
342:     PetscCall(BVMultVec(V,1.0,0.0,y,Zi));
343:     if (eps->ispositive) {
344:       PetscCall(BVCreateVec(V,&z));
345:       PetscCall(STApply(eps->st,y,z));
346:       PetscCall(VecNorm(z,NORM_2,&norm));
347:       PetscCall(VecScale(z,1.0/norm));
348:       PetscCall(VecCopy(z,y));
349:       PetscCall(VecDestroy(&z));
350:     }
351:     if (eps->balance!=EPS_BALANCE_NONE && eps->D) PetscCall(VecPointwiseDivide(y,y,eps->D));
352:   } else
353: #endif
354:     PetscCall(VecSet(y,0.0));

356:   /* normalize eigenvectors (when using balancing) */
357:   if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
358: #if !defined(PETSC_USE_COMPLEX)
359:     if (Zi) PetscCall(VecNormalizeComplex(x,y,PETSC_TRUE,NULL));
360:     else
361: #endif
362:     PetscCall(VecNormalize(x,NULL));
363:   }
364:   PetscCall(BVSetActiveColumns(V,l,k));
365:   PetscFunctionReturn(PETSC_SUCCESS);
366: }

368: /*
369:   EPSBuildBalance_Krylov - uses a Krylov subspace method to compute the
370:   diagonal matrix to be applied for balancing in non-Hermitian problems.
371: */
372: PetscErrorCode EPSBuildBalance_Krylov(EPS eps)
373: {
374:   Vec               z,p,r;
375:   PetscInt          i,j;
376:   PetscReal         norma;
377:   PetscScalar       *pz,*pD;
378:   const PetscScalar *pr,*pp;
379:   PetscRandom       rand;

381:   PetscFunctionBegin;
382:   PetscCall(EPSSetWorkVecs(eps,3));
383:   PetscCall(BVGetRandomContext(eps->V,&rand));
384:   r = eps->work[0];
385:   p = eps->work[1];
386:   z = eps->work[2];
387:   PetscCall(VecSet(eps->D,1.0));

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

391:     /* Build a random vector of +-1's */
392:     PetscCall(VecSetRandom(z,rand));
393:     PetscCall(VecGetArray(z,&pz));
394:     for (i=0;i<eps->nloc;i++) {
395:       if (PetscRealPart(pz[i])<0.5) pz[i]=-1.0;
396:       else pz[i]=1.0;
397:     }
398:     PetscCall(VecRestoreArray(z,&pz));

400:     /* Compute p=DA(D\z) */
401:     PetscCall(VecPointwiseDivide(r,z,eps->D));
402:     PetscCall(STApply(eps->st,r,p));
403:     PetscCall(VecPointwiseMult(p,p,eps->D));
404:     if (eps->balance == EPS_BALANCE_TWOSIDE) {
405:       if (j==0) {
406:         /* Estimate the matrix inf-norm */
407:         PetscCall(VecAbs(p));
408:         PetscCall(VecMax(p,NULL,&norma));
409:       }
410:       /* Compute r=D\(A'Dz) */
411:       PetscCall(VecPointwiseMult(z,z,eps->D));
412:       PetscCall(STApplyHermitianTranspose(eps->st,z,r));
413:       PetscCall(VecPointwiseDivide(r,r,eps->D));
414:     }

416:     /* Adjust values of D */
417:     PetscCall(VecGetArrayRead(r,&pr));
418:     PetscCall(VecGetArrayRead(p,&pp));
419:     PetscCall(VecGetArray(eps->D,&pD));
420:     for (i=0;i<eps->nloc;i++) {
421:       if (eps->balance == EPS_BALANCE_TWOSIDE) {
422:         if (PetscAbsScalar(pp[i])>eps->balance_cutoff*norma && pr[i]!=0.0)
423:           pD[i] *= PetscSqrtReal(PetscAbsScalar(pr[i]/pp[i]));
424:       } else {
425:         if (pp[i]!=0.0) pD[i] /= PetscAbsScalar(pp[i]);
426:       }
427:     }
428:     PetscCall(VecRestoreArrayRead(r,&pr));
429:     PetscCall(VecRestoreArrayRead(p,&pp));
430:     PetscCall(VecRestoreArray(eps->D,&pD));
431:   }
432:   PetscFunctionReturn(PETSC_SUCCESS);
433: }