Actual source code: pepdefault.c

slepc-3.10.2 2019-02-11
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:    Simple default routines for common PEP operations
 12: */

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

 16: /*@
 17:    PEPSetWorkVecs - Sets a number of work vectors into a PEP object.

 19:    Collective on PEP

 21:    Input Parameters:
 22: +  pep - polynomial eigensolver context
 23: -  nw  - number of work vectors to allocate

 25:    Developers Note:
 26:    This is PETSC_EXTERN because it may be required by user plugin PEP
 27:    implementations.

 29:    Level: developer
 30: @*/
 31: PetscErrorCode PEPSetWorkVecs(PEP pep,PetscInt nw)
 32: {
 34:   Vec            t;

 37:   if (pep->nwork < nw) {
 38:     VecDestroyVecs(pep->nwork,&pep->work);
 39:     pep->nwork = nw;
 40:     BVGetColumn(pep->V,0,&t);
 41:     VecDuplicateVecs(t,nw,&pep->work);
 42:     BVRestoreColumn(pep->V,0,&t);
 43:     PetscLogObjectParents(pep,nw,pep->work);
 44:   }
 45:   return(0);
 46: }

 48: /*
 49:   PEPConvergedRelative - Checks convergence relative to the eigenvalue.
 50: */
 51: PetscErrorCode PEPConvergedRelative(PEP pep,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
 52: {
 53:   PetscReal w;

 56:   w = SlepcAbsEigenvalue(eigr,eigi);
 57:   *errest = res/w;
 58:   return(0);
 59: }

 61: /*
 62:   PEPConvergedNorm - Checks convergence relative to the matrix norms.
 63: */
 64: PetscErrorCode PEPConvergedNorm(PEP pep,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
 65: {
 66:   PetscReal      w=0.0,t;
 67:   PetscInt       j;
 68:   PetscBool      flg;

 72:   /* initialization of matrix norms */
 73:   if (!pep->nrma[pep->nmat-1]) {
 74:     for (j=0;j<pep->nmat;j++) {
 75:       MatHasOperation(pep->A[j],MATOP_NORM,&flg);
 76:       if (!flg) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"The convergence test related to the matrix norms requires a matrix norm operation");
 77:       MatNorm(pep->A[j],NORM_INFINITY,&pep->nrma[j]);
 78:     }
 79:   }
 80:   t = SlepcAbsEigenvalue(eigr,eigi);
 81:   for (j=pep->nmat-1;j>=0;j--) {
 82:     w = w*t+pep->nrma[j];
 83:   }
 84:   *errest = res/w;
 85:   return(0);
 86: }

 88: /*
 89:   PEPConvergedAbsolute - Checks convergence absolutely.
 90: */
 91: PetscErrorCode PEPConvergedAbsolute(PEP pep,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
 92: {
 94:   *errest = res;
 95:   return(0);
 96: }

 98: /*@C
 99:    PEPStoppingBasic - Default routine to determine whether the outer eigensolver
100:    iteration must be stopped.

102:    Collective on PEP

104:    Input Parameters:
105: +  pep    - eigensolver context obtained from PEPCreate()
106: .  its    - current number of iterations
107: .  max_it - maximum number of iterations
108: .  nconv  - number of currently converged eigenpairs
109: .  nev    - number of requested eigenpairs
110: -  ctx    - context (not used here)

112:    Output Parameter:
113: .  reason - result of the stopping test

115:    Notes:
116:    A positive value of reason indicates that the iteration has finished successfully
117:    (converged), and a negative value indicates an error condition (diverged). If
118:    the iteration needs to be continued, reason must be set to PEP_CONVERGED_ITERATING
119:    (zero).

121:    PEPStoppingBasic() will stop if all requested eigenvalues are converged, or if
122:    the maximum number of iterations has been reached.

124:    Use PEPSetStoppingTest() to provide your own test instead of using this one.

126:    Level: advanced

128: .seealso: PEPSetStoppingTest(), PEPConvergedReason, PEPGetConvergedReason()
129: @*/
130: PetscErrorCode PEPStoppingBasic(PEP pep,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nev,PEPConvergedReason *reason,void *ctx)
131: {

135:   *reason = PEP_CONVERGED_ITERATING;
136:   if (nconv >= nev) {
137:     PetscInfo2(pep,"Polynomial eigensolver finished successfully: %D eigenpairs converged at iteration %D\n",nconv,its);
138:     *reason = PEP_CONVERGED_TOL;
139:   } else if (its >= max_it) {
140:     *reason = PEP_DIVERGED_ITS;
141:     PetscInfo1(pep,"Polynomial eigensolver iteration reached maximum number of iterations (%D)\n",its);
142:   }
143:   return(0);
144: }

146: PetscErrorCode PEPBackTransform_Default(PEP pep)
147: {

151:   STBackTransform(pep->st,pep->nconv,pep->eigr,pep->eigi);
152:   return(0);
153: }

155: PetscErrorCode PEPComputeVectors_Default(PEP pep)
156: {
158:   PetscInt       i;
159:   Vec            v;
160: #if !defined(PETSC_USE_COMPLEX)
161:   Vec            v1;
162: #endif

165:   PEPExtractVectors(pep);

167:   /* Fix eigenvectors if balancing was used */
168:   if ((pep->scale==PEP_SCALE_DIAGONAL || pep->scale==PEP_SCALE_BOTH) && pep->Dr && (pep->refine!=PEP_REFINE_MULTIPLE)) {
169:     for (i=0;i<pep->nconv;i++) {
170:       BVGetColumn(pep->V,i,&v);
171:       VecPointwiseMult(v,v,pep->Dr);
172:       BVRestoreColumn(pep->V,i,&v);
173:     }
174:   }

176:   /* normalization */
177:   for (i=0;i<pep->nconv;i++) {
178: #if !defined(PETSC_USE_COMPLEX)
179:     if (pep->eigi[i]!=0.0) {   /* first eigenvalue of a complex conjugate pair */
180:       BVGetColumn(pep->V,i,&v);
181:       BVGetColumn(pep->V,i+1,&v1);
182:       VecNormalizeComplex(v,v1,PETSC_TRUE,NULL);
183:       BVRestoreColumn(pep->V,i,&v);
184:       BVRestoreColumn(pep->V,i+1,&v1);
185:       i++;
186:     } else   /* real eigenvalue */
187: #endif
188:     {
189:       BVGetColumn(pep->V,i,&v);
190:       VecNormalizeComplex(v,NULL,PETSC_FALSE,NULL);
191:       BVRestoreColumn(pep->V,i,&v);
192:     }
193:   }
194:   return(0);
195: }

197: /*
198:    PEPKrylovConvergence - This is the analogue to EPSKrylovConvergence, but
199:    for polynomial Krylov methods.

201:    Differences:
202:    - Always non-symmetric
203:    - Does not check for STSHIFT
204:    - No correction factor
205:    - No support for true residual
206: */
207: PetscErrorCode PEPKrylovConvergence(PEP pep,PetscBool getall,PetscInt kini,PetscInt nits,PetscReal beta,PetscInt *kout)
208: {
210:   PetscInt       k,newk,marker,inside;
211:   PetscScalar    re,im;
212:   PetscReal      resnorm;
213:   PetscBool      istrivial;

216:   RGIsTrivial(pep->rg,&istrivial);
217:   marker = -1;
218:   if (pep->trackall) getall = PETSC_TRUE;
219:   for (k=kini;k<kini+nits;k++) {
220:     /* eigenvalue */
221:     re = pep->eigr[k];
222:     im = pep->eigi[k];
223:     if (!istrivial) {
224:       STBackTransform(pep->st,1,&re,&im);
225:       RGCheckInside(pep->rg,1,&re,&im,&inside);
226:       if (marker==-1 && inside<0) marker = k;
227:       re = pep->eigr[k];
228:       im = pep->eigi[k];
229:     }
230:     newk = k;
231:     DSVectors(pep->ds,DS_MAT_X,&newk,&resnorm);
232:     resnorm *= beta;
233:     /* error estimate */
234:     (*pep->converged)(pep,re,im,resnorm,&pep->errest[k],pep->convergedctx);
235:     if (marker==-1 && pep->errest[k] >= pep->tol) marker = k;
236:     if (newk==k+1) {
237:       pep->errest[k+1] = pep->errest[k];
238:       k++;
239:     }
240:     if (marker!=-1 && !getall) break;
241:   }
242:   if (marker!=-1) k = marker;
243:   *kout = k;
244:   return(0);
245: }

247: /*
248:   PEPBuildDiagonalScaling - compute two diagonal matrices to be applied for balancing
249:   in polynomial eigenproblems.
250: */
251: PetscErrorCode PEPBuildDiagonalScaling(PEP pep)
252: {
254:   PetscInt       it,i,j,k,nmat,nr,e,nz,lst,lend,nc=0,*cols,emax,emin,emaxl,eminl;
255:   const PetscInt *cidx,*ridx;
256:   Mat            M,*T,A;
257:   PetscMPIInt    n;
258:   PetscBool      cont=PETSC_TRUE,flg=PETSC_FALSE;
259:   PetscScalar    *array,*Dr,*Dl,t;
260:   PetscReal      l2,d,*rsum,*aux,*csum,w=1.0;
261:   MatStructure   str;
262:   MatInfo        info;

265:   l2 = 2*PetscLogReal(2.0);
266:   nmat = pep->nmat;
267:   PetscMPIIntCast(pep->n,&n);
268:   STGetMatStructure(pep->st,&str);
269:   PetscMalloc1(nmat,&T);
270:   for (k=0;k<nmat;k++) {
271:     STGetMatrixTransformed(pep->st,k,&T[k]);
272:   }
273:   /* Form local auxiliar matrix M */
274:   PetscObjectTypeCompareAny((PetscObject)T[0],&cont,MATMPIAIJ,MATSEQAIJ,"");
275:   if (!cont) SETERRQ(PetscObjectComm((PetscObject)T[0]),PETSC_ERR_SUP,"Only for MPIAIJ or SEQAIJ matrix types");
276:   PetscObjectTypeCompare((PetscObject)T[0],MATMPIAIJ,&cont);
277:   if (cont) {
278:     MatMPIAIJGetLocalMat(T[0],MAT_INITIAL_MATRIX,&M);
279:     flg = PETSC_TRUE;
280:   } else {
281:     MatDuplicate(T[0],MAT_COPY_VALUES,&M);
282:   }
283:   MatGetInfo(M,MAT_LOCAL,&info);
284:   nz = (PetscInt)info.nz_used;
285:   MatSeqAIJGetArray(M,&array);
286:   for (i=0;i<nz;i++) {
287:     t = PetscAbsScalar(array[i]);
288:     array[i] = t*t;
289:   }
290:   MatSeqAIJRestoreArray(M,&array);
291:   for (k=1;k<nmat;k++) {
292:     if (flg) {
293:       MatMPIAIJGetLocalMat(T[k],MAT_INITIAL_MATRIX,&A);
294:     } else {
295:       if (str==SAME_NONZERO_PATTERN) {
296:         MatCopy(T[k],A,SAME_NONZERO_PATTERN);
297:       } else {
298:         MatDuplicate(T[k],MAT_COPY_VALUES,&A);
299:       }
300:     }
301:     MatGetInfo(A,MAT_LOCAL,&info);
302:     nz = (PetscInt)info.nz_used;
303:     MatSeqAIJGetArray(A,&array);
304:     for (i=0;i<nz;i++) {
305:       t = PetscAbsScalar(array[i]);
306:       array[i] = t*t;
307:     }
308:     MatSeqAIJRestoreArray(A,&array);
309:     w *= pep->slambda*pep->slambda*pep->sfactor;
310:     MatAXPY(M,w,A,str);
311:     if (flg || str!=SAME_NONZERO_PATTERN || k==nmat-2) {
312:       MatDestroy(&A);
313:     }
314:   }
315:   MatGetRowIJ(M,0,PETSC_FALSE,PETSC_FALSE,&nr,&ridx,&cidx,&cont);
316:   if (!cont) SETERRQ(PetscObjectComm((PetscObject)T[0]),PETSC_ERR_SUP,"It is not possible to compute scaling diagonals for these PEP matrices");
317:   MatGetInfo(M,MAT_LOCAL,&info);
318:   nz = (PetscInt)info.nz_used;
319:   VecGetOwnershipRange(pep->Dl,&lst,&lend);
320:   PetscMalloc4(nr,&rsum,pep->n,&csum,pep->n,&aux,PetscMin(pep->n-lend+lst,nz),&cols);
321:   VecSet(pep->Dr,1.0);
322:   VecSet(pep->Dl,1.0);
323:   VecGetArray(pep->Dl,&Dl);
324:   VecGetArray(pep->Dr,&Dr);
325:   MatSeqAIJGetArray(M,&array);
326:   PetscMemzero(aux,pep->n*sizeof(PetscReal));
327:   for (j=0;j<nz;j++) {
328:     /* Search non-zero columns outsize lst-lend */
329:     if (aux[cidx[j]]==0 && (cidx[j]<lst || lend<=cidx[j])) cols[nc++] = cidx[j];
330:     /* Local column sums */
331:     aux[cidx[j]] += PetscAbsScalar(array[j]);
332:   }
333:   for (it=0;it<pep->sits && cont;it++) {
334:     emaxl = 0; eminl = 0;
335:     /* Column sum  */
336:     if (it>0) { /* it=0 has been already done*/
337:       MatSeqAIJGetArray(M,&array);
338:       PetscMemzero(aux,pep->n*sizeof(PetscReal));
339:       for (j=0;j<nz;j++) aux[cidx[j]] += PetscAbsScalar(array[j]);
340:       MatSeqAIJRestoreArray(M,&array);
341:     }
342:     MPI_Allreduce(aux,csum,n,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)pep->Dr));
343:     /* Update Dr */
344:     for (j=lst;j<lend;j++) {
345:       d = PetscLogReal(csum[j])/l2;
346:       e = -(PetscInt)((d < 0)?(d-0.5):(d+0.5));
347:       d = PetscPowReal(2.0,e);
348:       Dr[j-lst] *= d;
349:       aux[j] = d*d;
350:       emaxl = PetscMax(emaxl,e);
351:       eminl = PetscMin(eminl,e);
352:     }
353:     for (j=0;j<nc;j++) {
354:       d = PetscLogReal(csum[cols[j]])/l2;
355:       e = -(PetscInt)((d < 0)?(d-0.5):(d+0.5));
356:       d = PetscPowReal(2.0,e);
357:       aux[cols[j]] = d*d;
358:       emaxl = PetscMax(emaxl,e);
359:       eminl = PetscMin(eminl,e);
360:     }
361:     /* Scale M */
362:     MatSeqAIJGetArray(M,&array);
363:     for (j=0;j<nz;j++) {
364:       array[j] *= aux[cidx[j]];
365:     }
366:     MatSeqAIJRestoreArray(M,&array);
367:     /* Row sum */
368:     PetscMemzero(rsum,nr*sizeof(PetscReal));
369:     MatSeqAIJGetArray(M,&array);
370:     for (i=0;i<nr;i++) {
371:       for (j=ridx[i];j<ridx[i+1];j++) rsum[i] += PetscAbsScalar(array[j]);
372:       /* Update Dl */
373:       d = PetscLogReal(rsum[i])/l2;
374:       e = -(PetscInt)((d < 0)?(d-0.5):(d+0.5));
375:       d = PetscPowReal(2.0,e);
376:       Dl[i] *= d;
377:       /* Scale M */
378:       for (j=ridx[i];j<ridx[i+1];j++) array[j] *= d*d;
379:       emaxl = PetscMax(emaxl,e);
380:       eminl = PetscMin(eminl,e);
381:     }
382:     MatSeqAIJRestoreArray(M,&array);
383:     /* Compute global max and min */
384:     MPI_Allreduce(&emaxl,&emax,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)pep->Dl));
385:     MPI_Allreduce(&eminl,&emin,1,MPIU_INT,MPI_MIN,PetscObjectComm((PetscObject)pep->Dl));
386:     if (emax<=emin+2) cont = PETSC_FALSE;
387:   }
388:   VecRestoreArray(pep->Dr,&Dr);
389:   VecRestoreArray(pep->Dl,&Dl);
390:   /* Free memory*/
391:   MatDestroy(&M);
392:   PetscFree4(rsum,csum,aux,cols);
393:   PetscFree(T);
394:   return(0);
395: }

397: /*
398:    PEPComputeScaleFactor - compute sfactor as described in [Betcke 2008].
399: */
400: PetscErrorCode PEPComputeScaleFactor(PEP pep)
401: {
403:   PetscBool      has0,has1,flg;
404:   PetscReal      norm0,norm1;
405:   Mat            T[2];
406:   PEPBasis       basis;
407:   PetscInt       i;

410:   if (pep->scale==PEP_SCALE_NONE || pep->scale==PEP_SCALE_DIAGONAL) {  /* no scalar scaling */
411:     pep->sfactor = 1.0;
412:     pep->dsfactor = 1.0;
413:     return(0);
414:   }
415:   if (pep->sfactor_set) return(0);  /* user provided value */
416:   pep->sfactor = 1.0;
417:   pep->dsfactor = 1.0;
418:   PEPGetBasis(pep,&basis);
419:   if (basis==PEP_BASIS_MONOMIAL) {
420:     STGetTransform(pep->st,&flg);
421:     if (flg) {
422:       STGetMatrixTransformed(pep->st,0,&T[0]);
423:       STGetMatrixTransformed(pep->st,pep->nmat-1,&T[1]);
424:     } else {
425:       T[0] = pep->A[0];
426:       T[1] = pep->A[pep->nmat-1];
427:     }
428:     if (pep->nmat>2) {
429:       MatHasOperation(T[0],MATOP_NORM,&has0);
430:       MatHasOperation(T[1],MATOP_NORM,&has1);
431:       if (has0 && has1) {
432:         MatNorm(T[0],NORM_INFINITY,&norm0);
433:         MatNorm(T[1],NORM_INFINITY,&norm1);
434:         pep->sfactor = PetscPowReal(norm0/norm1,1.0/(pep->nmat-1));
435:         pep->dsfactor = norm1;
436:         for (i=pep->nmat-2;i>0;i--) {
437:           STGetMatrixTransformed(pep->st,i,&T[1]);
438:           MatHasOperation(T[1],MATOP_NORM,&has1);
439:           if (has1) {
440:             MatNorm(T[1],NORM_INFINITY,&norm1);
441:             pep->dsfactor = pep->dsfactor*pep->sfactor+norm1;
442:           } else break;
443:         }
444:         if (has1) {
445:           pep->dsfactor = pep->dsfactor*pep->sfactor+norm0;
446:           pep->dsfactor = pep->nmat/pep->dsfactor;
447:         } else pep->dsfactor = 1.0;
448:       }
449:     }
450:   }
451:   return(0);
452: }

454: /*
455:    PEPBasisCoefficients - compute polynomial basis coefficients
456: */
457: PetscErrorCode PEPBasisCoefficients(PEP pep,PetscReal *pbc)
458: {
459:   PetscReal *ca,*cb,*cg;
460:   PetscInt  k,nmat=pep->nmat;

463:   ca = pbc;
464:   cb = pbc+nmat;
465:   cg = pbc+2*nmat;
466:   switch (pep->basis) {
467:   case PEP_BASIS_MONOMIAL:
468:     for (k=0;k<nmat;k++) {
469:       ca[k] = 1.0; cb[k] = 0.0; cg[k] = 0.0;
470:     }
471:     break;
472:   case PEP_BASIS_CHEBYSHEV1:
473:     ca[0] = 1.0; cb[0] = 0.0; cg[0] = 0.0;
474:     for (k=1;k<nmat;k++) {
475:       ca[k] = .5; cb[k] = 0.0; cg[k] = .5;
476:     }
477:     break;
478:   case PEP_BASIS_CHEBYSHEV2:
479:     ca[0] = .5; cb[0] = 0.0; cg[0] = 0.0;
480:     for (k=1;k<nmat;k++) {
481:       ca[k] = .5; cb[k] = 0.0; cg[k] = .5;
482:     }
483:     break;
484:   case PEP_BASIS_LEGENDRE:
485:     ca[0] = 1.0; cb[0] = 0.0; cg[0] = 0.0;
486:     for (k=1;k<nmat;k++) {
487:       ca[k] = k+1; cb[k] = -2*k; cg[k] = k;
488:     }
489:     break;
490:   case PEP_BASIS_LAGUERRE:
491:     ca[0] = -1.0; cb[0] = 0.0; cg[0] = 0.0;
492:     for (k=1;k<nmat;k++) {
493:       ca[k] = -(k+1); cb[k] = 2*k+1; cg[k] = -k;
494:     }
495:     break;
496:   case PEP_BASIS_HERMITE:
497:     ca[0] = .5; cb[0] = 0.0; cg[0] = 0.0;
498:     for (k=1;k<nmat;k++) {
499:       ca[k] = .5; cb[k] = 0.0; cg[k] = -k;
500:     }
501:     break;
502:   }
503:   return(0);
504: }