Actual source code: pepsolve.c

slepc-3.10.1 2018-10-23
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2018, 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
 12: */

 14: #include <slepc/private/pepimpl.h>       /*I "slepcpep.h" I*/
 15: #include <petscdraw.h>

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

 30: PetscErrorCode PEPComputeVectors(PEP pep)
 31: {

 35:   PEPCheckSolved(pep,1);
 36:   if (pep->state==PEP_STATE_SOLVED && pep->ops->computevectors) {
 37:     (*pep->ops->computevectors)(pep);
 38:   }
 39:   pep->state = PEP_STATE_EIGENVECTORS;
 40:   return(0);
 41: }

 43: PetscErrorCode PEPExtractVectors(PEP pep)
 44: {

 48:   PEPCheckSolved(pep,1);
 49:   if (pep->state==PEP_STATE_SOLVED && pep->ops->extractvectors) {
 50:     (*pep->ops->extractvectors)(pep);
 51:   }
 52:   return(0);
 53: }

 55: /*@
 56:    PEPSolve - Solves the polynomial eigensystem.

 58:    Collective on PEP

 60:    Input Parameter:
 61: .  pep - eigensolver context obtained from PEPCreate()

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

 74:    Level: beginner

 76: .seealso: PEPCreate(), PEPSetUp(), PEPDestroy(), PEPSetTolerances()
 77: @*/
 78: PetscErrorCode PEPSolve(PEP pep)
 79: {
 81:   PetscInt       i,k;
 82:   PetscBool      flg,islinear;
 83: #define OPTLEN 16
 84:   char           str[OPTLEN];

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

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

105:   (*pep->ops->solve)(pep);

107:   if (!pep->reason) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");

109:   PetscObjectTypeCompare((PetscObject)pep,PEPLINEAR,&islinear);
110:   if (!islinear) {
111:     STPostSolve(pep->st);
112:     /* Map eigenvalues back to the original problem */
113:     STGetTransform(pep->st,&flg);
114:     if (flg && pep->ops->backtransform) {
115:       (*pep->ops->backtransform)(pep);
116:     }
117:   }

119:   pep->state = PEP_STATE_SOLVED;

121: #if !defined(PETSC_USE_COMPLEX)
122:   /* reorder conjugate eigenvalues (positive imaginary first) */
123:   for (i=0;i<pep->nconv-1;i++) {
124:     if (pep->eigi[i] != 0) {
125:       if (pep->eigi[i] < 0) {
126:         pep->eigi[i] = -pep->eigi[i];
127:         pep->eigi[i+1] = -pep->eigi[i+1];
128:         /* the next correction only works with eigenvectors */
129:         PEPComputeVectors(pep);
130:         BVScaleColumn(pep->V,i+1,-1.0);
131:       }
132:       i++;
133:     }
134:   }
135: #endif

137:   if (pep->refine!=PEP_REFINE_NONE) {
138:     PetscCitationsRegister(citation,&cited);
139:   }

141:   if (pep->refine==PEP_REFINE_SIMPLE && pep->rits>0 && pep->nconv>0) {
142:     PEPComputeVectors(pep);
143:     PEPNewtonRefinementSimple(pep,&pep->rits,pep->rtol,pep->nconv);
144:   }

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

150:   /* various viewers */
151:   PEPViewFromOptions(pep,NULL,"-pep_view");
152:   PEPReasonViewFromOptions(pep);
153:   PEPErrorViewFromOptions(pep);
154:   PEPValuesViewFromOptions(pep);
155:   PEPVectorsViewFromOptions(pep);
156:   for (i=0;i<pep->nmat;i++) {
157:     PetscSNPrintf(str,OPTLEN,"-pep_view_mat%d",(int)i);
158:     MatViewFromOptions(pep->A[i],(PetscObject)pep,str);
159:   }

161:   /* Remove the initial subspace */
162:   pep->nini = 0;
163:   return(0);
164: }

166: /*@
167:    PEPGetIterationNumber - Gets the current iteration number. If the
168:    call to PEPSolve() is complete, then it returns the number of iterations
169:    carried out by the solution method.

171:    Not Collective

173:    Input Parameter:
174: .  pep - the polynomial eigensolver context

176:    Output Parameter:
177: .  its - number of iterations

179:    Note:
180:    During the i-th iteration this call returns i-1. If PEPSolve() is
181:    complete, then parameter "its" contains either the iteration number at
182:    which convergence was successfully reached, or failure was detected.
183:    Call PEPGetConvergedReason() to determine if the solver converged or
184:    failed and why.

186:    Level: intermediate

188: .seealso: PEPGetConvergedReason(), PEPSetTolerances()
189: @*/
190: PetscErrorCode PEPGetIterationNumber(PEP pep,PetscInt *its)
191: {
195:   *its = pep->its;
196:   return(0);
197: }

199: /*@
200:    PEPGetConverged - Gets the number of converged eigenpairs.

202:    Not Collective

204:    Input Parameter:
205: .  pep - the polynomial eigensolver context

207:    Output Parameter:
208: .  nconv - number of converged eigenpairs

210:    Note:
211:    This function should be called after PEPSolve() has finished.

213:    Level: beginner

215: .seealso: PEPSetDimensions(), PEPSolve()
216: @*/
217: PetscErrorCode PEPGetConverged(PEP pep,PetscInt *nconv)
218: {
222:   PEPCheckSolved(pep,1);
223:   *nconv = pep->nconv;
224:   return(0);
225: }

227: /*@
228:    PEPGetConvergedReason - Gets the reason why the PEPSolve() iteration was
229:    stopped.

231:    Not Collective

233:    Input Parameter:
234: .  pep - the polynomial eigensolver context

236:    Output Parameter:
237: .  reason - negative value indicates diverged, positive value converged

239:    Notes:

241:    Possible values for reason are
242: +  PEP_CONVERGED_TOL - converged up to tolerance
243: .  PEP_CONVERGED_USER - converged due to a user-defined condition
244: .  PEP_DIVERGED_ITS - required more than max_it iterations to reach convergence
245: .  PEP_DIVERGED_BREAKDOWN - generic breakdown in method
246: -  PEP_DIVERGED_SYMMETRY_LOST - pseudo-Lanczos was not able to keep symmetry

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

250:    Level: intermediate

252: .seealso: PEPSetTolerances(), PEPSolve(), PEPConvergedReason
253: @*/
254: PetscErrorCode PEPGetConvergedReason(PEP pep,PEPConvergedReason *reason)
255: {
259:   PEPCheckSolved(pep,1);
260:   *reason = pep->reason;
261:   return(0);
262: }

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

268:    Logically Collective on EPS

270:    Input Parameters:
271: +  pep - polynomial eigensolver context
272: -  i   - index of the solution

274:    Output Parameters:
275: +  eigr - real part of eigenvalue
276: .  eigi - imaginary part of eigenvalue
277: .  Vr   - real part of eigenvector
278: -  Vi   - imaginary part of eigenvector

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

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

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

294:    Level: beginner

296: .seealso: PEPSolve(), PEPGetConverged(), PEPSetWhichEigenpairs()
297: @*/
298: PetscErrorCode PEPGetEigenpair(PEP pep,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
299: {
300:   PetscInt       k;

308:   PEPCheckSolved(pep,1);
309:   if (i<0 || i>=pep->nconv) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");

311:   PEPComputeVectors(pep);
312:   k = pep->perm[i];

314:   /* eigenvalue */
315: #if defined(PETSC_USE_COMPLEX)
316:   if (eigr) *eigr = pep->eigr[k];
317:   if (eigi) *eigi = 0;
318: #else
319:   if (eigr) *eigr = pep->eigr[k];
320:   if (eigi) *eigi = pep->eigi[k];
321: #endif

323:   /* eigenvector */
324: #if defined(PETSC_USE_COMPLEX)
325:   if (Vr) { BVCopyVec(pep->V,k,Vr); }
326:   if (Vi) { VecSet(Vi,0.0); }
327: #else
328:   if (pep->eigi[k]>0) { /* first value of conjugate pair */
329:     if (Vr) { BVCopyVec(pep->V,k,Vr); }
330:     if (Vi) { BVCopyVec(pep->V,k+1,Vi); }
331:   } else if (pep->eigi[k]<0) { /* second value of conjugate pair */
332:     if (Vr) { BVCopyVec(pep->V,k-1,Vr); }
333:     if (Vi) {
334:       BVCopyVec(pep->V,k,Vi);
335:       VecScale(Vi,-1.0);
336:     }
337:   } else { /* real eigenvalue */
338:     if (Vr) { BVCopyVec(pep->V,k,Vr); }
339:     if (Vi) { VecSet(Vi,0.0); }
340:   }
341: #endif
342:   return(0);
343: }

345: /*@
346:    PEPGetErrorEstimate - Returns the error estimate associated to the i-th
347:    computed eigenpair.

349:    Not Collective

351:    Input Parameter:
352: +  pep - polynomial eigensolver context
353: -  i   - index of eigenpair

355:    Output Parameter:
356: .  errest - the error estimate

358:    Notes:
359:    This is the error estimate used internally by the eigensolver. The actual
360:    error bound can be computed with PEPComputeError(). See also the users
361:    manual for details.

363:    Level: advanced

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

378: /*
379:    PEPComputeResidualNorm_Private - Computes the norm of the residual vector
380:    associated with an eigenpair.

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

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

462: /*@
463:    PEPComputeError - Computes the error (based on the residual norm) associated
464:    with the i-th computed eigenpair.

466:    Collective on PEP

468:    Input Parameter:
469: +  pep  - the polynomial eigensolver context
470: .  i    - the solution index
471: -  type - the type of error to compute

473:    Output Parameter:
474: .  error - the error

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

481:    Level: beginner

483: .seealso: PEPErrorType, PEPSolve(), PEPGetErrorEstimate()
484: @*/
485: PetscErrorCode PEPComputeError(PEP pep,PetscInt i,PEPErrorType type,PetscReal *error)
486: {
488:   Vec            xr,xi,w[4];
489:   PetscScalar    kr,ki;
490:   PetscReal      t,z=0.0;
491:   PetscInt       j;
492:   PetscBool      flg;

499:   PEPCheckSolved(pep,1);

501:   /* allocate work vectors */
502: #if defined(PETSC_USE_COMPLEX)
503:   PEPSetWorkVecs(pep,3);
504:   xi   = NULL;
505:   w[2] = NULL;
506:   w[3] = NULL;
507: #else
508:   PEPSetWorkVecs(pep,6);
509:   xi   = pep->work[3];
510:   w[2] = pep->work[4];
511:   w[3] = pep->work[5];
512: #endif
513:   xr   = pep->work[0];
514:   w[0] = pep->work[1];
515:   w[1] = pep->work[2];

517:   /* compute residual norms */
518:   PEPGetEigenpair(pep,i,&kr,&ki,xr,xi);
519:   PEPComputeResidualNorm_Private(pep,kr,ki,xr,xi,w,error);

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