Actual source code: epssolve.c

slepc-main 2024-11-09
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:       }
172:       i++;
173:     }
174:   }
175: #endif

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

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

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

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

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

215:    Not Collective

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

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

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

230:    Level: intermediate

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

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

246:    Not Collective

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

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

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

257:    Level: beginner

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

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

275:    Not Collective

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

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

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

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

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

296:    Level: intermediate

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

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

314:    Collective

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

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

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

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

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

335:    Level: intermediate

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

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

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

371:    Collective

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

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

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

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

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

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

401:    Level: beginner

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

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

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

425:    Not Collective

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

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

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

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

444:    Level: beginner

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

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

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

485:    Collective

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

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

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

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

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

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

512:    Level: beginner

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

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

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

537:    Collective

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

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

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

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

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

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

563:    Level: intermediate

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

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

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

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

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

617:    Not Collective

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

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

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

631:    Level: advanced

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

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

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

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

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

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

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

724:    Collective

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

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

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

738:    Level: beginner

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

750:   PetscFunctionBegin;
754:   PetscAssertPointer(error,4);
755:   EPSCheckSolved(eps,1);

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

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

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

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

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

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

821:    Collective

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

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

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

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

847:   PetscFunctionBegin;

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

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

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

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

884:   PetscFunctionBegin;

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

891:   /* Orthonormalize the vector with respect to previous vectors */
892:   PetscCall(BVOrthogonalizeColumn(eps->W,i,NULL,&norm,&lindep));
893:   if (breakdown) *breakdown = lindep;
894:   else if (lindep || norm == 0.0) {
895:     PetscCheck(i,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Left initial vector is zero");
896:     SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Unable to generate more left start vectors");
897:   }
898:   PetscCall(BVScaleColumn(eps->W,i,1.0/norm));
899:   PetscFunctionReturn(PETSC_SUCCESS);
900: }