Actual source code: epssolve.c

slepc-3.18.0 2022-10-01
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:   EPSCheckSolved(eps,1);
 21:   if (eps->state==EPS_STATE_SOLVED) PetscTryTypeMethod(eps,computevectors);
 22:   eps->state = EPS_STATE_EIGENVECTORS;
 23:   return 0;
 24: }

 26: #define SWAP(a,b,t) {t=a;a=b;b=t;}

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

 34:   switch (eps->categ) {
 35:     case EPS_CATEGORY_KRYLOV:
 36:     case EPS_CATEGORY_OTHER:
 37:       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:         DSGetDimensions(eps->ds,&n,NULL,NULL,NULL);
 44:         if (!n) break;
 45:         EPSGetOperators(eps,&A,&B);
 46:         BVSetActiveColumns(eps->V,0,n);
 47:         DSGetCompact(eps->ds,&iscomp);
 48:         DSSetCompact(eps->ds,PETSC_FALSE);
 49:         DSGetMat(eps->ds,DS_MAT_A,&G);
 50:         BVMatProject(eps->V,A,eps->V,G);
 51:         DSRestoreMat(eps->ds,DS_MAT_A,&G);
 52:         if (B) {
 53:           DSGetMat(eps->ds,DS_MAT_B,&G);
 54:           BVMatProject(eps->V,B,eps->V,G);
 55:           DSRestoreMat(eps->ds,DS_MAT_A,&G);
 56:         }
 57:         DSSolve(eps->ds,eps->eigr,eps->eigi);
 58:         DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL);
 59:         DSSynchronize(eps->ds,eps->eigr,eps->eigi);
 60:         DSSetCompact(eps->ds,iscomp);
 61:         if (eps->ishermitian && (!eps->isgeneralized || eps->ispositive)) { /* V = V * Z */
 62:           DSVectors(eps->ds,DS_MAT_X,NULL,NULL);
 63:           DSGetMat(eps->ds,DS_MAT_X,&Z);
 64:           BVMultInPlace(eps->V,Z,0,n);
 65:           DSRestoreMat(eps->ds,DS_MAT_X,&Z);
 66:         }
 67:         /* in case of STFILTER discard computed eigenvalues that lie outside the wanted interval */
 68:         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) { SWAP(eps->perm[i],eps->perm[eps->nconv],aux); i--; }
 75:             }
 76:           }
 77:           if (nconv0>eps->nconv) 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:   return 0;
 87: }

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

 92:    Collective on eps

 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:   STMatMode      matmode;
123:   Mat            A,B;

126:   if (eps->state>=EPS_STATE_SOLVED) return 0;
127:   PetscLogEventBegin(EPS_Solve,eps,0,0,0);

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

142:   /* Call solver */
143:   PetscUseTypeMethod(eps,solve);
145:   eps->state = EPS_STATE_SOLVED;

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

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

156:   /* Map eigenvalues back to the original problem if appropriate */
157:   EPSComputeValues(eps);

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

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

179:   /* Various viewers */
180:   EPSViewFromOptions(eps,NULL,"-eps_view");
181:   EPSConvergedReasonViewFromOptions(eps);
182:   EPSErrorViewFromOptions(eps);
183:   EPSValuesViewFromOptions(eps);
184:   EPSVectorsViewFromOptions(eps);
185:   EPSGetOperators(eps,&A,&B);
186:   MatViewFromOptions(A,(PetscObject)eps,"-eps_view_mat0");
187:   if (eps->isgeneralized) MatViewFromOptions(B,(PetscObject)eps,"-eps_view_mat1");

189:   /* Remove deflation and initial subspaces */
190:   if (eps->nds) {
191:     BVSetNumConstraints(eps->V,0);
192:     eps->nds = 0;
193:   }
194:   eps->nini = 0;
195:   return 0;
196: }

198: /*@
199:    EPSGetIterationNumber - Gets the current iteration number. If the
200:    call to EPSSolve() is complete, then it returns the number of iterations
201:    carried out by the solution method.

203:    Not Collective

205:    Input Parameter:
206: .  eps - the eigensolver context

208:    Output Parameter:
209: .  its - number of iterations

211:    Note:
212:    During the i-th iteration this call returns i-1. If EPSSolve() is
213:    complete, then parameter "its" contains either the iteration number at
214:    which convergence was successfully reached, or failure was detected.
215:    Call EPSGetConvergedReason() to determine if the solver converged or
216:    failed and why.

218:    Level: intermediate

220: .seealso: EPSGetConvergedReason(), EPSSetTolerances()
221: @*/
222: PetscErrorCode EPSGetIterationNumber(EPS eps,PetscInt *its)
223: {
226:   *its = eps->its;
227:   return 0;
228: }

230: /*@
231:    EPSGetConverged - Gets the number of converged eigenpairs.

233:    Not Collective

235:    Input Parameter:
236: .  eps - the eigensolver context

238:    Output Parameter:
239: .  nconv - number of converged eigenpairs

241:    Note:
242:    This function should be called after EPSSolve() has finished.

244:    Level: beginner

246: .seealso: EPSSetDimensions(), EPSSolve(), EPSGetEigenpair()
247: @*/
248: PetscErrorCode EPSGetConverged(EPS eps,PetscInt *nconv)
249: {
252:   EPSCheckSolved(eps,1);
253:   *nconv = eps->nconv;
254:   return 0;
255: }

257: /*@
258:    EPSGetConvergedReason - Gets the reason why the EPSSolve() iteration was
259:    stopped.

261:    Not Collective

263:    Input Parameter:
264: .  eps - the eigensolver context

266:    Output Parameter:
267: .  reason - negative value indicates diverged, positive value converged

269:    Options Database Key:
270: .  -eps_converged_reason - print the reason to a viewer

272:    Notes:
273:    Possible values for reason are
274: +  EPS_CONVERGED_TOL - converged up to tolerance
275: .  EPS_CONVERGED_USER - converged due to a user-defined condition
276: .  EPS_DIVERGED_ITS - required more than max_it iterations to reach convergence
277: .  EPS_DIVERGED_BREAKDOWN - generic breakdown in method
278: -  EPS_DIVERGED_SYMMETRY_LOST - pseudo-Lanczos was not able to keep symmetry

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

282:    Level: intermediate

284: .seealso: EPSSetTolerances(), EPSSolve(), EPSConvergedReason
285: @*/
286: PetscErrorCode EPSGetConvergedReason(EPS eps,EPSConvergedReason *reason)
287: {
290:   EPSCheckSolved(eps,1);
291:   *reason = eps->reason;
292:   return 0;
293: }

295: /*@
296:    EPSGetInvariantSubspace - Gets an orthonormal basis of the computed invariant
297:    subspace.

299:    Not Collective, but vectors are shared by all processors that share the EPS

301:    Input Parameter:
302: .  eps - the eigensolver context

304:    Output Parameter:
305: .  v - an array of vectors

307:    Notes:
308:    This function should be called after EPSSolve() has finished.

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

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

320:    Level: intermediate

322: .seealso: EPSGetEigenpair(), EPSGetConverged(), EPSSolve()
323: @*/
324: PetscErrorCode EPSGetInvariantSubspace(EPS eps,Vec v[])
325: {
326:   PetscInt       i;
327:   BV             V=eps->V;
328:   Vec            w;

333:   EPSCheckSolved(eps,1);
335:   if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
336:     BVDuplicateResize(eps->V,eps->nconv,&V);
337:     BVSetActiveColumns(eps->V,0,eps->nconv);
338:     BVCopy(eps->V,V);
339:     for (i=0;i<eps->nconv;i++) {
340:       BVGetColumn(V,i,&w);
341:       VecPointwiseDivide(w,w,eps->D);
342:       BVRestoreColumn(V,i,&w);
343:     }
344:     BVOrthogonalize(V,NULL);
345:   }
346:   for (i=0;i<eps->nconv;i++) BVCopyVec(V,i,v[i]);
347:   if (eps->balance!=EPS_BALANCE_NONE && eps->D) BVDestroy(&V);
348:   return 0;
349: }

351: /*@C
352:    EPSGetEigenpair - Gets the i-th solution of the eigenproblem as computed by
353:    EPSSolve(). The solution consists in both the eigenvalue and the eigenvector.

355:    Logically Collective on eps

357:    Input Parameters:
358: +  eps - eigensolver context
359: -  i   - index of the solution

361:    Output Parameters:
362: +  eigr - real part of eigenvalue
363: .  eigi - imaginary part of eigenvalue
364: .  Vr   - real part of eigenvector
365: -  Vi   - imaginary part of eigenvector

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

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

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

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

385:    Level: beginner

387: .seealso: EPSGetEigenvalue(), EPSGetEigenvector(), EPSGetLeftEigenvector(), EPSSolve(),
388:           EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetInvariantSubspace()
389: @*/
390: PetscErrorCode EPSGetEigenpair(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
391: {
394:   EPSCheckSolved(eps,1);
397:   EPSGetEigenvalue(eps,i,eigr,eigi);
398:   if (Vr || Vi) EPSGetEigenvector(eps,i,Vr,Vi);
399:   return 0;
400: }

402: /*@C
403:    EPSGetEigenvalue - Gets the i-th eigenvalue as computed by EPSSolve().

405:    Not Collective

407:    Input Parameters:
408: +  eps - eigensolver context
409: -  i   - index of the solution

411:    Output Parameters:
412: +  eigr - real part of eigenvalue
413: -  eigi - imaginary part of eigenvalue

415:    Notes:
416:    If the eigenvalue is real, then eigi is set to zero. If PETSc is
417:    configured with complex scalars the eigenvalue is stored
418:    directly in eigr (eigi is set to zero).

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

424:    Level: beginner

426: .seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetEigenpair()
427: @*/
428: PetscErrorCode EPSGetEigenvalue(EPS eps,PetscInt i,PetscScalar *eigr,PetscScalar *eigi)
429: {
430:   PetscInt k;

433:   EPSCheckSolved(eps,1);
436:   k = eps->perm[i];
437: #if defined(PETSC_USE_COMPLEX)
438:   if (eigr) *eigr = eps->eigr[k];
439:   if (eigi) *eigi = 0;
440: #else
441:   if (eigr) *eigr = eps->eigr[k];
442:   if (eigi) *eigi = eps->eigi[k];
443: #endif
444:   return 0;
445: }

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

450:    Logically Collective on eps

452:    Input Parameters:
453: +  eps - eigensolver context
454: -  i   - index of the solution

456:    Output Parameters:
457: +  Vr   - real part of eigenvector
458: -  Vi   - imaginary part of eigenvector

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

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

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

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

477:    Level: beginner

479: .seealso: EPSSolve(), EPSGetConverged(), EPSSetWhichEigenpairs(), EPSGetEigenpair(), EPSGetLeftEigenvector()
480: @*/
481: PetscErrorCode EPSGetEigenvector(EPS eps,PetscInt i,Vec Vr,Vec Vi)
482: {
483:   PetscInt       k;

489:   EPSCheckSolved(eps,1);
492:   EPSComputeVectors(eps);
493:   k = eps->perm[i];
494:   BV_GetEigenvector(eps->V,k,eps->eigi[k],Vr,Vi);
495:   return 0;
496: }

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

501:    Logically Collective on eps

503:    Input Parameters:
504: +  eps - eigensolver context
505: -  i   - index of the solution

507:    Output Parameters:
508: +  Wr   - real part of left eigenvector
509: -  Wi   - imaginary part of left eigenvector

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

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

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

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

527:    Level: intermediate

529: .seealso: EPSGetEigenvector(), EPSGetConverged(), EPSSetWhichEigenpairs(), EPSSetTwoSided()
530: @*/
531: PetscErrorCode EPSGetLeftEigenvector(EPS eps,PetscInt i,Vec Wr,Vec Wi)
532: {
533:   PetscInt       k;

539:   EPSCheckSolved(eps,1);
543:   EPSComputeVectors(eps);
544:   k = eps->perm[i];
545:   BV_GetEigenvector(eps->W,k,eps->eigi[k],Wr,Wi);
546:   return 0;
547: }

549: /*@
550:    EPSGetErrorEstimate - Returns the error estimate associated to the i-th
551:    computed eigenpair.

553:    Not Collective

555:    Input Parameters:
556: +  eps - eigensolver context
557: -  i   - index of eigenpair

559:    Output Parameter:
560: .  errest - the error estimate

562:    Notes:
563:    This is the error estimate used internally by the eigensolver. The actual
564:    error bound can be computed with EPSComputeError(). See also the users
565:    manual for details.

567:    Level: advanced

569: .seealso: EPSComputeError()
570: @*/
571: PetscErrorCode EPSGetErrorEstimate(EPS eps,PetscInt i,PetscReal *errest)
572: {
575:   EPSCheckSolved(eps,1);
578:   *errest = eps->errest[eps->perm[i]];
579:   return 0;
580: }

582: /*
583:    EPSComputeResidualNorm_Private - Computes the norm of the residual vector
584:    associated with an eigenpair.

586:    Input Parameters:
587:      trans - whether A' must be used instead of A
588:      kr,ki - eigenvalue
589:      xr,xi - eigenvector
590:      z     - three work vectors (the second one not referenced in complex scalars)
591: */
592: PetscErrorCode EPSComputeResidualNorm_Private(EPS eps,PetscBool trans,PetscScalar kr,PetscScalar ki,Vec xr,Vec xi,Vec *z,PetscReal *norm)
593: {
594:   PetscInt       nmat;
595:   Mat            A,B;
596:   Vec            u,w;
597:   PetscScalar    alpha;
598: #if !defined(PETSC_USE_COMPLEX)
599:   Vec            v;
600:   PetscReal      ni,nr;
601: #endif
602:   PetscErrorCode (*matmult)(Mat,Vec,Vec) = trans? MatMultHermitianTranspose: MatMult;

604:   u = z[0]; w = z[2];
605:   STGetNumMatrices(eps->st,&nmat);
606:   STGetMatrix(eps->st,0,&A);
607:   if (nmat>1) STGetMatrix(eps->st,1,&B);

609: #if !defined(PETSC_USE_COMPLEX)
610:   v = z[1];
611:   if (ki == 0 || PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON)) {
612: #endif
613:     (*matmult)(A,xr,u);                          /* u=A*x */
614:     if (PetscAbsScalar(kr) > PETSC_MACHINE_EPSILON) {
615:       if (nmat>1) (*matmult)(B,xr,w);
616:       else VecCopy(xr,w);                        /* w=B*x */
617:       alpha = trans? -PetscConj(kr): -kr;
618:       VecAXPY(u,alpha,w);                        /* u=A*x-k*B*x */
619:     }
620:     VecNorm(u,NORM_2,norm);
621: #if !defined(PETSC_USE_COMPLEX)
622:   } else {
623:     (*matmult)(A,xr,u);                          /* u=A*xr */
624:     if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
625:       if (nmat>1) (*matmult)(B,xr,v);
626:       else VecCopy(xr,v);                        /* v=B*xr */
627:       VecAXPY(u,-kr,v);                          /* u=A*xr-kr*B*xr */
628:       if (nmat>1) (*matmult)(B,xi,w);
629:       else VecCopy(xi,w);                        /* w=B*xi */
630:       VecAXPY(u,trans?-ki:ki,w);                 /* u=A*xr-kr*B*xr+ki*B*xi */
631:     }
632:     VecNorm(u,NORM_2,&nr);
633:     (*matmult)(A,xi,u);                          /* u=A*xi */
634:     if (SlepcAbsEigenvalue(kr,ki) > PETSC_MACHINE_EPSILON) {
635:       VecAXPY(u,-kr,w);                          /* u=A*xi-kr*B*xi */
636:       VecAXPY(u,trans?ki:-ki,v);                 /* u=A*xi-kr*B*xi-ki*B*xr */
637:     }
638:     VecNorm(u,NORM_2,&ni);
639:     *norm = SlepcAbsEigenvalue(nr,ni);
640:   }
641: #endif
642:   return 0;
643: }

645: /*@
646:    EPSComputeError - Computes the error (based on the residual norm) associated
647:    with the i-th computed eigenpair.

649:    Collective on eps

651:    Input Parameters:
652: +  eps  - the eigensolver context
653: .  i    - the solution index
654: -  type - the type of error to compute

656:    Output Parameter:
657: .  error - the error

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

663:    Level: beginner

665: .seealso: EPSErrorType, EPSSolve(), EPSGetErrorEstimate()
666: @*/
667: PetscErrorCode EPSComputeError(EPS eps,PetscInt i,EPSErrorType type,PetscReal *error)
668: {
669:   Mat            A,B;
670:   Vec            xr,xi,w[3];
671:   PetscReal      t,vecnorm=1.0,errorl;
672:   PetscScalar    kr,ki;
673:   PetscBool      flg;

679:   EPSCheckSolved(eps,1);

681:   /* allocate work vectors */
682: #if defined(PETSC_USE_COMPLEX)
683:   EPSSetWorkVecs(eps,3);
684:   xi   = NULL;
685:   w[1] = NULL;
686: #else
687:   EPSSetWorkVecs(eps,5);
688:   xi   = eps->work[3];
689:   w[1] = eps->work[4];
690: #endif
691:   xr   = eps->work[0];
692:   w[0] = eps->work[1];
693:   w[2] = eps->work[2];

695:   /* compute residual norm */
696:   EPSGetEigenpair(eps,i,&kr,&ki,xr,xi);
697:   EPSComputeResidualNorm_Private(eps,PETSC_FALSE,kr,ki,xr,xi,w,error);

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

702:   /* if two-sided, compute left residual norm and take the maximum */
703:   if (eps->twosided) {
704:     EPSGetLeftEigenvector(eps,i,xr,xi);
705:     EPSComputeResidualNorm_Private(eps,PETSC_TRUE,kr,ki,xr,xi,w,&errorl);
706:     *error = PetscMax(*error,errorl);
707:   }

709:   /* compute error */
710:   switch (type) {
711:     case EPS_ERROR_ABSOLUTE:
712:       break;
713:     case EPS_ERROR_RELATIVE:
714:       *error /= SlepcAbsEigenvalue(kr,ki)*vecnorm;
715:       break;
716:     case EPS_ERROR_BACKWARD:
717:       /* initialization of matrix norms */
718:       if (!eps->nrma) {
719:         STGetMatrix(eps->st,0,&A);
720:         MatHasOperation(A,MATOP_NORM,&flg);
722:         MatNorm(A,NORM_INFINITY,&eps->nrma);
723:       }
724:       if (eps->isgeneralized) {
725:         if (!eps->nrmb) {
726:           STGetMatrix(eps->st,1,&B);
727:           MatHasOperation(B,MATOP_NORM,&flg);
729:           MatNorm(B,NORM_INFINITY,&eps->nrmb);
730:         }
731:       } else eps->nrmb = 1.0;
732:       t = SlepcAbsEigenvalue(kr,ki);
733:       *error /= (eps->nrma+t*eps->nrmb)*vecnorm;
734:       break;
735:     default:
736:       SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
737:   }
738:   return 0;
739: }

741: /*
742:    EPSGetStartVector - Generate a suitable vector to be used as the starting vector
743:    for the recurrence that builds the right subspace.

745:    Collective on eps

747:    Input Parameters:
748: +  eps - the eigensolver context
749: -  i   - iteration number

751:    Output Parameters:
752: .  breakdown - flag indicating that a breakdown has occurred

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

761:    The flag breakdown is set to true if either i=0 and the vector belongs to the
762:    deflation space, or i>0 and the vector is linearly dependent with respect
763:    to the V-vectors.
764: */
765: PetscErrorCode EPSGetStartVector(EPS eps,PetscInt i,PetscBool *breakdown)
766: {
767:   PetscReal      norm;
768:   PetscBool      lindep;
769:   Vec            w,z;


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

777:   /* Force the vector to be in the range of OP for definite generalized problems */
778:   if (eps->ispositive || (eps->isgeneralized && eps->ishermitian)) {
779:     BVCreateVec(eps->V,&w);
780:     BVCopyVec(eps->V,i,w);
781:     BVGetColumn(eps->V,i,&z);
782:     STApply(eps->st,w,z);
783:     BVRestoreColumn(eps->V,i,&z);
784:     VecDestroy(&w);
785:   }

787:   /* Orthonormalize the vector with respect to previous vectors */
788:   BVOrthogonalizeColumn(eps->V,i,NULL,&norm,&lindep);
789:   if (breakdown) *breakdown = lindep;
790:   else if (lindep || norm == 0.0) {
793:   }
794:   BVScaleColumn(eps->V,i,1.0/norm);
795:   return 0;
796: }

798: /*
799:    EPSGetLeftStartVector - Generate a suitable vector to be used as the left starting
800:    vector for the recurrence that builds the left subspace. See EPSGetStartVector().
801: */
802: PetscErrorCode EPSGetLeftStartVector(EPS eps,PetscInt i,PetscBool *breakdown)
803: {
804:   PetscReal      norm;
805:   PetscBool      lindep;


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

813:   /* Orthonormalize the vector with respect to previous vectors */
814:   BVOrthogonalizeColumn(eps->W,i,NULL,&norm,&lindep);
815:   if (breakdown) *breakdown = lindep;
816:   else if (lindep || norm == 0.0) {
818:     SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_CONV_FAILED,"Unable to generate more left start vectors");
819:   }
820:   BVScaleColumn(eps->W,i,1.0/norm);
821:   return 0;
822: }