LCOV - code coverage report
Current view: top level - nep/interface - nepsetup.c (source / functions) Hit Total Coverage
Test: SLEPc Lines: 142 179 79.3 %
Date: 2021-08-02 00:35:43 Functions: 4 4 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-2021, 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             :    NEP routines related to problem setup
      12             : */
      13             : 
      14             : #include <slepc/private/nepimpl.h>       /*I "slepcnep.h" I*/
      15             : 
      16             : /*@
      17             :    NEPSetUp - Sets up all the internal data structures necessary for the
      18             :    execution of the NEP solver.
      19             : 
      20             :    Collective on nep
      21             : 
      22             :    Input Parameter:
      23             : .  nep   - solver context
      24             : 
      25             :    Notes:
      26             :    This function need not be called explicitly in most cases, since NEPSolve()
      27             :    calls it. It can be useful when one wants to measure the set-up time
      28             :    separately from the solve time.
      29             : 
      30             :    Level: developer
      31             : 
      32             : .seealso: NEPCreate(), NEPSolve(), NEPDestroy()
      33             : @*/
      34         124 : PetscErrorCode NEPSetUp(NEP nep)
      35             : {
      36         124 :   PetscErrorCode ierr;
      37         124 :   PetscInt       k;
      38         124 :   SlepcSC        sc;
      39         124 :   Mat            T;
      40         124 :   PetscBool      flg;
      41         124 :   KSP            ksp;
      42         124 :   PC             pc;
      43         124 :   PetscMPIInt    size;
      44         124 :   MatSolverType  stype;
      45             : 
      46         124 :   PetscFunctionBegin;
      47         124 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
      48         124 :   NEPCheckProblem(nep,1);
      49         124 :   if (nep->state) PetscFunctionReturn(0);
      50         124 :   ierr = PetscLogEventBegin(NEP_SetUp,nep,0,0,0);CHKERRQ(ierr);
      51             : 
      52             :   /* reset the convergence flag from the previous solves */
      53         124 :   nep->reason = NEP_CONVERGED_ITERATING;
      54             : 
      55             :   /* set default solver type (NEPSetFromOptions was not called) */
      56         124 :   if (!((PetscObject)nep)->type_name) {
      57           0 :     ierr = NEPSetType(nep,NEPRII);CHKERRQ(ierr);
      58             :   }
      59         124 :   if (nep->useds && !nep->ds) { ierr = NEPGetDS(nep,&nep->ds);CHKERRQ(ierr); }
      60         124 :   if (!nep->rg) { ierr = NEPGetRG(nep,&nep->rg);CHKERRQ(ierr); }
      61         124 :   if (!((PetscObject)nep->rg)->type_name) {
      62           0 :     ierr = RGSetType(nep->rg,RGINTERVAL);CHKERRQ(ierr);
      63             :   }
      64             : 
      65             :   /* set problem dimensions */
      66         124 :   switch (nep->fui) {
      67          19 :   case NEP_USER_INTERFACE_CALLBACK:
      68          19 :     ierr = NEPGetFunction(nep,&T,NULL,NULL,NULL);CHKERRQ(ierr);
      69          19 :     ierr = MatGetSize(T,&nep->n,NULL);CHKERRQ(ierr);
      70          19 :     ierr = MatGetLocalSize(T,&nep->nloc,NULL);CHKERRQ(ierr);
      71             :     break;
      72         105 :   case NEP_USER_INTERFACE_SPLIT:
      73         105 :     ierr = MatDuplicate(nep->A[0],MAT_DO_NOT_COPY_VALUES,&nep->function);CHKERRQ(ierr);
      74         105 :     ierr = MatDuplicate(nep->A[0],MAT_DO_NOT_COPY_VALUES,&nep->jacobian);CHKERRQ(ierr);
      75         105 :     ierr = PetscLogObjectParent((PetscObject)nep,(PetscObject)nep->function);CHKERRQ(ierr);
      76         105 :     ierr = PetscLogObjectParent((PetscObject)nep,(PetscObject)nep->jacobian);CHKERRQ(ierr);
      77         105 :     ierr = MatGetSize(nep->A[0],&nep->n,NULL);CHKERRQ(ierr);
      78         105 :     ierr = MatGetLocalSize(nep->A[0],&nep->nloc,NULL);CHKERRQ(ierr);
      79             :     break;
      80             :   }
      81             : 
      82             :   /* set default problem type */
      83         124 :   if (!nep->problem_type) {
      84          85 :     ierr = NEPSetProblemType(nep,NEP_GENERAL);CHKERRQ(ierr);
      85             :   }
      86             : 
      87             :   /* check consistency of refinement options */
      88         124 :   if (nep->refine) {
      89          13 :     if (nep->fui!=NEP_USER_INTERFACE_SPLIT) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Iterative refinement only implemented in split form");
      90          13 :     if (!nep->scheme) {  /* set default scheme */
      91           1 :       ierr = NEPRefineGetKSP(nep,&ksp);CHKERRQ(ierr);
      92           1 :       ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
      93           1 :       ierr = PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&flg);CHKERRQ(ierr);
      94           1 :       if (flg) {
      95           0 :         ierr = PetscObjectTypeCompareAny((PetscObject)pc,&flg,PCLU,PCCHOLESKY,"");CHKERRQ(ierr);
      96             :       }
      97           2 :       nep->scheme = flg? NEP_REFINE_SCHEME_MBE: NEP_REFINE_SCHEME_SCHUR;
      98             :     }
      99          13 :     if (nep->scheme==NEP_REFINE_SCHEME_MBE) {
     100           3 :       ierr = NEPRefineGetKSP(nep,&ksp);CHKERRQ(ierr);
     101           3 :       ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
     102           3 :       ierr = PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&flg);CHKERRQ(ierr);
     103           3 :       if (flg) {
     104           3 :         ierr = PetscObjectTypeCompareAny((PetscObject)pc,&flg,PCLU,PCCHOLESKY,"");CHKERRQ(ierr);
     105             :       }
     106           3 :       if (!flg) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"The MBE scheme for refinement requires a direct solver in KSP");
     107           3 :       ierr = MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);CHKERRMPI(ierr);
     108           3 :       if (size>1) {   /* currently selected PC is a factorization */
     109           0 :         ierr = PCFactorGetMatSolverType(pc,&stype);CHKERRQ(ierr);
     110           0 :         ierr = PetscStrcmp(stype,MATSOLVERPETSC,&flg);CHKERRQ(ierr);
     111           0 :         if (flg) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"For Newton refinement, you chose to solve linear systems with a factorization, but in parallel runs you need to select an external package");
     112             :       }
     113             :     }
     114          13 :     if (nep->scheme==NEP_REFINE_SCHEME_SCHUR) {
     115           2 :       if (nep->npart>1) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"The Schur scheme for refinement does not support subcommunicators");
     116             :     }
     117             :   }
     118             :   /* call specific solver setup */
     119         124 :   ierr = (*nep->ops->setup)(nep);CHKERRQ(ierr);
     120             : 
     121             :   /* set tolerance if not yet set */
     122         124 :   if (nep->tol==PETSC_DEFAULT) nep->tol = SLEPC_DEFAULT_TOL;
     123         124 :   if (nep->refine) {
     124          25 :     if (nep->rtol==PETSC_DEFAULT) nep->rtol = PetscMax(nep->tol/1000,PETSC_MACHINE_EPSILON);
     125          13 :     if (nep->rits==PETSC_DEFAULT) nep->rits = (nep->refine==NEP_REFINE_SIMPLE)? 10: 1;
     126             :   }
     127             : 
     128             :   /* fill sorting criterion context */
     129         124 :   switch (nep->which) {
     130           0 :     case NEP_LARGEST_MAGNITUDE:
     131           0 :       nep->sc->comparison    = SlepcCompareLargestMagnitude;
     132           0 :       nep->sc->comparisonctx = NULL;
     133           0 :       break;
     134           0 :     case NEP_SMALLEST_MAGNITUDE:
     135           0 :       nep->sc->comparison    = SlepcCompareSmallestMagnitude;
     136           0 :       nep->sc->comparisonctx = NULL;
     137           0 :       break;
     138           0 :     case NEP_LARGEST_REAL:
     139           0 :       nep->sc->comparison    = SlepcCompareLargestReal;
     140           0 :       nep->sc->comparisonctx = NULL;
     141           0 :       break;
     142           0 :     case NEP_SMALLEST_REAL:
     143           0 :       nep->sc->comparison    = SlepcCompareSmallestReal;
     144           0 :       nep->sc->comparisonctx = NULL;
     145           0 :       break;
     146           0 :     case NEP_LARGEST_IMAGINARY:
     147           0 :       nep->sc->comparison    = SlepcCompareLargestImaginary;
     148           0 :       nep->sc->comparisonctx = NULL;
     149           0 :       break;
     150           0 :     case NEP_SMALLEST_IMAGINARY:
     151           0 :       nep->sc->comparison    = SlepcCompareSmallestImaginary;
     152           0 :       nep->sc->comparisonctx = NULL;
     153           0 :       break;
     154          98 :     case NEP_TARGET_MAGNITUDE:
     155          98 :       nep->sc->comparison    = SlepcCompareTargetMagnitude;
     156          98 :       nep->sc->comparisonctx = &nep->target;
     157          98 :       break;
     158           1 :     case NEP_TARGET_REAL:
     159           1 :       nep->sc->comparison    = SlepcCompareTargetReal;
     160           1 :       nep->sc->comparisonctx = &nep->target;
     161           1 :       break;
     162           0 :     case NEP_TARGET_IMAGINARY:
     163             : #if defined(PETSC_USE_COMPLEX)
     164           0 :       nep->sc->comparison    = SlepcCompareTargetImaginary;
     165           0 :       nep->sc->comparisonctx = &nep->target;
     166             : #endif
     167           0 :       break;
     168          25 :     case NEP_ALL:
     169          25 :       nep->sc->comparison    = SlepcCompareSmallestReal;
     170          25 :       nep->sc->comparisonctx = NULL;
     171          25 :       break;
     172             :     case NEP_WHICH_USER:
     173             :       break;
     174             :   }
     175             : 
     176         124 :   nep->sc->map    = NULL;
     177         124 :   nep->sc->mapobj = NULL;
     178             : 
     179             :   /* fill sorting criterion for DS */
     180         124 :   if (nep->useds) {
     181          97 :     ierr = DSGetSlepcSC(nep->ds,&sc);CHKERRQ(ierr);
     182          97 :     sc->comparison    = nep->sc->comparison;
     183          97 :     sc->comparisonctx = nep->sc->comparisonctx;
     184          97 :     ierr = PetscObjectTypeCompare((PetscObject)nep,NEPNLEIGS,&flg);CHKERRQ(ierr);
     185          97 :     if (!flg) {
     186          81 :       sc->map    = NULL;
     187          81 :       sc->mapobj = NULL;
     188             :     }
     189             :   }
     190         124 :   if (nep->nev > nep->ncv) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"nev bigger than ncv");
     191             : 
     192             :   /* process initial vectors */
     193         124 :   if (nep->nini<0) {
     194           4 :     k = -nep->nini;
     195           4 :     if (k>nep->ncv) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_USER_INPUT,"The number of initial vectors is larger than ncv");
     196           4 :     ierr = BVInsertVecs(nep->V,0,&k,nep->IS,PETSC_TRUE);CHKERRQ(ierr);
     197           4 :     ierr = SlepcBasisDestroy_Private(&nep->nini,&nep->IS);CHKERRQ(ierr);
     198           4 :     nep->nini = k;
     199             :   }
     200         124 :   ierr = PetscLogEventEnd(NEP_SetUp,nep,0,0,0);CHKERRQ(ierr);
     201         124 :   nep->state = NEP_STATE_SETUP;
     202         124 :   PetscFunctionReturn(0);
     203             : }
     204             : 
     205             : /*@C
     206             :    NEPSetInitialSpace - Specify a basis of vectors that constitute the initial
     207             :    space, that is, the subspace from which the solver starts to iterate.
     208             : 
     209             :    Collective on nep
     210             : 
     211             :    Input Parameters:
     212             : +  nep   - the nonlinear eigensolver context
     213             : .  n     - number of vectors
     214             : -  is    - set of basis vectors of the initial space
     215             : 
     216             :    Notes:
     217             :    Some solvers start to iterate on a single vector (initial vector). In that case,
     218             :    the other vectors are ignored.
     219             : 
     220             :    These vectors do not persist from one NEPSolve() call to the other, so the
     221             :    initial space should be set every time.
     222             : 
     223             :    The vectors do not need to be mutually orthonormal, since they are explicitly
     224             :    orthonormalized internally.
     225             : 
     226             :    Common usage of this function is when the user can provide a rough approximation
     227             :    of the wanted eigenspace. Then, convergence may be faster.
     228             : 
     229             :    Level: intermediate
     230             : @*/
     231           6 : PetscErrorCode NEPSetInitialSpace(NEP nep,PetscInt n,Vec is[])
     232             : {
     233           6 :   PetscErrorCode ierr;
     234             : 
     235           6 :   PetscFunctionBegin;
     236           6 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
     237          12 :   PetscValidLogicalCollectiveInt(nep,n,2);
     238           6 :   if (n<0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Argument n cannot be negative");
     239           6 :   if (n>0) {
     240           6 :     PetscValidPointer(is,3);
     241           6 :     PetscValidHeaderSpecific(*is,VEC_CLASSID,3);
     242             :   }
     243           6 :   ierr = SlepcBasisReference_Private(n,is,&nep->nini,&nep->IS);CHKERRQ(ierr);
     244           6 :   if (n>0) nep->state = NEP_STATE_INITIAL;
     245           6 :   PetscFunctionReturn(0);
     246             : }
     247             : 
     248             : /*
     249             :   NEPSetDimensions_Default - Set reasonable values for ncv, mpd if not set
     250             :   by the user. This is called at setup.
     251             :  */
     252          61 : PetscErrorCode NEPSetDimensions_Default(NEP nep,PetscInt nev,PetscInt *ncv,PetscInt *mpd)
     253             : {
     254          61 :   PetscFunctionBegin;
     255          61 :   if (*ncv!=PETSC_DEFAULT) { /* ncv set */
     256          15 :     if (*ncv<nev) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_USER_INPUT,"The value of ncv must be at least nev");
     257          46 :   } else if (*mpd!=PETSC_DEFAULT) { /* mpd set */
     258           0 :     *ncv = PetscMin(nep->n,nev+(*mpd));
     259             :   } else { /* neither set: defaults depend on nev being small or large */
     260          46 :     if (nev<500) *ncv = PetscMin(nep->n,PetscMax(2*nev,nev+15));
     261             :     else {
     262           0 :       *mpd = 500;
     263           0 :       *ncv = PetscMin(nep->n,nev+(*mpd));
     264             :     }
     265             :   }
     266          61 :   if (*mpd==PETSC_DEFAULT) *mpd = *ncv;
     267          61 :   PetscFunctionReturn(0);
     268             : }
     269             : 
     270             : /*@
     271             :    NEPAllocateSolution - Allocate memory storage for common variables such
     272             :    as eigenvalues and eigenvectors.
     273             : 
     274             :    Collective on nep
     275             : 
     276             :    Input Parameters:
     277             : +  nep   - eigensolver context
     278             : -  extra - number of additional positions, used for methods that require a
     279             :            working basis slightly larger than ncv
     280             : 
     281             :    Developers Note:
     282             :    This is SLEPC_EXTERN because it may be required by user plugin NEP
     283             :    implementations.
     284             : 
     285             :    Level: developer
     286             : @*/
     287         127 : PetscErrorCode NEPAllocateSolution(NEP nep,PetscInt extra)
     288             : {
     289         127 :   PetscErrorCode ierr;
     290         127 :   PetscInt       oldsize,newc,requested;
     291         127 :   PetscLogDouble cnt;
     292         127 :   PetscRandom    rand;
     293         127 :   Mat            T;
     294         127 :   Vec            t;
     295             : 
     296         127 :   PetscFunctionBegin;
     297         127 :   requested = nep->ncv + extra;
     298             : 
     299             :   /* oldsize is zero if this is the first time setup is called */
     300         127 :   ierr = BVGetSizes(nep->V,NULL,NULL,&oldsize);CHKERRQ(ierr);
     301         127 :   newc = PetscMax(0,requested-oldsize);
     302             : 
     303             :   /* allocate space for eigenvalues and friends */
     304         127 :   if (requested != oldsize || !nep->eigr) {
     305         127 :     ierr = PetscFree4(nep->eigr,nep->eigi,nep->errest,nep->perm);CHKERRQ(ierr);
     306         127 :     ierr = PetscMalloc4(requested,&nep->eigr,requested,&nep->eigi,requested,&nep->errest,requested,&nep->perm);CHKERRQ(ierr);
     307         127 :     cnt = newc*sizeof(PetscScalar) + newc*sizeof(PetscReal) + newc*sizeof(PetscInt);
     308         127 :     ierr = PetscLogObjectMemory((PetscObject)nep,cnt);CHKERRQ(ierr);
     309             :   }
     310             : 
     311             :   /* allocate V */
     312         127 :   if (!nep->V) { ierr = NEPGetBV(nep,&nep->V);CHKERRQ(ierr); }
     313         127 :   if (!oldsize) {
     314         124 :     if (!((PetscObject)(nep->V))->type_name) {
     315          14 :       ierr = BVSetType(nep->V,BVSVEC);CHKERRQ(ierr);
     316             :     }
     317         124 :     if (nep->fui==NEP_USER_INTERFACE_SPLIT) T = nep->A[0];
     318             :     else {
     319          19 :       ierr = NEPGetFunction(nep,&T,NULL,NULL,NULL);CHKERRQ(ierr);
     320             :     }
     321         124 :     ierr = MatCreateVecsEmpty(T,&t,NULL);CHKERRQ(ierr);
     322         124 :     ierr = BVSetSizesFromVec(nep->V,t,requested);CHKERRQ(ierr);
     323         124 :     ierr = VecDestroy(&t);CHKERRQ(ierr);
     324             :   } else {
     325           3 :     ierr = BVResize(nep->V,requested,PETSC_FALSE);CHKERRQ(ierr);
     326             :   }
     327             : 
     328             :   /* allocate W */
     329         127 :   if (nep->twosided) {
     330           7 :     ierr = BVGetRandomContext(nep->V,&rand);CHKERRQ(ierr);  /* make sure the random context is available when duplicating */
     331           7 :     ierr = BVDestroy(&nep->W);CHKERRQ(ierr);
     332           7 :     ierr = BVDuplicate(nep->V,&nep->W);CHKERRQ(ierr);
     333             :   }
     334         127 :   PetscFunctionReturn(0);
     335             : }
     336             : 

Generated by: LCOV version 1.14