LCOV - code coverage report
Current view: top level - eps/impls/subspace - subspace.c (source / functions) Hit Total Coverage
Test: SLEPc Lines: 218 223 97.8 %
Date: 2024-03-29 00:55:19 Functions: 8 8 100.0 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.14