Actual source code: pepdefault.c

slepc-3.14.2 2021-02-01
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2020, 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>

 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 SLEPC_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:   PEPSetWhichEigenpairs_Default - Sets the default value for which,
 90:   depending on the ST.
 91:  */
 92: PetscErrorCode PEPSetWhichEigenpairs_Default(PEP pep)
 93: {
 94:   PetscBool      target;

 98:   PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&target);
 99:   if (target) pep->which = PEP_TARGET_MAGNITUDE;
100:   else pep->which = PEP_LARGEST_MAGNITUDE;
101:   return(0);
102: }

104: /*
105:   PEPConvergedAbsolute - Checks convergence absolutely.
106: */
107: PetscErrorCode PEPConvergedAbsolute(PEP pep,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
108: {
110:   *errest = res;
111:   return(0);
112: }

114: /*@C
115:    PEPStoppingBasic - Default routine to determine whether the outer eigensolver
116:    iteration must be stopped.

118:    Collective on pep

120:    Input Parameters:
121: +  pep    - eigensolver context obtained from PEPCreate()
122: .  its    - current number of iterations
123: .  max_it - maximum number of iterations
124: .  nconv  - number of currently converged eigenpairs
125: .  nev    - number of requested eigenpairs
126: -  ctx    - context (not used here)

128:    Output Parameter:
129: .  reason - result of the stopping test

131:    Notes:
132:    A positive value of reason indicates that the iteration has finished successfully
133:    (converged), and a negative value indicates an error condition (diverged). If
134:    the iteration needs to be continued, reason must be set to PEP_CONVERGED_ITERATING
135:    (zero).

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

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

142:    Level: advanced

144: .seealso: PEPSetStoppingTest(), PEPConvergedReason, PEPGetConvergedReason()
145: @*/
146: PetscErrorCode PEPStoppingBasic(PEP pep,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nev,PEPConvergedReason *reason,void *ctx)
147: {

151:   *reason = PEP_CONVERGED_ITERATING;
152:   if (nconv >= nev) {
153:     PetscInfo2(pep,"Polynomial eigensolver finished successfully: %D eigenpairs converged at iteration %D\n",nconv,its);
154:     *reason = PEP_CONVERGED_TOL;
155:   } else if (its >= max_it) {
156:     *reason = PEP_DIVERGED_ITS;
157:     PetscInfo1(pep,"Polynomial eigensolver iteration reached maximum number of iterations (%D)\n",its);
158:   }
159:   return(0);
160: }

162: PetscErrorCode PEPBackTransform_Default(PEP pep)
163: {

167:   STBackTransform(pep->st,pep->nconv,pep->eigr,pep->eigi);
168:   return(0);
169: }

171: PetscErrorCode PEPComputeVectors_Default(PEP pep)
172: {
174:   PetscInt       i;
175:   Vec            v;

178:   PEPExtractVectors(pep);

180:   /* Fix eigenvectors if balancing was used */
181:   if ((pep->scale==PEP_SCALE_DIAGONAL || pep->scale==PEP_SCALE_BOTH) && pep->Dr && (pep->refine!=PEP_REFINE_MULTIPLE)) {
182:     for (i=0;i<pep->nconv;i++) {
183:       BVGetColumn(pep->V,i,&v);
184:       VecPointwiseMult(v,v,pep->Dr);
185:       BVRestoreColumn(pep->V,i,&v);
186:     }
187:   }

189:   /* normalization */
190:   BVNormalize(pep->V,pep->eigi);
191:   return(0);
192: }

194: /*
195:   PEPBuildDiagonalScaling - compute two diagonal matrices to be applied for balancing
196:   in polynomial eigenproblems.
197: */
198: PetscErrorCode PEPBuildDiagonalScaling(PEP pep)
199: {
201:   PetscInt       it,i,j,k,nmat,nr,e,nz,lst,lend,nc=0,*cols,emax,emin,emaxl,eminl;
202:   const PetscInt *cidx,*ridx;
203:   Mat            M,*T,A;
204:   PetscMPIInt    n;
205:   PetscBool      cont=PETSC_TRUE,flg=PETSC_FALSE;
206:   PetscScalar    *array,*Dr,*Dl,t;
207:   PetscReal      l2,d,*rsum,*aux,*csum,w=1.0;
208:   MatStructure   str;
209:   MatInfo        info;

212:   l2 = 2*PetscLogReal(2.0);
213:   nmat = pep->nmat;
214:   PetscMPIIntCast(pep->n,&n);
215:   STGetMatStructure(pep->st,&str);
216:   PetscMalloc1(nmat,&T);
217:   for (k=0;k<nmat;k++) {
218:     STGetMatrixTransformed(pep->st,k,&T[k]);
219:   }
220:   /* Form local auxiliar matrix M */
221:   PetscObjectTypeCompareAny((PetscObject)T[0],&cont,MATMPIAIJ,MATSEQAIJ,"");
222:   if (!cont) SETERRQ(PetscObjectComm((PetscObject)T[0]),PETSC_ERR_SUP,"Only for MPIAIJ or SEQAIJ matrix types");
223:   PetscObjectTypeCompare((PetscObject)T[0],MATMPIAIJ,&cont);
224:   if (cont) {
225:     MatMPIAIJGetLocalMat(T[0],MAT_INITIAL_MATRIX,&M);
226:     flg = PETSC_TRUE;
227:   } else {
228:     MatDuplicate(T[0],MAT_COPY_VALUES,&M);
229:   }
230:   MatGetInfo(M,MAT_LOCAL,&info);
231:   nz = (PetscInt)info.nz_used;
232:   MatSeqAIJGetArray(M,&array);
233:   for (i=0;i<nz;i++) {
234:     t = PetscAbsScalar(array[i]);
235:     array[i] = t*t;
236:   }
237:   MatSeqAIJRestoreArray(M,&array);
238:   for (k=1;k<nmat;k++) {
239:     if (flg) {
240:       MatMPIAIJGetLocalMat(T[k],MAT_INITIAL_MATRIX,&A);
241:     } else {
242:       if (str==SAME_NONZERO_PATTERN) {
243:         MatCopy(T[k],A,SAME_NONZERO_PATTERN);
244:       } else {
245:         MatDuplicate(T[k],MAT_COPY_VALUES,&A);
246:       }
247:     }
248:     MatGetInfo(A,MAT_LOCAL,&info);
249:     nz = (PetscInt)info.nz_used;
250:     MatSeqAIJGetArray(A,&array);
251:     for (i=0;i<nz;i++) {
252:       t = PetscAbsScalar(array[i]);
253:       array[i] = t*t;
254:     }
255:     MatSeqAIJRestoreArray(A,&array);
256:     w *= pep->slambda*pep->slambda*pep->sfactor;
257:     MatAXPY(M,w,A,str);
258:     if (flg || str!=SAME_NONZERO_PATTERN || k==nmat-2) {
259:       MatDestroy(&A);
260:     }
261:   }
262:   MatGetRowIJ(M,0,PETSC_FALSE,PETSC_FALSE,&nr,&ridx,&cidx,&cont);
263:   if (!cont) SETERRQ(PetscObjectComm((PetscObject)T[0]),PETSC_ERR_SUP,"It is not possible to compute scaling diagonals for these PEP matrices");
264:   MatGetInfo(M,MAT_LOCAL,&info);
265:   nz = (PetscInt)info.nz_used;
266:   VecGetOwnershipRange(pep->Dl,&lst,&lend);
267:   PetscMalloc4(nr,&rsum,pep->n,&csum,pep->n,&aux,PetscMin(pep->n-lend+lst,nz),&cols);
268:   VecSet(pep->Dr,1.0);
269:   VecSet(pep->Dl,1.0);
270:   VecGetArray(pep->Dl,&Dl);
271:   VecGetArray(pep->Dr,&Dr);
272:   MatSeqAIJGetArray(M,&array);
273:   PetscArrayzero(aux,pep->n);
274:   for (j=0;j<nz;j++) {
275:     /* Search non-zero columns outsize lst-lend */
276:     if (aux[cidx[j]]==0 && (cidx[j]<lst || lend<=cidx[j])) cols[nc++] = cidx[j];
277:     /* Local column sums */
278:     aux[cidx[j]] += PetscAbsScalar(array[j]);
279:   }
280:   for (it=0;it<pep->sits && cont;it++) {
281:     emaxl = 0; eminl = 0;
282:     /* Column sum  */
283:     if (it>0) { /* it=0 has been already done*/
284:       MatSeqAIJGetArray(M,&array);
285:       PetscArrayzero(aux,pep->n);
286:       for (j=0;j<nz;j++) aux[cidx[j]] += PetscAbsScalar(array[j]);
287:       MatSeqAIJRestoreArray(M,&array);
288:     }
289:     MPIU_Allreduce(aux,csum,n,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)pep->Dr));
290:     /* Update Dr */
291:     for (j=lst;j<lend;j++) {
292:       d = PetscLogReal(csum[j])/l2;
293:       e = -(PetscInt)((d < 0)?(d-0.5):(d+0.5));
294:       d = PetscPowReal(2.0,e);
295:       Dr[j-lst] *= d;
296:       aux[j] = d*d;
297:       emaxl = PetscMax(emaxl,e);
298:       eminl = PetscMin(eminl,e);
299:     }
300:     for (j=0;j<nc;j++) {
301:       d = PetscLogReal(csum[cols[j]])/l2;
302:       e = -(PetscInt)((d < 0)?(d-0.5):(d+0.5));
303:       d = PetscPowReal(2.0,e);
304:       aux[cols[j]] = d*d;
305:       emaxl = PetscMax(emaxl,e);
306:       eminl = PetscMin(eminl,e);
307:     }
308:     /* Scale M */
309:     MatSeqAIJGetArray(M,&array);
310:     for (j=0;j<nz;j++) {
311:       array[j] *= aux[cidx[j]];
312:     }
313:     MatSeqAIJRestoreArray(M,&array);
314:     /* Row sum */
315:     PetscArrayzero(rsum,nr);
316:     MatSeqAIJGetArray(M,&array);
317:     for (i=0;i<nr;i++) {
318:       for (j=ridx[i];j<ridx[i+1];j++) rsum[i] += PetscAbsScalar(array[j]);
319:       /* Update Dl */
320:       d = PetscLogReal(rsum[i])/l2;
321:       e = -(PetscInt)((d < 0)?(d-0.5):(d+0.5));
322:       d = PetscPowReal(2.0,e);
323:       Dl[i] *= d;
324:       /* Scale M */
325:       for (j=ridx[i];j<ridx[i+1];j++) array[j] *= d*d;
326:       emaxl = PetscMax(emaxl,e);
327:       eminl = PetscMin(eminl,e);
328:     }
329:     MatSeqAIJRestoreArray(M,&array);
330:     /* Compute global max and min */
331:     MPIU_Allreduce(&emaxl,&emax,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)pep->Dl));
332:     MPIU_Allreduce(&eminl,&emin,1,MPIU_INT,MPI_MIN,PetscObjectComm((PetscObject)pep->Dl));
333:     if (emax<=emin+2) cont = PETSC_FALSE;
334:   }
335:   VecRestoreArray(pep->Dr,&Dr);
336:   VecRestoreArray(pep->Dl,&Dl);
337:   /* Free memory*/
338:   MatDestroy(&M);
339:   PetscFree4(rsum,csum,aux,cols);
340:   PetscFree(T);
341:   return(0);
342: }

344: /*
345:    PEPComputeScaleFactor - compute sfactor as described in [Betcke 2008].
346: */
347: PetscErrorCode PEPComputeScaleFactor(PEP pep)
348: {
350:   PetscBool      has0,has1,flg;
351:   PetscReal      norm0,norm1;
352:   Mat            T[2];
353:   PEPBasis       basis;
354:   PetscInt       i;

357:   if (pep->scale==PEP_SCALE_NONE || pep->scale==PEP_SCALE_DIAGONAL) {  /* no scalar scaling */
358:     pep->sfactor = 1.0;
359:     pep->dsfactor = 1.0;
360:     return(0);
361:   }
362:   if (pep->sfactor_set) return(0);  /* user provided value */
363:   pep->sfactor = 1.0;
364:   pep->dsfactor = 1.0;
365:   PEPGetBasis(pep,&basis);
366:   if (basis==PEP_BASIS_MONOMIAL) {
367:     STGetTransform(pep->st,&flg);
368:     if (flg) {
369:       STGetMatrixTransformed(pep->st,0,&T[0]);
370:       STGetMatrixTransformed(pep->st,pep->nmat-1,&T[1]);
371:     } else {
372:       T[0] = pep->A[0];
373:       T[1] = pep->A[pep->nmat-1];
374:     }
375:     if (pep->nmat>2) {
376:       MatHasOperation(T[0],MATOP_NORM,&has0);
377:       MatHasOperation(T[1],MATOP_NORM,&has1);
378:       if (has0 && has1) {
379:         MatNorm(T[0],NORM_INFINITY,&norm0);
380:         MatNorm(T[1],NORM_INFINITY,&norm1);
381:         pep->sfactor = PetscPowReal(norm0/norm1,1.0/(pep->nmat-1));
382:         pep->dsfactor = norm1;
383:         for (i=pep->nmat-2;i>0;i--) {
384:           STGetMatrixTransformed(pep->st,i,&T[1]);
385:           MatHasOperation(T[1],MATOP_NORM,&has1);
386:           if (has1) {
387:             MatNorm(T[1],NORM_INFINITY,&norm1);
388:             pep->dsfactor = pep->dsfactor*pep->sfactor+norm1;
389:           } else break;
390:         }
391:         if (has1) {
392:           pep->dsfactor = pep->dsfactor*pep->sfactor+norm0;
393:           pep->dsfactor = pep->nmat/pep->dsfactor;
394:         } else pep->dsfactor = 1.0;
395:       }
396:     }
397:   }
398:   return(0);
399: }

401: /*
402:    PEPBasisCoefficients - compute polynomial basis coefficients
403: */
404: PetscErrorCode PEPBasisCoefficients(PEP pep,PetscReal *pbc)
405: {
406:   PetscReal *ca,*cb,*cg;
407:   PetscInt  k,nmat=pep->nmat;

410:   ca = pbc;
411:   cb = pbc+nmat;
412:   cg = pbc+2*nmat;
413:   switch (pep->basis) {
414:   case PEP_BASIS_MONOMIAL:
415:     for (k=0;k<nmat;k++) {
416:       ca[k] = 1.0; cb[k] = 0.0; cg[k] = 0.0;
417:     }
418:     break;
419:   case PEP_BASIS_CHEBYSHEV1:
420:     ca[0] = 1.0; cb[0] = 0.0; cg[0] = 0.0;
421:     for (k=1;k<nmat;k++) {
422:       ca[k] = .5; cb[k] = 0.0; cg[k] = .5;
423:     }
424:     break;
425:   case PEP_BASIS_CHEBYSHEV2:
426:     ca[0] = .5; cb[0] = 0.0; cg[0] = 0.0;
427:     for (k=1;k<nmat;k++) {
428:       ca[k] = .5; cb[k] = 0.0; cg[k] = .5;
429:     }
430:     break;
431:   case PEP_BASIS_LEGENDRE:
432:     ca[0] = 1.0; cb[0] = 0.0; cg[0] = 0.0;
433:     for (k=1;k<nmat;k++) {
434:       ca[k] = k+1; cb[k] = -2*k; cg[k] = k;
435:     }
436:     break;
437:   case PEP_BASIS_LAGUERRE:
438:     ca[0] = -1.0; cb[0] = 0.0; cg[0] = 0.0;
439:     for (k=1;k<nmat;k++) {
440:       ca[k] = -(k+1); cb[k] = 2*k+1; cg[k] = -k;
441:     }
442:     break;
443:   case PEP_BASIS_HERMITE:
444:     ca[0] = .5; cb[0] = 0.0; cg[0] = 0.0;
445:     for (k=1;k<nmat;k++) {
446:       ca[k] = .5; cb[k] = 0.0; cg[k] = -k;
447:     }
448:     break;
449:   }
450:   return(0);
451: }