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 :     PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
      73           1 :     PetscCheck(isfilt,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Spectrum slicing not supported in this solver");
      74           1 :     PetscCall(EPSSetUp_Subspace_Filter(eps));
      75             :   } else {
      76          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");
      77          17 :     PetscCall(EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd));
      78             :   }
      79          18 :   EPSCheckUnsupported(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_EXTRACTION | EPS_FEATURE_TWOSIDED);
      80          18 :   PetscCheck(eps->converged==EPSConvergedRelative,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver only supports relative convergence test");
      81             : 
      82          18 :   PetscCall(EPSAllocateSolution(eps,0));
      83          18 :   PetscCall(EPS_SetInnerProduct(eps));
      84          18 :   if (eps->ishermitian) PetscCall(DSSetType(eps->ds,DSHEP));
      85           6 :   else PetscCall(DSSetType(eps->ds,DSNHEP));
      86          18 :   PetscCall(DSAllocate(eps->ds,eps->ncv));
      87          18 :   PetscFunctionReturn(PETSC_SUCCESS);
      88             : }
      89             : 
      90          18 : static PetscErrorCode EPSSetUpSort_Subspace(EPS eps)
      91             : {
      92          18 :   SlepcSC sc;
      93             : 
      94          18 :   PetscFunctionBegin;
      95          18 :   PetscCall(EPSSetUpSort_Default(eps));
      96          18 :   if (eps->which==EPS_ALL) {
      97           1 :     PetscCall(DSGetSlepcSC(eps->ds,&sc));
      98           1 :     sc->rg            = NULL;
      99           1 :     sc->comparison    = SlepcCompareLargestReal;
     100           1 :     sc->comparisonctx = NULL;
     101           1 :     sc->map           = NULL;
     102           1 :     sc->mapobj        = NULL;
     103             :   }
     104          18 :   PetscFunctionReturn(PETSC_SUCCESS);
     105             : }
     106             : 
     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:
     112             : 
     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         544 : 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         544 :   PetscInt  i;
     121         544 :   PetscReal rmod,rmod1;
     122             : 
     123         544 :   PetscFunctionBegin;
     124         544 :   *ngrp = 0;
     125         544 :   *ctr = 0;
     126         544 :   rmod = SlepcAbsEigenvalue(wr[l],wi[l]);
     127             : 
     128        1921 :   for (i=l;i<m;) {
     129        1357 :     rmod1 = SlepcAbsEigenvalue(wr[i],wi[i]);
     130        1357 :     if (PetscAbsReal(rmod-rmod1) > grptol*(rmod+rmod1)) break;
     131         833 :     *ctr = (rmod+rmod1)/2.0;
     132         833 :     if (wi[i] == 0.0) {
     133         833 :       (*ngrp)++;
     134         833 :       i++;
     135             :     } else {
     136           0 :       (*ngrp)+=2;
     137           0 :       i+=2;
     138             :     }
     139             :   }
     140             : 
     141         544 :   *ae = 0;
     142         544 :   *arsd = 0;
     143         544 :   if (*ngrp) {
     144        1377 :     for (i=l;i<l+*ngrp;i++) {
     145         833 :       (*ae) += PetscRealPart(wr[i]);
     146         833 :       (*arsd) += rsd[i]*rsd[i];
     147             :     }
     148         544 :     *ae = *ae / *ngrp;
     149         544 :     *arsd = PetscSqrtReal(*arsd / *ngrp);
     150             :   }
     151         544 :   PetscFunctionReturn(PETSC_SUCCESS);
     152             : }
     153             : 
     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         158 : static PetscErrorCode EPSSubspaceResidualNorms(BV R,BV V,Mat T,PetscInt l,PetscInt m,PetscScalar *eigi,PetscReal *rsd)
     160             : {
     161         158 :   PetscInt       i;
     162             : 
     163         158 :   PetscFunctionBegin;
     164         158 :   PetscCall(BVMult(R,-1.0,1.0,V,T));
     165        2650 :   for (i=l;i<m;i++) PetscCall(BVNormColumnBegin(R,i,NORM_2,rsd+i));
     166        2650 :   for (i=l;i<m;i++) PetscCall(BVNormColumnEnd(R,i,NORM_2,rsd+i));
     167             : #if !defined(PETSC_USE_COMPLEX)
     168        2492 :   for (i=l;i<m-1;i++) {
     169        2334 :     if (eigi[i]!=0.0) {
     170           0 :       rsd[i]   = SlepcAbs(rsd[i],rsd[i+1]);
     171           0 :       rsd[i+1] = rsd[i];
     172           0 :       i++;
     173             :     }
     174             :   }
     175             : #endif
     176         158 :   PetscFunctionReturn(PETSC_SUCCESS);
     177             : }
     178             : 
     179          18 : static PetscErrorCode EPSSolve_Subspace(EPS eps)
     180             : {
     181          18 :   Mat            H,Q,S,T,B;
     182          18 :   BV             AV,R;
     183          18 :   PetscBool      indef;
     184          18 :   PetscInt       i,k,ld,ngrp,nogrp,*itrsd,*itrsdold;
     185          18 :   PetscInt       nxtsrr,idsrr,idort,nxtort,nv,ncv = eps->ncv,its,ninside;
     186          18 :   PetscReal      arsd,oarsd,ctr,octr,ae,oae,*rsd,*orsd,tcond=1.0,gamma;
     187          18 :   PetscScalar    *oeigr,*oeigi;
     188             :   /* Parameters */
     189          18 :   PetscInt       init = 5;        /* Number of initial iterations */
     190          18 :   PetscReal      stpfac = 1.5;    /* Max num of iter before next SRR step */
     191          18 :   PetscReal      alpha = 1.0;     /* Used to predict convergence of next residual */
     192          18 :   PetscReal      beta = 1.1;      /* Used to predict convergence of next residual */
     193          18 :   PetscReal      grptol = SLEPC_DEFAULT_TOL;   /* Tolerance for EPSSubspaceFindGroup */
     194          18 :   PetscReal      cnvtol = 1e-6;   /* Convergence criterion for cnv */
     195          18 :   PetscInt       orttol = 2;      /* Number of decimal digits whose loss
     196             :                                      can be tolerated in orthogonalization */
     197             : 
     198          18 :   PetscFunctionBegin;
     199          18 :   its = 0;
     200          18 :   PetscCall(PetscMalloc6(ncv,&rsd,ncv,&orsd,ncv,&oeigr,ncv,&oeigi,ncv,&itrsd,ncv,&itrsdold));
     201          18 :   PetscCall(DSGetLeadingDimension(eps->ds,&ld));
     202          18 :   PetscCall(BVDuplicate(eps->V,&AV));
     203          18 :   PetscCall(BVDuplicate(eps->V,&R));
     204          18 :   PetscCall(STGetOperator(eps->st,&S));
     205             : 
     206         331 :   for (i=0;i<ncv;i++) {
     207         313 :     rsd[i] = 0.0;
     208         313 :     itrsd[i] = -1;
     209             :   }
     210             : 
     211             :   /* Complete the initial basis with random vectors and orthonormalize them */
     212         329 :   for (k=eps->nini;k<ncv;k++) {
     213         311 :     PetscCall(BVSetRandomColumn(eps->V,k));
     214         311 :     PetscCall(BVOrthonormalizeColumn(eps->V,k,PETSC_TRUE,NULL,NULL));
     215             :   }
     216             : 
     217         158 :   while (eps->reason == EPS_CONVERGED_ITERATING) {
     218         158 :     eps->its++;
     219         158 :     nv = PetscMin(eps->nconv+eps->mpd,ncv);
     220         158 :     PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,0));
     221             : 
     222        2650 :     for (i=eps->nconv;i<nv;i++) {
     223        2492 :       oeigr[i] = eps->eigr[i];
     224        2492 :       oeigi[i] = eps->eigi[i];
     225        2492 :       orsd[i]  = rsd[i];
     226             :     }
     227             : 
     228             :     /* AV(:,idx) = OP * V(:,idx) */
     229         158 :     PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
     230         158 :     PetscCall(BVSetActiveColumns(AV,eps->nconv,nv));
     231         158 :     PetscCall(BVMatMult(eps->V,S,AV));
     232             : 
     233             :     /* T(:,idx) = V' * AV(:,idx) */
     234         158 :     PetscCall(BVSetActiveColumns(eps->V,0,nv));
     235         158 :     PetscCall(DSGetMat(eps->ds,DS_MAT_A,&H));
     236         158 :     PetscCall(BVDot(AV,eps->V,H));
     237         158 :     PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&H));
     238         158 :     PetscCall(DSSetState(eps->ds,DS_STATE_RAW));
     239             : 
     240             :     /* Solve projected problem */
     241         158 :     PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
     242         158 :     PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
     243         158 :     PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
     244             : 
     245             :     /* Update vectors V(:,idx) = V * U(:,idx) */
     246         158 :     PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&Q));
     247         158 :     PetscCall(BVSetActiveColumns(AV,0,nv));
     248         158 :     PetscCall(BVSetActiveColumns(R,0,nv));
     249         158 :     PetscCall(BVMultInPlace(eps->V,Q,eps->nconv,nv));
     250         158 :     PetscCall(BVMultInPlace(AV,Q,eps->nconv,nv));
     251         158 :     PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&Q));
     252         158 :     PetscCall(BVCopy(AV,R));
     253             : 
     254             :     /* Convergence check */
     255         158 :     PetscCall(DSGetMat(eps->ds,DS_MAT_A,&T));
     256         158 :     PetscCall(EPSSubspaceResidualNorms(R,eps->V,T,eps->nconv,nv,eps->eigi,rsd));
     257         158 :     PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&T));
     258             : 
     259         158 :     if (eps->which==EPS_ALL && eps->its>1) {   /* adjust eigenvalue count */
     260          10 :       ninside = 0;
     261          10 :       PetscCall(STFilterGetThreshold(eps->st,&gamma));
     262          91 :       for (i=eps->nconv;i<nv;i++) {
     263          91 :         if (PetscRealPart(eps->eigr[i]) < gamma) break;
     264          81 :         ninside++;
     265             :       }
     266          10 :       eps->nev = eps->nconv+ninside;
     267             :     }
     268        2650 :     for (i=eps->nconv;i<nv;i++) {
     269        2492 :       itrsdold[i] = itrsd[i];
     270        2492 :       itrsd[i] = its;
     271        2492 :       eps->errest[i] = rsd[i];
     272             :     }
     273             : 
     274         272 :     for (;;) {
     275             :       /* Find clusters of computed eigenvalues */
     276         272 :       PetscCall(EPSSubspaceFindGroup(eps->nconv,nv,eps->eigr,eps->eigi,eps->errest,grptol,&ngrp,&ctr,&ae,&arsd));
     277         272 :       PetscCall(EPSSubspaceFindGroup(eps->nconv,nv,oeigr,oeigi,orsd,grptol,&nogrp,&octr,&oae,&oarsd));
     278         272 :       if (ngrp!=nogrp) break;
     279         254 :       if (ngrp==0) break;
     280         254 :       if (PetscAbsReal(ae-oae)>ctr*cnvtol*(itrsd[eps->nconv]-itrsdold[eps->nconv])) break;
     281         184 :       if (arsd>ctr*eps->tol) break;
     282         115 :       eps->nconv = eps->nconv + ngrp;
     283         115 :       if (eps->nconv>=nv) break;
     284             :     }
     285             : 
     286         158 :     PetscCall(EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,nv));
     287         158 :     PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx));
     288         158 :     if (eps->reason != EPS_CONVERGED_ITERATING) break;
     289             : 
     290             :     /* Compute nxtsrr (iteration of next projection step) */
     291         140 :     nxtsrr = PetscMin(eps->max_it,PetscMax((PetscInt)PetscFloorReal(stpfac*its),init));
     292             : 
     293         140 :     if (ngrp!=nogrp || ngrp==0 || arsd>=oarsd) {
     294          21 :       idsrr = nxtsrr - its;
     295             :     } else {
     296         119 :       idsrr = (PetscInt)PetscFloorReal(alpha+beta*(itrsdold[eps->nconv]-itrsd[eps->nconv])*PetscLogReal(arsd/eps->tol)/PetscLogReal(arsd/oarsd));
     297         119 :       idsrr = PetscMax(1,idsrr);
     298             :     }
     299         140 :     nxtsrr = PetscMin(nxtsrr,its+idsrr);
     300             : 
     301             :     /* Compute nxtort (iteration of next orthogonalization step) */
     302         140 :     PetscCall(DSCond(eps->ds,&tcond));
     303         140 :     idort = PetscMax(1,(PetscInt)PetscFloorReal(orttol/PetscMax(1,PetscLog10Real(tcond))));
     304         140 :     nxtort = PetscMin(its+idort,nxtsrr);
     305         140 :     PetscCall(PetscInfo(eps,"Updated iteration counts: nxtort=%" PetscInt_FMT ", nxtsrr=%" PetscInt_FMT "\n",nxtort,nxtsrr));
     306             : 
     307             :     /* V(:,idx) = AV(:,idx) */
     308         140 :     PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
     309         140 :     PetscCall(BVSetActiveColumns(AV,eps->nconv,nv));
     310         140 :     PetscCall(BVCopy(AV,eps->V));
     311         140 :     its++;
     312             : 
     313             :     /* Orthogonalization loop */
     314        1742 :     do {
     315        1742 :       PetscCall(BVGetMatrix(eps->V,&B,&indef));
     316        1742 :       PetscCall(BVSetMatrix(eps->V,NULL,PETSC_FALSE));
     317        4905 :       while (its<nxtort) {
     318             :         /* A(:,idx) = OP*V(:,idx) with normalization */
     319        3163 :         PetscCall(BVMatMult(eps->V,S,AV));
     320        3163 :         PetscCall(BVCopy(AV,eps->V));
     321        3163 :         PetscCall(BVNormalize(eps->V,NULL));
     322        3163 :         its++;
     323             :       }
     324        1742 :       PetscCall(BVSetMatrix(eps->V,B,indef));
     325             :       /* Orthonormalize vectors */
     326        1742 :       PetscCall(BVOrthogonalize(eps->V,NULL));
     327        1742 :       nxtort = PetscMin(its+idort,nxtsrr);
     328        1742 :     } while (its<nxtsrr);
     329             :   }
     330             : 
     331          18 :   PetscCall(PetscFree6(rsd,orsd,oeigr,oeigi,itrsd,itrsdold));
     332          18 :   PetscCall(BVDestroy(&AV));
     333          18 :   PetscCall(BVDestroy(&R));
     334          18 :   PetscCall(STRestoreOperator(eps->st,&S));
     335          18 :   PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
     336          18 :   PetscFunctionReturn(PETSC_SUCCESS);
     337             : }
     338             : 
     339          16 : static PetscErrorCode EPSDestroy_Subspace(EPS eps)
     340             : {
     341          16 :   PetscFunctionBegin;
     342          16 :   PetscCall(PetscFree(eps->data));
     343          16 :   PetscFunctionReturn(PETSC_SUCCESS);
     344             : }
     345             : 
     346          16 : SLEPC_EXTERN PetscErrorCode EPSCreate_Subspace(EPS eps)
     347             : {
     348          16 :   EPS_SUBSPACE *ctx;
     349             : 
     350          16 :   PetscFunctionBegin;
     351          16 :   PetscCall(PetscNew(&ctx));
     352          16 :   eps->data  = (void*)ctx;
     353             : 
     354          16 :   eps->useds = PETSC_TRUE;
     355          16 :   eps->categ = EPS_CATEGORY_OTHER;
     356             : 
     357          16 :   eps->ops->solve          = EPSSolve_Subspace;
     358          16 :   eps->ops->setup          = EPSSetUp_Subspace;
     359          16 :   eps->ops->setupsort      = EPSSetUpSort_Subspace;
     360          16 :   eps->ops->destroy        = EPSDestroy_Subspace;
     361          16 :   eps->ops->backtransform  = EPSBackTransform_Default;
     362          16 :   eps->ops->computevectors = EPSComputeVectors_Schur;
     363          16 :   PetscFunctionReturn(PETSC_SUCCESS);
     364             : }