LCOV - code coverage report
Current view: top level - eps/impls/external/trlan - trlan.c (source / functions) Hit Total Coverage
Test: coverage.info Lines: 96 98 98.0 %
Date: 2020-05-28 03:14:05 Functions: 6 6 100.0 %

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

Generated by: LCOV version 1.13