LCOV - code coverage report
Current view: top level - eps/impls/davidson - davidson.c (source / functions) Hit Total Coverage
Test: coverage.info Lines: 179 188 95.2 %
Date: 2019-07-19 01:26:51 Functions: 14 14 100.0 %

          Line data    Source code
       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-2019, 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             :    Skeleton of Davidson solver. Actual solvers are GD and JD.
      12             : */
      13             : 
      14             : #include "davidson.h"
      15             : 
      16             : static PetscBool  cited = PETSC_FALSE;
      17             : static const char citation[] =
      18             :   "@Article{slepc-davidson,\n"
      19             :   "   author = \"E. Romero and J. E. Roman\",\n"
      20             :   "   title = \"A parallel implementation of {Davidson} methods for large-scale eigenvalue problems in {SLEPc}\",\n"
      21             :   "   journal = \"{ACM} Trans. Math. Software\",\n"
      22             :   "   volume = \"40\",\n"
      23             :   "   number = \"2\",\n"
      24             :   "   pages = \"13:1--13:29\",\n"
      25             :   "   year = \"2014,\"\n"
      26             :   "   doi = \"https://doi.org/10.1145/2543696\"\n"
      27             :   "}\n";
      28             : 
      29          90 : PetscErrorCode EPSSetUp_XD(EPS eps)
      30             : {
      31             :   PetscErrorCode ierr;
      32          90 :   EPS_DAVIDSON   *data = (EPS_DAVIDSON*)eps->data;
      33          90 :   dvdDashboard   *dvd = &data->ddb;
      34             :   dvdBlackboard  b;
      35             :   PetscInt       min_size_V,bs,initv,nmat;
      36             :   Mat            A,B;
      37             :   KSP            ksp;
      38             :   PetscBool      ipB,ispositive;
      39             :   HarmType_t     harm;
      40             :   InitType_t     init;
      41             :   PetscScalar    target;
      42             : 
      43          90 :   PetscFunctionBegin;
      44             :   /* Setup EPS options and get the problem specification */
      45          90 :   bs = data->blocksize;
      46          90 :   if (bs <= 0) bs = 1;
      47          90 :   if (eps->ncv) {
      48          39 :     if (eps->ncv<eps->nev) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The value of ncv must be at least nev");
      49          51 :   } else if (eps->mpd) eps->ncv = eps->mpd + eps->nev + bs;
      50          50 :   else if (eps->nev<500) eps->ncv = PetscMin(eps->n-bs,PetscMax(2*eps->nev,eps->nev+15))+bs;
      51           0 :   else eps->ncv = PetscMin(eps->n-bs,eps->nev+500)+bs;
      52          90 :   if (!eps->mpd) eps->mpd = eps->ncv;
      53          90 :   if (eps->mpd > eps->ncv) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The mpd has to be less or equal than ncv");
      54          90 :   if (eps->mpd < 2) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The mpd has to be greater than 2");
      55          90 :   if (!eps->max_it) eps->max_it = PetscMax(100*eps->ncv,2*eps->n);
      56          90 :   if (!eps->which) eps->which = EPS_LARGEST_MAGNITUDE;
      57          90 :   if (eps->ishermitian && (eps->which==EPS_LARGEST_IMAGINARY || eps->which==EPS_SMALLEST_IMAGINARY)) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Wrong value of eps->which");
      58          90 :   if (!(eps->nev + bs <= eps->ncv)) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The ncv has to be greater than nev plus blocksize");
      59          90 :   if (eps->trueres) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"-eps_true_residual is temporally disable in this solver.");
      60             : 
      61          90 :   ierr = EPSXDSetRestart_XD(eps,data->minv,data->plusk);CHKERRQ(ierr);
      62          90 :   min_size_V = data->minv;
      63          90 :   if (!min_size_V) min_size_V = PetscMin(PetscMax(bs,5),eps->mpd/2);
      64          90 :   if (!(min_size_V+bs <= eps->mpd)) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The value of minv must be less than mpd minus blocksize");
      65          90 :   initv = data->initialsize;
      66          90 :   if (eps->mpd < initv) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The initv has to be less or equal than mpd");
      67             : 
      68             :   /* Change the default sigma to inf if necessary */
      69          90 :   if (eps->which == EPS_LARGEST_MAGNITUDE || eps->which == EPS_LARGEST_REAL || eps->which == EPS_LARGEST_IMAGINARY) {
      70          44 :     ierr = STSetDefaultShift(eps->st,PETSC_MAX_REAL);CHKERRQ(ierr);
      71             :   }
      72             : 
      73             :   /* Set up preconditioner */
      74          90 :   ierr = STSetUp(eps->st);CHKERRQ(ierr);
      75             : 
      76             :   /* Setup problem specification in dvd */
      77          90 :   ierr = STGetNumMatrices(eps->st,&nmat);CHKERRQ(ierr);
      78          90 :   ierr = STGetMatrix(eps->st,0,&A);CHKERRQ(ierr);
      79          90 :   if (nmat>1) { ierr = STGetMatrix(eps->st,1,&B);CHKERRQ(ierr); }
      80          90 :   ierr = EPSReset_XD(eps);CHKERRQ(ierr);
      81          90 :   ierr = PetscMemzero(dvd,sizeof(dvdDashboard));CHKERRQ(ierr);
      82          90 :   dvd->A = A; dvd->B = eps->isgeneralized? B: NULL;
      83          90 :   ispositive = eps->ispositive;
      84          90 :   dvd->sA = DVD_MAT_IMPLICIT | (eps->ishermitian? DVD_MAT_HERMITIAN: 0) | ((ispositive && !eps->isgeneralized) ? DVD_MAT_POS_DEF: 0);
      85             :   /* Asume -eps_hermitian means hermitian-definite in generalized problems */
      86          90 :   if (!ispositive && !eps->isgeneralized && eps->ishermitian) ispositive = PETSC_TRUE;
      87          90 :   if (!eps->isgeneralized) dvd->sB = DVD_MAT_IMPLICIT | DVD_MAT_HERMITIAN | DVD_MAT_IDENTITY | DVD_MAT_UNITARY | DVD_MAT_POS_DEF;
      88          12 :   else dvd->sB = DVD_MAT_IMPLICIT | (eps->ishermitian? DVD_MAT_HERMITIAN: 0) | (ispositive? DVD_MAT_POS_DEF: 0);
      89          90 :   ipB = (dvd->B && data->ipB && DVD_IS(dvd->sB,DVD_MAT_HERMITIAN))?PETSC_TRUE:PETSC_FALSE;
      90          90 :   dvd->sEP = ((!eps->isgeneralized || (eps->isgeneralized && ipB))? DVD_EP_STD: 0) | (ispositive? DVD_EP_HERMITIAN: 0) | ((eps->problem_type == EPS_GHIEP && ipB) ? DVD_EP_INDEFINITE : 0);
      91          90 :   if (data->ipB && !ipB) data->ipB = PETSC_FALSE;
      92          90 :   dvd->correctXnorm = (dvd->B && (DVD_IS(dvd->sB,DVD_MAT_HERMITIAN)||DVD_IS(dvd->sEP,DVD_EP_INDEFINITE)))?PETSC_TRUE:PETSC_FALSE;
      93          90 :   dvd->nev        = eps->nev;
      94          90 :   dvd->which      = eps->which;
      95          90 :   dvd->withTarget = PETSC_TRUE;
      96          90 :   switch (eps->which) {
      97             :     case EPS_TARGET_MAGNITUDE:
      98             :     case EPS_TARGET_IMAGINARY:
      99          14 :       dvd->target[0] = target = eps->target;
     100          14 :       dvd->target[1] = 1.0;
     101          14 :       break;
     102             :     case EPS_TARGET_REAL:
     103           0 :       dvd->target[0] = PetscRealPart(target = eps->target);
     104           0 :       dvd->target[1] = 1.0;
     105           0 :       break;
     106             :     case EPS_LARGEST_REAL:
     107             :     case EPS_LARGEST_MAGNITUDE:
     108             :     case EPS_LARGEST_IMAGINARY: /* TODO: think about this case */
     109          44 :       dvd->target[0] = 1.0;
     110          44 :       dvd->target[1] = target = 0.0;
     111          44 :       break;
     112             :     case EPS_SMALLEST_MAGNITUDE:
     113             :     case EPS_SMALLEST_REAL:
     114             :     case EPS_SMALLEST_IMAGINARY: /* TODO: think about this case */
     115          28 :       dvd->target[0] = target = 0.0;
     116          28 :       dvd->target[1] = 1.0;
     117          28 :       break;
     118             :     case EPS_WHICH_USER:
     119           4 :       ierr = STGetShift(eps->st,&target);CHKERRQ(ierr);
     120           4 :       dvd->target[0] = target;
     121           4 :       dvd->target[1] = 1.0;
     122           4 :       break;
     123             :     case EPS_ALL:
     124           0 :       SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported option: which == EPS_ALL");
     125             :       break;
     126             :     default:
     127           0 :       SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported value of option 'which'");
     128             :   }
     129          90 :   dvd->tol = (eps->tol==PETSC_DEFAULT)? SLEPC_DEFAULT_TOL: eps->tol;
     130          90 :   dvd->eps = eps;
     131             : 
     132             :   /* Setup the extraction technique */
     133          90 :   if (!eps->extraction) {
     134          82 :     if (ipB || ispositive) eps->extraction = EPS_RITZ;
     135             :     else {
     136          16 :       switch (eps->which) {
     137             :         case EPS_TARGET_REAL:
     138             :         case EPS_TARGET_MAGNITUDE:
     139             :         case EPS_TARGET_IMAGINARY:
     140             :         case EPS_SMALLEST_MAGNITUDE:
     141             :         case EPS_SMALLEST_REAL:
     142             :         case EPS_SMALLEST_IMAGINARY:
     143           6 :           eps->extraction = EPS_HARMONIC;
     144           6 :           break;
     145             :         case EPS_LARGEST_REAL:
     146             :         case EPS_LARGEST_MAGNITUDE:
     147             :         case EPS_LARGEST_IMAGINARY:
     148           7 :           eps->extraction = EPS_HARMONIC_LARGEST;
     149           7 :           break;
     150             :         default:
     151           3 :           eps->extraction = EPS_RITZ;
     152             :       }
     153             :     }
     154             :   }
     155          90 :   switch (eps->extraction) {
     156             :     case EPS_RITZ:              harm = DVD_HARM_NONE; break;
     157          13 :     case EPS_HARMONIC:          harm = DVD_HARM_RR; break;
     158           0 :     case EPS_HARMONIC_RELATIVE: harm = DVD_HARM_RRR; break;
     159           0 :     case EPS_HARMONIC_RIGHT:    harm = DVD_HARM_REIGS; break;
     160           8 :     case EPS_HARMONIC_LARGEST:  harm = DVD_HARM_LEIGS; break;
     161           0 :     default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");
     162             :   }
     163             : 
     164             :   /* Setup the type of starting subspace */
     165          90 :   init = data->krylovstart? DVD_INITV_KRYLOV: DVD_INITV_CLASSIC;
     166             : 
     167             :   /* Preconfigure dvd */
     168          90 :   ierr = STGetKSP(eps->st,&ksp);CHKERRQ(ierr);
     169          90 :   ierr = dvd_schm_basic_preconf(dvd,&b,eps->mpd,min_size_V,bs,initv,PetscAbs(eps->nini),data->plusk,harm,ksp,init,eps->trackall,data->ipB,data->doubleexp);CHKERRQ(ierr);
     170             : 
     171             :   /* Allocate memory */
     172          90 :   ierr = EPSAllocateSolution(eps,0);CHKERRQ(ierr);
     173             : 
     174             :   /* Setup orthogonalization */
     175          90 :   ierr = EPS_SetInnerProduct(eps);CHKERRQ(ierr);
     176          90 :   if (!(ipB && dvd->B)) {
     177          82 :     ierr = BVSetMatrix(eps->V,NULL,PETSC_FALSE);CHKERRQ(ierr);
     178             :   }
     179             : 
     180             :   /* Configure dvd for a basic GD */
     181          90 :   ierr = dvd_schm_basic_conf(dvd,&b,eps->mpd,min_size_V,bs,initv,PetscAbs(eps->nini),data->plusk,harm,dvd->withTarget,target,ksp,data->fix,init,eps->trackall,data->ipB,data->dynamic,data->doubleexp);CHKERRQ(ierr);
     182          90 :   PetscFunctionReturn(0);
     183             : }
     184             : 
     185          90 : PetscErrorCode EPSSolve_XD(EPS eps)
     186             : {
     187          90 :   EPS_DAVIDSON   *data = (EPS_DAVIDSON*)eps->data;
     188          90 :   dvdDashboard   *d = &data->ddb;
     189             :   PetscInt       l,k;
     190             :   PetscErrorCode ierr;
     191             : 
     192          90 :   PetscFunctionBegin;
     193          90 :   ierr = PetscCitationsRegister(citation,&cited);CHKERRQ(ierr);
     194             :   /* Call the starting routines */
     195          90 :   ierr = EPSDavidsonFLCall(d->startList,d);CHKERRQ(ierr);
     196             : 
     197        9949 :   while (eps->reason == EPS_CONVERGED_ITERATING) {
     198             : 
     199             :     /* Initialize V, if it is needed */
     200        9949 :     ierr = BVGetActiveColumns(d->eps->V,&l,&k);CHKERRQ(ierr);
     201        9949 :     if (l == k) { ierr = d->initV(d);CHKERRQ(ierr); }
     202             : 
     203             :     /* Find the best approximated eigenpairs in V, X */
     204        9949 :     ierr = d->calcPairs(d);CHKERRQ(ierr);
     205             : 
     206             :     /* Test for convergence */
     207        9949 :     ierr = (*eps->stopping)(eps,eps->its,eps->max_it,eps->nconv,eps->nev,&eps->reason,eps->stoppingctx);CHKERRQ(ierr);
     208        9949 :     if (eps->reason != EPS_CONVERGED_ITERATING) break;
     209             : 
     210             :     /* Expand the subspace */
     211        9859 :     ierr = d->updateV(d);CHKERRQ(ierr);
     212             : 
     213             :     /* Monitor progress */
     214        9859 :     eps->nconv = d->nconv;
     215        9859 :     eps->its++;
     216        9859 :     ierr = BVGetActiveColumns(d->eps->V,NULL,&k);CHKERRQ(ierr);
     217        9859 :     ierr = EPSMonitor(eps,eps->its,eps->nconv+d->npreconv,eps->eigr,eps->eigi,eps->errest,PetscMin(k,eps->nev));CHKERRQ(ierr);
     218             :   }
     219             : 
     220             :   /* Call the ending routines */
     221          90 :   ierr = EPSDavidsonFLCall(d->endList,d);CHKERRQ(ierr);
     222          90 :   PetscFunctionReturn(0);
     223             : }
     224             : 
     225         155 : PetscErrorCode EPSReset_XD(EPS eps)
     226             : {
     227         155 :   EPS_DAVIDSON   *data = (EPS_DAVIDSON*)eps->data;
     228         155 :   dvdDashboard   *dvd = &data->ddb;
     229             :   PetscErrorCode ierr;
     230             : 
     231         155 :   PetscFunctionBegin;
     232             :   /* Call step destructors and destroys the list */
     233         155 :   ierr = EPSDavidsonFLCall(dvd->destroyList,dvd);CHKERRQ(ierr);
     234         155 :   ierr = EPSDavidsonFLDestroy(&dvd->destroyList);CHKERRQ(ierr);
     235         155 :   ierr = EPSDavidsonFLDestroy(&dvd->startList);CHKERRQ(ierr);
     236         155 :   ierr = EPSDavidsonFLDestroy(&dvd->endList);CHKERRQ(ierr);
     237         155 :   PetscFunctionReturn(0);
     238             : }
     239             : 
     240           4 : PetscErrorCode EPSXDSetKrylovStart_XD(EPS eps,PetscBool krylovstart)
     241             : {
     242           4 :   EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
     243             : 
     244           4 :   PetscFunctionBegin;
     245           4 :   data->krylovstart = krylovstart;
     246           4 :   PetscFunctionReturn(0);
     247             : }
     248             : 
     249          55 : PetscErrorCode EPSXDGetKrylovStart_XD(EPS eps,PetscBool *krylovstart)
     250             : {
     251          55 :   EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
     252             : 
     253          55 :   PetscFunctionBegin;
     254          55 :   *krylovstart = data->krylovstart;
     255          55 :   PetscFunctionReturn(0);
     256             : }
     257             : 
     258           2 : PetscErrorCode EPSXDSetBlockSize_XD(EPS eps,PetscInt blocksize)
     259             : {
     260           2 :   EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
     261             : 
     262           2 :   PetscFunctionBegin;
     263           2 :   if (blocksize == PETSC_DEFAULT || blocksize == PETSC_DECIDE) blocksize = 1;
     264           2 :   if (blocksize <= 0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid blocksize value");
     265           2 :   data->blocksize = blocksize;
     266           2 :   PetscFunctionReturn(0);
     267             : }
     268             : 
     269          55 : PetscErrorCode EPSXDGetBlockSize_XD(EPS eps,PetscInt *blocksize)
     270             : {
     271          55 :   EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
     272             : 
     273          55 :   PetscFunctionBegin;
     274          55 :   *blocksize = data->blocksize;
     275          55 :   PetscFunctionReturn(0);
     276             : }
     277             : 
     278          93 : PetscErrorCode EPSXDSetRestart_XD(EPS eps,PetscInt minv,PetscInt plusk)
     279             : {
     280          93 :   EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
     281             : 
     282          93 :   PetscFunctionBegin;
     283          93 :   if (minv == PETSC_DEFAULT || minv == PETSC_DECIDE) minv = 5;
     284          93 :   if (minv <= 0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid minv value");
     285          93 :   if (plusk == PETSC_DEFAULT || plusk == PETSC_DECIDE) plusk = eps->problem_type == EPS_GHIEP?0:1;
     286          93 :   if (plusk < 0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid plusk value");
     287          93 :   data->minv = minv;
     288          93 :   data->plusk = plusk;
     289          93 :   PetscFunctionReturn(0);
     290             : }
     291             : 
     292          55 : PetscErrorCode EPSXDGetRestart_XD(EPS eps,PetscInt *minv,PetscInt *plusk)
     293             : {
     294          55 :   EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
     295             : 
     296          55 :   PetscFunctionBegin;
     297          55 :   if (minv) *minv = data->minv;
     298          55 :   if (plusk) *plusk = data->plusk;
     299          55 :   PetscFunctionReturn(0);
     300             : }
     301             : 
     302          54 : PetscErrorCode EPSXDGetInitialSize_XD(EPS eps,PetscInt *initialsize)
     303             : {
     304          54 :   EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
     305             : 
     306          54 :   PetscFunctionBegin;
     307          54 :   *initialsize = data->initialsize;
     308          54 :   PetscFunctionReturn(0);
     309             : }
     310             : 
     311           2 : PetscErrorCode EPSXDSetInitialSize_XD(EPS eps,PetscInt initialsize)
     312             : {
     313           2 :   EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
     314             : 
     315           2 :   PetscFunctionBegin;
     316           2 :   if (initialsize == PETSC_DEFAULT || initialsize == PETSC_DECIDE) initialsize = 5;
     317           2 :   if (initialsize <= 0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid initial size value");
     318           2 :   data->initialsize = initialsize;
     319           2 :   PetscFunctionReturn(0);
     320             : }
     321             : 
     322           2 : PetscErrorCode EPSXDSetBOrth_XD(EPS eps,PetscBool borth)
     323             : {
     324           2 :   EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
     325             : 
     326           2 :   PetscFunctionBegin;
     327           2 :   data->ipB = borth;
     328           2 :   PetscFunctionReturn(0);
     329             : }
     330             : 
     331          55 : PetscErrorCode EPSXDGetBOrth_XD(EPS eps,PetscBool *borth)
     332             : {
     333          55 :   EPS_DAVIDSON *data = (EPS_DAVIDSON*)eps->data;
     334             : 
     335          55 :   PetscFunctionBegin;
     336          55 :   *borth = data->ipB;
     337          55 :   PetscFunctionReturn(0);
     338             : }
     339             : 
     340             : /*
     341             :   EPSComputeVectors_XD - Compute eigenvectors from the vectors
     342             :   provided by the eigensolver. This version is intended for solvers
     343             :   that provide Schur vectors from the QZ decomposition. Given the partial
     344             :   Schur decomposition OP*V=V*T, the following steps are performed:
     345             :       1) compute eigenvectors of (S,T): S*Z=T*Z*D
     346             :       2) compute eigenvectors of OP: X=V*Z
     347             :  */
     348          90 : PetscErrorCode EPSComputeVectors_XD(EPS eps)
     349             : {
     350             :   PetscErrorCode ierr;
     351             :   Mat            X;
     352             :   PetscBool      symm;
     353             : 
     354          90 :   PetscFunctionBegin;
     355          90 :   ierr = PetscObjectTypeCompare((PetscObject)eps->ds,DSHEP,&symm);CHKERRQ(ierr);
     356          90 :   if (symm) PetscFunctionReturn(0);
     357          24 :   ierr = DSVectors(eps->ds,DS_MAT_X,NULL,NULL);CHKERRQ(ierr);
     358             : 
     359             :   /* V <- V * X */
     360          24 :   ierr = DSGetMat(eps->ds,DS_MAT_X,&X);CHKERRQ(ierr);
     361          24 :   ierr = BVSetActiveColumns(eps->V,0,eps->nconv);CHKERRQ(ierr);
     362          24 :   ierr = BVMultInPlace(eps->V,X,0,eps->nconv);CHKERRQ(ierr);
     363          24 :   ierr = MatDestroy(&X);CHKERRQ(ierr);
     364          24 :   PetscFunctionReturn(0);
     365             : }

Generated by: LCOV version 1.13