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

Generated by: LCOV version 1.13