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: "rqcg"
      12             : 
      13             :    Method: Rayleigh Quotient Conjugate Gradient
      14             : 
      15             :    Algorithm:
      16             : 
      17             :        Conjugate Gradient minimization of the Rayleigh quotient with
      18             :        periodic Rayleigh-Ritz acceleration.
      19             : 
      20             :    References:
      21             : 
      22             :        [1] L. Bergamaschi et al., "Parallel preconditioned conjugate gradient
      23             :            optimization of the Rayleigh quotient for the solution of sparse
      24             :            eigenproblems", Appl. Math. Comput. 175(2):1694-1715, 2006.
      25             : */
      26             : 
      27             : #include <slepc/private/epsimpl.h>                /*I "slepceps.h" I*/
      28             : 
      29             : static PetscErrorCode EPSSolve_RQCG(EPS);
      30             : 
      31             : typedef struct {
      32             :   PetscInt nrest;         /* user-provided reset parameter */
      33             :   PetscInt allocsize;     /* number of columns of work BV's allocated at setup */
      34             :   BV       AV,W,P,G;
      35             : } EPS_RQCG;
      36             : 
      37          14 : static PetscErrorCode EPSSetUp_RQCG(EPS eps)
      38             : {
      39          14 :   PetscInt       nmat;
      40          14 :   EPS_RQCG       *ctx = (EPS_RQCG*)eps->data;
      41             : 
      42          14 :   PetscFunctionBegin;
      43          14 :   EPSCheckHermitianDefinite(eps);
      44          14 :   EPSCheckNotStructured(eps);
      45          14 :   PetscCall(EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd));
      46          14 :   if (eps->max_it==PETSC_DETERMINE) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
      47          14 :   if (!eps->which) eps->which = EPS_SMALLEST_REAL;
      48          14 :   PetscCheck(eps->which==EPS_SMALLEST_REAL,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver supports only smallest real eigenvalues");
      49          14 :   EPSCheckUnsupported(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION);
      50          14 :   EPSCheckIgnored(eps,EPS_FEATURE_BALANCE);
      51             : 
      52          14 :   if (!ctx->nrest) ctx->nrest = 20;
      53             : 
      54          14 :   PetscCall(EPSAllocateSolution(eps,0));
      55          14 :   PetscCall(EPS_SetInnerProduct(eps));
      56             : 
      57          14 :   PetscCall(STGetNumMatrices(eps->st,&nmat));
      58          14 :   if (!ctx->allocsize) {
      59          11 :     ctx->allocsize = eps->mpd;
      60          11 :     PetscCall(BVDuplicateResize(eps->V,eps->mpd,&ctx->AV));
      61          11 :     if (nmat>1) PetscCall(BVDuplicate(ctx->AV,&ctx->W));
      62          11 :     PetscCall(BVDuplicate(ctx->AV,&ctx->P));
      63          11 :     PetscCall(BVDuplicate(ctx->AV,&ctx->G));
      64           3 :   } else if (ctx->allocsize!=eps->mpd) {
      65           1 :     ctx->allocsize = eps->mpd;
      66           1 :     PetscCall(BVResize(ctx->AV,eps->mpd,PETSC_FALSE));
      67           1 :     if (nmat>1) PetscCall(BVResize(ctx->W,eps->mpd,PETSC_FALSE));
      68           1 :     PetscCall(BVResize(ctx->P,eps->mpd,PETSC_FALSE));
      69           1 :     PetscCall(BVResize(ctx->G,eps->mpd,PETSC_FALSE));
      70             :   }
      71          14 :   PetscCall(DSSetType(eps->ds,DSHEP));
      72          14 :   PetscCall(DSAllocate(eps->ds,eps->ncv));
      73          14 :   PetscCall(EPSSetWorkVecs(eps,1));
      74          14 :   PetscFunctionReturn(PETSC_SUCCESS);
      75             : }
      76             : 
      77          14 : static PetscErrorCode EPSSolve_RQCG(EPS eps)
      78             : {
      79          14 :   EPS_RQCG       *ctx = (EPS_RQCG*)eps->data;
      80          14 :   PetscInt       i,j,k,ld,nv,ncv = eps->ncv,kini,nmat;
      81          14 :   PetscScalar    *C,*gamma,g,pap,pbp,pbx,pax,nu,mu,alpha,beta;
      82          14 :   PetscReal      resnorm,a,b,c,d,disc,t;
      83          14 :   PetscBool      reset;
      84          14 :   Mat            A,B,Q,Q1;
      85          14 :   Vec            v,av,bv,p,w=eps->work[0];
      86             : 
      87          14 :   PetscFunctionBegin;
      88          14 :   PetscCall(DSGetLeadingDimension(eps->ds,&ld));
      89          14 :   PetscCall(STGetNumMatrices(eps->st,&nmat));
      90          14 :   PetscCall(STGetMatrix(eps->st,0,&A));
      91          14 :   if (nmat>1) PetscCall(STGetMatrix(eps->st,1,&B));
      92          12 :   else B = NULL;
      93          14 :   PetscCall(PetscMalloc1(eps->mpd,&gamma));
      94             : 
      95          14 :   kini = eps->nini;
      96         593 :   while (eps->reason == EPS_CONVERGED_ITERATING) {
      97         579 :     eps->its++;
      98         579 :     nv = PetscMin(eps->nconv+eps->mpd,ncv);
      99         579 :     PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,0));
     100         857 :     for (;kini<nv;kini++) { /* Generate more initial vectors if necessary */
     101         278 :       PetscCall(BVSetRandomColumn(eps->V,kini));
     102         278 :       PetscCall(BVOrthonormalizeColumn(eps->V,kini,PETSC_TRUE,NULL,NULL));
     103             :     }
     104         579 :     reset = (eps->its>1 && (eps->its-1)%ctx->nrest==0)? PETSC_TRUE: PETSC_FALSE;
     105             : 
     106         579 :     if (reset) {
     107             :       /* Prevent BVDotVec below to use B-product, restored at the end */
     108          36 :       PetscCall(BVSetMatrix(eps->V,NULL,PETSC_FALSE));
     109             : 
     110             :       /* Compute Rayleigh quotient */
     111          36 :       PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
     112          36 :       PetscCall(BVSetActiveColumns(ctx->AV,0,nv-eps->nconv));
     113          36 :       PetscCall(BVMatMult(eps->V,A,ctx->AV));
     114          36 :       PetscCall(DSGetArray(eps->ds,DS_MAT_A,&C));
     115         804 :       for (i=eps->nconv;i<nv;i++) {
     116         768 :         PetscCall(BVSetActiveColumns(eps->V,eps->nconv,i+1));
     117         768 :         PetscCall(BVGetColumn(ctx->AV,i-eps->nconv,&av));
     118         768 :         PetscCall(BVDotVec(eps->V,av,C+eps->nconv+i*ld));
     119         768 :         PetscCall(BVRestoreColumn(ctx->AV,i-eps->nconv,&av));
     120        8962 :         for (j=eps->nconv;j<i-1;j++) C[i+j*ld] = PetscConj(C[j+i*ld]);
     121             :       }
     122          36 :       PetscCall(DSRestoreArray(eps->ds,DS_MAT_A,&C));
     123          36 :       PetscCall(DSSetState(eps->ds,DS_STATE_RAW));
     124             : 
     125             :       /* Solve projected problem */
     126          36 :       PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
     127          36 :       PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,NULL,NULL,NULL));
     128          36 :       PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
     129             : 
     130             :       /* Update vectors V(:,idx) = V * Y(:,idx) */
     131          36 :       PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&Q));
     132          36 :       PetscCall(BVMultInPlace(eps->V,Q,eps->nconv,nv));
     133          36 :       PetscCall(MatDenseGetSubMatrix(Q,eps->nconv,PETSC_DECIDE,eps->nconv,PETSC_DECIDE,&Q1));
     134          36 :       PetscCall(BVMultInPlace(ctx->AV,Q1,0,nv-eps->nconv));
     135          36 :       PetscCall(MatDenseRestoreSubMatrix(Q,&Q1));
     136          36 :       PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&Q));
     137          36 :       if (B) PetscCall(BVSetMatrix(eps->V,B,PETSC_FALSE));
     138             :     } else {
     139             :       /* No need to do Rayleigh-Ritz, just take diag(V'*A*V) */
     140        9938 :       for (i=eps->nconv;i<nv;i++) {
     141        9395 :         PetscCall(BVGetColumn(eps->V,i,&v));
     142        9395 :         PetscCall(BVGetColumn(ctx->AV,i-eps->nconv,&av));
     143        9395 :         PetscCall(MatMult(A,v,av));
     144        9395 :         PetscCall(VecDot(av,v,eps->eigr+i));
     145        9395 :         PetscCall(BVRestoreColumn(eps->V,i,&v));
     146        9395 :         PetscCall(BVRestoreColumn(ctx->AV,i-eps->nconv,&av));
     147             :       }
     148             :     }
     149             : 
     150             :     /* Compute gradient and check convergence */
     151         579 :     k = -1;
     152       10742 :     for (i=eps->nconv;i<nv;i++) {
     153       10163 :       PetscCall(BVGetColumn(eps->V,i,&v));
     154       10163 :       PetscCall(BVGetColumn(ctx->AV,i-eps->nconv,&av));
     155       10163 :       PetscCall(BVGetColumn(ctx->G,i-eps->nconv,&p));
     156       10163 :       if (B) {
     157        1415 :         PetscCall(BVGetColumn(ctx->W,i-eps->nconv,&bv));
     158        1415 :         PetscCall(MatMult(B,v,bv));
     159        1415 :         PetscCall(VecWAXPY(p,-eps->eigr[i],bv,av));
     160        1415 :         PetscCall(BVRestoreColumn(ctx->W,i-eps->nconv,&bv));
     161        8748 :       } else PetscCall(VecWAXPY(p,-eps->eigr[i],v,av));
     162       10163 :       PetscCall(BVRestoreColumn(eps->V,i,&v));
     163       10163 :       PetscCall(BVRestoreColumn(ctx->AV,i-eps->nconv,&av));
     164       10163 :       PetscCall(VecNorm(p,NORM_2,&resnorm));
     165       10163 :       PetscCall(BVRestoreColumn(ctx->G,i-eps->nconv,&p));
     166       10163 :       PetscCall((*eps->converged)(eps,eps->eigr[i],0.0,resnorm,&eps->errest[i],eps->convergedctx));
     167       10163 :       if (k==-1 && eps->errest[i] >= eps->tol) k = i;
     168             :     }
     169         579 :     if (k==-1) k = nv;
     170         579 :     PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
     171             : 
     172             :     /* The next lines are necessary to avoid DS zeroing eigr */
     173         579 :     PetscCall(DSGetArray(eps->ds,DS_MAT_A,&C));
     174         719 :     for (i=eps->nconv;i<k;i++) C[i+i*ld] = eps->eigr[i];
     175         579 :     PetscCall(DSRestoreArray(eps->ds,DS_MAT_A,&C));
     176             : 
     177         579 :     if (eps->reason == EPS_CONVERGED_ITERATING) {
     178             : 
     179             :       /* Search direction */
     180       10471 :       for (i=0;i<nv-eps->nconv;i++) {
     181        9906 :         PetscCall(BVGetColumn(ctx->G,i,&v));
     182        9906 :         PetscCall(STApply(eps->st,v,w));
     183        9906 :         PetscCall(VecDot(w,v,&g));
     184        9906 :         PetscCall(BVRestoreColumn(ctx->G,i,&v));
     185        9906 :         beta = (!reset && eps->its>1)? g/gamma[i]: 0.0;
     186        9906 :         gamma[i] = g;
     187        9906 :         PetscCall(BVGetColumn(ctx->P,i,&v));
     188        9906 :         PetscCall(VecAXPBY(v,1.0,beta,w));
     189        9906 :         if (i+eps->nconv>0) {
     190        9690 :           PetscCall(BVSetActiveColumns(eps->V,0,i+eps->nconv));
     191        9690 :           PetscCall(BVOrthogonalizeVec(eps->V,v,NULL,NULL,NULL));
     192             :         }
     193        9906 :         PetscCall(BVRestoreColumn(ctx->P,i,&v));
     194             :       }
     195             : 
     196             :       /* Minimization problem */
     197       10471 :       for (i=eps->nconv;i<nv;i++) {
     198        9906 :         PetscCall(BVGetColumn(eps->V,i,&v));
     199        9906 :         PetscCall(BVGetColumn(ctx->AV,i-eps->nconv,&av));
     200        9906 :         PetscCall(BVGetColumn(ctx->P,i-eps->nconv,&p));
     201        9906 :         PetscCall(VecDot(av,v,&nu));
     202        9906 :         PetscCall(VecDot(av,p,&pax));
     203        9906 :         PetscCall(MatMult(A,p,w));
     204        9906 :         PetscCall(VecDot(w,p,&pap));
     205        9906 :         if (B) {
     206        1381 :           PetscCall(BVGetColumn(ctx->W,i-eps->nconv,&bv));
     207        1381 :           PetscCall(VecDot(bv,v,&mu));
     208        1381 :           PetscCall(VecDot(bv,p,&pbx));
     209        1381 :           PetscCall(BVRestoreColumn(ctx->W,i-eps->nconv,&bv));
     210        1381 :           PetscCall(MatMult(B,p,w));
     211        1381 :           PetscCall(VecDot(w,p,&pbp));
     212             :         } else {
     213        8525 :           PetscCall(VecDot(v,v,&mu));
     214        8525 :           PetscCall(VecDot(v,p,&pbx));
     215        8525 :           PetscCall(VecDot(p,p,&pbp));
     216             :         }
     217        9906 :         PetscCall(BVRestoreColumn(ctx->AV,i-eps->nconv,&av));
     218        9906 :         a = PetscRealPart(pap*pbx-pax*pbp);
     219        9906 :         b = PetscRealPart(nu*pbp-mu*pap);
     220        9906 :         c = PetscRealPart(mu*pax-nu*pbx);
     221       10493 :         t = PetscMax(PetscMax(PetscAbsReal(a),PetscAbsReal(b)),PetscAbsReal(c));
     222        9906 :         if (t!=0.0) { a /= t; b /= t; c /= t; }
     223        9906 :         disc = b*b-4.0*a*c;
     224        9906 :         d = PetscSqrtReal(PetscAbsReal(disc));
     225        9906 :         if (b>=0.0 && a!=0.0) alpha = (b+d)/(2.0*a);
     226        9066 :         else if (b!=d) alpha = 2.0*c/(b-d);
     227             :         else alpha = 0;
     228             :         /* Next iterate */
     229        9906 :         if (alpha!=0.0) PetscCall(VecAXPY(v,alpha,p));
     230        9906 :         PetscCall(BVRestoreColumn(eps->V,i,&v));
     231        9906 :         PetscCall(BVRestoreColumn(ctx->P,i-eps->nconv,&p));
     232        9906 :         PetscCall(BVOrthonormalizeColumn(eps->V,i,PETSC_TRUE,NULL,NULL));
     233             :       }
     234             :     }
     235             : 
     236         579 :     PetscCall(EPSMonitor(eps,eps->its,k,eps->eigr,eps->eigi,eps->errest,nv));
     237         579 :     eps->nconv = k;
     238             :   }
     239             : 
     240          14 :   PetscCall(PetscFree(gamma));
     241          14 :   PetscFunctionReturn(PETSC_SUCCESS);
     242             : }
     243             : 
     244           2 : static PetscErrorCode EPSRQCGSetReset_RQCG(EPS eps,PetscInt nrest)
     245             : {
     246           2 :   EPS_RQCG *ctx = (EPS_RQCG*)eps->data;
     247             : 
     248           2 :   PetscFunctionBegin;
     249           2 :   if (nrest==PETSC_DEFAULT || nrest==PETSC_DECIDE) {
     250           0 :     ctx->nrest = 0;
     251           0 :     eps->state = EPS_STATE_INITIAL;
     252             :   } else {
     253           2 :     PetscCheck(nrest>0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Reset parameter must be >0");
     254           2 :     ctx->nrest = nrest;
     255             :   }
     256           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     257             : }
     258             : 
     259             : /*@
     260             :    EPSRQCGSetReset - Sets the reset parameter of the RQCG iteration. Every
     261             :    nrest iterations, the solver performs a Rayleigh-Ritz projection step.
     262             : 
     263             :    Logically Collective
     264             : 
     265             :    Input Parameters:
     266             : +  eps - the eigenproblem solver context
     267             : -  nrest - the number of iterations between resets
     268             : 
     269             :    Options Database Key:
     270             : .  -eps_rqcg_reset - Sets the reset parameter
     271             : 
     272             :    Level: advanced
     273             : 
     274             : .seealso: EPSRQCGGetReset()
     275             : @*/
     276           2 : PetscErrorCode EPSRQCGSetReset(EPS eps,PetscInt nrest)
     277             : {
     278           2 :   PetscFunctionBegin;
     279           2 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
     280           6 :   PetscValidLogicalCollectiveInt(eps,nrest,2);
     281           2 :   PetscTryMethod(eps,"EPSRQCGSetReset_C",(EPS,PetscInt),(eps,nrest));
     282           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     283             : }
     284             : 
     285           1 : static PetscErrorCode EPSRQCGGetReset_RQCG(EPS eps,PetscInt *nrest)
     286             : {
     287           1 :   EPS_RQCG *ctx = (EPS_RQCG*)eps->data;
     288             : 
     289           1 :   PetscFunctionBegin;
     290           1 :   *nrest = ctx->nrest;
     291           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     292             : }
     293             : 
     294             : /*@
     295             :    EPSRQCGGetReset - Gets the reset parameter used in the RQCG method.
     296             : 
     297             :    Not Collective
     298             : 
     299             :    Input Parameter:
     300             : .  eps - the eigenproblem solver context
     301             : 
     302             :    Output Parameter:
     303             : .  nrest - the reset parameter
     304             : 
     305             :    Level: advanced
     306             : 
     307             : .seealso: EPSRQCGSetReset()
     308             : @*/
     309           1 : PetscErrorCode EPSRQCGGetReset(EPS eps,PetscInt *nrest)
     310             : {
     311           1 :   PetscFunctionBegin;
     312           1 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
     313           1 :   PetscAssertPointer(nrest,2);
     314           1 :   PetscUseMethod(eps,"EPSRQCGGetReset_C",(EPS,PetscInt*),(eps,nrest));
     315           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     316             : }
     317             : 
     318          11 : static PetscErrorCode EPSReset_RQCG(EPS eps)
     319             : {
     320          11 :   EPS_RQCG       *ctx = (EPS_RQCG*)eps->data;
     321             : 
     322          11 :   PetscFunctionBegin;
     323          11 :   PetscCall(BVDestroy(&ctx->AV));
     324          11 :   PetscCall(BVDestroy(&ctx->W));
     325          11 :   PetscCall(BVDestroy(&ctx->P));
     326          11 :   PetscCall(BVDestroy(&ctx->G));
     327          11 :   ctx->allocsize = 0;
     328          11 :   PetscFunctionReturn(PETSC_SUCCESS);
     329             : }
     330             : 
     331          10 : static PetscErrorCode EPSSetFromOptions_RQCG(EPS eps,PetscOptionItems *PetscOptionsObject)
     332             : {
     333          10 :   PetscBool      flg;
     334          10 :   PetscInt       nrest;
     335             : 
     336          10 :   PetscFunctionBegin;
     337          10 :   PetscOptionsHeadBegin(PetscOptionsObject,"EPS RQCG Options");
     338             : 
     339          10 :     PetscCall(PetscOptionsInt("-eps_rqcg_reset","Reset parameter","EPSRQCGSetReset",20,&nrest,&flg));
     340          10 :     if (flg) PetscCall(EPSRQCGSetReset(eps,nrest));
     341             : 
     342          10 :   PetscOptionsHeadEnd();
     343          10 :   PetscFunctionReturn(PETSC_SUCCESS);
     344             : }
     345             : 
     346          11 : static PetscErrorCode EPSDestroy_RQCG(EPS eps)
     347             : {
     348          11 :   PetscFunctionBegin;
     349          11 :   PetscCall(PetscFree(eps->data));
     350          11 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSRQCGSetReset_C",NULL));
     351          11 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSRQCGGetReset_C",NULL));
     352          11 :   PetscFunctionReturn(PETSC_SUCCESS);
     353             : }
     354             : 
     355           1 : static PetscErrorCode EPSView_RQCG(EPS eps,PetscViewer viewer)
     356             : {
     357           1 :   EPS_RQCG       *ctx = (EPS_RQCG*)eps->data;
     358           1 :   PetscBool      isascii;
     359             : 
     360           1 :   PetscFunctionBegin;
     361           1 :   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
     362           1 :   if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer,"  reset every %" PetscInt_FMT " iterations\n",ctx->nrest));
     363           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     364             : }
     365             : 
     366          11 : SLEPC_EXTERN PetscErrorCode EPSCreate_RQCG(EPS eps)
     367             : {
     368          11 :   EPS_RQCG       *rqcg;
     369             : 
     370          11 :   PetscFunctionBegin;
     371          11 :   PetscCall(PetscNew(&rqcg));
     372          11 :   eps->data = (void*)rqcg;
     373             : 
     374          11 :   eps->useds = PETSC_TRUE;
     375          11 :   eps->categ = EPS_CATEGORY_PRECOND;
     376             : 
     377          11 :   eps->ops->solve          = EPSSolve_RQCG;
     378          11 :   eps->ops->setup          = EPSSetUp_RQCG;
     379          11 :   eps->ops->setupsort      = EPSSetUpSort_Default;
     380          11 :   eps->ops->setfromoptions = EPSSetFromOptions_RQCG;
     381          11 :   eps->ops->destroy        = EPSDestroy_RQCG;
     382          11 :   eps->ops->reset          = EPSReset_RQCG;
     383          11 :   eps->ops->view           = EPSView_RQCG;
     384          11 :   eps->ops->backtransform  = EPSBackTransform_Default;
     385          11 :   eps->ops->setdefaultst   = EPSSetDefaultST_GMRES;
     386             : 
     387          11 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSRQCGSetReset_C",EPSRQCGSetReset_RQCG));
     388          11 :   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSRQCGGetReset_C",EPSRQCGGetReset_RQCG));
     389          11 :   PetscFunctionReturn(PETSC_SUCCESS);
     390             : }