Actual source code: nepbasic.c

      1: /*
      2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
      4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

      6:    This file is part of SLEPc.
      7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
      8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      9: */
     10: /*
     11:    Basic NEP routines
     12: */

     14: #include <slepc/private/nepimpl.h>

     16: /* Logging support */
     17: PetscClassId      NEP_CLASSID = 0;
     18: PetscLogEvent     NEP_SetUp = 0,NEP_Solve = 0,NEP_Refine = 0,NEP_FunctionEval = 0,NEP_JacobianEval = 0,NEP_Resolvent = 0,NEP_CISS_SVD = 0;

     20: /* List of registered NEP routines */
     21: PetscFunctionList NEPList = NULL;
     22: PetscBool         NEPRegisterAllCalled = PETSC_FALSE;

     24: /* List of registered NEP monitors */
     25: PetscFunctionList NEPMonitorList              = NULL;
     26: PetscFunctionList NEPMonitorCreateList        = NULL;
     27: PetscFunctionList NEPMonitorDestroyList       = NULL;
     28: PetscBool         NEPMonitorRegisterAllCalled = PETSC_FALSE;

     30: /*@
     31:    NEPCreate - Creates the `NEP` context.

     33:    Collective

     35:    Input Parameter:
     36: .  comm - MPI communicator

     38:    Output Parameter:
     39: .  outnep - location to put the `NEP` context

     41:    Note:
     42:    The default `NEP` type is `NEPRII`.

     44:    Level: beginner

     46: .seealso: [](ch:nep), `NEPSetUp()`, `NEPSolve()`, `NEPDestroy()`, `NEP`
     47: @*/
     48: PetscErrorCode NEPCreate(MPI_Comm comm,NEP *outnep)
     49: {
     50:   NEP            nep;

     52:   PetscFunctionBegin;
     53:   PetscAssertPointer(outnep,2);
     54:   PetscCall(NEPInitializePackage());
     55:   PetscCall(SlepcHeaderCreate(nep,NEP_CLASSID,"NEP","Nonlinear Eigenvalue Problem","NEP",comm,NEPDestroy,NEPView));

     57:   nep->max_it          = PETSC_DETERMINE;
     58:   nep->nev             = 1;
     59:   nep->ncv             = PETSC_DETERMINE;
     60:   nep->mpd             = PETSC_DETERMINE;
     61:   nep->nini            = 0;
     62:   nep->target          = 0.0;
     63:   nep->tol             = PETSC_DETERMINE;
     64:   nep->conv            = NEP_CONV_REL;
     65:   nep->stop            = NEP_STOP_BASIC;
     66:   nep->which           = (NEPWhich)0;
     67:   nep->problem_type    = (NEPProblemType)0;
     68:   nep->refine          = NEP_REFINE_NONE;
     69:   nep->npart           = 1;
     70:   nep->rtol            = PETSC_DETERMINE;
     71:   nep->rits            = PETSC_DETERMINE;
     72:   nep->scheme          = (NEPRefineScheme)0;
     73:   nep->trackall        = PETSC_FALSE;
     74:   nep->twosided        = PETSC_FALSE;

     76:   nep->computefunction = NULL;
     77:   nep->computejacobian = NULL;
     78:   nep->functionctx     = NULL;
     79:   nep->jacobianctx     = NULL;
     80:   nep->converged       = NEPConvergedRelative;
     81:   nep->convergeduser   = NULL;
     82:   nep->convergeddestroy= NULL;
     83:   nep->stopping        = NEPStoppingBasic;
     84:   nep->stoppinguser    = NULL;
     85:   nep->stoppingdestroy = NULL;
     86:   nep->convergedctx    = NULL;
     87:   nep->stoppingctx     = NULL;
     88:   nep->numbermonitors  = 0;

     90:   nep->ds              = NULL;
     91:   nep->V               = NULL;
     92:   nep->W               = NULL;
     93:   nep->rg              = NULL;
     94:   nep->function        = NULL;
     95:   nep->function_pre    = NULL;
     96:   nep->jacobian        = NULL;
     97:   nep->A               = NULL;
     98:   nep->f               = NULL;
     99:   nep->nt              = 0;
    100:   nep->mstr            = UNKNOWN_NONZERO_PATTERN;
    101:   nep->P               = NULL;
    102:   nep->mstrp           = UNKNOWN_NONZERO_PATTERN;
    103:   nep->IS              = NULL;
    104:   nep->eigr            = NULL;
    105:   nep->eigi            = NULL;
    106:   nep->errest          = NULL;
    107:   nep->perm            = NULL;
    108:   nep->nwork           = 0;
    109:   nep->work            = NULL;
    110:   nep->data            = NULL;

    112:   nep->state           = NEP_STATE_INITIAL;
    113:   nep->nconv           = 0;
    114:   nep->its             = 0;
    115:   nep->n               = 0;
    116:   nep->nloc            = 0;
    117:   nep->nrma            = NULL;
    118:   nep->fui             = (NEPUserInterface)0;
    119:   nep->useds           = PETSC_FALSE;
    120:   nep->resolvent       = NULL;
    121:   nep->reason          = NEP_CONVERGED_ITERATING;

    123:   PetscCall(PetscNew(&nep->sc));
    124:   *outnep = nep;
    125:   PetscFunctionReturn(PETSC_SUCCESS);
    126: }

    128: /*@
    129:    NEPSetType - Selects the particular solver to be used in the `NEP` object.

    131:    Logically Collective

    133:    Input Parameters:
    134: +  nep  - the nonlinear eigensolver context
    135: -  type - a known method

    137:    Options Database Key:
    138: .  -nep_type type - sets the method; use `-help` for a list of available methods

    140:    Notes:
    141:    See `NEPType` for available methods. The default is `NEPRII`.

    143:    Normally, it is best to use the `NEPSetFromOptions()` command and
    144:    then set the `NEP` type from the options database rather than by using
    145:    this routine.  Using the options database provides the user with
    146:    maximum flexibility in evaluating the different available methods.
    147:    The `NEPSetType()` routine is provided for those situations where it
    148:    is necessary to set the iterative solver independently of the command
    149:    line or options database.

    151:    Level: intermediate

    153: .seealso: [](ch:nep), `NEPType`
    154: @*/
    155: PetscErrorCode NEPSetType(NEP nep,NEPType type)
    156: {
    157:   PetscErrorCode (*r)(NEP);
    158:   PetscBool      match;

    160:   PetscFunctionBegin;
    162:   PetscAssertPointer(type,2);

    164:   PetscCall(PetscObjectTypeCompare((PetscObject)nep,type,&match));
    165:   if (match) PetscFunctionReturn(PETSC_SUCCESS);

    167:   PetscCall(PetscFunctionListFind(NEPList,type,&r));
    168:   PetscCheck(r,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown NEP type given: %s",type);

    170:   PetscTryTypeMethod(nep,destroy);
    171:   PetscCall(PetscMemzero(nep->ops,sizeof(struct _NEPOps)));

    173:   nep->state = NEP_STATE_INITIAL;
    174:   PetscCall(PetscObjectChangeTypeName((PetscObject)nep,type));
    175:   PetscCall((*r)(nep));
    176:   PetscFunctionReturn(PETSC_SUCCESS);
    177: }

    179: /*@
    180:    NEPGetType - Gets the `NEP` type as a string from the `NEP` object.

    182:    Not Collective

    184:    Input Parameter:
    185: .  nep - the nonlinear eigensolver context

    187:    Output Parameter:
    188: .  type - name of `NEP` method

    190:    Level: intermediate

    192: .seealso: [](ch:nep), `NEPSetType()`
    193: @*/
    194: PetscErrorCode NEPGetType(NEP nep,NEPType *type)
    195: {
    196:   PetscFunctionBegin;
    198:   PetscAssertPointer(type,2);
    199:   *type = ((PetscObject)nep)->type_name;
    200:   PetscFunctionReturn(PETSC_SUCCESS);
    201: }

    203: /*@C
    204:    NEPRegister - Adds a method to the nonlinear eigenproblem solver package.

    206:    Not Collective

    208:    Input Parameters:
    209: +  name - name of a new user-defined solver
    210: -  function - routine to create the solver context

    212:    Note:
    213:    `NEPRegister()` may be called multiple times to add several user-defined solvers.

    215:    Example Usage:
    216: .vb
    217:    NEPRegister("my_solver",MySolverCreate);
    218: .ve

    220:    Then, your solver can be chosen with the procedural interface via
    221: .vb
    222:    NEPSetType(nep,"my_solver")
    223: .ve
    224:    or at runtime via the option `-nep_type my_solver`.

    226:    Level: advanced

    228: .seealso: [](ch:nep), `NEPRegisterAll()`
    229: @*/
    230: PetscErrorCode NEPRegister(const char *name,PetscErrorCode (*function)(NEP))
    231: {
    232:   PetscFunctionBegin;
    233:   PetscCall(NEPInitializePackage());
    234:   PetscCall(PetscFunctionListAdd(&NEPList,name,function));
    235:   PetscFunctionReturn(PETSC_SUCCESS);
    236: }

    238: /*@C
    239:    NEPMonitorRegister - Registers a `NEP` monitor routine that may be accessed with
    240:    `NEPMonitorSetFromOptions()`.

    242:    Not Collective

    244:    Input Parameters:
    245: +  name    - name of a new monitor routine
    246: .  vtype   - a `PetscViewerType` for the output
    247: .  format  - a `PetscViewerFormat` for the output
    248: .  monitor - monitor routine, see `NEPMonitorRegisterFn`
    249: .  create  - creation routine, or `NULL`
    250: -  destroy - destruction routine, or `NULL`

    252:    Notes:
    253:    `NEPMonitorRegister()` may be called multiple times to add several user-defined monitors.

    255:    The calling sequence for the given function matches the calling sequence of `NEPMonitorFn`
    256:    functions passed to `NEPMonitorSet()` with the additional requirement that its final argument
    257:    be a `PetscViewerAndFormat`.

    259:    Example Usage:
    260: .vb
    261:    NEPMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
    262: .ve

    264:    Then, your monitor can be chosen with the procedural interface via
    265: .vb
    266:    NEPMonitorSetFromOptions(nep,"-nep_monitor_my_monitor","my_monitor",NULL);
    267: .ve
    268:    or at runtime via the option `-nep_monitor_my_monitor`.

    270:    Level: advanced

    272: .seealso: [](ch:nep), `NEPMonitorSet()`, `NEPMonitorRegisterAll()`, `NEPMonitorSetFromOptions()`
    273: @*/
    274: PetscErrorCode NEPMonitorRegister(const char name[],PetscViewerType vtype,PetscViewerFormat format,NEPMonitorRegisterFn *monitor,NEPMonitorRegisterCreateFn *create,NEPMonitorRegisterDestroyFn *destroy)
    275: {
    276:   char           key[PETSC_MAX_PATH_LEN];

    278:   PetscFunctionBegin;
    279:   PetscCall(NEPInitializePackage());
    280:   PetscCall(SlepcMonitorMakeKey_Internal(name,vtype,format,key));
    281:   PetscCall(PetscFunctionListAdd(&NEPMonitorList,key,monitor));
    282:   if (create)  PetscCall(PetscFunctionListAdd(&NEPMonitorCreateList,key,create));
    283:   if (destroy) PetscCall(PetscFunctionListAdd(&NEPMonitorDestroyList,key,destroy));
    284:   PetscFunctionReturn(PETSC_SUCCESS);
    285: }

    287: /*
    288:    NEPReset_Problem - Destroys the problem matrices.
    289: */
    290: PetscErrorCode NEPReset_Problem(NEP nep)
    291: {
    292:   PetscInt       i;

    294:   PetscFunctionBegin;
    296:   PetscCall(MatDestroy(&nep->function));
    297:   PetscCall(MatDestroy(&nep->function_pre));
    298:   PetscCall(MatDestroy(&nep->jacobian));
    299:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
    300:     PetscCall(MatDestroyMatrices(nep->nt,&nep->A));
    301:     for (i=0;i<nep->nt;i++) PetscCall(FNDestroy(&nep->f[i]));
    302:     PetscCall(PetscFree(nep->f));
    303:     PetscCall(PetscFree(nep->nrma));
    304:     if (nep->P) PetscCall(MatDestroyMatrices(nep->nt,&nep->P));
    305:     nep->nt = 0;
    306:   }
    307:   PetscFunctionReturn(PETSC_SUCCESS);
    308: }
    309: /*@
    310:    NEPReset - Resets the `NEP` context to the initial state (prior to setup)
    311:    and destroys any allocated `Vec`s and `Mat`s.

    313:    Collective

    315:    Input Parameter:
    316: .  nep - the nonlinear eigensolver context

    318:    Level: advanced

    320: .seealso: [](ch:nep), `NEPDestroy()`
    321: @*/
    322: PetscErrorCode NEPReset(NEP nep)
    323: {
    324:   PetscFunctionBegin;
    326:   if (!nep) PetscFunctionReturn(PETSC_SUCCESS);
    327:   PetscTryTypeMethod(nep,reset);
    328:   if (nep->refineksp) PetscCall(KSPReset(nep->refineksp));
    329:   PetscCall(NEPReset_Problem(nep));
    330:   PetscCall(BVDestroy(&nep->V));
    331:   PetscCall(BVDestroy(&nep->W));
    332:   PetscCall(VecDestroyVecs(nep->nwork,&nep->work));
    333:   PetscCall(MatDestroy(&nep->resolvent));
    334:   nep->nwork = 0;
    335:   nep->state = NEP_STATE_INITIAL;
    336:   PetscFunctionReturn(PETSC_SUCCESS);
    337: }

    339: /*@
    340:    NEPDestroy - Destroys the `NEP` context.

    342:    Collective

    344:    Input Parameter:
    345: .  nep - the nonlinear eigensolver context

    347:    Level: beginner

    349: .seealso: [](ch:nep), `NEPCreate()`, `NEPSetUp()`, `NEPSolve()`
    350: @*/
    351: PetscErrorCode NEPDestroy(NEP *nep)
    352: {
    353:   PetscFunctionBegin;
    354:   if (!*nep) PetscFunctionReturn(PETSC_SUCCESS);
    356:   if (--((PetscObject)*nep)->refct > 0) { *nep = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
    357:   PetscCall(NEPReset(*nep));
    358:   PetscTryTypeMethod(*nep,destroy);
    359:   if ((*nep)->eigr) PetscCall(PetscFree4((*nep)->eigr,(*nep)->eigi,(*nep)->errest,(*nep)->perm));
    360:   PetscCall(RGDestroy(&(*nep)->rg));
    361:   PetscCall(DSDestroy(&(*nep)->ds));
    362:   PetscCall(KSPDestroy(&(*nep)->refineksp));
    363:   PetscCall(PetscSubcommDestroy(&(*nep)->refinesubc));
    364:   PetscCall(PetscFree((*nep)->sc));
    365:   /* just in case the initial vectors have not been used */
    366:   PetscCall(SlepcBasisDestroy_Private(&(*nep)->nini,&(*nep)->IS));
    367:   if ((*nep)->convergeddestroy) PetscCall((*(*nep)->convergeddestroy)(&(*nep)->convergedctx));
    368:   if ((*nep)->stoppingdestroy) PetscCall((*(*nep)->stoppingdestroy)(&(*nep)->stoppingctx));
    369:   PetscCall(NEPMonitorCancel(*nep));
    370:   PetscCall(PetscHeaderDestroy(nep));
    371:   PetscFunctionReturn(PETSC_SUCCESS);
    372: }

    374: /*@
    375:    NEPSetBV - Associates a basis vectors object to the nonlinear eigensolver.

    377:    Collective

    379:    Input Parameters:
    380: +  nep - the nonlinear eigensolver context
    381: -  bv  - the basis vectors object

    383:    Note:
    384:    Use `NEPGetBV()` to retrieve the basis vectors context (for example,
    385:    to free it at the end of the computations).

    387:    Level: advanced

    389: .seealso: [](ch:nep), `NEPGetBV()`
    390: @*/
    391: PetscErrorCode NEPSetBV(NEP nep,BV bv)
    392: {
    393:   PetscFunctionBegin;
    396:   PetscCheckSameComm(nep,1,bv,2);
    397:   PetscCall(PetscObjectReference((PetscObject)bv));
    398:   PetscCall(BVDestroy(&nep->V));
    399:   nep->V = bv;
    400:   PetscFunctionReturn(PETSC_SUCCESS);
    401: }

    403: /*@
    404:    NEPGetBV - Obtain the basis vectors object associated to the nonlinear
    405:    eigensolver object.

    407:    Not Collective

    409:    Input Parameter:
    410: .  nep - the nonlinear eigensolver context

    412:    Output Parameter:
    413: .  bv - basis vectors context

    415:    Level: advanced

    417: .seealso: [](ch:nep), `NEPSetBV()`
    418: @*/
    419: PetscErrorCode NEPGetBV(NEP nep,BV *bv)
    420: {
    421:   PetscFunctionBegin;
    423:   PetscAssertPointer(bv,2);
    424:   if (!nep->V) {
    425:     PetscCall(BVCreate(PetscObjectComm((PetscObject)nep),&nep->V));
    426:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)nep->V,(PetscObject)nep,0));
    427:     PetscCall(PetscObjectSetOptions((PetscObject)nep->V,((PetscObject)nep)->options));
    428:   }
    429:   *bv = nep->V;
    430:   PetscFunctionReturn(PETSC_SUCCESS);
    431: }

    433: /*@
    434:    NEPSetRG - Associates a region object to the nonlinear eigensolver.

    436:    Collective

    438:    Input Parameters:
    439: +  nep - the nonlinear eigensolver context
    440: -  rg  - the region object

    442:    Note:
    443:    Use `NEPGetRG()` to retrieve the region context (for example,
    444:    to free it at the end of the computations).

    446:    Level: advanced

    448: .seealso: [](ch:nep), `NEPGetRG()`
    449: @*/
    450: PetscErrorCode NEPSetRG(NEP nep,RG rg)
    451: {
    452:   PetscFunctionBegin;
    454:   if (rg) {
    456:     PetscCheckSameComm(nep,1,rg,2);
    457:   }
    458:   PetscCall(PetscObjectReference((PetscObject)rg));
    459:   PetscCall(RGDestroy(&nep->rg));
    460:   nep->rg = rg;
    461:   PetscFunctionReturn(PETSC_SUCCESS);
    462: }

    464: /*@
    465:    NEPGetRG - Obtain the region object associated to the
    466:    nonlinear eigensolver object.

    468:    Not Collective

    470:    Input Parameter:
    471: .  nep - the nonlinear eigensolver context

    473:    Output Parameter:
    474: .  rg - region context

    476:    Level: advanced

    478: .seealso: [](ch:nep), `NEPSetRG()`
    479: @*/
    480: PetscErrorCode NEPGetRG(NEP nep,RG *rg)
    481: {
    482:   PetscFunctionBegin;
    484:   PetscAssertPointer(rg,2);
    485:   if (!nep->rg) {
    486:     PetscCall(RGCreate(PetscObjectComm((PetscObject)nep),&nep->rg));
    487:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)nep->rg,(PetscObject)nep,0));
    488:     PetscCall(PetscObjectSetOptions((PetscObject)nep->rg,((PetscObject)nep)->options));
    489:   }
    490:   *rg = nep->rg;
    491:   PetscFunctionReturn(PETSC_SUCCESS);
    492: }

    494: /*@
    495:    NEPSetDS - Associates a direct solver object to the nonlinear eigensolver.

    497:    Collective

    499:    Input Parameters:
    500: +  nep - the nonlinear eigensolver context
    501: -  ds  - the direct solver object

    503:    Note:
    504:    Use `NEPGetDS()` to retrieve the direct solver context (for example,
    505:    to free it at the end of the computations).

    507:    Level: advanced

    509: .seealso: [](ch:nep), `NEPGetDS()`
    510: @*/
    511: PetscErrorCode NEPSetDS(NEP nep,DS ds)
    512: {
    513:   PetscFunctionBegin;
    516:   PetscCheckSameComm(nep,1,ds,2);
    517:   PetscCall(PetscObjectReference((PetscObject)ds));
    518:   PetscCall(DSDestroy(&nep->ds));
    519:   nep->ds = ds;
    520:   PetscFunctionReturn(PETSC_SUCCESS);
    521: }

    523: /*@
    524:    NEPGetDS - Obtain the direct solver object associated to the
    525:    nonlinear eigensolver object.

    527:    Not Collective

    529:    Input Parameter:
    530: .  nep - the nonlinear eigensolver context

    532:    Output Parameter:
    533: .  ds - direct solver context

    535:    Level: advanced

    537: .seealso: [](ch:nep), `NEPSetDS()`
    538: @*/
    539: PetscErrorCode NEPGetDS(NEP nep,DS *ds)
    540: {
    541:   PetscFunctionBegin;
    543:   PetscAssertPointer(ds,2);
    544:   if (!nep->ds) {
    545:     PetscCall(DSCreate(PetscObjectComm((PetscObject)nep),&nep->ds));
    546:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)nep->ds,(PetscObject)nep,0));
    547:     PetscCall(PetscObjectSetOptions((PetscObject)nep->ds,((PetscObject)nep)->options));
    548:   }
    549:   *ds = nep->ds;
    550:   PetscFunctionReturn(PETSC_SUCCESS);
    551: }

    553: /*@
    554:    NEPRefineGetKSP - Obtain the `KSP` object used by the eigensolver
    555:    object in the refinement phase.

    557:    Collective

    559:    Input Parameter:
    560: .  nep - the nonlinear eigensolver context

    562:    Output Parameter:
    563: .  ksp - the linear solver context

    565:    Level: advanced

    567: .seealso: [](ch:nep), `NEPSetRefine()`
    568: @*/
    569: PetscErrorCode NEPRefineGetKSP(NEP nep,KSP *ksp)
    570: {
    571:   MPI_Comm       comm;

    573:   PetscFunctionBegin;
    575:   PetscAssertPointer(ksp,2);
    576:   if (!nep->refineksp) {
    577:     if (nep->npart>1) {
    578:       /* Split in subcomunicators */
    579:       PetscCall(PetscSubcommCreate(PetscObjectComm((PetscObject)nep),&nep->refinesubc));
    580:       PetscCall(PetscSubcommSetNumber(nep->refinesubc,nep->npart));
    581:       PetscCall(PetscSubcommSetType(nep->refinesubc,PETSC_SUBCOMM_CONTIGUOUS));
    582:       PetscCall(PetscSubcommGetChild(nep->refinesubc,&comm));
    583:     } else PetscCall(PetscObjectGetComm((PetscObject)nep,&comm));
    584:     PetscCall(KSPCreate(comm,&nep->refineksp));
    585:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)nep->refineksp,(PetscObject)nep,0));
    586:     PetscCall(PetscObjectSetOptions((PetscObject)nep->refineksp,((PetscObject)nep)->options));
    587:     PetscCall(KSPSetOptionsPrefix(*ksp,((PetscObject)nep)->prefix));
    588:     PetscCall(KSPAppendOptionsPrefix(*ksp,"nep_refine_"));
    589:     PetscCall(KSPSetTolerances(nep->refineksp,SlepcDefaultTol(nep->rtol),PETSC_CURRENT,PETSC_CURRENT,PETSC_CURRENT));
    590:   }
    591:   *ksp = nep->refineksp;
    592:   PetscFunctionReturn(PETSC_SUCCESS);
    593: }

    595: /*@
    596:    NEPSetTarget - Sets the value of the target.

    598:    Logically Collective

    600:    Input Parameters:
    601: +  nep    - the nonlinear eigensolver context
    602: -  target - the value of the target

    604:    Options Database Key:
    605: .  -nep_target target - the value of the target

    607:    Notes:
    608:    The target is a scalar value used to determine the portion of the spectrum
    609:    of interest. It is used in combination with `NEPSetWhichEigenpairs()`.

    611:    When PETSc is built with real scalars, it is not possible to specify a
    612:    complex target.

    614:    In the case of complex scalars, a complex value can be provided in the
    615:    command line with `[+/-][realnumber][+/-]realnumberi` with no spaces, e.g.
    616:    `-nep_target 1.0+2.0i`.

    618:    Level: intermediate

    620: .seealso: [](ch:nep), `NEPGetTarget()`, `NEPSetWhichEigenpairs()`
    621: @*/
    622: PetscErrorCode NEPSetTarget(NEP nep,PetscScalar target)
    623: {
    624:   PetscFunctionBegin;
    627:   nep->target = target;
    628:   PetscFunctionReturn(PETSC_SUCCESS);
    629: }

    631: /*@
    632:    NEPGetTarget - Gets the value of the target.

    634:    Not Collective

    636:    Input Parameter:
    637: .  nep - the nonlinear eigensolver context

    639:    Output Parameter:
    640: .  target - the value of the target

    642:    Note:
    643:    If the target was not set by the user, then zero is returned.

    645:    Level: intermediate

    647: .seealso: [](ch:nep), `NEPSetTarget()`
    648: @*/
    649: PetscErrorCode NEPGetTarget(NEP nep,PetscScalar* target)
    650: {
    651:   PetscFunctionBegin;
    653:   PetscAssertPointer(target,2);
    654:   *target = nep->target;
    655:   PetscFunctionReturn(PETSC_SUCCESS);
    656: }

    658: /*@C
    659:    NEPSetFunction - Sets the function to compute the nonlinear Function $T(\lambda)$
    660:    as well as the location to store the matrix.

    662:    Collective

    664:    Input Parameters:
    665: +  nep - the nonlinear eigensolver context
    666: .  F   - Function matrix
    667: .  P   - preconditioner matrix (usually the same as `F`)
    668: .  fun - Function evaluation routine (if `NULL` then `NEP` retains any
    669:          previously set value), see `NEPFunctionFn` for the calling sequence
    670: -  ctx - [optional] user-defined context for private data for the Function
    671:          evaluation routine (may be `NULL`) (if `NULL` then `NEP` retains any
    672:          previously set value)

    674:    Level: beginner

    676: .seealso: [](ch:nep), `NEPGetFunction()`, `NEPSetJacobian()`
    677: @*/
    678: PetscErrorCode NEPSetFunction(NEP nep,Mat F,Mat P,NEPFunctionFn *fun,PetscCtx ctx)
    679: {
    680:   PetscFunctionBegin;
    684:   if (F) PetscCheckSameComm(nep,1,F,2);
    685:   if (P) PetscCheckSameComm(nep,1,P,3);

    687:   if (nep->state) PetscCall(NEPReset(nep));
    688:   else if (nep->fui && nep->fui!=NEP_USER_INTERFACE_CALLBACK) PetscCall(NEPReset_Problem(nep));

    690:   if (fun) nep->computefunction = fun;
    691:   if (ctx) nep->functionctx     = ctx;
    692:   if (F) {
    693:     PetscCall(PetscObjectReference((PetscObject)F));
    694:     PetscCall(MatDestroy(&nep->function));
    695:     nep->function = F;
    696:   }
    697:   if (P) {
    698:     PetscCall(PetscObjectReference((PetscObject)P));
    699:     PetscCall(MatDestroy(&nep->function_pre));
    700:     nep->function_pre = P;
    701:   }
    702:   nep->fui   = NEP_USER_INTERFACE_CALLBACK;
    703:   nep->state = NEP_STATE_INITIAL;
    704:   PetscFunctionReturn(PETSC_SUCCESS);
    705: }

    707: /*@C
    708:    NEPGetFunction - Returns the Function matrix and optionally the user
    709:    provided context for evaluating the Function.

    711:    Collective

    713:    Input Parameter:
    714: .  nep - the nonlinear eigensolver context

    716:    Output Parameters:
    717: +  F   - location to stash Function matrix (or `NULL`)
    718: .  P   - location to stash preconditioner matrix (or `NULL`)
    719: .  fun - location to put Function function (or `NULL`)
    720: -  ctx - location to stash Function context (or `NULL`)

    722:    Level: advanced

    724: .seealso: [](ch:nep), `NEPSetFunction()`
    725: @*/
    726: PetscErrorCode NEPGetFunction(NEP nep,Mat *F,Mat *P,NEPFunctionFn **fun,PetscCtxRt ctx)
    727: {
    728:   PetscFunctionBegin;
    730:   NEPCheckCallback(nep,1);
    731:   if (F)   *F   = nep->function;
    732:   if (P)   *P   = nep->function_pre;
    733:   if (fun) *fun = nep->computefunction;
    734:   if (ctx) *(void**)ctx = nep->functionctx;
    735:   PetscFunctionReturn(PETSC_SUCCESS);
    736: }

    738: /*@C
    739:    NEPSetJacobian - Sets the function to compute the Jacobian $T'(\lambda)$ as well
    740:    as the location to store the matrix.

    742:    Collective

    744:    Input Parameters:
    745: +  nep - the nonlinear eigensolver context
    746: .  J   - Jacobian matrix
    747: .  jac - Jacobian evaluation routine (if `NULL` then `NEP` retains any
    748:          previously set value), see `NEPJacobianFn` for the calling sequence
    749: -  ctx - [optional] user-defined context for private data for the Jacobian
    750:          evaluation routine (may be `NULL`) (if `NULL` then `NEP` retains any
    751:          previously set value)

    753:    Level: beginner

    755: .seealso: [](ch:nep), `NEPSetFunction()`, `NEPGetJacobian()`
    756: @*/
    757: PetscErrorCode NEPSetJacobian(NEP nep,Mat J,NEPJacobianFn *jac,PetscCtx ctx)
    758: {
    759:   PetscFunctionBegin;
    762:   if (J) PetscCheckSameComm(nep,1,J,2);

    764:   if (nep->state) PetscCall(NEPReset(nep));
    765:   else if (nep->fui && nep->fui!=NEP_USER_INTERFACE_CALLBACK) PetscCall(NEPReset_Problem(nep));

    767:   if (jac) nep->computejacobian = jac;
    768:   if (ctx) nep->jacobianctx     = ctx;
    769:   if (J) {
    770:     PetscCall(PetscObjectReference((PetscObject)J));
    771:     PetscCall(MatDestroy(&nep->jacobian));
    772:     nep->jacobian = J;
    773:   }
    774:   nep->fui   = NEP_USER_INTERFACE_CALLBACK;
    775:   nep->state = NEP_STATE_INITIAL;
    776:   PetscFunctionReturn(PETSC_SUCCESS);
    777: }

    779: /*@C
    780:    NEPGetJacobian - Returns the Jacobian matrix and optionally the user
    781:    provided routine and context for evaluating the Jacobian.

    783:    Collective

    785:    Input Parameter:
    786: .  nep - the nonlinear eigensolver context

    788:    Output Parameters:
    789: +  J   - location to stash Jacobian matrix (or `NULL`)
    790: .  jac - location to put Jacobian function (or `NULL`)
    791: -  ctx - location to stash Jacobian context (or `NULL`)

    793:    Level: advanced

    795: .seealso: [](ch:nep), `NEPSetJacobian()`
    796: @*/
    797: PetscErrorCode NEPGetJacobian(NEP nep,Mat *J,NEPJacobianFn **jac,PetscCtxRt ctx)
    798: {
    799:   PetscFunctionBegin;
    801:   NEPCheckCallback(nep,1);
    802:   if (J)   *J   = nep->jacobian;
    803:   if (jac) *jac = nep->computejacobian;
    804:   if (ctx) *(void**)ctx = nep->jacobianctx;
    805:   PetscFunctionReturn(PETSC_SUCCESS);
    806: }

    808: /*@
    809:    NEPSetSplitOperator - Sets the operator of the nonlinear eigenvalue problem
    810:    in split form.

    812:    Collective

    814:    Input Parameters:
    815: +  nep - the nonlinear eigensolver context
    816: .  nt  - number of terms in the split form
    817: .  A   - array of matrices
    818: .  f   - array of functions
    819: -  str - structure flag for matrices

    821:    Notes:
    822:    The nonlinear operator is written as $T(\lambda) = \sum_i A_i f_i(\lambda)$,
    823:    for $i=1,\dots,n$. The derivative $T'(\lambda)$ can be obtained using the
    824:    derivatives of $f_i$.

    826:    The structure flag provides information about $A_i$'s nonzero pattern
    827:    (see `MatStructure`). If all matrices have the same pattern, then
    828:    use `SAME_NONZERO_PATTERN`. If the patterns are different but contained
    829:    in the pattern of the first one, then use `SUBSET_NONZERO_PATTERN`. If
    830:    patterns are known to be different, use `DIFFERENT_NONZERO_PATTERN`.
    831:    If set to `UNKNOWN_NONZERO_PATTERN`, the patterns will be compared to
    832:    determine if they are equal.

    834:    This function must be called before `NEPSetUp()`. If it is called again
    835:    after `NEPSetUp()` then the `NEP` object is reset.

    837:    Level: beginner

    839: .seealso: [](ch:nep), `NEPGetSplitOperatorTerm()`, `NEPGetSplitOperatorInfo()`, `NEPSetSplitPreconditioner()`
    840: @*/
    841: PetscErrorCode NEPSetSplitOperator(NEP nep,PetscInt nt,Mat A[],FN f[],MatStructure str)
    842: {
    843:   PetscInt       i,n=0,m,m0=0,mloc,nloc,mloc0=0;

    845:   PetscFunctionBegin;
    848:   PetscCheck(nt>0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Must have one or more terms, you have %" PetscInt_FMT,nt);
    849:   PetscAssertPointer(A,3);
    850:   PetscAssertPointer(f,4);

    853:   for (i=0;i<nt;i++) {
    855:     PetscCheckSameComm(nep,1,A[i],3);
    857:     PetscCheckSameComm(nep,1,f[i],4);
    858:     PetscCall(MatGetSize(A[i],&m,&n));
    859:     PetscCall(MatGetLocalSize(A[i],&mloc,&nloc));
    860:     PetscCheck(m==n,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"A[%" PetscInt_FMT "] is a non-square matrix (%" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",i,m,n);
    861:     PetscCheck(mloc==nloc,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"A[%" PetscInt_FMT "] does not have equal row and column local sizes (%" PetscInt_FMT ", %" PetscInt_FMT ")",i,mloc,nloc);
    862:     if (!i) { m0 = m; mloc0 = mloc; }
    863:     PetscCheck(m==m0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_INCOMP,"Dimensions of A[%" PetscInt_FMT "] do not match with previous matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")",i,m,m0);
    864:     PetscCheck(mloc==mloc0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_INCOMP,"Local dimensions of A[%" PetscInt_FMT "] do not match with previous matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")",i,mloc,mloc0);
    865:     PetscCall(PetscObjectReference((PetscObject)A[i]));
    866:     PetscCall(PetscObjectReference((PetscObject)f[i]));
    867:   }

    869:   if (nep->state && (n!=nep->n || nloc!=nep->nloc)) PetscCall(NEPReset(nep));
    870:   else PetscCall(NEPReset_Problem(nep));

    872:   /* allocate space and copy matrices and functions */
    873:   PetscCall(PetscMalloc1(nt,&nep->A));
    874:   for (i=0;i<nt;i++) nep->A[i] = A[i];
    875:   PetscCall(PetscMalloc1(nt,&nep->f));
    876:   for (i=0;i<nt;i++) nep->f[i] = f[i];
    877:   PetscCall(PetscCalloc1(nt,&nep->nrma));
    878:   nep->nt    = nt;
    879:   nep->mstr  = str;
    880:   nep->fui   = NEP_USER_INTERFACE_SPLIT;
    881:   nep->state = NEP_STATE_INITIAL;
    882:   PetscFunctionReturn(PETSC_SUCCESS);
    883: }

    885: /*@
    886:    NEPGetSplitOperatorTerm - Gets the matrices and functions associated with
    887:    the nonlinear operator in split form.

    889:    Collective

    891:    Input Parameters:
    892: +  nep - the nonlinear eigensolver context
    893: -  k   - the index of the requested term (starting in 0)

    895:    Output Parameters:
    896: +  A - the matrix of the requested term
    897: -  f - the function of the requested term

    899:    Level: intermediate

    901: .seealso: [](ch:nep), `NEPSetSplitOperator()`, `NEPGetSplitOperatorInfo()`
    902: @*/
    903: PetscErrorCode NEPGetSplitOperatorTerm(NEP nep,PetscInt k,Mat *A,FN *f)
    904: {
    905:   PetscFunctionBegin;
    908:   NEPCheckSplit(nep,1);
    909:   PetscCheck(k>=0 && k<nep->nt,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,nep->nt-1);
    910:   if (A) *A = nep->A[k];
    911:   if (f) *f = nep->f[k];
    912:   PetscFunctionReturn(PETSC_SUCCESS);
    913: }

    915: /*@
    916:    NEPGetSplitOperatorInfo - Returns the number of terms of the split form of
    917:    the nonlinear operator, as well as the structure flag for matrices.

    919:    Not Collective

    921:    Input Parameter:
    922: .  nep - the nonlinear eigensolver context

    924:    Output Parameters:
    925: +  n   - the number of terms passed in `NEPSetSplitOperator()`
    926: -  str - the matrix structure flag passed in `NEPSetSplitOperator()`

    928:    Level: intermediate

    930: .seealso: [](ch:nep), `NEPSetSplitOperator()`, `NEPGetSplitOperatorTerm()`
    931: @*/
    932: PetscErrorCode NEPGetSplitOperatorInfo(NEP nep,PetscInt *n,MatStructure *str)
    933: {
    934:   PetscFunctionBegin;
    936:   NEPCheckSplit(nep,1);
    937:   if (n)   *n = nep->nt;
    938:   if (str) *str = nep->mstr;
    939:   PetscFunctionReturn(PETSC_SUCCESS);
    940: }

    942: /*@
    943:    NEPSetSplitPreconditioner - Sets an operator in split form from which
    944:    to build the preconditioner to be used when solving the nonlinear
    945:    eigenvalue problem in split form.

    947:    Collective

    949:    Input Parameters:
    950: +  nep  - the nonlinear eigensolver context
    951: .  ntp  - number of terms in the split preconditioner
    952: .  P    - array of matrices
    953: -  strp - structure flag for matrices

    955:    Notes:
    956:    The matrix for the preconditioner is expressed as $P(\lambda) =
    957:    \sum_i P_i f_i(\lambda)$, for $i=1,\dots,n$, where the $f_i$ functions
    958:    are the same as in `NEPSetSplitOperator()`. It is not necessary to call
    959:    this function. If it is not invoked, then the preconditioner is
    960:    built from $T(\lambda)$, i.e., both matrices and functions passed in
    961:    `NEPSetSplitOperator()`.

    963:    The structure flag provides information about $P_i$'s nonzero pattern
    964:    in the same way as in `NEPSetSplitOperator()`.

    966:    If the functions defining the preconditioner operator were different
    967:    from the ones given in `NEPSetSplitOperator()`, then the split form
    968:    cannot be used. Use the callback interface instead.

    970:    Use `ntp=0` to reset a previously set split preconditioner.

    972:    Level: advanced

    974: .seealso: [](ch:nep), `NEPGetSplitPreconditionerTerm()`, `NEPGetSplitPreconditionerInfo()`, `NEPSetSplitOperator()`
    975: @*/
    976: PetscErrorCode NEPSetSplitPreconditioner(NEP nep,PetscInt ntp,Mat P[],MatStructure strp)
    977: {
    978:   PetscInt       i,n=0,m,m0=0,mloc,nloc,mloc0=0;

    980:   PetscFunctionBegin;
    983:   PetscCheck(ntp>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Negative value of ntp = %" PetscInt_FMT,ntp);
    984:   PetscCheck(nep->fui==NEP_USER_INTERFACE_SPLIT,PetscObjectComm((PetscObject)nep),PETSC_ERR_ORDER,"Must call NEPSetSplitOperator first");
    985:   PetscCheck(ntp==0 || nep->nt==ntp,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"The number of terms must be the same as in NEPSetSplitOperator()");
    986:   if (ntp) PetscAssertPointer(P,3);

    989:   for (i=0;i<ntp;i++) {
    991:     PetscCheckSameComm(nep,1,P[i],3);
    992:     PetscCall(MatGetSize(P[i],&m,&n));
    993:     PetscCall(MatGetLocalSize(P[i],&mloc,&nloc));
    994:     PetscCheck(m==n,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"P[%" PetscInt_FMT "] is a non-square matrix (%" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",i,m,n);
    995:     PetscCheck(mloc==nloc,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"P[%" PetscInt_FMT "] does not have equal row and column local sizes (%" PetscInt_FMT ", %" PetscInt_FMT ")",i,mloc,nloc);
    996:     if (!i) { m0 = m; mloc0 = mloc; }
    997:     PetscCheck(m==m0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_INCOMP,"Dimensions of P[%" PetscInt_FMT "] do not match with previous matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")",i,m,m0);
    998:     PetscCheck(mloc==mloc0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_INCOMP,"Local dimensions of P[%" PetscInt_FMT "] do not match with previous matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")",i,mloc,mloc0);
    999:     PetscCall(PetscObjectReference((PetscObject)P[i]));
    1000:   }

    1002:   PetscCheck(!nep->state,PetscObjectComm((PetscObject)nep),PETSC_ERR_ORDER,"To call this function after NEPSetUp(), you must call NEPSetSplitOperator() again");
    1003:   if (nep->P) PetscCall(MatDestroyMatrices(nep->nt,&nep->P));

    1005:   /* allocate space and copy matrices */
    1006:   if (ntp) {
    1007:     PetscCall(PetscMalloc1(ntp,&nep->P));
    1008:     for (i=0;i<ntp;i++) nep->P[i] = P[i];
    1009:   }
    1010:   nep->mstrp = strp;
    1011:   nep->state = NEP_STATE_INITIAL;
    1012:   PetscFunctionReturn(PETSC_SUCCESS);
    1013: }

    1015: /*@
    1016:    NEPGetSplitPreconditionerTerm - Gets the matrices associated with
    1017:    the split preconditioner.

    1019:    Not Collective

    1021:    Input Parameters:
    1022: +  nep - the nonlinear eigensolver context
    1023: -  k   - the index of the requested term (starting in 0)

    1025:    Output Parameter:
    1026: .  P  - the matrix of the requested term

    1028:    Level: advanced

    1030: .seealso: [](ch:nep), `NEPSetSplitPreconditioner()`, `NEPGetSplitPreconditionerInfo()`
    1031: @*/
    1032: PetscErrorCode NEPGetSplitPreconditionerTerm(NEP nep,PetscInt k,Mat *P)
    1033: {
    1034:   PetscFunctionBegin;
    1037:   PetscAssertPointer(P,3);
    1038:   NEPCheckSplit(nep,1);
    1039:   PetscCheck(k>=0 && k<nep->nt,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,nep->nt-1);
    1040:   PetscCheck(nep->P,PetscObjectComm((PetscObject)nep),PETSC_ERR_ORDER,"You have not called NEPSetSplitPreconditioner()");
    1041:   *P = nep->P[k];
    1042:   PetscFunctionReturn(PETSC_SUCCESS);
    1043: }

    1045: /*@
    1046:    NEPGetSplitPreconditionerInfo - Returns the number of terms of the split
    1047:    preconditioner, as well as the structure flag for matrices.

    1049:    Not Collective

    1051:    Input Parameter:
    1052: .  nep - the nonlinear eigensolver context

    1054:    Output Parameters:
    1055: +  n    - the number of terms passed in `NEPSetSplitPreconditioner()`
    1056: -  strp - the matrix structure flag passed in `NEPSetSplitPreconditioner()`

    1058:    Level: advanced

    1060: .seealso: [](ch:nep), `NEPSetSplitPreconditioner()`, `NEPGetSplitPreconditionerTerm()`
    1061: @*/
    1062: PetscErrorCode NEPGetSplitPreconditionerInfo(NEP nep,PetscInt *n,MatStructure *strp)
    1063: {
    1064:   PetscFunctionBegin;
    1066:   NEPCheckSplit(nep,1);
    1067:   if (n)    *n    = nep->P? nep->nt: 0;
    1068:   if (strp) *strp = nep->mstrp;
    1069:   PetscFunctionReturn(PETSC_SUCCESS);
    1070: }