Actual source code: pepsolve.c

slepc-3.16.1 2021-11-17
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2021, 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:    PEP routines related to the solution process

 13:    References:

 15:        [1] C. Campos and J.E. Roman, "Parallel iterative refinement in
 16:            polynomial eigenvalue problems", Numer. Linear Algebra Appl.
 17:            23(4):730-745, 2016.
 18: */

 20: #include <slepc/private/pepimpl.h>
 21: #include <slepc/private/bvimpl.h>
 22: #include <petscdraw.h>

 24: static PetscBool  cited = PETSC_FALSE;
 25: static const char citation[] =
 26:   "@Article{slepc-pep-refine,\n"
 27:   "   author = \"C. Campos and J. E. Roman\",\n"
 28:   "   title = \"Parallel iterative refinement in polynomial eigenvalue problems\",\n"
 29:   "   journal = \"Numer. Linear Algebra Appl.\",\n"
 30:   "   volume = \"23\",\n"
 31:   "   number = \"4\",\n"
 32:   "   pages = \"730--745\",\n"
 33:   "   year = \"2016,\"\n"
 34:   "   doi = \"https://doi.org/10.1002/nla.2052\"\n"
 35:   "}\n";

 37: PetscErrorCode PEPComputeVectors(PEP pep)
 38: {

 42:   PEPCheckSolved(pep,1);
 43:   if (pep->state==PEP_STATE_SOLVED && pep->ops->computevectors) {
 44:     (*pep->ops->computevectors)(pep);
 45:   }
 46:   pep->state = PEP_STATE_EIGENVECTORS;
 47:   return(0);
 48: }

 50: PetscErrorCode PEPExtractVectors(PEP pep)
 51: {

 55:   PEPCheckSolved(pep,1);
 56:   if (pep->state==PEP_STATE_SOLVED && pep->ops->extractvectors) {
 57:     (*pep->ops->extractvectors)(pep);
 58:   }
 59:   return(0);
 60: }

 62: /*@
 63:    PEPSolve - Solves the polynomial eigensystem.

 65:    Collective on pep

 67:    Input Parameter:
 68: .  pep - eigensolver context obtained from PEPCreate()

 70:    Options Database Keys:
 71: +  -pep_view - print information about the solver used
 72: .  -pep_view_matk binary - save any of the coefficient matrices (Ak) to the
 73:                 default binary viewer (replace k by an integer from 0 to nmat-1)
 74: .  -pep_view_vectors binary - save the computed eigenvectors to the default binary viewer
 75: .  -pep_view_values - print computed eigenvalues
 76: .  -pep_converged_reason - print reason for convergence, and number of iterations
 77: .  -pep_error_absolute - print absolute errors of each eigenpair
 78: .  -pep_error_relative - print relative errors of each eigenpair
 79: -  -pep_error_backward - print backward errors of each eigenpair

 81:    Level: beginner

 83: .seealso: PEPCreate(), PEPSetUp(), PEPDestroy(), PEPSetTolerances()
 84: @*/
 85: PetscErrorCode PEPSolve(PEP pep)
 86: {
 88:   PetscInt       i,k;
 89:   PetscBool      flg,islinear;
 90:   char           str[16];

 94:   if (pep->state>=PEP_STATE_SOLVED) return(0);
 95:   PetscLogEventBegin(PEP_Solve,pep,0,0,0);

 97:   /* call setup */
 98:   PEPSetUp(pep);
 99:   pep->nconv = 0;
100:   pep->its   = 0;
101:   k = pep->lineariz? pep->ncv: pep->ncv*(pep->nmat-1);
102:   for (i=0;i<k;i++) {
103:     pep->eigr[i]   = 0.0;
104:     pep->eigi[i]   = 0.0;
105:     pep->errest[i] = 0.0;
106:     pep->perm[i]   = i;
107:   }
108:   PEPViewFromOptions(pep,NULL,"-pep_view_pre");
109:   RGViewFromOptions(pep->rg,NULL,"-rg_view");

111:   /* Call solver */
112:   (*pep->ops->solve)(pep);
113:   if (!pep->reason) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
114:   pep->state = PEP_STATE_SOLVED;

116:   /* Only the first nconv columns contain useful information */
117:   BVSetActiveColumns(pep->V,0,pep->nconv);

119:   PetscObjectTypeCompare((PetscObject)pep,PEPLINEAR,&islinear);
120:   if (!islinear) {
121:     STPostSolve(pep->st);
122:     /* Map eigenvalues back to the original problem */
123:     STGetTransform(pep->st,&flg);
124:     if (flg && pep->ops->backtransform) {
125:       (*pep->ops->backtransform)(pep);
126:     }
127:   }

129: #if !defined(PETSC_USE_COMPLEX)
130:   /* reorder conjugate eigenvalues (positive imaginary first) */
131:   for (i=0;i<pep->nconv-1;i++) {
132:     if (pep->eigi[i] != 0) {
133:       if (pep->eigi[i] < 0) {
134:         pep->eigi[i] = -pep->eigi[i];
135:         pep->eigi[i+1] = -pep->eigi[i+1];
136:         /* the next correction only works with eigenvectors */
137:         PEPComputeVectors(pep);
138:         BVScaleColumn(pep->V,i+1,-1.0);
139:       }
140:       i++;
141:     }
142:   }
143: #endif

145:   if (pep->refine!=PEP_REFINE_NONE) {
146:     PetscCitationsRegister(citation,&cited);
147:   }

149:   if (pep->refine==PEP_REFINE_SIMPLE && pep->rits>0 && pep->nconv>0) {
150:     PEPComputeVectors(pep);
151:     PEPNewtonRefinementSimple(pep,&pep->rits,pep->rtol,pep->nconv);
152:   }

154:   /* sort eigenvalues according to pep->which parameter */
155:   SlepcSortEigenvalues(pep->sc,pep->nconv,pep->eigr,pep->eigi,pep->perm);
156:   PetscLogEventEnd(PEP_Solve,pep,0,0,0);

158:   /* various viewers */
159:   PEPViewFromOptions(pep,NULL,"-pep_view");
160:   PEPConvergedReasonViewFromOptions(pep);
161:   PEPErrorViewFromOptions(pep);
162:   PEPValuesViewFromOptions(pep);
163:   PEPVectorsViewFromOptions(pep);
164:   for (i=0;i<pep->nmat;i++) {
165:     PetscSNPrintf(str,sizeof(str),"-pep_view_mat%d",(int)i);
166:     MatViewFromOptions(pep->A[i],(PetscObject)pep,str);
167:   }

169:   /* Remove the initial subspace */
170:   pep->nini = 0;
171:   return(0);
172: }

174: /*@
175:    PEPGetIterationNumber - Gets the current iteration number. If the
176:    call to PEPSolve() is complete, then it returns the number of iterations
177:    carried out by the solution method.

179:    Not Collective

181:    Input Parameter:
182: .  pep - the polynomial eigensolver context

184:    Output Parameter:
185: .  its - number of iterations

187:    Note:
188:    During the i-th iteration this call returns i-1. If PEPSolve() is
189:    complete, then parameter "its" contains either the iteration number at
190:    which convergence was successfully reached, or failure was detected.
191:    Call PEPGetConvergedReason() to determine if the solver converged or
192:    failed and why.

194:    Level: intermediate

196: .seealso: PEPGetConvergedReason(), PEPSetTolerances()
197: @*/
198: PetscErrorCode PEPGetIterationNumber(PEP pep,PetscInt *its)
199: {
203:   *its = pep->its;
204:   return(0);
205: }

207: /*@
208:    PEPGetConverged - Gets the number of converged eigenpairs.

210:    Not Collective

212:    Input Parameter:
213: .  pep - the polynomial eigensolver context

215:    Output Parameter:
216: .  nconv - number of converged eigenpairs

218:    Note:
219:    This function should be called after PEPSolve() has finished.

221:    Level: beginner

223: .seealso: PEPSetDimensions(), PEPSolve()
224: @*/
225: PetscErrorCode PEPGetConverged(PEP pep,PetscInt *nconv)
226: {
230:   PEPCheckSolved(pep,1);
231:   *nconv = pep->nconv;
232:   return(0);
233: }

235: /*@
236:    PEPGetConvergedReason - Gets the reason why the PEPSolve() iteration was
237:    stopped.

239:    Not Collective

241:    Input Parameter:
242: .  pep - the polynomial eigensolver context

244:    Output Parameter:
245: .  reason - negative value indicates diverged, positive value converged

247:    Notes:

249:    Possible values for reason are
250: +  PEP_CONVERGED_TOL - converged up to tolerance
251: .  PEP_CONVERGED_USER - converged due to a user-defined condition
252: .  PEP_DIVERGED_ITS - required more than max_it iterations to reach convergence
253: .  PEP_DIVERGED_BREAKDOWN - generic breakdown in method
254: -  PEP_DIVERGED_SYMMETRY_LOST - pseudo-Lanczos was not able to keep symmetry

256:    Can only be called after the call to PEPSolve() is complete.

258:    Level: intermediate

260: .seealso: PEPSetTolerances(), PEPSolve(), PEPConvergedReason
261: @*/
262: PetscErrorCode PEPGetConvergedReason(PEP pep,PEPConvergedReason *reason)
263: {
267:   PEPCheckSolved(pep,1);
268:   *reason = pep->reason;
269:   return(0);
270: }

272: /*@C
273:    PEPGetEigenpair - Gets the i-th solution of the eigenproblem as computed by
274:    PEPSolve(). The solution consists in both the eigenvalue and the eigenvector.

276:    Logically Collective on pep

278:    Input Parameters:
279: +  pep - polynomial eigensolver context
280: -  i   - index of the solution

282:    Output Parameters:
283: +  eigr - real part of eigenvalue
284: .  eigi - imaginary part of eigenvalue
285: .  Vr   - real part of eigenvector
286: -  Vi   - imaginary part of eigenvector

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

293:    If the eigenvalue is real, then eigi and Vi are set to zero. If PETSc is
294:    configured with complex scalars the eigenvalue is stored
295:    directly in eigr (eigi is set to zero) and the eigenvector in Vr (Vi is
296:    set to zero). In any case, the user can pass NULL in Vr or Vi if one of
297:    them is not required.

299:    The index i should be a value between 0 and nconv-1 (see PEPGetConverged()).
300:    Eigenpairs are indexed according to the ordering criterion established
301:    with PEPSetWhichEigenpairs().

303:    Level: beginner

305: .seealso: PEPSolve(), PEPGetConverged(), PEPSetWhichEigenpairs()
306: @*/
307: PetscErrorCode PEPGetEigenpair(PEP pep,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
308: {
309:   PetscInt       k;

317:   PEPCheckSolved(pep,1);
318:   if (i<0) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
319:   if (i>=pep->nconv) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see PEPGetConverged()");

321:   PEPComputeVectors(pep);
322:   k = pep->perm[i];

324:   /* eigenvalue */
325: #if defined(PETSC_USE_COMPLEX)
326:   if (eigr) *eigr = pep->eigr[k];
327:   if (eigi) *eigi = 0;
328: #else
329:   if (eigr) *eigr = pep->eigr[k];
330:   if (eigi) *eigi = pep->eigi[k];
331: #endif

333:   /* eigenvector */
334:   BV_GetEigenvector(pep->V,k,pep->eigi[k],Vr,Vi);
335:   return(0);
336: }

338: /*@
339:    PEPGetErrorEstimate - Returns the error estimate associated to the i-th
340:    computed eigenpair.

342:    Not Collective

344:    Input Parameters:
345: +  pep - polynomial eigensolver context
346: -  i   - index of eigenpair

348:    Output Parameter:
349: .  errest - the error estimate

351:    Notes:
352:    This is the error estimate used internally by the eigensolver. The actual
353:    error bound can be computed with PEPComputeError(). See also the users
354:    manual for details.

356:    Level: advanced

358: .seealso: PEPComputeError()
359: @*/
360: PetscErrorCode PEPGetErrorEstimate(PEP pep,PetscInt i,PetscReal *errest)
361: {
365:   PEPCheckSolved(pep,1);
366:   if (i<0 || i>=pep->nconv) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
367:   *errest = pep->errest[pep->perm[i]];
368:   return(0);
369: }

371: /*
372:    PEPComputeResidualNorm_Private - Computes the norm of the residual vector
373:    associated with an eigenpair.

375:    Input Parameters:
376:      kr,ki - eigenvalue
377:      xr,xi - eigenvector
378:      z     - array of 4 work vectors (z[2],z[3] not referenced in complex scalars)
379: */
380: PetscErrorCode PEPComputeResidualNorm_Private(PEP pep,PetscScalar kr,PetscScalar ki,Vec xr,Vec xi,Vec *z,PetscReal *norm)
381: {
383:   Mat            *A=pep->A;
384:   PetscInt       i,nmat=pep->nmat;
385:   PetscScalar    t[20],*vals=t,*ivals=NULL;
386:   Vec            u,w;
387: #if !defined(PETSC_USE_COMPLEX)
388:   Vec            ui,wi;
389:   PetscReal      ni;
390:   PetscBool      imag;
391:   PetscScalar    it[20];
392: #endif

395:   u = z[0]; w = z[1];
396:   VecSet(u,0.0);
397: #if !defined(PETSC_USE_COMPLEX)
398:   ui = z[2]; wi = z[3];
399:   ivals = it;
400: #endif
401:   if (nmat>20) {
402:     PetscMalloc1(nmat,&vals);
403: #if !defined(PETSC_USE_COMPLEX)
404:     PetscMalloc1(nmat,&ivals);
405: #endif
406:   }
407:   PEPEvaluateBasis(pep,kr,ki,vals,ivals);
408: #if !defined(PETSC_USE_COMPLEX)
409:   if (ki == 0 || PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON))
410:     imag = PETSC_FALSE;
411:   else {
412:     imag = PETSC_TRUE;
413:     VecSet(ui,0.0);
414:   }
415: #endif
416:   for (i=0;i<nmat;i++) {
417:     if (vals[i]!=0.0) {
418:       MatMult(A[i],xr,w);
419:       VecAXPY(u,vals[i],w);
420:     }
421: #if !defined(PETSC_USE_COMPLEX)
422:     if (imag) {
423:       if (ivals[i]!=0 || vals[i]!=0) {
424:         MatMult(A[i],xi,wi);
425:         if (vals[i]==0) {
426:           MatMult(A[i],xr,w);
427:         }
428:       }
429:       if (ivals[i]!=0) {
430:         VecAXPY(u,-ivals[i],wi);
431:         VecAXPY(ui,ivals[i],w);
432:       }
433:       if (vals[i]!=0) {
434:         VecAXPY(ui,vals[i],wi);
435:       }
436:     }
437: #endif
438:   }
439:   VecNorm(u,NORM_2,norm);
440: #if !defined(PETSC_USE_COMPLEX)
441:   if (imag) {
442:     VecNorm(ui,NORM_2,&ni);
443:     *norm = SlepcAbsEigenvalue(*norm,ni);
444:   }
445: #endif
446:   if (nmat>20) {
447:     PetscFree(vals);
448: #if !defined(PETSC_USE_COMPLEX)
449:     PetscFree(ivals);
450: #endif
451:   }
452:   return(0);
453: }

455: /*@
456:    PEPComputeError - Computes the error (based on the residual norm) associated
457:    with the i-th computed eigenpair.

459:    Collective on pep

461:    Input Parameters:
462: +  pep  - the polynomial eigensolver context
463: .  i    - the solution index
464: -  type - the type of error to compute

466:    Output Parameter:
467: .  error - the error

469:    Notes:
470:    The error can be computed in various ways, all of them based on the residual
471:    norm ||P(l)x||_2 where l is the eigenvalue and x is the eigenvector.
472:    See the users guide for additional details.

474:    Level: beginner

476: .seealso: PEPErrorType, PEPSolve(), PEPGetErrorEstimate()
477: @*/
478: PetscErrorCode PEPComputeError(PEP pep,PetscInt i,PEPErrorType type,PetscReal *error)
479: {
481:   Vec            xr,xi,w[4];
482:   PetscScalar    kr,ki;
483:   PetscReal      t,z=0.0;
484:   PetscInt       j;
485:   PetscBool      flg;

492:   PEPCheckSolved(pep,1);

494:   /* allocate work vectors */
495: #if defined(PETSC_USE_COMPLEX)
496:   PEPSetWorkVecs(pep,3);
497:   xi   = NULL;
498:   w[2] = NULL;
499:   w[3] = NULL;
500: #else
501:   PEPSetWorkVecs(pep,6);
502:   xi   = pep->work[3];
503:   w[2] = pep->work[4];
504:   w[3] = pep->work[5];
505: #endif
506:   xr   = pep->work[0];
507:   w[0] = pep->work[1];
508:   w[1] = pep->work[2];

510:   /* compute residual norms */
511:   PEPGetEigenpair(pep,i,&kr,&ki,xr,xi);
512:   PEPComputeResidualNorm_Private(pep,kr,ki,xr,xi,w,error);

514:   /* compute error */
515:   switch (type) {
516:     case PEP_ERROR_ABSOLUTE:
517:       break;
518:     case PEP_ERROR_RELATIVE:
519:       *error /= SlepcAbsEigenvalue(kr,ki);
520:       break;
521:     case PEP_ERROR_BACKWARD:
522:       /* initialization of matrix norms */
523:       if (!pep->nrma[pep->nmat-1]) {
524:         for (j=0;j<pep->nmat;j++) {
525:           MatHasOperation(pep->A[j],MATOP_NORM,&flg);
526:           if (!flg) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
527:           MatNorm(pep->A[j],NORM_INFINITY,&pep->nrma[j]);
528:         }
529:       }
530:       t = SlepcAbsEigenvalue(kr,ki);
531:       for (j=pep->nmat-1;j>=0;j--) {
532:         z = z*t+pep->nrma[j];
533:       }
534:       *error /= z;
535:       break;
536:     default:
537:       SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
538:   }
539:   return(0);
540: }