Actual source code: subspace.c

      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:    SLEPc eigensolver: "subspace"

     13:    Method: Subspace Iteration

     15:    Algorithm:

     17:        Subspace iteration with Rayleigh-Ritz projection and locking,
     18:        based on the SRRIT implementation.

     20:    References:

     22:        [1] "Subspace Iteration in SLEPc", SLEPc Technical Report STR-3,
     23:            available at https://slepc.upv.es.
     24: */

     26: #include <slepc/private/epsimpl.h>

     28: typedef struct {
     29:   PetscBool estimatedrange;     /* the filter range was not set by the user */
     30: } EPS_SUBSPACE;

     32: static PetscErrorCode EPSSetUp_Subspace_Filter(EPS eps)
     33: {
     34:   EPS_SUBSPACE   *ctx = (EPS_SUBSPACE*)eps->data;
     35:   PetscBool      estimaterange=PETSC_TRUE;
     36:   PetscReal      rleft,rright;
     37:   Mat            A;

     39:   PetscFunctionBegin;
     40:   EPSCheckHermitianCondition(eps,PETSC_TRUE," with polynomial filter");
     41:   EPSCheckStandardCondition(eps,PETSC_TRUE," with polynomial filter");
     42:   PetscCheck(eps->intb<PETSC_MAX_REAL || eps->inta>PETSC_MIN_REAL,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The defined computational interval should have at least one of their sides bounded");
     43:   EPSCheckUnsupportedCondition(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION,PETSC_TRUE," with polynomial filter");
     44:   PetscCall(STFilterSetInterval(eps->st,eps->inta,eps->intb));
     45:   if (!ctx->estimatedrange) {
     46:     PetscCall(STFilterGetRange(eps->st,&rleft,&rright));
     47:     estimaterange = (!rleft && !rright)? PETSC_TRUE: PETSC_FALSE;
     48:   }
     49:   if (estimaterange) { /* user did not set a range */
     50:     PetscCall(STGetMatrix(eps->st,0,&A));
     51:     PetscCall(MatEstimateSpectralRange_EPS(A,&rleft,&rright));
     52:     PetscCall(PetscInfo(eps,"Setting eigenvalue range to [%g,%g]\n",(double)rleft,(double)rright));
     53:     PetscCall(STFilterSetRange(eps->st,rleft,rright));
     54:     ctx->estimatedrange = PETSC_TRUE;
     55:   }
     56:   PetscCall(EPSSetStoppingTest(eps,EPS_STOP_THRESHOLD));
     57:   PetscCall(EPSSetDimensions_Default(eps,&eps->nev,&eps->ncv,&eps->mpd));
     58:   PetscCheck(eps->ncv<=eps->nev+eps->mpd,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must not be larger than nev+mpd");
     59:   PetscFunctionReturn(PETSC_SUCCESS);
     60: }

     62: static PetscErrorCode EPSSetUp_Subspace(EPS eps)
     63: {
     64:   PetscBool isfilt;

     66:   PetscFunctionBegin;
     67:   EPSCheckDefinite(eps);
     68:   EPSCheckNotStructured(eps);
     69:   if (eps->max_it==PETSC_DETERMINE) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
     70:   if (!eps->which) PetscCall(EPSSetWhichEigenpairs_Default(eps));
     71:   if (eps->which==EPS_ALL) {
     72:     PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
     73:     PetscCheck(isfilt,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Spectrum slicing not supported in this solver");
     74:     PetscCall(EPSSetUp_Subspace_Filter(eps));
     75:   } else {
     76:     PetscCheck(eps->which==EPS_LARGEST_MAGNITUDE || eps->which==EPS_TARGET_MAGNITUDE,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver supports only largest magnitude or target magnitude eigenvalues");
     77:     PetscCall(EPSSetDimensions_Default(eps,&eps->nev,&eps->ncv,&eps->mpd));
     78:   }
     79:   EPSCheckUnsupported(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_EXTRACTION | EPS_FEATURE_TWOSIDED);
     80:   PetscCheck(eps->converged==EPSConvergedRelative,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver only supports relative convergence test");

     82:   PetscCall(EPSAllocateSolution(eps,0));
     83:   PetscCall(EPS_SetInnerProduct(eps));
     84:   if (eps->ishermitian) PetscCall(DSSetType(eps->ds,DSHEP));
     85:   else PetscCall(DSSetType(eps->ds,DSNHEP));
     86:   PetscCall(DSAllocate(eps->ds,eps->ncv));
     87:   PetscFunctionReturn(PETSC_SUCCESS);
     88: }

     90: static PetscErrorCode EPSSetUpSort_Subspace(EPS eps)
     91: {
     92:   SlepcSC sc;

     94:   PetscFunctionBegin;
     95:   PetscCall(EPSSetUpSort_Default(eps));
     96:   if (eps->which==EPS_ALL) {
     97:     PetscCall(DSGetSlepcSC(eps->ds,&sc));
     98:     sc->rg            = NULL;
     99:     sc->comparison    = SlepcCompareLargestReal;
    100:     sc->comparisonctx = NULL;
    101:     sc->map           = NULL;
    102:     sc->mapobj        = NULL;
    103:   }
    104:   PetscFunctionReturn(PETSC_SUCCESS);
    105: }

    107: /*
    108:    EPSSubspaceFindGroup - Find a group of nearly equimodular eigenvalues, provided
    109:    in arrays wr and wi, according to the tolerance grptol. Also the 2-norms
    110:    of the residuals must be passed in (rsd). Arrays are processed from index
    111:    l to index m only. The output information is:

    113:    ngrp - number of entries of the group
    114:    ctr  - (w(l)+w(l+ngrp-1))/2
    115:    ae   - average of wr(l),...,wr(l+ngrp-1)
    116:    arsd - average of rsd(l),...,rsd(l+ngrp-1)
    117: */
    118: static PetscErrorCode EPSSubspaceFindGroup(PetscInt l,PetscInt m,PetscScalar *wr,PetscScalar *wi,PetscReal *rsd,PetscReal grptol,PetscInt *ngrp,PetscReal *ctr,PetscReal *ae,PetscReal *arsd)
    119: {
    120:   PetscInt  i;
    121:   PetscReal rmod,rmod1;

    123:   PetscFunctionBegin;
    124:   *ngrp = 0;
    125:   *ctr = 0;
    126:   rmod = SlepcAbsEigenvalue(wr[l],wi[l]);

    128:   for (i=l;i<m;) {
    129:     rmod1 = SlepcAbsEigenvalue(wr[i],wi[i]);
    130:     if (PetscAbsReal(rmod-rmod1) > grptol*(rmod+rmod1)) break;
    131:     *ctr = (rmod+rmod1)/2.0;
    132:     if (wi[i] == 0.0) {
    133:       (*ngrp)++;
    134:       i++;
    135:     } else {
    136:       (*ngrp)+=2;
    137:       i+=2;
    138:     }
    139:   }

    141:   *ae = 0;
    142:   *arsd = 0;
    143:   if (*ngrp) {
    144:     for (i=l;i<l+*ngrp;i++) {
    145:       (*ae) += PetscRealPart(wr[i]);
    146:       (*arsd) += rsd[i]*rsd[i];
    147:     }
    148:     *ae = *ae / *ngrp;
    149:     *arsd = PetscSqrtReal(*arsd / *ngrp);
    150:   }
    151:   PetscFunctionReturn(PETSC_SUCCESS);
    152: }

    154: /*
    155:    EPSSubspaceResidualNorms - Computes the column norms of residual vectors
    156:    OP*V(1:n,l:m) - V*T(1:m,l:m), where, on entry, OP*V has been computed and
    157:    stored in R. On exit, rsd(l) to rsd(m) contain the computed norms.
    158: */
    159: static PetscErrorCode EPSSubspaceResidualNorms(BV R,BV V,Mat T,PetscInt l,PetscInt m,PetscScalar *eigi,PetscReal *rsd)
    160: {
    161:   PetscInt       i;

    163:   PetscFunctionBegin;
    164:   PetscCall(BVMult(R,-1.0,1.0,V,T));
    165:   for (i=l;i<m;i++) PetscCall(BVNormColumnBegin(R,i,NORM_2,rsd+i));
    166:   for (i=l;i<m;i++) PetscCall(BVNormColumnEnd(R,i,NORM_2,rsd+i));
    167: #if !defined(PETSC_USE_COMPLEX)
    168:   for (i=l;i<m-1;i++) {
    169:     if (eigi[i]!=0.0) {
    170:       rsd[i]   = SlepcAbs(rsd[i],rsd[i+1]);
    171:       rsd[i+1] = rsd[i];
    172:       i++;
    173:     }
    174:   }
    175: #endif
    176:   PetscFunctionReturn(PETSC_SUCCESS);
    177: }

    179: static PetscErrorCode EPSSolve_Subspace(EPS eps)
    180: {
    181:   Mat            H,Q,S,T,B;
    182:   BV             AV,R;
    183:   PetscBool      indef;
    184:   PetscInt       i,k,ngrp,nogrp,*itrsd,*itrsdold,*p_itrsd,*p_itrsdold;
    185:   PetscInt       nxtsrr,idsrr,idort,nxtort,nv,its;
    186:   PetscReal      arsd,oarsd,ctr,octr,ae,oae,*rsd,*orsd,*p_rsd,*p_orsd,tcond=1.0;
    187:   PetscScalar    *oeigr,*oeigi,*p_oeigr,*p_oeigi;
    188:   /* Parameters */
    189:   PetscInt       init = 5;        /* Number of initial iterations */
    190:   PetscReal      stpfac = 1.5;    /* Max num of iter before next SRR step */
    191:   PetscReal      alpha = 1.0;     /* Used to predict convergence of next residual */
    192:   PetscReal      beta = 1.1;      /* Used to predict convergence of next residual */
    193:   PetscReal      grptol = SLEPC_DEFAULT_TOL;   /* Tolerance for EPSSubspaceFindGroup */
    194:   PetscReal      cnvtol = 1e-6;   /* Convergence criterion for cnv */
    195:   PetscInt       orttol = 2;      /* Number of decimal digits whose loss
    196:                                      can be tolerated in orthogonalization */

    198:   PetscFunctionBegin;
    199:   its = 0;
    200:   PetscCall(PetscMalloc6(eps->ncv,&rsd,eps->ncv,&orsd,eps->ncv,&oeigr,eps->ncv,&oeigi,eps->ncv,&itrsd,eps->ncv,&itrsdold));
    201:   PetscCall(BVDuplicate(eps->V,&AV));
    202:   PetscCall(BVDuplicate(eps->V,&R));
    203:   PetscCall(STGetOperator(eps->st,&S));

    205:   for (i=0;i<eps->ncv;i++) {
    206:     rsd[i] = 0.0;
    207:     itrsd[i] = -1;
    208:   }

    210:   /* Complete the initial basis with random vectors and orthonormalize them */
    211:   for (k=eps->nini;k<eps->ncv;k++) {
    212:     PetscCall(BVSetRandomColumn(eps->V,k));
    213:     PetscCall(BVOrthonormalizeColumn(eps->V,k,PETSC_TRUE,NULL,NULL));
    214:   }

    216:   while (eps->reason == EPS_CONVERGED_ITERATING) {
    217:     eps->its++;
    218:     nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
    219:     PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,0));

    221:     for (i=eps->nconv;i<nv;i++) {
    222:       oeigr[i] = eps->eigr[i];
    223:       oeigi[i] = eps->eigi[i];
    224:       orsd[i]  = rsd[i];
    225:     }

    227:     /* AV(:,idx) = OP * V(:,idx) */
    228:     PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
    229:     PetscCall(BVSetActiveColumns(AV,eps->nconv,nv));
    230:     PetscCall(BVMatMult(eps->V,S,AV));

    232:     /* T(:,idx) = V' * AV(:,idx) */
    233:     PetscCall(BVSetActiveColumns(eps->V,0,nv));
    234:     PetscCall(DSGetMat(eps->ds,DS_MAT_A,&H));
    235:     PetscCall(BVDot(AV,eps->V,H));
    236:     PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&H));
    237:     PetscCall(DSSetState(eps->ds,DS_STATE_RAW));

    239:     /* Solve projected problem */
    240:     PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
    241:     PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
    242:     PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));

    244:     /* Update vectors V(:,idx) = V * U(:,idx) */
    245:     PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&Q));
    246:     PetscCall(BVSetActiveColumns(AV,0,nv));
    247:     PetscCall(BVSetActiveColumns(R,0,nv));
    248:     PetscCall(BVMultInPlace(eps->V,Q,eps->nconv,nv));
    249:     PetscCall(BVMultInPlace(AV,Q,eps->nconv,nv));
    250:     PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&Q));
    251:     PetscCall(BVCopy(AV,R));

    253:     /* Convergence check */
    254:     PetscCall(DSGetMat(eps->ds,DS_MAT_A,&T));
    255:     PetscCall(EPSSubspaceResidualNorms(R,eps->V,T,eps->nconv,nv,eps->eigi,rsd));
    256:     PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&T));

    258:     for (i=eps->nconv;i<nv;i++) {
    259:       itrsdold[i] = itrsd[i];
    260:       itrsd[i] = its;
    261:       eps->errest[i] = rsd[i];
    262:     }

    264:     for (;;) {
    265:       /* Find clusters of computed eigenvalues */
    266:       PetscCall(EPSSubspaceFindGroup(eps->nconv,nv,eps->eigr,eps->eigi,eps->errest,grptol,&ngrp,&ctr,&ae,&arsd));
    267:       PetscCall(EPSSubspaceFindGroup(eps->nconv,nv,oeigr,oeigi,orsd,grptol,&nogrp,&octr,&oae,&oarsd));
    268:       if (ngrp!=nogrp) break;
    269:       if (ngrp==0) break;
    270:       if (PetscAbsReal(ae-oae)>ctr*cnvtol*(itrsd[eps->nconv]-itrsdold[eps->nconv])) break;
    271:       if (arsd>ctr*eps->tol) break;
    272:       eps->nconv = eps->nconv + ngrp;
    273:       if (eps->nconv>=nv) break;
    274:     }

    276:     PetscCall(EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,nv));
    277:     EPSSetCtxThreshold(eps,eps->eigr,eps->eigi,eps->errest,eps->nconv,nv);
    278:     PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx));
    279:     if (eps->reason != EPS_CONVERGED_ITERATING) break;

    281:     /* Compute nxtsrr (iteration of next projection step) */
    282:     nxtsrr = PetscMax((PetscInt)PetscFloorReal(stpfac*its),init);

    284:     if (ngrp!=nogrp || ngrp==0 || arsd>=oarsd) {
    285:       idsrr = nxtsrr - its;
    286:     } else {
    287:       idsrr = (PetscInt)PetscFloorReal(alpha+beta*(itrsdold[eps->nconv]-itrsd[eps->nconv])*PetscLogReal(arsd/eps->tol)/PetscLogReal(arsd/oarsd));
    288:       idsrr = PetscMax(1,idsrr);
    289:     }
    290:     nxtsrr = PetscMin(nxtsrr,its+idsrr);

    292:     /* Compute nxtort (iteration of next orthogonalization step) */
    293:     PetscCall(DSCond(eps->ds,&tcond));
    294:     idort = PetscMax(1,(PetscInt)PetscFloorReal(orttol/PetscMax(1,PetscLog10Real(tcond))));
    295:     nxtort = PetscMin(its+idort,nxtsrr);
    296:     PetscCall(PetscInfo(eps,"Updated iteration counts: nxtort=%" PetscInt_FMT ", nxtsrr=%" PetscInt_FMT "\n",nxtort,nxtsrr));

    298:     /* V(:,idx) = AV(:,idx) */
    299:     PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
    300:     PetscCall(BVSetActiveColumns(AV,eps->nconv,nv));
    301:     PetscCall(BVCopy(AV,eps->V));
    302:     its++;

    304:     /* Orthogonalization loop */
    305:     do {
    306:       PetscCall(BVGetMatrix(eps->V,&B,&indef));
    307:       PetscCall(BVSetMatrix(eps->V,NULL,PETSC_FALSE));
    308:       while (its<nxtort) {
    309:         /* A(:,idx) = OP*V(:,idx) with normalization */
    310:         PetscCall(BVMatMult(eps->V,S,AV));
    311:         PetscCall(BVCopy(AV,eps->V));
    312:         PetscCall(BVNormalize(eps->V,NULL));
    313:         its++;
    314:       }
    315:       PetscCall(BVSetMatrix(eps->V,B,indef));
    316:       /* Orthonormalize vectors */
    317:       PetscCall(BVOrthogonalize(eps->V,NULL));
    318:       nxtort = PetscMin(its+idort,nxtsrr);
    319:     } while (its<nxtsrr);

    321:     if (eps->stop==EPS_STOP_THRESHOLD && nv/(PetscReal)eps->nconv<1.3) {  /* reallocate */
    322:       eps->ncv = 1.3*eps->nconv+1;
    323:       PetscCall(EPSReallocateSolution(eps,eps->ncv));
    324:       PetscCall(BVResize(AV,eps->ncv,PETSC_TRUE));
    325:       PetscCall(BVResize(R,eps->ncv,PETSC_TRUE));
    326:       for (i=nv;i<eps->ncv;i++) eps->perm[i] = i;
    327:       PetscCall(DSReallocate(eps->ds,eps->ncv));
    328:       for (k=nv;k<eps->ncv;k++) {
    329:         PetscCall(BVSetRandomColumn(eps->V,k));
    330:         PetscCall(BVOrthonormalizeColumn(eps->V,k,PETSC_TRUE,NULL,NULL));
    331:       }
    332:       p_rsd      = rsd;
    333:       p_orsd     = orsd;
    334:       p_oeigr    = oeigr;
    335:       p_oeigi    = oeigi;
    336:       p_itrsd    = itrsd;
    337:       p_itrsdold = itrsdold;
    338:       PetscCall(PetscMalloc6(eps->ncv,&rsd,eps->ncv,&orsd,eps->ncv,&oeigr,eps->ncv,&oeigi,eps->ncv,&itrsd,eps->ncv,&itrsdold));
    339:       for (i=0;i<nv;i++) {
    340:         rsd[i]      = p_rsd[i];
    341:         orsd[i]     = p_orsd[i];
    342:         oeigr[i]    = p_oeigr[i];
    343:         oeigi[i]    = p_oeigi[i];
    344:         itrsd[i]    = p_itrsd[i];
    345:         itrsdold[i] = p_itrsdold[i];
    346:       }
    347:       PetscCall(PetscFree6(p_rsd,p_orsd,p_oeigr,p_oeigi,p_itrsd,p_itrsdold));
    348:     }

    350:   }

    352:   PetscCall(PetscFree6(rsd,orsd,oeigr,oeigi,itrsd,itrsdold));
    353:   PetscCall(BVDestroy(&AV));
    354:   PetscCall(BVDestroy(&R));
    355:   PetscCall(STRestoreOperator(eps->st,&S));
    356:   PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
    357:   PetscFunctionReturn(PETSC_SUCCESS);
    358: }

    360: static PetscErrorCode EPSDestroy_Subspace(EPS eps)
    361: {
    362:   PetscFunctionBegin;
    363:   PetscCall(PetscFree(eps->data));
    364:   PetscFunctionReturn(PETSC_SUCCESS);
    365: }

    367: /*MC
    368:    EPSSUBSPACE - EPSSUBSPACE = "subspace" - The simple subspace iteration method.

    370:    Notes:
    371:    This solver is very basic and is not recommended in general, since it
    372:    will not be competitive with respect to other solvers.

    374:    The implemented method is subspace iteration with Rayleigh-Ritz projection and
    375:    locking, based on the SRRIT code {cite:p}`Bai97`.

    377:    Level: beginner

    379: .seealso: [](ch:eps), `EPS`, `EPSType`, `EPSSetType()`
    380: M*/
    381: SLEPC_EXTERN PetscErrorCode EPSCreate_Subspace(EPS eps)
    382: {
    383:   EPS_SUBSPACE *ctx;

    385:   PetscFunctionBegin;
    386:   PetscCall(PetscNew(&ctx));
    387:   eps->data  = (void*)ctx;

    389:   eps->useds = PETSC_TRUE;
    390:   eps->categ = EPS_CATEGORY_OTHER;

    392:   eps->ops->solve          = EPSSolve_Subspace;
    393:   eps->ops->setup          = EPSSetUp_Subspace;
    394:   eps->ops->setupsort      = EPSSetUpSort_Subspace;
    395:   eps->ops->destroy        = EPSDestroy_Subspace;
    396:   eps->ops->backtransform  = EPSBackTransform_Default;
    397:   eps->ops->computevectors = EPSComputeVectors_Schur;
    398:   PetscFunctionReturn(PETSC_SUCCESS);
    399: }