LCOV - code coverage report
Current view: top level - eps/impls/subspace - subspace.c (source / functions) Hit Total Coverage
Test: SLEPc Lines: 160 165 97.0 %
Date: 2021-08-02 00:32:28 Functions: 5 5 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-2021, 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          17 : PetscErrorCode EPSSetUp_Subspace(EPS eps)
      29             : {
      30          17 :   PetscErrorCode ierr;
      31             : 
      32          17 :   PetscFunctionBegin;
      33          17 :   EPSCheckDefinite(eps);
      34          17 :   ierr = EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd);CHKERRQ(ierr);
      35          17 :   if (eps->max_it==PETSC_DEFAULT) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
      36          17 :   if (!eps->which) { ierr = EPSSetWhichEigenpairs_Default(eps);CHKERRQ(ierr); }
      37          17 :   if (eps->which!=EPS_LARGEST_MAGNITUDE && eps->which!=EPS_TARGET_MAGNITUDE) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver supports only largest magnitude or target magnitude eigenvalues");
      38          17 :   EPSCheckUnsupported(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_EXTRACTION | EPS_FEATURE_TWOSIDED);
      39          17 :   if (eps->converged != EPSConvergedRelative) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver only supports relative convergence test");
      40             : 
      41          17 :   ierr = EPSAllocateSolution(eps,0);CHKERRQ(ierr);
      42          17 :   ierr = EPS_SetInnerProduct(eps);CHKERRQ(ierr);
      43          17 :   if (eps->ishermitian) {
      44          11 :     ierr = DSSetType(eps->ds,DSHEP);CHKERRQ(ierr);
      45             :   } else {
      46           6 :     ierr = DSSetType(eps->ds,DSNHEP);CHKERRQ(ierr);
      47             :   }
      48          17 :   ierr = DSAllocate(eps->ds,eps->ncv);CHKERRQ(ierr);
      49          17 :   PetscFunctionReturn(0);
      50             : }
      51             : 
      52             : /*
      53             :    EPSSubspaceFindGroup - Find a group of nearly equimodular eigenvalues, provided
      54             :    in arrays wr and wi, according to the tolerance grptol. Also the 2-norms
      55             :    of the residuals must be passed in (rsd). Arrays are processed from index
      56             :    l to index m only. The output information is:
      57             : 
      58             :    ngrp - number of entries of the group
      59             :    ctr  - (w(l)+w(l+ngrp-1))/2
      60             :    ae   - average of wr(l),...,wr(l+ngrp-1)
      61             :    arsd - average of rsd(l),...,rsd(l+ngrp-1)
      62             : */
      63         843 : static PetscErrorCode EPSSubspaceFindGroup(PetscInt l,PetscInt m,PetscScalar *wr,PetscScalar *wi,PetscReal *rsd,PetscReal grptol,PetscInt *ngrp,PetscReal *ctr,PetscReal *ae,PetscReal *arsd)
      64             : {
      65         843 :   PetscInt  i;
      66         843 :   PetscReal rmod,rmod1;
      67             : 
      68         843 :   PetscFunctionBegin;
      69         843 :   *ngrp = 0;
      70         843 :   *ctr = 0;
      71         843 :   rmod = SlepcAbsEigenvalue(wr[l],wi[l]);
      72             : 
      73        2793 :   for (i=l;i<m;) {
      74        1933 :     rmod1 = SlepcAbsEigenvalue(wr[i],wi[i]);
      75        1933 :     if (PetscAbsReal(rmod-rmod1) > grptol*(rmod+rmod1)) break;
      76        1107 :     *ctr = (rmod+rmod1)/2.0;
      77        1107 :     if (wi[i] != 0.0) {
      78           0 :       (*ngrp)+=2;
      79           0 :       i+=2;
      80             :     } else {
      81        1107 :       (*ngrp)++;
      82        1107 :       i++;
      83             :     }
      84             :   }
      85             : 
      86         843 :   *ae = 0;
      87         843 :   *arsd = 0;
      88         843 :   if (*ngrp) {
      89        1950 :     for (i=l;i<l+*ngrp;i++) {
      90        1107 :       (*ae) += PetscRealPart(wr[i]);
      91        1107 :       (*arsd) += rsd[i]*rsd[i];
      92             :     }
      93         843 :     *ae = *ae / *ngrp;
      94         843 :     *arsd = PetscSqrtScalar(*arsd / *ngrp);
      95             :   }
      96         843 :   PetscFunctionReturn(0);
      97             : }
      98             : 
      99             : /*
     100             :    EPSSubspaceResidualNorms - Computes the column norms of residual vectors
     101             :    OP*V(1:n,l:m) - V*T(1:m,l:m), where, on entry, OP*V has been computed and
     102             :    stored in R. On exit, rsd(l) to rsd(m) contain the computed norms.
     103             : */
     104         387 : static PetscErrorCode EPSSubspaceResidualNorms(BV R,BV V,Mat T,PetscInt l,PetscInt m,PetscScalar *eigi,PetscReal *rsd)
     105             : {
     106         387 :   PetscErrorCode ierr;
     107         387 :   PetscInt       i;
     108             : 
     109         387 :   PetscFunctionBegin;
     110         387 :   ierr = BVMult(R,-1.0,1.0,V,T);CHKERRQ(ierr);
     111        4783 :   for (i=l;i<m;i++) { ierr = BVNormColumnBegin(R,i,NORM_2,rsd+i);CHKERRQ(ierr); }
     112        4783 :   for (i=l;i<m;i++) { ierr = BVNormColumnEnd(R,i,NORM_2,rsd+i);CHKERRQ(ierr); }
     113             : #if !defined(PETSC_USE_COMPLEX)
     114        4396 :   for (i=l;i<m-1;i++) {
     115        4009 :     if (eigi[i]!=0.0) {
     116           0 :       rsd[i]   = SlepcAbs(rsd[i],rsd[i+1]);
     117           0 :       rsd[i+1] = rsd[i];
     118           0 :       i++;
     119             :     }
     120             :   }
     121             : #endif
     122         387 :   PetscFunctionReturn(0);
     123             : }
     124             : 
     125          17 : PetscErrorCode EPSSolve_Subspace(EPS eps)
     126             : {
     127          17 :   PetscErrorCode ierr;
     128          17 :   Mat            H,Q,S,T,B;
     129          17 :   BV             AV,R;
     130          17 :   PetscBool      indef;
     131          17 :   PetscInt       i,k,ld,ngrp,nogrp,*itrsd,*itrsdold;
     132          17 :   PetscInt       nxtsrr,idsrr,idort,nxtort,nv,ncv = eps->ncv,its;
     133          17 :   PetscReal      arsd,oarsd,ctr,octr,ae,oae,*rsd,tcond=1.0;
     134             :   /* Parameters */
     135          17 :   PetscInt       init = 5;        /* Number of initial iterations */
     136          17 :   PetscReal      stpfac = 1.5;    /* Max num of iter before next SRR step */
     137          17 :   PetscReal      alpha = 1.0;     /* Used to predict convergence of next residual */
     138          17 :   PetscReal      beta = 1.1;      /* Used to predict convergence of next residual */
     139          17 :   PetscReal      grptol = SLEPC_DEFAULT_TOL;   /* Tolerance for EPSSubspaceFindGroup */
     140          17 :   PetscReal      cnvtol = 1e-6;   /* Convergence criterion for cnv */
     141          17 :   PetscInt       orttol = 2;      /* Number of decimal digits whose loss
     142             :                                      can be tolerated in orthogonalization */
     143             : 
     144          17 :   PetscFunctionBegin;
     145          17 :   its = 0;
     146          17 :   ierr = PetscMalloc3(ncv,&rsd,ncv,&itrsd,ncv,&itrsdold);CHKERRQ(ierr);
     147          17 :   ierr = DSGetLeadingDimension(eps->ds,&ld);CHKERRQ(ierr);
     148          17 :   ierr = BVDuplicate(eps->V,&AV);CHKERRQ(ierr);
     149          17 :   ierr = BVDuplicate(eps->V,&R);CHKERRQ(ierr);
     150          17 :   ierr = STGetOperator(eps->st,&S);CHKERRQ(ierr);
     151             : 
     152         306 :   for (i=0;i<ncv;i++) {
     153         289 :     rsd[i] = 0.0;
     154         289 :     itrsd[i] = -1;
     155             :   }
     156             : 
     157             :   /* Complete the initial basis with random vectors and orthonormalize them */
     158         304 :   for (k=eps->nini;k<ncv;k++) {
     159         287 :     ierr = BVSetRandomColumn(eps->V,k);CHKERRQ(ierr);
     160         287 :     ierr = BVOrthonormalizeColumn(eps->V,k,PETSC_TRUE,NULL,NULL);CHKERRQ(ierr);
     161             :   }
     162             : 
     163         387 :   while (eps->reason == EPS_CONVERGED_ITERATING) {
     164         387 :     eps->its++;
     165         387 :     nv = PetscMin(eps->nconv+eps->mpd,ncv);
     166         387 :     ierr = DSSetDimensions(eps->ds,nv,eps->nconv,0);CHKERRQ(ierr);
     167             : 
     168             :     /* Find group in previously computed eigenvalues */
     169         387 :     ierr = EPSSubspaceFindGroup(eps->nconv,nv,eps->eigr,eps->eigi,rsd,grptol,&nogrp,&octr,&oae,&oarsd);CHKERRQ(ierr);
     170             : 
     171             :     /* AV(:,idx) = OP * V(:,idx) */
     172         387 :     ierr = BVSetActiveColumns(eps->V,eps->nconv,nv);CHKERRQ(ierr);
     173         387 :     ierr = BVSetActiveColumns(AV,eps->nconv,nv);CHKERRQ(ierr);
     174         387 :     ierr = BVMatMult(eps->V,S,AV);CHKERRQ(ierr);
     175             : 
     176             :     /* T(:,idx) = V' * AV(:,idx) */
     177         387 :     ierr = BVSetActiveColumns(eps->V,0,nv);CHKERRQ(ierr);
     178         387 :     ierr = DSGetMat(eps->ds,DS_MAT_A,&H);CHKERRQ(ierr);
     179         387 :     ierr = BVDot(AV,eps->V,H);CHKERRQ(ierr);
     180         387 :     ierr = DSRestoreMat(eps->ds,DS_MAT_A,&H);CHKERRQ(ierr);
     181         387 :     ierr = DSSetState(eps->ds,DS_STATE_RAW);CHKERRQ(ierr);
     182             : 
     183             :     /* Solve projected problem */
     184         387 :     ierr = DSSolve(eps->ds,eps->eigr,eps->eigi);CHKERRQ(ierr);
     185         387 :     ierr = DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL);CHKERRQ(ierr);
     186         387 :     ierr = DSSynchronize(eps->ds,eps->eigr,eps->eigi);CHKERRQ(ierr);
     187             : 
     188             :     /* Update vectors V(:,idx) = V * U(:,idx) */
     189         387 :     ierr = DSGetMat(eps->ds,DS_MAT_Q,&Q);CHKERRQ(ierr);
     190         387 :     ierr = BVSetActiveColumns(AV,0,nv);CHKERRQ(ierr);
     191         387 :     ierr = BVSetActiveColumns(R,0,nv);CHKERRQ(ierr);
     192         387 :     ierr = BVMultInPlace(eps->V,Q,eps->nconv,nv);CHKERRQ(ierr);
     193         387 :     ierr = BVMultInPlace(AV,Q,eps->nconv,nv);CHKERRQ(ierr);
     194         387 :     ierr = BVCopy(AV,R);CHKERRQ(ierr);
     195         387 :     ierr = MatDestroy(&Q);CHKERRQ(ierr);
     196             : 
     197             :     /* Convergence check */
     198         387 :     ierr = DSGetMat(eps->ds,DS_MAT_A,&T);CHKERRQ(ierr);
     199         387 :     ierr = EPSSubspaceResidualNorms(R,eps->V,T,eps->nconv,nv,eps->eigi,rsd);CHKERRQ(ierr);
     200         387 :     ierr = DSRestoreMat(eps->ds,DS_MAT_A,&T);CHKERRQ(ierr);
     201             : 
     202        4783 :     for (i=eps->nconv;i<nv;i++) {
     203        4396 :       itrsdold[i] = itrsd[i];
     204        4396 :       itrsd[i] = its;
     205        4396 :       eps->errest[i] = rsd[i];
     206             :     }
     207             : 
     208         456 :     for (;;) {
     209             :       /* Find group in currently computed eigenvalues */
     210         456 :       ierr = EPSSubspaceFindGroup(eps->nconv,nv,eps->eigr,eps->eigi,eps->errest,grptol,&ngrp,&ctr,&ae,&arsd);CHKERRQ(ierr);
     211         456 :       if (ngrp!=nogrp) break;
     212         436 :       if (ngrp==0) break;
     213         436 :       if (PetscAbsReal(ae-oae)>ctr*cnvtol*(itrsd[eps->nconv]-itrsdold[eps->nconv])) break;
     214         290 :       if (arsd>ctr*eps->tol) break;
     215          69 :       eps->nconv = eps->nconv + ngrp;
     216          69 :       if (eps->nconv>=nv) break;
     217             :     }
     218             : 
     219         387 :     ierr = EPSMonitor(eps,eps->its,eps->nconv,eps->eigr,eps->eigi,eps->errest,nv);CHKERRQ(ierr);
     220         387 :     ierr = (*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx);CHKERRQ(ierr);
     221         387 :     if (eps->reason != EPS_CONVERGED_ITERATING) break;
     222             : 
     223             :     /* Compute nxtsrr (iteration of next projection step) */
     224         370 :     nxtsrr = PetscMin(eps->max_it,PetscMax((PetscInt)PetscFloorReal(stpfac*its),init));
     225             : 
     226         370 :     if (ngrp!=nogrp || ngrp==0 || arsd>=oarsd) {
     227          55 :       idsrr = nxtsrr - its;
     228             :     } else {
     229         315 :       idsrr = (PetscInt)PetscFloorReal(alpha+beta*(itrsdold[eps->nconv]-itrsd[eps->nconv])*PetscLogReal(arsd/eps->tol)/PetscLogReal(arsd/oarsd));
     230         315 :       idsrr = PetscMax(1,idsrr);
     231             :     }
     232         370 :     nxtsrr = PetscMin(nxtsrr,its+idsrr);
     233             : 
     234             :     /* Compute nxtort (iteration of next orthogonalization step) */
     235         370 :     ierr = DSCond(eps->ds,&tcond);CHKERRQ(ierr);
     236         370 :     idort = PetscMax(1,(PetscInt)PetscFloorReal(orttol/PetscMax(1,PetscLog10Real(tcond))));
     237         370 :     nxtort = PetscMin(its+idort,nxtsrr);
     238             : 
     239             :     /* V(:,idx) = AV(:,idx) */
     240         370 :     ierr = BVSetActiveColumns(eps->V,eps->nconv,nv);CHKERRQ(ierr);
     241         370 :     ierr = BVSetActiveColumns(AV,eps->nconv,nv);CHKERRQ(ierr);
     242         370 :     ierr = BVCopy(AV,eps->V);CHKERRQ(ierr);
     243         370 :     its++;
     244             : 
     245             :     /* Orthogonalization loop */
     246        2920 :     do {
     247        2920 :       ierr = BVGetMatrix(eps->V,&B,&indef);CHKERRQ(ierr);
     248        2920 :       ierr = BVSetMatrix(eps->V,NULL,PETSC_FALSE);CHKERRQ(ierr);
     249        7979 :       while (its<nxtort) {
     250             :         /* A(:,idx) = OP*V(:,idx) with normalization */
     251        5059 :         ierr = BVMatMult(eps->V,S,AV);CHKERRQ(ierr);
     252        5059 :         ierr = BVCopy(AV,eps->V);CHKERRQ(ierr);
     253        5059 :         ierr = BVNormalize(eps->V,NULL);CHKERRQ(ierr);
     254        5059 :         its++;
     255             :       }
     256        2920 :       ierr = BVSetMatrix(eps->V,B,indef);CHKERRQ(ierr);
     257             :       /* Orthonormalize vectors */
     258        2920 :       ierr = BVOrthogonalize(eps->V,NULL);CHKERRQ(ierr);
     259        2920 :       nxtort = PetscMin(its+idort,nxtsrr);
     260        2920 :     } while (its<nxtsrr);
     261             :   }
     262             : 
     263          17 :   ierr = PetscFree3(rsd,itrsd,itrsdold);CHKERRQ(ierr);
     264          17 :   ierr = BVDestroy(&AV);CHKERRQ(ierr);
     265          17 :   ierr = BVDestroy(&R);CHKERRQ(ierr);
     266          17 :   ierr = STRestoreOperator(eps->st,&S);CHKERRQ(ierr);
     267          17 :   ierr = DSTruncate(eps->ds,eps->nconv,PETSC_TRUE);CHKERRQ(ierr);
     268          17 :   PetscFunctionReturn(0);
     269             : }
     270             : 
     271          15 : SLEPC_EXTERN PetscErrorCode EPSCreate_Subspace(EPS eps)
     272             : {
     273          15 :   PetscFunctionBegin;
     274          15 :   eps->useds = PETSC_TRUE;
     275          15 :   eps->categ = EPS_CATEGORY_OTHER;
     276             : 
     277          15 :   eps->ops->solve          = EPSSolve_Subspace;
     278          15 :   eps->ops->setup          = EPSSetUp_Subspace;
     279          15 :   eps->ops->setupsort      = EPSSetUpSort_Default;
     280          15 :   eps->ops->backtransform  = EPSBackTransform_Default;
     281          15 :   eps->ops->computevectors = EPSComputeVectors_Schur;
     282          15 :   PetscFunctionReturn(0);
     283             : }
     284             : 

Generated by: LCOV version 1.14