Actual source code: epssolve.c

  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));

134:   /* Safeguard for matrices of size 0 */
135:   if (eps->n == 0) {
136:     eps->nconv  = 0;
137:     eps->reason = EPS_CONVERGED_TOL;
138:     eps->state  = EPS_STATE_SOLVED;
139:     PetscFunctionReturn(PETSC_SUCCESS);
140:   }

142:   eps->nconv = 0;
143:   eps->its   = 0;
144:   for (i=0;i<eps->ncv;i++) {
145:     eps->eigr[i]   = 0.0;
146:     eps->eigi[i]   = 0.0;
147:     eps->errest[i] = 0.0;
148:     eps->perm[i]   = i;
149:   }
150:   PetscCall(EPSViewFromOptions(eps,NULL,"-eps_view_pre"));
151:   PetscCall(RGViewFromOptions(eps->rg,NULL,"-rg_view"));

153:   /* Call solver */
154:   PetscUseTypeMethod(eps,solve);
155:   PetscCheck(eps->reason,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
156:   eps->state = EPS_STATE_SOLVED;

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

162:   /* If inplace, purify eigenvectors before reverting operator */
163:   PetscCall(STGetMatMode(eps->st,&matmode));
164:   if (matmode == ST_MATMODE_INPLACE && eps->ispositive) PetscCall(EPSComputeVectors(eps));
165:   PetscCall(STPostSolve(eps->st));

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

170: #if !defined(PETSC_USE_COMPLEX)
171:   /* Reorder conjugate eigenvalues (positive imaginary first) */
172:   for (i=0;i<eps->nconv-1;i++) {
173:     if (eps->eigi[i] != 0 && (eps->problem_type!=EPS_HAMILT || eps->eigr[i]!=0)) {
174:       /* conjugate eigenvalues */
175:       if (eps->eigi[i] < 0) {
176:         eps->eigi[i] = -eps->eigi[i];
177:         eps->eigi[i+1] = -eps->eigi[i+1];
178:         /* the next correction only works with eigenvectors */
179:         PetscCall(EPSComputeVectors(eps));
180:         PetscCall(BVScaleColumn(eps->V,i+1,-1.0));
181:         if (eps->W) PetscCall(BVScaleColumn(eps->W,i+1,-1.0));
182:       }
183:       i++;
184:     }
185:   }
186: #endif

188:   /* Sort eigenvalues according to eps->which parameter */
189:   if (eps->problem_type==EPS_HAMILT) PetscCall(SlepcSortEigenvaluesSpecial(eps->sc,eps->nconv,eps->eigr,eps->eigi,eps->perm));
190:   else PetscCall(SlepcSortEigenvalues(eps->sc,eps->nconv,eps->eigr,eps->eigi,eps->perm));
191:   PetscCall(PetscLogEventEnd(EPS_Solve,eps,0,0,0));

193:   /* Various viewers */
194:   PetscCall(EPSViewFromOptions(eps,NULL,"-eps_view"));
195:   PetscCall(EPSConvergedReasonViewFromOptions(eps));
196:   PetscCall(EPSErrorViewFromOptions(eps));
197:   PetscCall(EPSValuesViewFromOptions(eps));
198:   PetscCall(EPSVectorsViewFromOptions(eps));

200:   PetscCall(PetscOptionsHasName(((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_view_mat0",&hasname));
201:   if (hasname) {
202:     PetscCall(EPSGetOperators(eps,&A,NULL));
203:     PetscCall(MatViewFromOptions(A,(PetscObject)eps,"-eps_view_mat0"));
204:   }
205:   if (eps->isgeneralized) {
206:     PetscCall(PetscOptionsHasName(((PetscObject)eps)->options,((PetscObject)eps)->prefix,"-eps_view_mat1",&hasname));
207:     if (hasname) {
208:       PetscCall(EPSGetOperators(eps,NULL,&B));
209:       PetscCall(MatViewFromOptions(B,(PetscObject)eps,"-eps_view_mat1"));
210:     }
211:   }

213:   /* Remove deflation and initial subspaces */
214:   if (eps->nds) {
215:     PetscCall(BVSetNumConstraints(eps->V,0));
216:     eps->nds = 0;
217:   }
218:   eps->nini = 0;
219:   PetscFunctionReturn(PETSC_SUCCESS);
220: }

222: /*@
223:    EPSGetIterationNumber - Gets the current iteration number. If the
224:    call to EPSSolve() is complete, then it returns the number of iterations
225:    carried out by the solution method.

227:    Not Collective

229:    Input Parameter:
230: .  eps - the eigensolver context

232:    Output Parameter:
233: .  its - number of iterations

235:    Note:
236:    During the i-th iteration this call returns i-1. If EPSSolve() is
237:    complete, then parameter "its" contains either the iteration number at
238:    which convergence was successfully reached, or failure was detected.
239:    Call EPSGetConvergedReason() to determine if the solver converged or
240:    failed and why.

242:    Level: intermediate

244: .seealso: `EPSGetConvergedReason()`, `EPSSetTolerances()`
245: @*/
246: PetscErrorCode EPSGetIterationNumber(EPS eps,PetscInt *its)
247: {
248:   PetscFunctionBegin;
250:   PetscAssertPointer(its,2);
251:   *its = eps->its;
252:   PetscFunctionReturn(PETSC_SUCCESS);
253: }

255: /*@
256:    EPSGetConverged - Gets the number of converged eigenpairs.

258:    Not Collective

260:    Input Parameter:
261: .  eps - the eigensolver context

263:    Output Parameter:
264: .  nconv - number of converged eigenpairs

266:    Note:
267:    This function should be called after EPSSolve() has finished.

269:    Level: beginner

271: .seealso: `EPSSetDimensions()`, `EPSSolve()`, `EPSGetEigenpair()`
272: @*/
273: PetscErrorCode EPSGetConverged(EPS eps,PetscInt *nconv)
274: {
275:   PetscFunctionBegin;
277:   PetscAssertPointer(nconv,2);
278:   EPSCheckSolved(eps,1);
279:   PetscCall(EPS_GetActualConverged(eps,nconv));
280:   PetscFunctionReturn(PETSC_SUCCESS);
281: }

283: /*@
284:    EPSGetConvergedReason - Gets the reason why the EPSSolve() iteration was
285:    stopped.

287:    Not Collective

289:    Input Parameter:
290: .  eps - the eigensolver context

292:    Output Parameter:
293: .  reason - negative value indicates diverged, positive value converged

295:    Options Database Key:
296: .  -eps_converged_reason - print the reason to a viewer

298:    Notes:
299:    Possible values for reason are
300: +  EPS_CONVERGED_TOL - converged up to tolerance
301: .  EPS_CONVERGED_USER - converged due to a user-defined condition
302: .  EPS_DIVERGED_ITS - required more than max_it iterations to reach convergence
303: .  EPS_DIVERGED_BREAKDOWN - generic breakdown in method
304: -  EPS_DIVERGED_SYMMETRY_LOST - pseudo-Lanczos was not able to keep symmetry

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

308:    Level: intermediate

310: .seealso: `EPSSetTolerances()`, `EPSSolve()`, `EPSConvergedReason`
311: @*/
312: PetscErrorCode EPSGetConvergedReason(EPS eps,EPSConvergedReason *reason)
313: {
314:   PetscFunctionBegin;
316:   PetscAssertPointer(reason,2);
317:   EPSCheckSolved(eps,1);
318:   *reason = eps->reason;
319:   PetscFunctionReturn(PETSC_SUCCESS);
320: }

322: /*@
323:    EPSGetInvariantSubspace - Gets an orthonormal basis of the computed invariant
324:    subspace.

326:    Collective

328:    Input Parameter:
329: .  eps - the eigensolver context

331:    Output Parameter:
332: .  v - an array of vectors

334:    Notes:
335:    This function should be called after EPSSolve() has finished.

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

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

347:    Level: intermediate

349: .seealso: `EPSGetEigenpair()`, `EPSGetConverged()`, `EPSSolve()`
350: @*/
351: PetscErrorCode EPSGetInvariantSubspace(EPS eps,Vec v[])
352: {
353:   PetscInt       i;
354:   BV             V=eps->V;
355:   Vec            w;

357:   PetscFunctionBegin;
359:   PetscAssertPointer(v,2);
361:   EPSCheckSolved(eps,1);
362:   PetscCheck(eps->ishermitian || eps->state!=EPS_STATE_EIGENVECTORS,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"EPSGetInvariantSubspace must be called before EPSGetEigenpair,EPSGetEigenvector or EPSComputeError");
363:   if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
364:     PetscCall(BVDuplicateResize(eps->V,eps->nconv,&V));
365:     PetscCall(BVSetActiveColumns(eps->V,0,eps->nconv));
366:     PetscCall(BVCopy(eps->V,V));
367:     for (i=0;i<eps->nconv;i++) {
368:       PetscCall(BVGetColumn(V,i,&w));
369:       PetscCall(VecPointwiseDivide(w,w,eps->D));
370:       PetscCall(BVRestoreColumn(V,i,&w));
371:     }
372:     PetscCall(BVOrthogonalize(V,NULL));
373:   }
374:   for (i=0;i<eps->nconv;i++) PetscCall(BVCopyVec(V,i,v[i]));
375:   if (eps->balance!=EPS_BALANCE_NONE && eps->D) PetscCall(BVDestroy(&V));
376:   PetscFunctionReturn(PETSC_SUCCESS);
377: }

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

383:    Collective

385:    Input Parameters:
386: +  eps - eigensolver context
387: -  i   - index of the solution

389:    Output Parameters:
390: +  eigr - real part of eigenvalue
391: .  eigi - imaginary part of eigenvalue
392: .  Vr   - real part of eigenvector
393: -  Vi   - imaginary part of eigenvector

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

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

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

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

413:    Level: beginner

415: .seealso: `EPSGetEigenvalue()`, `EPSGetEigenvector()`, `EPSGetLeftEigenvector()`, `EPSSolve()`,
416:           `EPSGetConverged()`, `EPSSetWhichEigenpairs()`, `EPSGetInvariantSubspace()`
417: @*/
418: PetscErrorCode EPSGetEigenpair(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
419: {
420:   PetscInt nconv;

422:   PetscFunctionBegin;
425:   EPSCheckSolved(eps,1);
426:   PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
427:   PetscCall(EPS_GetActualConverged(eps,&nconv));
428:   PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
429:   PetscCall(EPSGetEigenvalue(eps,i,eigr,eigi));
430:   if (Vr || Vi) PetscCall(EPSGetEigenvector(eps,i,Vr,Vi));
431:   PetscFunctionReturn(PETSC_SUCCESS);
432: }

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

437:    Not Collective

439:    Input Parameters:
440: +  eps - eigensolver context
441: -  i   - index of the solution

443:    Output Parameters:
444: +  eigr - real part of eigenvalue
445: -  eigi - imaginary part of eigenvalue

447:    Notes:
448:    If the eigenvalue is real, then eigi is set to zero. If PETSc is
449:    configured with complex scalars the eigenvalue is stored
450:    directly in eigr (eigi is set to zero).

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

456:    Level: beginner

458: .seealso: `EPSSolve()`, `EPSGetConverged()`, `EPSSetWhichEigenpairs()`, `EPSGetEigenpair()`
459: @*/
460: PetscErrorCode EPSGetEigenvalue(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi)
461: {
462:   PetscInt  k,nconv;
463: #if !defined(PETSC_USE_COMPLEX)
464:   PetscInt  k2, iquad;
465: #endif

467:   PetscFunctionBegin;
469:   EPSCheckSolved(eps,1);
470:   PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
471:   PetscCall(EPS_GetActualConverged(eps,&nconv));
472:   PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
473:   if (nconv==eps->nconv) {
474:     k = eps->perm[i];
475: #if defined(PETSC_USE_COMPLEX)
476:     if (eigr) *eigr = eps->eigr[k];
477:     if (eigi) *eigi = 0;
478: #else
479:     if (eigr) *eigr = eps->eigr[k];
480:     if (eigi) *eigi = eps->eigi[k];
481: #endif
482:   } else {
483:     PetscCheck(eps->problem_type==EPS_BSE || eps->problem_type==EPS_HAMILT,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Problem type should be BSE or Hamiltonian");
484:     if (eps->problem_type==EPS_BSE) {
485:       /* BSE problem, even index is +lambda, odd index is -lambda */
486:       k = eps->perm[i/2];
487: #if defined(PETSC_USE_COMPLEX)
488:       if (eigr) *eigr = (i%2)? -eps->eigr[k]: eps->eigr[k];
489:       if (eigi) *eigi = 0;
490: #else
491:       if (eigr) *eigr = (i%2)? -eps->eigr[k]: eps->eigr[k];
492:       if (eigi) *eigi = eps->eigi[k];
493: #endif
494:     } else if (eps->problem_type==EPS_HAMILT) {
495:       /* Hamiltonian eigenproblem */
496:       k = eps->perm[i/2];
497: #if defined(PETSC_USE_COMPLEX)
498:       if (eigr) *eigr = (i%2)? -eps->eigr[k]: eps->eigr[k];
499:       if (eigi) *eigi = 0;
500: #else
501:       if (eps->eigi[k]==0.0) { /* real eigenvalue */
502:         if (eigr) *eigr = (i%2)? -eps->eigr[k]: eps->eigr[k];
503:         if (eigi) *eigi = 0.0;
504:       } else if (eps->eigr[k]==0.0) { /* purely imaginary eigenvalue */
505:         if (eigr) *eigr = 0.0;
506:         if (eigi) *eigi = (i%2)? -eps->eigi[k]: eps->eigi[k];
507:       } else { /* quadruple eigenvalue (-conj(lambda),-lambda,lambda,conj(lambda)) */
508:         iquad = i%2;  /* index within the 4 values */
509:         if (i>1) {
510:           k2 = eps->perm[(i-2)/2];
511:           if (eps->eigr[k]==eps->eigr[k2] && eps->eigi[k]==-eps->eigi[k2]) iquad += 2;
512:         }
513:         if (eigr) *eigr = (iquad<2)? -eps->eigr[k]: eps->eigr[k];
514:         if (eigi) *eigi = (iquad%3)? -eps->eigi[k]: eps->eigi[k];
515:       }
516: #endif
517:     }
518:   }
519:   PetscFunctionReturn(PETSC_SUCCESS);
520: }

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

525:    Collective

527:    Input Parameters:
528: +  eps - eigensolver context
529: -  i   - index of the solution

531:    Output Parameters:
532: +  Vr   - real part of eigenvector
533: -  Vi   - imaginary part of eigenvector

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

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

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

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

552:    Level: beginner

554: .seealso: `EPSSolve()`, `EPSGetConverged()`, `EPSSetWhichEigenpairs()`, `EPSGetEigenpair()`, `EPSGetLeftEigenvector()`
555: @*/
556: PetscErrorCode EPSGetEigenvector(EPS eps,PetscInt i,Vec Vr,Vec Vi)
557: {
558:   PetscInt nconv;

560:   PetscFunctionBegin;
565:   EPSCheckSolved(eps,1);
566:   PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
567:   PetscCall(EPS_GetActualConverged(eps,&nconv));
568:   PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
569:   PetscCall(EPSComputeVectors(eps));
570:   PetscCall(EPS_GetEigenvector(eps,eps->V,i,Vr,Vi));
571:   PetscFunctionReturn(PETSC_SUCCESS);
572: }

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

577:    Collective

579:    Input Parameters:
580: +  eps - eigensolver context
581: -  i   - index of the solution

583:    Output Parameters:
584: +  Wr   - real part of left eigenvector
585: -  Wi   - imaginary part of left eigenvector

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

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

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

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

603:    Level: intermediate

605: .seealso: `EPSGetEigenvector()`, `EPSGetConverged()`, `EPSSetWhichEigenpairs()`, `EPSSetTwoSided()`
606: @*/
607: PetscErrorCode EPSGetLeftEigenvector(EPS eps,PetscInt i,Vec Wr,Vec Wi)
608: {
609:   PetscInt    nconv;
610:   PetscBool   trivial;
611:   Mat         H;
612:   IS          is[2];
613:   Vec         v;

615:   PetscFunctionBegin;
620:   EPSCheckSolved(eps,1);
621:   PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
622:   PetscCall(EPS_GetActualConverged(eps,&nconv));
623:   PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");

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

628:   PetscCall(EPSComputeVectors(eps));
629:   if (trivial) {
630:     PetscCall(EPS_GetEigenvector(eps,eps->V,i,Wr,Wi));
631:     if (eps->problem_type==EPS_BSE) {   /* change sign of bottom part of the vector */
632:       PetscCall(STGetMatrix(eps->st,0,&H));
633:       PetscCall(MatNestGetISs(H,is,NULL));
634:       if (Wr) {
635:         PetscCall(VecGetSubVector(Wr,is[1],&v));
636:         PetscCall(VecScale(v,-1.0));
637:         PetscCall(VecRestoreSubVector(Wr,is[1],&v));
638:       }
639: #if !defined(PETSC_USE_COMPLEX)
640:       if (Wi) {
641:         PetscCall(VecGetSubVector(Wi,is[1],&v));
642:         PetscCall(VecScale(v,-1.0));
643:         PetscCall(VecRestoreSubVector(Wi,is[1],&v));
644:       }
645: #endif
646:     }
647:   } else {
648:     PetscCall(EPS_GetEigenvector(eps,eps->W,i,Wr,Wi));
649:   }
650:   PetscFunctionReturn(PETSC_SUCCESS);
651: }

653: /*@
654:    EPSGetErrorEstimate - Returns the error estimate associated to the i-th
655:    computed eigenpair.

657:    Not Collective

659:    Input Parameters:
660: +  eps - eigensolver context
661: -  i   - index of eigenpair

663:    Output Parameter:
664: .  errest - the error estimate

666:    Notes:
667:    This is the error estimate used internally by the eigensolver. The actual
668:    error bound can be computed with EPSComputeError(). See also the users
669:    manual for details.

671:    Level: advanced

673: .seealso: `EPSComputeError()`
674: @*/
675: PetscErrorCode EPSGetErrorEstimate(EPS eps,PetscInt i,PetscReal *errest)
676: {
677:   PetscInt nconv;

679:   PetscFunctionBegin;
681:   PetscAssertPointer(errest,3);
682:   EPSCheckSolved(eps,1);
683:   PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
684:   PetscCall(EPS_GetActualConverged(eps,&nconv));
685:   PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
686:   if (nconv==eps->nconv) {
687:     *errest = eps->errest[eps->perm[i]];
688:   } else {
689:     PetscCheck(eps->problem_type==EPS_BSE,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Problem type should be BSE");
690:     /* BSE problem, even index is +lambda, odd index is -lambda, assume both have same error */
691:     *errest = eps->errest[eps->perm[i/2]];
692:   }
693:   PetscFunctionReturn(PETSC_SUCCESS);
694: }

696: /*
697:    EPSComputeResidualNorm_Private - Computes the norm of the residual vector
698:    associated with an eigenpair.

700:    Input Parameters:
701:      trans - whether A' must be used instead of A
702:      kr,ki - eigenvalue
703:      xr,xi - eigenvector
704:      z     - three work vectors (the second one not referenced in complex scalars)
705: */
706: PetscErrorCode EPSComputeResidualNorm_Private(EPS eps,PetscBool trans,PetscScalar kr,PetscScalar ki,Vec xr,Vec xi,Vec *z,PetscReal *norm)
707: {
708:   PetscInt       nmat;
709:   Mat            A,B;
710:   Vec            u,w;
711:   PetscScalar    alpha;
712: #if !defined(PETSC_USE_COMPLEX)
713:   Vec            v;
714:   PetscReal      ni,nr;
715: #endif
716:   PetscErrorCode (*matmult)(Mat,Vec,Vec) = trans? MatMultHermitianTranspose: MatMult;

718:   PetscFunctionBegin;
719:   u = z[0]; w = z[2];
720:   PetscCall(STGetNumMatrices(eps->st,&nmat));
721:   PetscCall(STGetMatrix(eps->st,0,&A));
722:   if (nmat>1) PetscCall(STGetMatrix(eps->st,1,&B));

724: #if !defined(PETSC_USE_COMPLEX)
725:   v = z[1];
726:   if (ki == 0 || PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
727: #endif
728:     PetscCall((*matmult)(A,xr,u));                          /* u=A*x */
729:     if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
730:       if (nmat>1) PetscCall((*matmult)(B,xr,w));
731:       else PetscCall(VecCopy(xr,w));                        /* w=B*x */
732:       alpha = trans? -PetscConj(kr): -kr;
733:       PetscCall(VecAXPY(u,alpha,w));                        /* u=A*x-k*B*x */
734:     }
735:     PetscCall(VecNorm(u,NORM_2,norm));
736: #if !defined(PETSC_USE_COMPLEX)
737:   } else {
738:     PetscCall((*matmult)(A,xr,u));                          /* u=A*xr */
739:     if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
740:       if (nmat>1) PetscCall((*matmult)(B,xr,v));
741:       else PetscCall(VecCopy(xr,v));                        /* v=B*xr */
742:       PetscCall(VecAXPY(u,-kr,v));                          /* u=A*xr-kr*B*xr */
743:       if (nmat>1) PetscCall((*matmult)(B,xi,w));
744:       else PetscCall(VecCopy(xi,w));                        /* w=B*xi */
745:       PetscCall(VecAXPY(u,trans?-ki:ki,w));                 /* u=A*xr-kr*B*xr+ki*B*xi */
746:     }
747:     PetscCall(VecNorm(u,NORM_2,&nr));
748:     PetscCall((*matmult)(A,xi,u));                          /* u=A*xi */
749:     if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
750:       PetscCall(VecAXPY(u,-kr,w));                          /* u=A*xi-kr*B*xi */
751:       PetscCall(VecAXPY(u,trans?ki:-ki,v));                 /* u=A*xi-kr*B*xi-ki*B*xr */
752:     }
753:     PetscCall(VecNorm(u,NORM_2,&ni));
754:     *norm = SlepcAbsEigenvalue(nr,ni);
755:   }
756: #endif
757:   PetscFunctionReturn(PETSC_SUCCESS);
758: }

760: /*@
761:    EPSComputeError - Computes the error (based on the residual norm) associated
762:    with the i-th computed eigenpair.

764:    Collective

766:    Input Parameters:
767: +  eps  - the eigensolver context
768: .  i    - the solution index
769: -  type - the type of error to compute

771:    Output Parameter:
772: .  error - the error

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

778:    Level: beginner

780: .seealso: `EPSErrorType`, `EPSSolve()`, `EPSGetErrorEstimate()`
781: @*/
782: PetscErrorCode EPSComputeError(EPS eps,PetscInt i,EPSErrorType type,PetscReal *error)
783: {
784:   Mat            A,B;
785:   Vec            xr,xi,w[3];
786:   PetscReal      t,vecnorm=1.0,errorl;
787:   PetscScalar    kr,ki;
788:   PetscBool      flg;

790:   PetscFunctionBegin;
794:   PetscAssertPointer(error,4);
795:   EPSCheckSolved(eps,1);

797:   /* allocate work vectors */
798: #if defined(PETSC_USE_COMPLEX)
799:   PetscCall(EPSSetWorkVecs(eps,3));
800:   xi   = NULL;
801:   w[1] = NULL;
802: #else
803:   PetscCall(EPSSetWorkVecs(eps,5));
804:   xi   = eps->work[3];
805:   w[1] = eps->work[4];
806: #endif
807:   xr   = eps->work[0];
808:   w[0] = eps->work[1];
809:   w[2] = eps->work[2];

811:   /* compute residual norm */
812:   PetscCall(EPSGetEigenpair(eps,i,&kr,&ki,xr,xi));
813:   PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_FALSE,kr,ki,xr,xi,w,error));

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

818:   /* if two-sided, compute left residual norm and take the maximum */
819:   if (eps->twosided) {
820:     PetscCall(EPSGetLeftEigenvector(eps,i,xr,xi));
821:     PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_TRUE,kr,ki,xr,xi,w,&errorl));
822:     *error = PetscMax(*error,errorl);
823:   }

825:   /* compute error */
826:   switch (type) {
827:     case EPS_ERROR_ABSOLUTE:
828:       break;
829:     case EPS_ERROR_RELATIVE:
830:       *error /= SlepcAbsEigenvalue(kr,ki)*vecnorm;
831:       break;
832:     case EPS_ERROR_BACKWARD:
833:       /* initialization of matrix norms */
834:       if (!eps->nrma) {
835:         PetscCall(STGetMatrix(eps->st,0,&A));
836:         PetscCall(MatHasOperation(A,MATOP_NORM,&flg));
837:         PetscCheck(flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
838:         PetscCall(MatNorm(A,NORM_INFINITY,&eps->nrma));
839:       }
840:       if (eps->isgeneralized) {
841:         if (!eps->nrmb) {
842:           PetscCall(STGetMatrix(eps->st,1,&B));
843:           PetscCall(MatHasOperation(B,MATOP_NORM,&flg));
844:           PetscCheck(flg,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
845:           PetscCall(MatNorm(B,NORM_INFINITY,&eps->nrmb));
846:         }
847:       } else eps->nrmb = 1.0;
848:       t = SlepcAbsEigenvalue(kr,ki);
849:       *error /= (eps->nrma+t*eps->nrmb)*vecnorm;
850:       break;
851:     default:
852:       SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
853:   }
854:   PetscFunctionReturn(PETSC_SUCCESS);
855: }

857: /*
858:    EPSGetStartVector - Generate a suitable vector to be used as the starting vector
859:    for the recurrence that builds the right subspace.

861:    Collective

863:    Input Parameters:
864: +  eps - the eigensolver context
865: -  i   - iteration number

867:    Output Parameters:
868: .  breakdown - flag indicating that a breakdown has occurred

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

877:    The flag breakdown is set to true if either i=0 and the vector belongs to the
878:    deflation space, or i>0 and the vector is linearly dependent with respect
879:    to the V-vectors.
880: */
881: PetscErrorCode EPSGetStartVector(EPS eps,PetscInt i,PetscBool *breakdown)
882: {
883:   PetscReal      norm;
884:   PetscBool      lindep;
885:   Vec            w,z;

887:   PetscFunctionBegin;

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

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

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

915: /*
916:    EPSGetLeftStartVector - Generate a suitable vector to be used as the left starting
917:    vector for the recurrence that builds the left subspace. See EPSGetStartVector().
918: */
919: PetscErrorCode EPSGetLeftStartVector(EPS eps,PetscInt i,PetscBool *breakdown)
920: {
921:   PetscReal      norm;
922:   PetscBool      lindep;
923:   Vec            w,z;

925:   PetscFunctionBegin;

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

932:   /* Force the vector to be in the range of OP' for generalized problems with B-inner product */
933:   if (eps->ispositive || (eps->isgeneralized && eps->ishermitian)) {
934:     PetscCall(BVCreateVec(eps->W,&w));
935:     PetscCall(BVCopyVec(eps->W,i,w));
936:     PetscCall(BVGetColumn(eps->W,i,&z));
937:     PetscCall(STApplyHermitianTranspose(eps->st,w,z));
938:     PetscCall(BVRestoreColumn(eps->W,i,&z));
939:     PetscCall(VecDestroy(&w));
940:   }

942:   /* Orthonormalize the vector with respect to previous vectors */
943:   PetscCall(BVOrthogonalizeColumn(eps->W,i,NULL,&norm,&lindep));
944:   if (breakdown) *breakdown = lindep;
945:   else if (lindep || norm == 0.0) {
946:     PetscCheck(i,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Left initial vector is zero");
947:     PetscCheck(!i,PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Unable to generate more left start vectors");
948:   }
949:   PetscCall(BVScaleColumn(eps->W,i,1.0/norm));
950:   PetscFunctionReturn(PETSC_SUCCESS);
951: }