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