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 - the linear eigensolver context

 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: [](ch:eps), `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 linear 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: [](ch:eps), `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 linear 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: [](ch:eps), `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 linear eigensolver context

292:    Output Parameter:
293: .  reason - negative value indicates diverged, positive value converged, see
294:    `EPSConvergedReason` for the possible values

296:    Options Database Key:
297: .  -eps_converged_reason - print reason for convergence/divergence, and number of iterations

299:    Note:
300:    If this routine is called before or doing the `EPSSolve()` the value of
301:    `EPS_CONVERGED_ITERATING` is returned.

303:    Level: intermediate

305: .seealso: [](ch:eps), `EPSSetTolerances()`, `EPSSolve()`, `EPSConvergedReason`
306: @*/
307: PetscErrorCode EPSGetConvergedReason(EPS eps,EPSConvergedReason *reason)
308: {
309:   PetscFunctionBegin;
311:   PetscAssertPointer(reason,2);
312:   EPSCheckSolved(eps,1);
313:   *reason = eps->reason;
314:   PetscFunctionReturn(PETSC_SUCCESS);
315: }

317: /*@
318:    EPSGetInvariantSubspace - Gets an orthonormal basis of the computed invariant
319:    subspace.

321:    Collective

323:    Input Parameter:
324: .  eps - the linear eigensolver context

326:    Output Parameter:
327: .  v - an array of vectors

329:    Notes:
330:    This function should be called after EPSSolve() has finished.

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

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

342:    Level: intermediate

344: .seealso: [](ch:eps), `EPSGetEigenpair()`, `EPSGetConverged()`, `EPSSolve()`
345: @*/
346: PetscErrorCode EPSGetInvariantSubspace(EPS eps,Vec v[])
347: {
348:   PetscInt       i;
349:   BV             V=eps->V;
350:   Vec            w;

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

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

378:    Collective

380:    Input Parameters:
381: +  eps - the linear eigensolver context
382: -  i   - index of the solution

384:    Output Parameters:
385: +  eigr - real part of eigenvalue
386: .  eigi - imaginary part of eigenvalue
387: .  Vr   - real part of eigenvector
388: -  Vi   - imaginary part of eigenvector

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

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

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

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

408:    Level: beginner

410: .seealso: [](ch:eps), `EPSGetEigenvalue()`, `EPSGetEigenvector()`, `EPSGetLeftEigenvector()`, `EPSSolve()`, `EPSGetConverged()`, `EPSSetWhichEigenpairs()`, `EPSGetInvariantSubspace()`
411: @*/
412: PetscErrorCode EPSGetEigenpair(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
413: {
414:   PetscInt nconv;

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

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

431:    Not Collective

433:    Input Parameters:
434: +  eps - the linear eigensolver context
435: -  i   - index of the solution

437:    Output Parameters:
438: +  eigr - real part of eigenvalue
439: -  eigi - imaginary part of eigenvalue

441:    Notes:
442:    If the eigenvalue is real, then eigi is set to zero. If PETSc is
443:    configured with complex scalars the eigenvalue is stored
444:    directly in eigr (eigi is set to zero).

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

450:    Level: beginner

452: .seealso: [](ch:eps), `EPSSolve()`, `EPSGetConverged()`, `EPSSetWhichEigenpairs()`, `EPSGetEigenpair()`
453: @*/
454: PetscErrorCode EPSGetEigenvalue(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi)
455: {
456:   PetscInt  k,nconv;
457: #if !defined(PETSC_USE_COMPLEX)
458:   PetscInt  k2, iquad;
459: #endif

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

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

519:    Collective

521:    Input Parameters:
522: +  eps - the linear eigensolver context
523: -  i   - index of the solution

525:    Output Parameters:
526: +  Vr   - real part of eigenvector
527: -  Vi   - imaginary part of eigenvector

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

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

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

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

546:    Level: beginner

548: .seealso: [](ch:eps), `EPSSolve()`, `EPSGetConverged()`, `EPSSetWhichEigenpairs()`, `EPSGetEigenpair()`, `EPSGetLeftEigenvector()`
549: @*/
550: PetscErrorCode EPSGetEigenvector(EPS eps,PetscInt i,Vec Vr,Vec Vi)
551: {
552:   PetscInt nconv;

554:   PetscFunctionBegin;
559:   EPSCheckSolved(eps,1);
560:   PetscCheck(i>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
561:   PetscCall(EPS_GetActualConverged(eps,&nconv));
562:   PetscCheck(i<nconv,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see EPSGetConverged()");
563:   PetscCall(EPSComputeVectors(eps));
564:   PetscCall(EPS_GetEigenvector(eps,eps->V,i,Vr,Vi));
565:   PetscFunctionReturn(PETSC_SUCCESS);
566: }

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

571:    Collective

573:    Input Parameters:
574: +  eps - the linear eigensolver context
575: -  i   - index of the solution

577:    Output Parameters:
578: +  Wr   - real part of left eigenvector
579: -  Wi   - imaginary part of left eigenvector

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

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

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

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

597:    Level: intermediate

599: .seealso: [](ch:eps), `EPSGetEigenvector()`, `EPSGetConverged()`, `EPSSetWhichEigenpairs()`, `EPSSetTwoSided()`
600: @*/
601: PetscErrorCode EPSGetLeftEigenvector(EPS eps,PetscInt i,Vec Wr,Vec Wi)
602: {
603:   PetscInt    nconv;
604:   PetscBool   trivial;
605:   Mat         H;
606:   IS          is[2];
607:   Vec         v;

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

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

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

647: /*@
648:    EPSGetErrorEstimate - Returns the error estimate associated to the i-th
649:    computed eigenpair.

651:    Not Collective

653:    Input Parameters:
654: +  eps - the linear eigensolver context
655: -  i   - index of eigenpair

657:    Output Parameter:
658: .  errest - the error estimate

660:    Notes:
661:    This is the error estimate used internally by the eigensolver. The actual
662:    error bound can be computed with EPSComputeError(). See also the users
663:    manual for details.

665:    Level: advanced

667: .seealso: [](ch:eps), `EPSComputeError()`
668: @*/
669: PetscErrorCode EPSGetErrorEstimate(EPS eps,PetscInt i,PetscReal *errest)
670: {
671:   PetscInt nconv;

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

690: /*
691:    EPSComputeResidualNorm_Private - Computes the norm of the residual vector
692:    associated with an eigenpair.

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

712:   PetscFunctionBegin;
713:   u = z[0]; w = z[2];
714:   PetscCall(STGetNumMatrices(eps->st,&nmat));
715:   PetscCall(STGetMatrix(eps->st,0,&A));
716:   if (nmat>1) PetscCall(STGetMatrix(eps->st,1,&B));

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

754: /*@
755:    EPSComputeError - Computes the error (based on the residual norm) associated
756:    with the i-th computed eigenpair.

758:    Collective

760:    Input Parameters:
761: +  eps  - the linear eigensolver context
762: .  i    - the solution index
763: -  type - the type of error to compute

765:    Output Parameter:
766: .  error - the error

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

772:    Level: beginner

774: .seealso: [](ch:eps), `EPSErrorType`, `EPSSolve()`, `EPSGetErrorEstimate()`
775: @*/
776: PetscErrorCode EPSComputeError(EPS eps,PetscInt i,EPSErrorType type,PetscReal *error)
777: {
778:   Mat            A,B;
779:   Vec            xr,xi,w[3];
780:   PetscReal      t,vecnorm=1.0,errorl;
781:   PetscScalar    kr,ki;
782:   PetscBool      flg;

784:   PetscFunctionBegin;
788:   PetscAssertPointer(error,4);
789:   EPSCheckSolved(eps,1);

791:   /* allocate work vectors */
792: #if defined(PETSC_USE_COMPLEX)
793:   PetscCall(EPSSetWorkVecs(eps,3));
794:   xi   = NULL;
795:   w[1] = NULL;
796: #else
797:   PetscCall(EPSSetWorkVecs(eps,5));
798:   xi   = eps->work[3];
799:   w[1] = eps->work[4];
800: #endif
801:   xr   = eps->work[0];
802:   w[0] = eps->work[1];
803:   w[2] = eps->work[2];

805:   /* compute residual norm */
806:   PetscCall(EPSGetEigenpair(eps,i,&kr,&ki,xr,xi));
807:   PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_FALSE,kr,ki,xr,xi,w,error));

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

812:   /* if two-sided, compute left residual norm and take the maximum */
813:   if (eps->twosided) {
814:     PetscCall(EPSGetLeftEigenvector(eps,i,xr,xi));
815:     PetscCall(EPSComputeResidualNorm_Private(eps,PETSC_TRUE,kr,ki,xr,xi,w,&errorl));
816:     *error = PetscMax(*error,errorl);
817:   }

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

851: /*
852:    EPSGetStartVector - Generate a suitable vector to be used as the starting vector
853:    for the recurrence that builds the right subspace.

855:    Collective

857:    Input Parameters:
858: +  eps - the linear eigensolver context
859: -  i   - iteration number

861:    Output Parameter:
862: .  breakdown - flag indicating that a breakdown has occurred

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

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

881:   PetscFunctionBegin;

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

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

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

909: /*
910:    EPSGetLeftStartVector - Generate a suitable vector to be used as the left starting
911:    vector for the recurrence that builds the left subspace. See EPSGetStartVector().
912: */
913: PetscErrorCode EPSGetLeftStartVector(EPS eps,PetscInt i,PetscBool *breakdown)
914: {
915:   PetscReal      norm;
916:   PetscBool      lindep;
917:   Vec            w,z;

919:   PetscFunctionBegin;

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

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

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