Line data    Source code

       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
       5             : 
       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"
      12             : 
      13             :    Method: Subspace Iteration
      14             : 
      15             :    Algorithm:
      16             : 
      17             :        Subspace iteration with Rayleigh-Ritz projection and locking,
      18             :        based on the SRRIT implementation.
      19             : 
      20             :    References:
      21             : 
      22             :        [1] "Subspace Iteration in SLEPc", SLEPc Technical Report STR-3,
      23             :            available at https://slepc.upv.es.
      24             : */
      25             : 
      26             : #include <slepc/private/epsimpl.h>
      27             : 
      28             : typedef struct {
      29             :   PetscBool estimatedrange;     /* the filter range was not set by the user */
      30             : } EPS_SUBSPACE;
      31             : 
      32           1 : static PetscErrorCode EPSSetUp_Subspace_Filter(EPS eps)
      33             : {
      34           1 :   EPS_SUBSPACE   *ctx = (EPS_SUBSPACE*)eps->data;
      35           1 :   PetscBool      estimaterange=PETSC_TRUE;
      36           1 :   PetscReal      rleft,rright;
      37           1 :   Mat            A;
      38             : 
      39           1 :   PetscFunctionBegin;
      40           1 :   EPSCheckHermitianCondition(eps,PETSC_TRUE," with polynomial filter");
      41           1 :   EPSCheckStandardCondition(eps,PETSC_TRUE," with polynomial filter");
      42           1 :   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           1 :   EPSCheckUnsupportedCondition(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION,PETSC_TRUE," with polynomial filter");
      44           1 :   PetscCall(STFilterSetInterval(eps->st,eps->inta,eps->intb));
      45           1 :   if (!ctx->estimatedrange) {
      46           1 :     PetscCall(STFilterGetRange(eps->st,&rleft,&rright));
      47           1 :     estimaterange = (!rleft && !rright)? PETSC_TRUE: PETSC_FALSE;
      48             :   }
      49           1 :   if (estimaterange) { /* user did not set a range */
      50           1 :     PetscCall(STGetMatrix(eps->st,0,&A));
      51           1 :     PetscCall(MatEstimateSpectralRange_EPS(A,&rleft,&rright));
      52           1 :     PetscCall(PetscInfo(eps,"Setting eigenvalue range to [%g,%g]\n",(double)rleft,(double)rright));
      53           1 :     PetscCall(STFilterSetRange(eps->st,rleft,rright));
      54           1 :     ctx->estimatedrange = PETSC_TRUE;
      55             :   }
      56           1 :   if (eps->ncv==PETSC_DETERMINE && eps->nev==1) eps->nev = 40;  /* user did not provide nev estimation */
      57           1 :   PetscCall(EPSSetDimensions_Default(eps,&eps->nev,&eps->ncv,&eps->mpd));
      58           1 :   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           1 :   PetscFunctionReturn(PETSC_SUCCESS);
      60             : }
      61             : 
      62          18 : static PetscErrorCode EPSSetUp_Subspace(EPS eps)
      63             : {
      64          18 :   PetscBool isfilt;
      65             : 
      66          18 :   PetscFunctionBegin;
      67          18 :   EPSCheckDefinite(eps);
      68          18 :   EPSCheckNotStructured(eps);
      69          18 :   if (eps->max_it==PETSC_DETERMINE) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
      70          18 :   if (!eps->which) PetscCall(EPSSetWhichEigenpairs_Default(eps));
      71          18 :   if (eps->which==EPS_ALL) {
      72           1 :     if (eps->nev==0) eps->nev = 1;
      73           1 :     PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
      74           1 :     PetscCheck(isfilt,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Spectrum slicing not supported in this solver");
      75           1 :     PetscCall(EPSSetUp_Subspace_Filter(eps));
      76             :   } else {
      77          17 :     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");
      78          17 :     PetscCall(EPSSetDimensions_Default(eps,&eps->nev,&eps->ncv,&eps->mpd));
      79             :   }
      80          18 :   EPSCheckUnsupported(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_EXTRACTION | EPS_FEATURE_THRESHOLD | EPS_FEATURE_TWOSIDED);
      81          18 :   PetscCheck(eps->converged==EPSConvergedRelative,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver only supports relative convergence test");
      82             : 
      83          18 :   PetscCall(EPSAllocateSolution(eps,0));
      84          18 :   PetscCall(EPS_SetInnerProduct(eps));
      85          18 :   if (eps->ishermitian) PetscCall(DSSetType(eps->ds,DSHEP));
      86           6 :   else PetscCall(DSSetType(eps->ds,DSNHEP));
      87          18 :   PetscCall(DSAllocate(eps->ds,eps->ncv));
      88          18 :   PetscFunctionReturn(PETSC_SUCCESS);
      89             : }
      90             : 
      91          18 : static PetscErrorCode EPSSetUpSort_Subspace(EPS eps)
      92             : {
      93          18 :   SlepcSC sc;
      94             : 
      95          18 :   PetscFunctionBegin;
      96          18 :   PetscCall(EPSSetUpSort_Default(eps));
      97          18 :   if (eps->which==EPS_ALL) {
      98           1 :     PetscCall(DSGetSlepcSC(eps->ds,&sc));
      99           1 :     sc->rg            = NULL;
     100           1 :     sc->comparison    = SlepcCompareLargestReal;
     101           1 :     sc->comparisonctx = NULL;
     102           1 :     sc->map           = NULL;
     103           1 :     sc->mapobj        = NULL;
     104             :   }
     105          18 :   PetscFunctionReturn(PETSC_SUCCESS);
     106             : }
     107             : 
     108             : /*
     109             :    EPSSubspaceFindGroup - Find a group of nearly equimodular eigenvalues, provided
     110             :    in arrays wr and wi, according to the tolerance grptol. Also the 2-norms
     111             :    of the residuals must be passed in (rsd). Arrays are processed from index
     112             :    l to index m only. The output information is:
     113             : 
     114             :    ngrp - number of entries of the group
     115             :    ctr  - (w(l)+w(l+ngrp-1))/2
     116             :    ae   - average of wr(l),...,wr(l+ngrp-1)
     117             :    arsd - average of rsd(l),...,rsd(l+ngrp-1)
     118             : */
     119         544 : static PetscErrorCode EPSSubspaceFindGroup(PetscInt l,PetscInt m,PetscScalar *wr,PetscScalar *wi,PetscReal *rsd,PetscReal grptol,PetscInt *ngrp,PetscReal *ctr,PetscReal *ae,PetscReal *arsd)
     120             : {
     121         544 :   PetscInt  i;
     122         544 :   PetscReal rmod,rmod1;
     123             : 
     124         544 :   PetscFunctionBegin;
     125         544 :   *ngrp = 0;
     126         544 :   *ctr = 0;
     127         544 :   rmod = SlepcAbsEigenvalue(wr[l],wi[l]);
     128             : 
     129        1921 :   for (i=l;i<m;) {
     130        1357 :     rmod1 = SlepcAbsEigenvalue(wr[i],wi[i]);
     131        1357 :     if (PetscAbsReal(rmod-rmod1) > grptol*(rmod+rmod1)) break;
     132         833 :     *ctr = (rmod+rmod1)/2.0;
     133         833 :     if (wi[i] == 0.0) {
     134         833 :       (*ngrp)++;
     135         833 :       i++;
     136             :     } else {
     137           0 :       (*ngrp)+=2;
     138           0 :       i+=2;
     139             :     }
     140             :   }
     141             : 
     142         544 :   *ae = 0;
     143         544 :   *arsd = 0;
     144         544 :   if (*ngrp) {
     145        1377 :     for (i=l;i<l+*ngrp;i++) {
     146         833 :       (*ae) += PetscRealPart(wr[i]);
     147         833 :       (*arsd) += rsd[i]*rsd[i];
     148             :     }
     149         544 :     *ae = *ae / *ngrp;
     150         544 :     *arsd = PetscSqrtReal(*arsd / *ngrp);
     151             :   }
     152         544 :   PetscFunctionReturn(PETSC_SUCCESS);
     153             : }
     154             : 
     155             : /*
     156             :    EPSSubspaceResidualNorms - Computes the column norms of residual vectors
     157             :    OP*V(1:n,l:m) - V*T(1:m,l:m), where, on entry, OP*V has been computed and
     158             :    stored in R. On exit, rsd(l) to rsd(m) contain the computed norms.
     159             : */
     160         158 : static PetscErrorCode EPSSubspaceResidualNorms(BV R,BV V,Mat T,PetscInt l,PetscInt m,PetscScalar *eigi,PetscReal *rsd)
     161             : {
     162         158 :   PetscInt       i;
     163             : 
     164         158 :   PetscFunctionBegin;
     165         158 :   PetscCall(BVMult(R,-1.0,1.0,V,T));
     166        2650 :   for (i=l;i<m;i++) PetscCall(BVNormColumnBegin(R,i,NORM_2,rsd+i));
     167        2650 :   for (i=l;i<m;i++) PetscCall(BVNormColumnEnd(R,i,NORM_2,rsd+i));
     168             : #if !defined(PETSC_USE_COMPLEX)
     169        2492 :   for (i=l;i<m-1;i++) {
     170        2334 :     if (eigi[i]!=0.0) {
     171           0 :       rsd[i]   = SlepcAbs(rsd[i],rsd[i+1]);
     172           0 :       rsd[i+1] = rsd[i];
     173           0 :       i++;
     174             :     }
     175             :   }
     176             : #endif
     177         158 :   PetscFunctionReturn(PETSC_SUCCESS);
     178             : }
     179             : 
     180          18 : static PetscErrorCode EPSSolve_Subspace(EPS eps)
     181             : {
     182          18 :   Mat            H,Q,S,T,B;
     183          18 :   BV             AV,R;
     184          18 :   PetscBool      indef;
     185          18 :   PetscInt       i,k,ld,ngrp,nogrp,*itrsd,*itrsdold;
     186          18 :   PetscInt       nxtsrr,idsrr,idort,nxtort,nv,ncv = eps->ncv,its,ninside;
     187          18 :   PetscReal      arsd,oarsd,ctr,octr,ae,oae,*rsd,*orsd,tcond=1.0,gamma;
     188          18 :   PetscScalar    *oeigr,*oeigi;
     189             :   /* Parameters */
     190          18 :   PetscInt       init = 5;        /* Number of initial iterations */
     191          18 :   PetscReal      stpfac = 1.5;    /* Max num of iter before next SRR step */
     192          18 :   PetscReal      alpha = 1.0;     /* Used to predict convergence of next residual */
     193          18 :   PetscReal      beta = 1.1;      /* Used to predict convergence of next residual */
     194          18 :   PetscReal      grptol = SLEPC_DEFAULT_TOL;   /* Tolerance for EPSSubspaceFindGroup */
     195          18 :   PetscReal      cnvtol = 1e-6;   /* Convergence criterion for cnv */
     196          18 :   PetscInt       orttol = 2;      /* Number of decimal digits whose loss
     197             :                                      can be tolerated in orthogonalization */
     198             : 
     199          18 :   PetscFunctionBegin;
     200          18 :   its = 0;
     201          18 :   PetscCall(PetscMalloc6(ncv,&rsd,ncv,&orsd,ncv,&oeigr,ncv,&oeigi,ncv,&itrsd,ncv,&itrsdold));
     202          18 :   PetscCall(DSGetLeadingDimension(eps->ds,&ld));
     203          18 :   PetscCall(BVDuplicate(eps->V,&AV));
     204          18 :   PetscCall(BVDuplicate(eps->V,&R));
     205          18 :   PetscCall(STGetOperator(eps->st,&S));
     206             : 
     207         331 :   for (i=0;i<ncv;i++) {
     208         313 :     rsd[i] = 0.0;
     209         313 :     itrsd[i] = -1;
     210             :   }
     211             : 
     212             :   /* Complete the initial basis with random vectors and orthonormalize them */
     213         329 :   for (k=eps->nini;k<ncv;k++) {
     214         311 :     PetscCall(BVSetRandomColumn(eps->V,k));
     215         311 :     PetscCall(BVOrthonormalizeColumn(eps->V,k,PETSC_TRUE,NULL,NULL));
     216             :   }
     217             : 
     218         158 :   while (eps->reason == EPS_CONVERGED_ITERATING) {
     219         158 :     eps->its++;
     220         158 :     nv = PetscMin(eps->nconv+eps->mpd,ncv);
     221         158 :     PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,0));
     222             : 
     223        2650 :     for (i=eps->nconv;i<nv;i++) {
     224        2492 :       oeigr[i] = eps->eigr[i];
     225        2492 :       oeigi[i] = eps->eigi[i];
     226        2492 :       orsd[i]  = rsd[i];
     227             :     }
     228             : 
     229             :     /* AV(:,idx) = OP * V(:,idx) */
     230         158 :     PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
     231         158 :     PetscCall(BVSetActiveColumns(AV,eps->nconv,nv));
     232         158 :     PetscCall(BVMatMult(eps->V,S,AV));
     233             : 
     234             :     /* T(:,idx) = V' * AV(:,idx) */
     235         158 :     PetscCall(BVSetActiveColumns(eps->V,0,nv));
     236         158 :     PetscCall(DSGetMat(eps->ds,DS_MAT_A,&H));
     237         158 :     PetscCall(BVDot(AV,eps->V,H));
     238         158 :     PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&H));
     239         158 :     PetscCall(DSSetState(eps->ds,DS_STATE_RAW));
     240             : 
     241             :     /* Solve projected problem */
     242         158 :     PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
     243         158 :     PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
     244         158 :     PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
     245             : 
     246             :     /* Update vectors V(:,idx) = V * U(:,idx) */
     247         158 :     PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&Q));
     248         158 :     PetscCall(BVSetActiveColumns(AV,0,nv));
     249         158 :     PetscCall(BVSetActiveColumns(R,0,nv));
     250         158 :     PetscCall(BVMultInPlace(eps->V,Q,eps->nconv,nv));
     251         158 :     PetscCall(BVMultInPlace(AV,Q,eps->nconv,nv));
     252         158 :     PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&Q));
     253         158 :     PetscCall(BVCopy(AV,R));
     254             : 
     255             :     /* Convergence check */
     256         158 :     PetscCall(DSGetMat(eps->ds,DS_MAT_A,&T));
     257         158 :     PetscCall(EPSSubspaceResidualNorms(R,eps->V,T,eps->nconv,nv,eps->eigi,rsd));
     258         158 :     PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&T));
     259             : 
     260         158 :     if (eps->which==EPS_ALL && eps->its>1) {   /* adjust eigenvalue count */
     261          10 :       ninside = 0;
     262          10 :       PetscCall(STFilterGetThreshold(eps->st,&gamma));
     263          91 :       for (i=eps->nconv;i<nv;i++) {
     264          91 :         if (PetscRealPart(eps->eigr[i]) < gamma) break;
     265          81 :         ninside++;
     266             :       }
     267          10 :       eps->nev = eps->nconv+ninside;
     268             :     }
     269        2650 :     for (i=eps->nconv;i<nv;i++) {
     270        2492 :       itrsdold[i] = itrsd[i];
     271        2492 :       itrsd[i] = its;
     272        2492 :       eps->errest[i] = rsd[i];
     273             :     }
     274             : 
     275         272 :     for (;;) {
     276             :       /* Find clusters of computed eigenvalues */
     277         272 :       PetscCall(EPSSubspaceFindGroup(eps->nconv,nv,eps->eigr,eps->eigi,eps->errest,grptol,&ngrp,&ctr,&ae,&arsd));
     278         272 :       PetscCall(EPSSubspaceFindGroup(eps->nconv,nv,oeigr,oeigi,orsd,grptol,&nogrp,&octr,&oae,&oarsd));
     279         272 :       if (ngrp!=nogrp) break;
     280         254 :       if (ngrp==0) break;
     281         254 :       if (PetscAbsReal(ae-oae)>ctr*cnvtol*(itrsd[eps->nconv]-itrsdold[eps->nconv])) break;
     282         184 :       if (arsd>ctr*eps->tol) break;
     283         115 :       eps->nconv = eps->nconv + ngrp;
     284         115 :       if (eps->nconv>=nv) break;
     285             :     }
     286             : 
     287         158 :     PetscCall(EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,nv));
     288         158 :     PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx));
     289         158 :     if (eps->reason != EPS_CONVERGED_ITERATING) break;
     290             : 
     291             :     /* Compute nxtsrr (iteration of next projection step) */
     292         140 :     nxtsrr = PetscMin(eps->max_it,PetscMax((PetscInt)PetscFloorReal(stpfac*its),init));
     293             : 
     294         140 :     if (ngrp!=nogrp || ngrp==0 || arsd>=oarsd) {
     295          21 :       idsrr = nxtsrr - its;
     296             :     } else {
     297         119 :       idsrr = (PetscInt)PetscFloorReal(alpha+beta*(itrsdold[eps->nconv]-itrsd[eps->nconv])*PetscLogReal(arsd/eps->tol)/PetscLogReal(arsd/oarsd));
     298         119 :       idsrr = PetscMax(1,idsrr);
     299             :     }
     300         140 :     nxtsrr = PetscMin(nxtsrr,its+idsrr);
     301             : 
     302             :     /* Compute nxtort (iteration of next orthogonalization step) */
     303         140 :     PetscCall(DSCond(eps->ds,&tcond));
     304         140 :     idort = PetscMax(1,(PetscInt)PetscFloorReal(orttol/PetscMax(1,PetscLog10Real(tcond))));
     305         140 :     nxtort = PetscMin(its+idort,nxtsrr);
     306         140 :     PetscCall(PetscInfo(eps,"Updated iteration counts: nxtort=%" PetscInt_FMT ", nxtsrr=%" PetscInt_FMT "\n",nxtort,nxtsrr));
     307             : 
     308             :     /* V(:,idx) = AV(:,idx) */
     309         140 :     PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
     310         140 :     PetscCall(BVSetActiveColumns(AV,eps->nconv,nv));
     311         140 :     PetscCall(BVCopy(AV,eps->V));
     312         140 :     its++;
     313             : 
     314             :     /* Orthogonalization loop */
     315        1742 :     do {
     316        1742 :       PetscCall(BVGetMatrix(eps->V,&B,&indef));
     317        1742 :       PetscCall(BVSetMatrix(eps->V,NULL,PETSC_FALSE));
     318        4905 :       while (its<nxtort) {
     319             :         /* A(:,idx) = OP*V(:,idx) with normalization */
     320        3163 :         PetscCall(BVMatMult(eps->V,S,AV));
     321        3163 :         PetscCall(BVCopy(AV,eps->V));
     322        3163 :         PetscCall(BVNormalize(eps->V,NULL));
     323        3163 :         its++;
     324             :       }
     325        1742 :       PetscCall(BVSetMatrix(eps->V,B,indef));
     326             :       /* Orthonormalize vectors */
     327        1742 :       PetscCall(BVOrthogonalize(eps->V,NULL));
     328        1742 :       nxtort = PetscMin(its+idort,nxtsrr);
     329        1742 :     } while (its<nxtsrr);
     330             :   }
     331             : 
     332          18 :   PetscCall(PetscFree6(rsd,orsd,oeigr,oeigi,itrsd,itrsdold));
     333          18 :   PetscCall(BVDestroy(&AV));
     334          18 :   PetscCall(BVDestroy(&R));
     335          18 :   PetscCall(STRestoreOperator(eps->st,&S));
     336          18 :   PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
     337          18 :   PetscFunctionReturn(PETSC_SUCCESS);
     338             : }
     339             : 
     340          16 : static PetscErrorCode EPSDestroy_Subspace(EPS eps)
     341             : {
     342          16 :   PetscFunctionBegin;
     343          16 :   PetscCall(PetscFree(eps->data));
     344          16 :   PetscFunctionReturn(PETSC_SUCCESS);
     345             : }
     346             : 
     347          16 : SLEPC_EXTERN PetscErrorCode EPSCreate_Subspace(EPS eps)
     348             : {
     349          16 :   EPS_SUBSPACE *ctx;
     350             : 
     351          16 :   PetscFunctionBegin;
     352          16 :   PetscCall(PetscNew(&ctx));
     353          16 :   eps->data  = (void*)ctx;
     354             : 
     355          16 :   eps->useds = PETSC_TRUE;
     356          16 :   eps->categ = EPS_CATEGORY_OTHER;
     357             : 
     358          16 :   eps->ops->solve          = EPSSolve_Subspace;
     359          16 :   eps->ops->setup          = EPSSetUp_Subspace;
     360          16 :   eps->ops->setupsort      = EPSSetUpSort_Subspace;
     361          16 :   eps->ops->destroy        = EPSDestroy_Subspace;
     362          16 :   eps->ops->backtransform  = EPSBackTransform_Default;
     363          16 :   eps->ops->computevectors = EPSComputeVectors_Schur;
     364          16 :   PetscFunctionReturn(PETSC_SUCCESS);
     365             : }