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