Actual source code: lmesetup.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:    LME routines related to problem setup
     12: */

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

     16: static inline PetscErrorCode LMESetUp_Lyapunov(LME lme)
     17: {
     18:   Mat            C1,C2,X1,X2;
     19:   Vec            dc,dx;

     21:   PetscFunctionBegin;
     22:   PetscCall(MatLRCGetMats(lme->C,NULL,&C1,&dc,&C2));
     23:   PetscCheck(C1==C2,PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"Lyapunov matrix equation requires symmetric right-hand side C");
     24:   PetscCheck(!dc,PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"Lyapunov solvers currently require positive-definite right-hand side C");
     25:   if (lme->X) {
     26:     PetscCall(MatLRCGetMats(lme->X,NULL,&X1,&dx,&X2));
     27:     PetscCheck(X1==X2,PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"Lyapunov matrix equation requires symmetric solution X");
     28:     PetscCheck(!dx,PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"Lyapunov solvers currently assume a positive-definite solution X");
     29:   }
     30:   PetscFunctionReturn(PETSC_SUCCESS);
     31: }

     33: /*@
     34:    LMESetUp - Sets up all the internal data structures necessary for the
     35:    execution of the linear matrix equation solver.

     37:    Collective

     39:    Input Parameter:
     40: .  lme - the linear matrix equation solver context

     42:    Notes:
     43:    This function need not be called explicitly in most cases, since `LMESolve()`
     44:    calls it. It can be useful when one wants to measure the set-up time
     45:    separately from the solve time.

     47:    Level: developer

     49: .seealso: [](ch:lme), `LMECreate()`, `LMESolve()`, `LMEDestroy()`
     50: @*/
     51: PetscErrorCode LMESetUp(LME lme)
     52: {
     53:   PetscInt       N;

     55:   PetscFunctionBegin;

     58:   /* reset the convergence flag from the previous solves */
     59:   lme->reason = LME_CONVERGED_ITERATING;

     61:   if (lme->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
     62:   PetscCall(PetscLogEventBegin(LME_SetUp,lme,0,0,0));

     64:   /* Set default solver type (LMESetFromOptions was not called) */
     65:   if (!((PetscObject)lme)->type_name) PetscCall(LMESetType(lme,LMEKRYLOV));

     67:   /* Check problem dimensions */
     68:   PetscCheck(lme->A,PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"LMESetCoefficients must be called first");
     69:   PetscCheck(lme->C,PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONGSTATE,"LMESetRHS must be called first");
     70:   PetscCall(MatGetSize(lme->A,&N,NULL));
     71:   if (lme->ncv > N) lme->ncv = N;

     73:   /* setup options for the particular equation type */
     74:   switch (lme->problem_type) {
     75:     case LME_LYAPUNOV:
     76:       PetscCall(LMESetUp_Lyapunov(lme));
     77:       break;
     78:     case LME_SYLVESTER:
     79:       LMECheckCoeff(lme,lme->B,"B","Sylvester");
     80:       break;
     81:     case LME_GEN_LYAPUNOV:
     82:       LMECheckCoeff(lme,lme->D,"D","Generalized Lyapunov");
     83:       break;
     84:     case LME_GEN_SYLVESTER:
     85:       LMECheckCoeff(lme,lme->B,"B","Generalized Sylvester");
     86:       LMECheckCoeff(lme,lme->D,"D","Generalized Sylvester");
     87:       LMECheckCoeff(lme,lme->E,"E","Generalized Sylvester");
     88:       break;
     89:     case LME_DT_LYAPUNOV:
     90:       break;
     91:     case LME_STEIN:
     92:       LMECheckCoeff(lme,lme->E,"E","Stein");
     93:       break;
     94:   }
     95:   PetscCheck(lme->problem_type==LME_LYAPUNOV,PetscObjectComm((PetscObject)lme),PETSC_ERR_SUP,"There is no solver yet for this matrix equation type");

     97:   /* call specific solver setup */
     98:   PetscUseTypeMethod(lme,setup);

    100:   /* set tolerance if not yet set */
    101:   if (lme->tol==(PetscReal)PETSC_DETERMINE) lme->tol = SLEPC_DEFAULT_TOL;

    103:   PetscCall(PetscLogEventEnd(LME_SetUp,lme,0,0,0));
    104:   lme->setupcalled = 1;
    105:   PetscFunctionReturn(PETSC_SUCCESS);
    106: }

    108: static inline PetscErrorCode LMESetCoefficients_Private(LME lme,Mat A,Mat *lmeA)
    109: {
    110:   PetscInt       m,n;

    112:   PetscFunctionBegin;
    113:   PetscCall(MatGetSize(A,&m,&n));
    114:   PetscCheck(m==n,PetscObjectComm((PetscObject)lme),PETSC_ERR_ARG_WRONG,"Matrix is non-square");
    115:   if (!lme->setupcalled) PetscCall(MatDestroy(lmeA));
    116:   PetscCall(PetscObjectReference((PetscObject)A));
    117:   *lmeA = A;
    118:   PetscFunctionReturn(PETSC_SUCCESS);
    119: }

    121: /*@
    122:    LMESetCoefficients - Sets the coefficient matrices that define the linear matrix
    123:    equation to be solved.

    125:    Collective

    127:    Input Parameters:
    128: +  lme - the linear matrix equation solver context
    129: .  A   - first coefficient matrix
    130: .  B   - second coefficient matrix
    131: .  D   - third coefficient matrix
    132: -  E   - fourth coefficient matrix

    134:    Notes:
    135:    The matrix equation takes the general form $AXE+DXB=C$, where matrix $C$ is not
    136:    provided here but with `LMESetRHS()`. Not all four matrices must be passed, some
    137:    can be `NULL` instead, see `LMESetProblemType()` for details.

    139:    It must be called before `LMESetUp()`. If it is called again after `LMESetUp()` then
    140:    the `LME` object is reset.

    142:    In order to delete a previously set matrix, pass `NULL` in the corresponding
    143:    argument.

    145:    Level: beginner

    147: .seealso: [](ch:lme), `LMESolve()`, `LMESetUp()`, `LMESetRHS()`, `LMESetProblemType()`
    148: @*/
    149: PetscErrorCode LMESetCoefficients(LME lme,Mat A,Mat B,Mat D,Mat E)
    150: {
    151:   PetscFunctionBegin;
    154:   PetscCheckSameComm(lme,1,A,2);
    155:   if (B) {
    157:     PetscCheckSameComm(lme,1,B,3);
    158:   }
    159:   if (D) {
    161:     PetscCheckSameComm(lme,1,D,4);
    162:   }
    163:   if (E) {
    165:     PetscCheckSameComm(lme,1,E,5);
    166:   }

    168:   if (lme->setupcalled) PetscCall(LMEReset(lme));

    170:   PetscCall(LMESetCoefficients_Private(lme,A,&lme->A));
    171:   if (B) PetscCall(LMESetCoefficients_Private(lme,B,&lme->B));
    172:   else if (!lme->setupcalled) PetscCall(MatDestroy(&lme->B));
    173:   if (D) PetscCall(LMESetCoefficients_Private(lme,D,&lme->D));
    174:   else if (!lme->setupcalled) PetscCall(MatDestroy(&lme->D));
    175:   if (E) PetscCall(LMESetCoefficients_Private(lme,E,&lme->E));
    176:   else if (!lme->setupcalled) PetscCall(MatDestroy(&lme->E));

    178:   lme->setupcalled = 0;
    179:   PetscFunctionReturn(PETSC_SUCCESS);
    180: }

    182: /*@
    183:    LMEGetCoefficients - Gets the coefficient matrices of the matrix equation.

    185:    Collective

    187:    Input Parameter:
    188: .  lme - the linear matrix equation solver context

    190:    Output Parameters:
    191: +  A   - first coefficient matrix
    192: .  B   - second coefficient matrix
    193: .  D   - third coefficient matrix
    194: -  E   - fourth coefficient matrix

    196:    Level: intermediate

    198: .seealso: [](ch:lme), `LMESolve()`, `LMESetCoefficients()`
    199: @*/
    200: PetscErrorCode LMEGetCoefficients(LME lme,Mat *A,Mat *B,Mat *D,Mat *E)
    201: {
    202:   PetscFunctionBegin;
    204:   if (A) *A = lme->A;
    205:   if (B) *B = lme->B;
    206:   if (D) *D = lme->D;
    207:   if (E) *E = lme->E;
    208:   PetscFunctionReturn(PETSC_SUCCESS);
    209: }

    211: /*@
    212:    LMESetRHS - Sets the right-hand side of the matrix equation, as a low-rank
    213:    matrix.

    215:    Collective

    217:    Input Parameters:
    218: +  lme - the linear matrix equation solver context
    219: -  C   - the right-hand side matrix

    221:    Notes:
    222:    The matrix equation takes the general form $AXE+DXB=C$, where matrix $C$ is
    223:    given with this function. $C$ must be a low-rank matrix of type `MATLRC`, that is,
    224:    $C = UDV^*$ where $D$ is diagonal of order $k$, and $U$, $V$ are dense tall-skinny
    225:    matrices with $k$ columns. No sparse matrix must be provided when creating the
    226:    `MATLRC` matrix.

    228:    In equation types that require $C$ to be symmetric, such as Lyapunov, `C` must be
    229:    created with $V=U$ (or `V=NULL`).

    231:    Level: beginner

    233: .seealso: [](ch:lme), `LMESetSolution()`, `LMESetProblemType()`, `LMEGetRHS()`
    234: @*/
    235: PetscErrorCode LMESetRHS(LME lme,Mat C)
    236: {
    237:   Mat            A;

    239:   PetscFunctionBegin;
    242:   PetscCheckSameComm(lme,1,C,2);
    243:   PetscCheckTypeName(C,MATLRC);

    245:   PetscCall(MatLRCGetMats(C,&A,NULL,NULL,NULL));
    246:   PetscCheck(!A,PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"The MatLRC must not have a sparse matrix term");

    248:   PetscCall(PetscObjectReference((PetscObject)C));
    249:   PetscCall(MatDestroy(&lme->C));
    250:   lme->C = C;
    251:   PetscFunctionReturn(PETSC_SUCCESS);
    252: }

    254: /*@
    255:    LMEGetRHS - Gets the right-hand side of the matrix equation.

    257:    Collective

    259:    Input Parameter:
    260: .  lme - the linear matrix equation solver context

    262:    Output Parameter:
    263: .  C   - the right-hand side matrix

    265:    Level: intermediate

    267: .seealso: [](ch:lme), `LMESolve()`, `LMESetRHS()`
    268: @*/
    269: PetscErrorCode LMEGetRHS(LME lme,Mat *C)
    270: {
    271:   PetscFunctionBegin;
    273:   PetscAssertPointer(C,2);
    274:   *C = lme->C;
    275:   PetscFunctionReturn(PETSC_SUCCESS);
    276: }

    278: /*@
    279:    LMESetSolution - Sets the placeholder for the solution of the matrix
    280:    equation, as a low-rank matrix.

    282:    Collective

    284:    Input Parameters:
    285: +  lme - the linear matrix equation solver context
    286: -  X   - the solution matrix

    288:    Notes:
    289:    The matrix equation takes the general form $AXE+DXB=C$, where the solution
    290:    matrix is of low rank and is written in factored form $X = UDV^*$. This function
    291:    provides a `Mat` object of type `MATLRC` that stores $U$, $V$ and (optionally) $D$.
    292:    These factors will be computed during `LMESolve()`.

    294:    In equation types whose solution $X$ is symmetric, such as Lyapunov, $X$ must be
    295:    created with $V=U$ (or `V=NULL`).

    297:    If the user provides `X` with this function, then the solver will
    298:    return a solution with rank at most the number of columns of $U$. Alternatively,
    299:    it is possible to let the solver choose the rank of the solution, by
    300:    setting `X` to `NULL` and then calling `LMEGetSolution()` after `LMESolve()`.

    302:    Level: intermediate

    304: .seealso: [](ch:lme), `LMEGetSolution()`, `LMESetRHS()`, `LMESetProblemType()`, `LMESolve()`
    305: @*/
    306: PetscErrorCode LMESetSolution(LME lme,Mat X)
    307: {
    308:   Mat            A;

    310:   PetscFunctionBegin;
    312:   if (X) {
    314:     PetscCheckSameComm(lme,1,X,2);
    315:     PetscCheckTypeName(X,MATLRC);
    316:     PetscCall(MatLRCGetMats(X,&A,NULL,NULL,NULL));
    317:     PetscCheck(!A,PetscObjectComm((PetscObject)X),PETSC_ERR_SUP,"The MatLRC must not have a sparse matrix term");
    318:     PetscCall(PetscObjectReference((PetscObject)X));
    319:   }
    320:   PetscCall(MatDestroy(&lme->X));
    321:   lme->X = X;
    322:   PetscFunctionReturn(PETSC_SUCCESS);
    323: }

    325: /*@
    326:    LMEGetSolution - Gets the solution of the matrix equation.

    328:    Collective

    330:    Input Parameter:
    331: .  lme - the linear matrix equation solver context

    333:    Output Parameter:
    334: .  X   - the low-rank matrix

    336:    Notes:
    337:    If called after `LMESolve()`, `X` will contain the solution of the equation.

    339:    The matrix `X` may have been passed by the user via `LMESetSolution()`,
    340:    although this is not required.

    342:    Level: beginner

    344: .seealso: [](ch:lme), `LMESolve()`, `LMESetSolution()`
    345: @*/
    346: PetscErrorCode LMEGetSolution(LME lme,Mat *X)
    347: {
    348:   PetscFunctionBegin;
    350:   PetscAssertPointer(X,2);
    351:   *X = lme->X;
    352:   PetscFunctionReturn(PETSC_SUCCESS);
    353: }

    355: /*@
    356:    LMEAllocateSolution - Allocate memory storage for common variables such
    357:    as the basis vectors.

    359:    Collective

    361:    Input Parameters:
    362: +  lme   - the linear matrix equation solver context
    363: -  extra - number of additional positions, used for methods that require a
    364:            working basis slightly larger than `ncv`

    366:    Developer Notes:
    367:    This is `SLEPC_EXTERN` because it may be required by user plugin `LME`
    368:    implementations.

    370:    Level: developer

    372: .seealso: [](ch:lme), `LMESetUp()`, `LMESetDimensions()`
    373: @*/
    374: PetscErrorCode LMEAllocateSolution(LME lme,PetscInt extra)
    375: {
    376:   PetscInt       oldsize,requested;
    377:   Vec            t;

    379:   PetscFunctionBegin;
    380:   requested = lme->ncv + extra;

    382:   /* oldsize is zero if this is the first time setup is called */
    383:   PetscCall(BVGetSizes(lme->V,NULL,NULL,&oldsize));

    385:   /* allocate basis vectors */
    386:   if (!lme->V) PetscCall(LMEGetBV(lme,&lme->V));
    387:   if (!oldsize) {
    388:     if (!((PetscObject)lme->V)->type_name) PetscCall(BVSetType(lme->V,BVMAT));
    389:     PetscCall(MatCreateVecsEmpty(lme->A,&t,NULL));
    390:     PetscCall(BVSetSizesFromVec(lme->V,t,requested));
    391:     PetscCall(VecDestroy(&t));
    392:   } else PetscCall(BVResize(lme->V,requested,PETSC_FALSE));
    393:   PetscFunctionReturn(PETSC_SUCCESS);
    394: }