Actual source code: epssolve.c

slepc-3.23.2 2025-06-30
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:    EPS routines related to the solution process
 12: */

 14: #include <slepc/private/epsimpl.h>
 15: #include <slepc/private/bvimpl.h>
 16: #include <petscdraw.h>

 18: PetscErrorCode EPSComputeVectors(EPS eps)
 19: {
 20:   PetscFunctionBegin;
 21:   EPSCheckSolved(eps,1);
 22:   if (eps->state==EPS_STATE_SOLVED) PetscTryTypeMethod(eps,computevectors);
 23:   eps->state = EPS_STATE_EIGENVECTORS;
 24:   PetscFunctionReturn(PETSC_SUCCESS);
 25: }

 27: static PetscErrorCode EPSComputeValues(EPS eps)
 28: {
 29:   PetscBool      injective,iscomp,isfilter;
 30:   PetscInt       i,n,aux,nconv0;
 31:   Mat            A,B=NULL,G,Z;

 33:   PetscFunctionBegin;
 34:   switch (eps->categ) {
 35:     case EPS_CATEGORY_KRYLOV:
 36:     case EPS_CATEGORY_OTHER:
 37:       PetscCall(STIsInjective(eps->st,&injective));
 38:       if (injective) {
 39:         /* one-to-one mapping: backtransform eigenvalues */
 40:         PetscUseTypeMethod(eps,backtransform);
 41:       } else {
 42:         /* compute eigenvalues from Rayleigh quotient */
 43:         PetscCall(DSGetDimensions(eps->ds,&n,NULL,NULL,NULL));
 44:         if (!n) break;
 45:         PetscCall(EPSGetOperators(eps,&A,&B));
 46:         PetscCall(BVSetActiveColumns(eps->V,0,n));
 47:         PetscCall(DSGetCompact(eps->ds,&iscomp));
 48:         PetscCall(DSSetCompact(eps->ds,PETSC_FALSE));
 49:         PetscCall(DSGetMat(eps->ds,DS_MAT_A,&G));
 50:         PetscCall(BVMatProject(eps->V,A,eps->V,G));
 51:         PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&G));
 52:         if (B) {
 53:           PetscCall(DSGetMat(eps->ds,DS_MAT_B,&G));
 54:           PetscCall(BVMatProject(eps->V,B,eps->V,G));
 55:           PetscCall(DSRestoreMat(eps->ds,DS_MAT_B,&G));
 56:         }
 57:         PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
 58:         PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
 59:         PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
 60:         PetscCall(DSSetCompact(eps->ds,iscomp));
 61:         if (eps->ishermitian && (!eps->isgeneralized || eps->ispositive)) { /* V = V * Z */
 62:           PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));
 63:           PetscCall(DSGetMat(eps->ds,DS_MAT_X,&Z));
 64:           PetscCall(BVMultInPlace(eps->V,Z,0,n));
 65:           PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&Z));
 66:         }
 67:         /* in case of STFILTER discard computed eigenvalues that lie outside the wanted interval */
 68:         PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilter));
 69:         if (isfilter) {
 70:           nconv0 = eps->nconv;
 71:           for (i=0;i<eps->nconv;i++) {
 72:             if (PetscRealPart(eps->eigr[eps->perm[i]])<eps->inta || PetscRealPart(eps->eigr[eps->perm[i]])>eps->intb) {
 73:               eps->nconv--;
 74:               if (i<eps->nconv) { SlepcSwap(eps->perm[i],eps->perm[eps->nconv],aux); i--; }
 75:             }
 76:           }
 77:           if (nconv0>eps->nconv) PetscCall(PetscInfo(eps,"Discarded %" PetscInt_FMT " computed eigenvalues lying outside the interval\n",nconv0-eps->nconv));
 78:         }
 79:       }
 80:       break;
 81:     case EPS_CATEGORY_PRECOND:
 82:     case EPS_CATEGORY_CONTOUR:
 83:       /* eigenvalues already available as an output of the solver */
 84:       break;
 85:   }
 86:   PetscFunctionReturn(PETSC_SUCCESS);
 87: }

 89: /*@
 90:    EPSSolve - Solves the eigensystem.

 92:    Collective

 94:    Input Parameter:
 95: .  eps - eigensolver context obtained from EPSCreate()

 97:    Options Database Keys:
 98: +  -eps_view - print information about the solver used
 99: .  -eps_view_mat0 - view the first matrix (A)
100: .  -eps_view_mat1 - view the second matrix (B)
101: .  -eps_view_vectors - view the computed eigenvectors
102: .  -eps_view_values - view the computed eigenvalues
103: .  -eps_converged_reason - print reason for convergence, and number of iterations
104: .  -eps_error_absolute - print absolute errors of each eigenpair
105: .  -eps_error_relative - print relative errors of each eigenpair
106: -  -eps_error_backward - print backward errors of each eigenpair

108:    Notes:
109:    All the command-line options listed above admit an optional argument specifying
110:    the viewer type and options. For instance, use '-eps_view_mat0 binary:amatrix.bin'
111:    to save the A matrix to a binary file, '-eps_view_values draw' to draw the computed
112:    eigenvalues graphically, or '-eps_error_relative :myerr.m:ascii_matlab' to save
113:    the errors in a file that can be executed in Matlab.

115:    Level: beginner

117: .seealso: EPSCreate(), EPSSetUp(), EPSDestroy(), EPSSetTolerances()
118: @*/
119: PetscErrorCode EPSSolve(EPS eps)
120: {
121:   PetscInt       i;
122:   PetscBool      hasname;
123:   STMatMode      matmode;
124:   Mat            A,B;

126:   PetscFunctionBegin;
128:   if (eps->state>=EPS_STATE_SOLVED) PetscFunctionReturn(PETSC_SUCCESS);
129:   PetscCall(PetscLogEventBegin(EPS_Solve,eps,0,0,0));

131:   /* Call setup */
132:   PetscCall(EPSSetUp(eps));
133:   eps->nconv = 0;
134:   eps->its   = 0;
135:   for (i=0;i<eps->ncv;i++) {
136:     eps->eigr[i]   = 0.0;
137:     eps->eigi[i]   = 0.0;
138:     eps->errest[i] = 0.0;
139:     eps->perm[i]   = i;
140:   }
141:   PetscCall(EPSViewFromOptions(eps,NULL,"-eps_view_pre"));
142:   PetscCall(RGViewFromOptions(eps->rg,NULL,"-rg_view"));

144:   /* Call solver */
145:   PetscUseTypeMethod(eps,solve);
146:   PetscCheck(eps->reason,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
147:   eps->state = EPS_STATE_SOLVED;

149:   /* Only the first nconv columns contain useful information (except in CISS) */
150:   PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
151:   if (eps->twosided) PetscCall(BVSetActiveColumns(eps->W,0,eps->nconv));

153:   /* If inplace, purify eigenvectors before reverting operator */
154:   PetscCall(STGetMatMode(eps->st,&matmode));
155:   if (matmode == ST_MATMODE_INPLACE && eps->ispositive) PetscCall(EPSComputeVectors(eps));
156:   PetscCall(STPostSolve(eps->st));

158:   /* Map eigenvalues back to the original problem if appropriate */
159:   PetscCall(EPSComputeValues(eps));

161: #if !defined(PETSC_USE_COMPLEX)
162:   /* Reorder conjugate eigenvalues (positive imaginary first) */
163:   for (i=0;i<eps->nconv-1;i++) {
164:     if (eps->eigi[i] != 0) {
165:       if (eps->eigi[i] < 0) {
166:         eps->eigi[i] = -eps->eigi[i];
167:         eps->eigi[i+1] = -eps->eigi[i+1];
168:         /* the next correction only works with eigenvectors */
169:         PetscCall(EPSComputeVectors(eps));
170:         PetscCall(BVScaleColumn(eps->V,i+1,-1.0));
171:         if (eps->W) PetscCall(BVScaleColumn(eps->W,i+1,-1.0));
172:       }
173:       i++;
174:     }
175:   }
176: #endif

178:   /* Sort eigenvalues according to eps->which parameter */
179:   PetscCall(SlepcSortEigenvalues(eps->sc,eps->nconv,eps->eigr,eps->eigi,eps->perm));
180:   PetscCall(PetscLogEventEnd(EPS_Solve,eps,0,0,0));

182:   /* Various viewers */
183:   PetscCall(EPSViewFromOptions(eps,NULL,"-eps_view"));
184:   PetscCall(EPSConvergedReasonViewFromOptions(eps));
185:   PetscCall(EPSErrorViewFromOptions(eps));
186:   PetscCall(EPSValuesViewFromOptions(eps));
187:   PetscCall(EPSVectorsViewFromOptions(eps));

189:   PetscCall(PetscOptionsHasName(((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_view_mat0",&hasname));
190:   if (hasname) {
191:     PetscCall(EPSGetOperators(eps,&A,NULL));
192:     PetscCall(MatViewFromOptions(A,(PetscObject)eps,"-eps_view_mat0"));
193:   }
194:   if (eps->isgeneralized) {
195:     PetscCall(PetscOptionsHasName(((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_view_mat1",&hasname));
196:     if (hasname) {
197:       PetscCall(EPSGetOperators(eps,NULL,&B));
198:       PetscCall(MatViewFromOptions(B,(PetscObject)eps,"-eps_view_mat1"));
199:     }
200:   }

202:   /* Remove deflation and initial subspaces */
203:   if (eps->nds) {
204:     PetscCall(BVSetNumConstraints(eps->V,0));
205:     eps->nds = 0;
206:   }
207:   eps->nini = 0;
208:   PetscFunctionReturn(PETSC_SUCCESS);
209: }

211: /*@
212:    EPSGetIterationNumber - Gets the current iteration number. If the
213:    call to EPSSolve() is complete, then it returns the number of iterations
214:    carried out by the solution method.

216:    Not Collective

218:    Input Parameter:
219: .  eps - the eigensolver context

221:    Output Parameter:
222: .  its - number of iterations

224:    Note:
225:    During the i-th iteration this call returns i-1. If EPSSolve() is
226:    complete, then parameter "its" contains either the iteration number at
227:    which convergence was successfully reached, or failure was detected.
228:    Call EPSGetConvergedReason() to determine if the solver converged or
229:    failed and why.

231:    Level: intermediate

233: .seealso: EPSGetConvergedReason(), EPSSetTolerances()
234: @*/
235: PetscErrorCode EPSGetIterationNumber(EPS eps,PetscInt *its)
236: {
237:   PetscFunctionBegin;
239:   PetscAssertPointer(its,2);
240:   *its = eps->its;
241:   PetscFunctionReturn(PETSC_SUCCESS);
242: }

244: /*@
245:    EPSGetConverged - Gets the number of converged eigenpairs.

247:    Not Collective

249:    Input Parameter:
250: .  eps - the eigensolver context

252:    Output Parameter:
253: .  nconv - number of converged eigenpairs

255:    Note:
256:    This function should be called after EPSSolve() has finished.

258:    Level: beginner

260: .seealso: EPSSetDimensions(), EPSSolve(), EPSGetEigenpair()
261: @*/
262: PetscErrorCode EPSGetConverged(EPS eps,PetscInt *nconv)
263: {
264:   PetscFunctionBegin;
266:   PetscAssertPointer(nconv,2);
267:   EPSCheckSolved(eps,1);
268:   PetscCall(EPS_GetActualConverged(eps,nconv));
269:   PetscFunctionReturn(PETSC_SUCCESS);
270: }

272: /*@
273:    EPSGetConvergedReason - Gets the reason why the EPSSolve() iteration was
274:    stopped.

276:    Not Collective

278:    Input Parameter:
279: .  eps - the eigensolver context

281:    Output Parameter:
282: .  reason - negative value indicates diverged, positive value converged

284:    Options Database Key:
285: .  -eps_converged_reason - print the reason to a viewer

287:    Notes:
288:    Possible values for reason are
289: +  EPS_CONVERGED_TOL - converged up to tolerance
290: .  EPS_CONVERGED_USER - converged due to a user-defined condition
291: .  EPS_DIVERGED_ITS - required more than max_it iterations to reach convergence
292: .  EPS_DIVERGED_BREAKDOWN - generic breakdown in method
293: -  EPS_DIVERGED_SYMMETRY_LOST - pseudo-Lanczos was not able to keep symmetry

295:    Can only be called after the call to EPSSolve() is complete.

297:    Level: intermediate

299: .seealso: EPSSetTolerances(), EPSSolve(), EPSConvergedReason
300: @*/
301: PetscErrorCode EPSGetConvergedReason(EPS eps,EPSConvergedReason *reason)
302: {
303:   PetscFunctionBegin;
305:   PetscAssertPointer(reason,2);
306:   EPSCheckSolved(eps,1);
307:   *reason = eps->reason;
308:   PetscFunctionReturn(PETSC_SUCCESS);
309: }

311: /*@
312:    EPSGetInvariantSubspace - Gets an orthonormal basis of the computed invariant
313:    subspace.

315:    Collective

317:    Input Parameter:
318: .  eps - the eigensolver context

320:    Output Parameter:
321: .  v - an array of vectors

323:    Notes:
324:    This function should be called after EPSSolve() has finished.

326:    The user should provide in v an array of nconv vectors, where nconv is
327:    the value returned by EPSGetConverged().

329:    The first k vectors returned in v span an invariant subspace associated
330:    with the first k computed eigenvalues (note that this is not true if the
331:    k-th eigenvalue is complex and matrix A is real; in this case the first
332:    k+1 vectors should be used). An invariant subspace X of A satisfies Ax
333:    in X for all x in X (a similar definition applies for generalized
334:    eigenproblems).

336:    Level: intermediate

338: .seealso: EPSGetEigenpair(), EPSGetConverged(), EPSSolve()
339: @*/
340: PetscErrorCode EPSGetInvariantSubspace(EPS eps,Vec v[])
341: {
342:   PetscInt       i;
343:   BV             V=eps->V;
344:   Vec            w;

346:   PetscFunctionBegin;
348:   PetscAssertPointer(v,2);
350:   EPSCheckSolved(eps,1);
351:   PetscCheck(eps->ishermitian || eps->state!=EPS_STATE_EIGENVECTORS,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"EPSGetInvariantSubspace must be called before EPSGetEigenpair,EPSGetEigenvector or EPSComputeError");
352:   if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
353:     PetscCall(BVDuplicateResize(eps->V,eps->nconv,&V));
354:     PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
355:     PetscCall(BVCopy(eps->V,V));
356:     for (i=0;i<eps->nconv;i++) {
357:       PetscCall(BVGetColumn(V,i,&w));
358:       PetscCall(VecPointwiseDivide(w,w,eps->D));
359:       PetscCall(BVRestoreColumn(V,i,&w));
360:     }
361:     PetscCall(BVOrthogonalize(V,NULL));
362:   }
363:   for (i=0;i<eps->nconv;i++) PetscCall(BVCopyVec(V,i,v[i]));
364:   if (eps->balance!=EPS_BALANCE_NONE && eps->D) PetscCall(BVDestroy(&V));
365:   PetscFunctionReturn(PETSC_SUCCESS);
366: }

368: /*@
369:    EPSGetEigenpair - Gets the i-th solution of the eigenproblem as computed by
370:    EPSSolve(). The solution consists in both the eigenvalue and the eigenvector.

372:    Collective

374:    Input Parameters:
375: +  eps - eigensolver context
376: -  i   - index of the solution

378:    Output Parameters:
379: +  eigr - real part of eigenvalue
380: .  eigi - imaginary part of eigenvalue
381: .  Vr   - real part of eigenvector
382: -  Vi   - imaginary part of eigenvector

384:    Notes:
385:    It is allowed to pass NULL for Vr and Vi, if the eigenvector is not
386:    required. Otherwise, the caller must provide valid Vec objects, i.e.,
387:    they must be created by the calling program with e.g. MatCreateVecs().

389:    If the eigenvalue is real, then eigi and Vi are set to zero. If PETSc is
390:    configured with complex scalars the eigenvalue is stored
391:    directly in eigr (eigi is set to zero) and the eigenvector in Vr (Vi is
392:    set to zero). In both cases, the user can pass NULL in eigi and Vi.

394:    The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
395:    Eigenpairs are indexed according to the ordering criterion established
396:    with EPSSetWhichEigenpairs().

398:    The 2-norm of the eigenvector is one unless the problem is generalized
399:    Hermitian. In this case the eigenvector is normalized with respect to the
400:    norm defined by the B matrix.

402:    Level: beginner

404: .seealso: EPSGetEigenvalue(), EPSGetEigenvector(), EPSGetLeftEigenvector(), EPSSolve(),
405:           EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetInvariantSubspace()
406: @*/
407: PetscErrorCode EPSGetEigenpair(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
408: {
409:   PetscInt nconv;

411:   PetscFunctionBegin;
414:   EPSCheckSolved(eps,1);
415:   PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
416:   PetscCall(EPS_GetActualConverged(eps,&nconv));
417:   PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
418:   PetscCall(EPSGetEigenvalue(eps,i,eigr,eigi));
419:   if (Vr || Vi) PetscCall(EPSGetEigenvector(eps,i,Vr,Vi));
420:   PetscFunctionReturn(PETSC_SUCCESS);
421: }

423: /*@
424:    EPSGetEigenvalue - Gets the i-th eigenvalue as computed by EPSSolve().

426:    Not Collective

428:    Input Parameters:
429: +  eps - eigensolver context
430: -  i   - index of the solution

432:    Output Parameters:
433: +  eigr - real part of eigenvalue
434: -  eigi - imaginary part of eigenvalue

436:    Notes:
437:    If the eigenvalue is real, then eigi is set to zero. If PETSc is
438:    configured with complex scalars the eigenvalue is stored
439:    directly in eigr (eigi is set to zero).

441:    The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
442:    Eigenpairs are indexed according to the ordering criterion established
443:    with EPSSetWhichEigenpairs().

445:    Level: beginner

447: .seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetEigenpair()
448: @*/
449: PetscErrorCode EPSGetEigenvalue(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi)
450: {
451:   PetscInt k,nconv;

453:   PetscFunctionBegin;
455:   EPSCheckSolved(eps,1);
456:   PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
457:   PetscCall(EPS_GetActualConverged(eps,&nconv));
458:   PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
459:   if (nconv==eps->nconv) {
460:     k = eps->perm[i];
461: #if defined(PETSC_USE_COMPLEX)
462:     if (eigr) *eigr = eps->eigr[k];
463:     if (eigi) *eigi = 0;
464: #else
465:     if (eigr) *eigr = eps->eigr[k];
466:     if (eigi) *eigi = eps->eigi[k];
467: #endif
468:   } else {
469:     PetscCheck(eps->problem_type==EPS_BSE,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Problem type should be BSE");
470:     /* BSE problem, even index is +lambda, odd index is -lambda */
471:     k = eps->perm[i/2];
472: #if defined(PETSC_USE_COMPLEX)
473:     if (eigr) *eigr = (i%2)? -eps->eigr[k]: eps->eigr[k];
474:     if (eigi) *eigi = 0;
475: #else
476:     if (eigr) *eigr = (i%2)? -eps->eigr[k]: eps->eigr[k];
477:     if (eigi) *eigi = eps->eigi[k];
478: #endif
479:   }
480:   PetscFunctionReturn(PETSC_SUCCESS);
481: }

483: /*@
484:    EPSGetEigenvector - Gets the i-th right eigenvector as computed by EPSSolve().

486:    Collective

488:    Input Parameters:
489: +  eps - eigensolver context
490: -  i   - index of the solution

492:    Output Parameters:
493: +  Vr   - real part of eigenvector
494: -  Vi   - imaginary part of eigenvector

496:    Notes:
497:    The caller must provide valid Vec objects, i.e., they must be created
498:    by the calling program with e.g. MatCreateVecs().

500:    If the corresponding eigenvalue is real, then Vi is set to zero. If PETSc is
501:    configured with complex scalars the eigenvector is stored
502:    directly in Vr (Vi is set to zero). In any case, the user can pass NULL in Vr
503:    or Vi if one of them is not required.

505:    The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
506:    Eigenpairs are indexed according to the ordering criterion established
507:    with EPSSetWhichEigenpairs().

509:    The 2-norm of the eigenvector is one unless the problem is generalized
510:    Hermitian. In this case the eigenvector is normalized with respect to the
511:    norm defined by the B matrix.

513:    Level: beginner

515: .seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetEigenpair(), EPSGetLeftEigenvector()
516: @*/
517: PetscErrorCode EPSGetEigenvector(EPS eps,PetscInt i,Vec Vr,Vec Vi)
518: {
519:   PetscInt nconv;

521:   PetscFunctionBegin;
526:   EPSCheckSolved(eps,1);
527:   PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
528:   PetscCall(EPS_GetActualConverged(eps,&nconv));
529:   PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
530:   PetscCall(EPSComputeVectors(eps));
531:   PetscCall(EPS_GetEigenvector(eps,eps->V,i,Vr,Vi));
532:   PetscFunctionReturn(PETSC_SUCCESS);
533: }

535: /*@
536:    EPSGetLeftEigenvector - Gets the i-th left eigenvector as computed by EPSSolve().

538:    Collective

540:    Input Parameters:
541: +  eps - eigensolver context
542: -  i   - index of the solution

544:    Output Parameters:
545: +  Wr   - real part of left eigenvector
546: -  Wi   - imaginary part of left eigenvector

548:    Notes:
549:    The caller must provide valid Vec objects, i.e., they must be created
550:    by the calling program with e.g. MatCreateVecs().

552:    If the corresponding eigenvalue is real, then Wi is set to zero. If PETSc is
553:    configured with complex scalars the eigenvector is stored directly in Wr
554:    (Wi is set to zero). In any case, the user can pass NULL in Wr or Wi if
555:    one of them is not required.

557:    The index i should be a value between 0 and nconv-1 (see EPSGetConverged()).
558:    Eigensolutions are indexed according to the ordering criterion established
559:    with EPSSetWhichEigenpairs().

561:    Left eigenvectors are available only if the twosided flag was set, see
562:    EPSSetTwoSided().

564:    Level: intermediate

566: .seealso: EPSGetEigenvector(), EPSGetConverged(), EPSSetWhichEigenpairs(), EPSSetTwoSided()
567: @*/
568: PetscErrorCode EPSGetLeftEigenvector(EPS eps,PetscInt i,Vec Wr,Vec Wi)
569: {
570:   PetscInt    nconv;
571:   PetscBool   trivial;
572:   Mat         H;
573:   IS          is[2];
574:   Vec         v;

576:   PetscFunctionBegin;
581:   EPSCheckSolved(eps,1);
582:   PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
583:   PetscCall(EPS_GetActualConverged(eps,&nconv));
584:   PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");

586:   trivial = (eps->problem_type==EPS_HEP || eps->problem_type==EPS_GHEP || eps->problem_type==EPS_BSE)? PETSC_TRUE: PETSC_FALSE;
587:   if (!trivial) PetscCheck(eps->twosided,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must request left vectors with EPSSetTwoSided");

589:   PetscCall(EPSComputeVectors(eps));
590:   if (trivial) {
591:     PetscCall(EPS_GetEigenvector(eps,eps->V,i,Wr,Wi));
592:     if (eps->problem_type==EPS_BSE) {   /* change sign of bottom part of the vector */
593:       PetscCall(STGetMatrix(eps->st,0,&H));
594:       PetscCall(MatNestGetISs(H,is,NULL));
595:       if (Wr) {
596:         PetscCall(VecGetSubVector(Wr,is[1],&v));
597:         PetscCall(VecScale(v,-1.0));
598:         PetscCall(VecRestoreSubVector(Wr,is[1],&v));
599:       }
600: #if !defined(PETSC_USE_COMPLEX)
601:       if (Wi) {
602:         PetscCall(VecGetSubVector(Wi,is[1],&v));
603:         PetscCall(VecScale(v,-1.0));
604:         PetscCall(VecRestoreSubVector(Wi,is[1],&v));
605:       }
606: #endif
607:     }
608:   } else {
609:     PetscCall(EPS_GetEigenvector(eps,eps->W,i,Wr,Wi));
610:   }
611:   PetscFunctionReturn(PETSC_SUCCESS);
612: }

614: /*@
615:    EPSGetErrorEstimate - Returns the error estimate associated to the i-th
616:    computed eigenpair.

618:    Not Collective

620:    Input Parameters:
621: +  eps - eigensolver context
622: -  i   - index of eigenpair

624:    Output Parameter:
625: .  errest - the error estimate

627:    Notes:
628:    This is the error estimate used internally by the eigensolver. The actual
629:    error bound can be computed with EPSComputeError(). See also the users
630:    manual for details.

632:    Level: advanced

634: .seealso: EPSComputeError()
635: @*/
636: PetscErrorCode EPSGetErrorEstimate(EPS eps,PetscInt i,PetscReal *errest)
637: {
638:   PetscInt nconv;

640:   PetscFunctionBegin;
642:   PetscAssertPointer(errest,3);
643:   EPSCheckSolved(eps,1);
644:   PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
645:   PetscCall(EPS_GetActualConverged(eps,&nconv));
646:   PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
647:   if (nconv==eps->nconv) {
648:     *errest = eps->errest[eps->perm[i]];
649:   } else {
650:     PetscCheck(eps->problem_type==EPS_BSE,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Problem type should be BSE");
651:     /* BSE problem, even index is +lambda, odd index is -lambda, assume both have same error */
652:     *errest = eps->errest[eps->perm[i/2]];
653:   }
654:   PetscFunctionReturn(PETSC_SUCCESS);
655: }

657: /*
658:    EPSComputeResidualNorm_Private - Computes the norm of the residual vector
659:    associated with an eigenpair.

661:    Input Parameters:
662:      trans - whether A' must be used instead of A
663:      kr,ki - eigenvalue
664:      xr,xi - eigenvector
665:      z     - three work vectors (the second one not referenced in complex scalars)
666: */
667: PetscErrorCode EPSComputeResidualNorm_Private(EPS eps,PetscBool trans,PetscScalar kr,PetscScalar ki,Vec xr,Vec xi,Vec *z,PetscReal *norm)
668: {
669:   PetscInt       nmat;
670:   Mat            A,B;
671:   Vec            u,w;
672:   PetscScalar    alpha;
673: #if !defined(PETSC_USE_COMPLEX)
674:   Vec            v;
675:   PetscReal      ni,nr;
676: #endif
677:   PetscErrorCode (*matmult)(Mat,Vec,Vec) = trans? MatMultHermitianTranspose: MatMult;

679:   PetscFunctionBegin;
680:   u = z[0]; w = z[2];
681:   PetscCall(STGetNumMatrices(eps->st,&nmat));
682:   PetscCall(STGetMatrix(eps->st,0,&A));
683:   if (nmat>1) PetscCall(STGetMatrix(eps->st,1,&B));

685: #if !defined(PETSC_USE_COMPLEX)
686:   v = z[1];
687:   if (ki == 0 || PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
688: #endif
689:     PetscCall((*matmult)(A,xr,u));                          /* u=A*x */
690:     if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
691:       if (nmat>1) PetscCall((*matmult)(B,xr,w));
692:       else PetscCall(VecCopy(xr,w));                        /* w=B*x */
693:       alpha = trans? -PetscConj(kr): -kr;
694:       PetscCall(VecAXPY(u,alpha,w));                        /* u=A*x-k*B*x */
695:     }
696:     PetscCall(VecNorm(u,NORM_2,norm));
697: #if !defined(PETSC_USE_COMPLEX)
698:   } else {
699:     PetscCall((*matmult)(A,xr,u));                          /* u=A*xr */
700:     if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
701:       if (nmat>1) PetscCall((*matmult)(B,xr,v));
702:       else PetscCall(VecCopy(xr,v));                        /* v=B*xr */
703:       PetscCall(VecAXPY(u,-kr,v));                          /* u=A*xr-kr*B*xr */
704:       if (nmat>1) PetscCall((*matmult)(B,xi,w));
705:       else PetscCall(VecCopy(xi,w));                        /* w=B*xi */
706:       PetscCall(VecAXPY(u,trans?-ki:ki,w));                 /* u=A*xr-kr*B*xr+ki*B*xi */
707:     }
708:     PetscCall(VecNorm(u,NORM_2,&nr));
709:     PetscCall((*matmult)(A,xi,u));                          /* u=A*xi */
710:     if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
711:       PetscCall(VecAXPY(u,-kr,w));                          /* u=A*xi-kr*B*xi */
712:       PetscCall(VecAXPY(u,trans?ki:-ki,v));                 /* u=A*xi-kr*B*xi-ki*B*xr */
713:     }
714:     PetscCall(VecNorm(u,NORM_2,&ni));
715:     *norm = SlepcAbsEigenvalue(nr,ni);
716:   }
717: #endif
718:   PetscFunctionReturn(PETSC_SUCCESS);
719: }

721: /*@
722:    EPSComputeError - Computes the error (based on the residual norm) associated
723:    with the i-th computed eigenpair.

725:    Collective

727:    Input Parameters:
728: +  eps  - the eigensolver context
729: .  i    - the solution index
730: -  type - the type of error to compute

732:    Output Parameter:
733: .  error - the error

735:    Notes:
736:    The error can be computed in various ways, all of them based on the residual
737:    norm ||Ax-kBx||_2 where k is the eigenvalue and x is the eigenvector.

739:    Level: beginner

741: .seealso: EPSErrorType, EPSSolve(), EPSGetErrorEstimate()
742: @*/
743: PetscErrorCode EPSComputeError(EPS eps,PetscInt i,EPSErrorType type,PetscReal *error)
744: {
745:   Mat            A,B;
746:   Vec            xr,xi,w[3];
747:   PetscReal      t,vecnorm=1.0,errorl;
748:   PetscScalar    kr,ki;
749:   PetscBool      flg;

751:   PetscFunctionBegin;
755:   PetscAssertPointer(error,4);
756:   EPSCheckSolved(eps,1);

758:   /* allocate work vectors */
759: #if defined(PETSC_USE_COMPLEX)
760:   PetscCall(EPSSetWorkVecs(eps,3));
761:   xi   = NULL;
762:   w[1] = NULL;
763: #else
764:   PetscCall(EPSSetWorkVecs(eps,5));
765:   xi   = eps->work[3];
766:   w[1] = eps->work[4];
767: #endif
768:   xr   = eps->work[0];
769:   w[0] = eps->work[1];
770:   w[2] = eps->work[2];

772:   /* compute residual norm */
773:   PetscCall(EPSGetEigenpair(eps,i,&kr,&ki,xr,xi));
774:   PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_FALSE,kr,ki,xr,xi,w,error));

776:   /* compute 2-norm of eigenvector */
777:   if (eps->problem_type==EPS_GHEP) PetscCall(VecNorm(xr,NORM_2,&vecnorm));

779:   /* if two-sided, compute left residual norm and take the maximum */
780:   if (eps->twosided) {
781:     PetscCall(EPSGetLeftEigenvector(eps,i,xr,xi));
782:     PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_TRUE,kr,ki,xr,xi,w,&errorl));
783:     *error = PetscMax(*error,errorl);
784:   }

786:   /* compute error */
787:   switch (type) {
788:     case EPS_ERROR_ABSOLUTE:
789:       break;
790:     case EPS_ERROR_RELATIVE:
791:       *error /= SlepcAbsEigenvalue(kr,ki)*vecnorm;
792:       break;
793:     case EPS_ERROR_BACKWARD:
794:       /* initialization of matrix norms */
795:       if (!eps->nrma) {
796:         PetscCall(STGetMatrix(eps->st,0,&A));
797:         PetscCall(MatHasOperation(A,MATOP_NORM,&flg));
798:         PetscCheck(flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
799:         PetscCall(MatNorm(A,NORM_INFINITY,&eps->nrma));
800:       }
801:       if (eps->isgeneralized) {
802:         if (!eps->nrmb) {
803:           PetscCall(STGetMatrix(eps->st,1,&B));
804:           PetscCall(MatHasOperation(B,MATOP_NORM,&flg));
805:           PetscCheck(flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
806:           PetscCall(MatNorm(B,NORM_INFINITY,&eps->nrmb));
807:         }
808:       } else eps->nrmb = 1.0;
809:       t = SlepcAbsEigenvalue(kr,ki);
810:       *error /= (eps->nrma+t*eps->nrmb)*vecnorm;
811:       break;
812:     default:
813:       SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
814:   }
815:   PetscFunctionReturn(PETSC_SUCCESS);
816: }

818: /*
819:    EPSGetStartVector - Generate a suitable vector to be used as the starting vector
820:    for the recurrence that builds the right subspace.

822:    Collective

824:    Input Parameters:
825: +  eps - the eigensolver context
826: -  i   - iteration number

828:    Output Parameters:
829: .  breakdown - flag indicating that a breakdown has occurred

831:    Notes:
832:    The start vector is computed from another vector: for the first step (i=0),
833:    the first initial vector is used (see EPSSetInitialSpace()); otherwise a random
834:    vector is created. Then this vector is forced to be in the range of OP (only
835:    for generalized definite problems) and orthonormalized with respect to all
836:    V-vectors up to i-1. The resulting vector is placed in V[i].

838:    The flag breakdown is set to true if either i=0 and the vector belongs to the
839:    deflation space, or i>0 and the vector is linearly dependent with respect
840:    to the V-vectors.
841: */
842: PetscErrorCode EPSGetStartVector(EPS eps,PetscInt i,PetscBool *breakdown)
843: {
844:   PetscReal      norm;
845:   PetscBool      lindep;
846:   Vec            w,z;

848:   PetscFunctionBegin;

852:   /* For the first step, use the first initial vector, otherwise a random one */
853:   if (i>0 || eps->nini==0) PetscCall(BVSetRandomColumn(eps->V,i));

855:   /* Force the vector to be in the range of OP for generalized problems with B-inner product */
856:   if (eps->ispositive || (eps->isgeneralized && eps->ishermitian)) {
857:     PetscCall(BVCreateVec(eps->V,&w));
858:     PetscCall(BVCopyVec(eps->V,i,w));
859:     PetscCall(BVGetColumn(eps->V,i,&z));
860:     PetscCall(STApply(eps->st,w,z));
861:     PetscCall(BVRestoreColumn(eps->V,i,&z));
862:     PetscCall(VecDestroy(&w));
863:   }

865:   /* Orthonormalize the vector with respect to previous vectors */
866:   PetscCall(BVOrthogonalizeColumn(eps->V,i,NULL,&norm,&lindep));
867:   if (breakdown) *breakdown = lindep;
868:   else if (lindep || norm == 0.0) {
869:     PetscCheck(i,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Initial vector is zero or belongs to the deflation space");
870:     PetscCheck(!i,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Unable to generate more start vectors");
871:   }
872:   PetscCall(BVScaleColumn(eps->V,i,1.0/norm));
873:   PetscFunctionReturn(PETSC_SUCCESS);
874: }

876: /*
877:    EPSGetLeftStartVector - Generate a suitable vector to be used as the left starting
878:    vector for the recurrence that builds the left subspace. See EPSGetStartVector().
879: */
880: PetscErrorCode EPSGetLeftStartVector(EPS eps,PetscInt i,PetscBool *breakdown)
881: {
882:   PetscReal      norm;
883:   PetscBool      lindep;
884:   Vec            w,z;

886:   PetscFunctionBegin;

890:   /* For the first step, use the first initial vector, otherwise a random one */
891:   if (i>0 || eps->ninil==0) PetscCall(BVSetRandomColumn(eps->W,i));

893:   /* Force the vector to be in the range of OP' for generalized problems with B-inner product */
894:   if (eps->ispositive || (eps->isgeneralized && eps->ishermitian)) {
895:     PetscCall(BVCreateVec(eps->W,&w));
896:     PetscCall(BVCopyVec(eps->W,i,w));
897:     PetscCall(BVGetColumn(eps->W,i,&z));
898:     PetscCall(STApplyHermitianTranspose(eps->st,w,z));
899:     PetscCall(BVRestoreColumn(eps->W,i,&z));
900:     PetscCall(VecDestroy(&w));
901:   }

903:   /* Orthonormalize the vector with respect to previous vectors */
904:   PetscCall(BVOrthogonalizeColumn(eps->W,i,NULL,&norm,&lindep));
905:   if (breakdown) *breakdown = lindep;
906:   else if (lindep || norm == 0.0) {
907:     PetscCheck(i,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Left initial vector is zero");
908:     PetscCheck(!i,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Unable to generate more left start vectors");
909:   }
910:   PetscCall(BVScaleColumn(eps->W,i,1.0/norm));
911:   PetscFunctionReturn(PETSC_SUCCESS);
912: }