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             :    This file implements a wrapper to the TRLAN package
      12             : */
      13             : 
      14             : #include <slepc/private/epsimpl.h>
      15             : #include "trlan.h"
      16             : 
      17             : /* Nasty global variable to access EPS data from TRLan_ */
      18             : static struct {
      19             :   EPS eps;
      20             :   Vec x,y;
      21             : } globaldata;
      22             : 
      23           7 : static PetscErrorCode EPSSetUp_TRLAN(EPS eps)
      24             : {
      25           7 :   EPS_TRLAN      *tr = (EPS_TRLAN*)eps->data;
      26             : 
      27           7 :   PetscFunctionBegin;
      28           7 :   EPSCheckHermitian(eps);
      29           7 :   EPSCheckStandard(eps);
      30           7 :   EPSCheckNotStructured(eps);
      31           7 :   if (eps->nev==0) eps->nev = 1;
      32           7 :   PetscCall(PetscBLASIntCast(PetscMax(7,eps->nev+PetscMin(eps->nev,6)),&tr->maxlan));
      33           7 :   if (eps->ncv!=PETSC_DETERMINE) {
      34           2 :     PetscCheck(eps->ncv>=eps->nev,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must be at least nev");
      35           5 :   } else eps->ncv = tr->maxlan;
      36           7 :   if (eps->mpd!=PETSC_DETERMINE) PetscCall(PetscInfo(eps,"Warning: parameter mpd ignored\n"));
      37           7 :   if (eps->max_it==PETSC_DETERMINE) eps->max_it = PetscMax(1000,eps->n);
      38             : 
      39           7 :   if (!eps->which) eps->which = EPS_LARGEST_REAL;
      40           7 :   PetscCheck(eps->which==EPS_SMALLEST_REAL || eps->which==EPS_LARGEST_REAL || eps->which==EPS_TARGET_REAL,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver supports only smallest, largest or target real eigenvalues");
      41           7 :   EPSCheckUnsupported(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_CONVERGENCE | EPS_FEATURE_STOPPING);
      42           7 :   EPSCheckIgnored(eps,EPS_FEATURE_BALANCE | EPS_FEATURE_EXTRACTION);
      43             : 
      44           7 :   tr->restart = 0;
      45           7 :   if (tr->maxlan+1-eps->ncv<=0) PetscCall(PetscBLASIntCast(tr->maxlan*(tr->maxlan+10),&tr->lwork));
      46           7 :   else PetscCall(PetscBLASIntCast(eps->nloc*(tr->maxlan+1-eps->ncv) + tr->maxlan*(tr->maxlan+10),&tr->lwork));
      47           7 :   if (tr->work) PetscCall(PetscFree(tr->work));
      48           7 :   PetscCall(PetscMalloc1(tr->lwork,&tr->work));
      49             : 
      50           7 :   PetscCall(EPSAllocateSolution(eps,0));
      51           7 :   PetscFunctionReturn(PETSC_SUCCESS);
      52             : }
      53             : 
      54         657 : static PetscBLASInt MatMult_TRLAN(PetscBLASInt *n,PetscBLASInt *m,PetscReal *xin,PetscBLASInt *ldx,PetscReal *yout,PetscBLASInt *ldy)
      55             : {
      56         657 :   Vec            x=globaldata.x,y=globaldata.y;
      57         657 :   EPS            eps=globaldata.eps;
      58         657 :   PetscBLASInt   i;
      59             : 
      60         657 :   PetscFunctionBegin;
      61        1314 :   for (i=0;i<*m;i++) {
      62         657 :     PetscCall(VecPlaceArray(x,(PetscScalar*)xin+i*(*ldx)));
      63         657 :     PetscCall(VecPlaceArray(y,(PetscScalar*)yout+i*(*ldy)));
      64         657 :     PetscCall(STApply(eps->st,x,y));
      65         657 :     PetscCall(BVOrthogonalizeVec(eps->V,y,NULL,NULL,NULL));
      66         657 :     PetscCall(VecResetArray(x));
      67         657 :     PetscCall(VecResetArray(y));
      68             :   }
      69         657 :   PetscFunctionReturn(PETSC_SUCCESS);
      70             : }
      71             : 
      72           7 : static PetscErrorCode EPSSolve_TRLAN(EPS eps)
      73             : {
      74           7 :   PetscInt       i;
      75           7 :   PetscBLASInt   ipar[32],n,lohi,stat,ncv;
      76           7 :   EPS_TRLAN      *tr = (EPS_TRLAN*)eps->data;
      77           7 :   PetscScalar    *pV;
      78           7 :   Vec            v0;
      79           7 :   Mat            A;
      80             : #if !defined(PETSC_HAVE_MPIUNI)
      81           7 :   MPI_Fint       fcomm;
      82             : #endif
      83             : 
      84           7 :   PetscFunctionBegin;
      85           7 :   PetscCall(PetscBLASIntCast(eps->ncv,&ncv));
      86           7 :   PetscCall(PetscBLASIntCast(eps->nloc,&n));
      87             : 
      88           7 :   PetscCheck(eps->which==EPS_LARGEST_REAL || eps->which==EPS_TARGET_REAL || eps->which==EPS_SMALLEST_REAL,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"Wrong value of eps->which");
      89           7 :   lohi = (eps->which==EPS_SMALLEST_REAL)? -1: 1;
      90             : 
      91           7 :   globaldata.eps = eps;
      92           7 :   PetscCall(STGetMatrix(eps->st,0,&A));
      93           7 :   PetscCall(MatCreateVecsEmpty(A,&globaldata.x,&globaldata.y));
      94             : 
      95           7 :   ipar[0]  = 0;            /* stat: error flag */
      96           7 :   ipar[1]  = lohi;         /* smallest (lohi<0) or largest eigenvalues (lohi>0) */
      97           7 :   PetscCall(PetscBLASIntCast(eps->nev,&ipar[2])); /* number of desired eigenpairs */
      98           7 :   ipar[3]  = 0;            /* number of eigenpairs already converged */
      99           7 :   ipar[4]  = tr->maxlan;   /* maximum Lanczos basis size */
     100           7 :   ipar[5]  = tr->restart;  /* restarting scheme */
     101           7 :   PetscCall(PetscBLASIntCast(eps->max_it,&ipar[6])); /* maximum number of MATVECs */
     102             : #if !defined(PETSC_HAVE_MPIUNI)
     103           7 :   fcomm    = MPI_Comm_c2f(PetscObjectComm((PetscObject)eps));
     104           7 :   ipar[7]  = fcomm;
     105             : #endif
     106           7 :   ipar[8]  = 0;            /* verboseness */
     107           7 :   ipar[9]  = 99;           /* Fortran IO unit number used to write log messages */
     108           7 :   ipar[10] = 1;            /* use supplied starting vector */
     109           7 :   ipar[11] = 0;            /* checkpointing flag */
     110           7 :   ipar[12] = 98;           /* Fortran IO unit number used to write checkpoint files */
     111           7 :   ipar[13] = 0;            /* number of flops per matvec per PE (not used) */
     112           7 :   tr->work[0] = eps->tol;  /* relative tolerance on residual norms */
     113             : 
     114          63 :   for (i=0;i<eps->ncv;i++) eps->eigr[i]=0.0;
     115           7 :   PetscCall(EPSGetStartVector(eps,0,NULL));
     116           7 :   PetscCall(BVSetActiveColumns(eps->V,0,0));  /* just for deflation space */
     117           7 :   PetscCall(BVGetColumn(eps->V,0,&v0));
     118           7 :   PetscCall(VecGetArray(v0,&pV));
     119             : 
     120           7 :   PetscStackCallExternalVoid("TRLan",TRLan_(MatMult_TRLAN,ipar,&n,&ncv,eps->eigr,pV,&n,tr->work,&tr->lwork));
     121             : 
     122           7 :   PetscCall(VecRestoreArray(v0,&pV));
     123           7 :   PetscCall(BVRestoreColumn(eps->V,0,&v0));
     124             : 
     125           7 :   stat        = ipar[0];
     126           7 :   eps->nconv  = ipar[3];
     127           7 :   eps->its    = ipar[25];
     128           7 :   eps->reason = EPS_CONVERGED_TOL;
     129             : 
     130           7 :   PetscCall(VecDestroy(&globaldata.x));
     131           7 :   PetscCall(VecDestroy(&globaldata.y));
     132           7 :   PetscCheck(stat==0,PetscObjectComm((PetscObject)eps),PETSC_ERR_LIB,"Error in TRLAN (code=%" PetscBLASInt_FMT ")",stat);
     133           7 :   PetscFunctionReturn(PETSC_SUCCESS);
     134             : }
     135             : 
     136           5 : static PetscErrorCode EPSReset_TRLAN(EPS eps)
     137             : {
     138           5 :   EPS_TRLAN      *tr = (EPS_TRLAN*)eps->data;
     139             : 
     140           5 :   PetscFunctionBegin;
     141           5 :   PetscCall(PetscFree(tr->work));
     142           5 :   PetscFunctionReturn(PETSC_SUCCESS);
     143             : }
     144             : 
     145           5 : static PetscErrorCode EPSDestroy_TRLAN(EPS eps)
     146             : {
     147           5 :   PetscFunctionBegin;
     148           5 :   PetscCall(PetscFree(eps->data));
     149           5 :   PetscFunctionReturn(PETSC_SUCCESS);
     150             : }
     151             : 
     152           5 : SLEPC_EXTERN PetscErrorCode EPSCreate_TRLAN(EPS eps)
     153             : {
     154           5 :   EPS_TRLAN      *ctx;
     155             : 
     156           5 :   PetscFunctionBegin;
     157           5 :   PetscCall(PetscNew(&ctx));
     158           5 :   eps->data = (void*)ctx;
     159             : 
     160           5 :   eps->ops->solve          = EPSSolve_TRLAN;
     161           5 :   eps->ops->setup          = EPSSetUp_TRLAN;
     162           5 :   eps->ops->setupsort      = EPSSetUpSort_Basic;
     163           5 :   eps->ops->destroy        = EPSDestroy_TRLAN;
     164           5 :   eps->ops->reset          = EPSReset_TRLAN;
     165           5 :   eps->ops->backtransform  = EPSBackTransform_Default;
     166           5 :   PetscFunctionReturn(PETSC_SUCCESS);
     167             : }