Actual source code: cayley.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:    Implements the Cayley spectral transform
     12: */

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

     16: typedef struct {
     17:   PetscScalar nu;
     18:   PetscBool   nu_set;
     19: } ST_CAYLEY;

     21: static PetscErrorCode MatMult_Cayley(Mat B,Vec x,Vec y)
     22: {
     23:   ST             st;
     24:   ST_CAYLEY      *ctx;
     25:   PetscScalar    nu;

     27:   PetscFunctionBegin;
     28:   PetscCall(MatShellGetContext(B,&st));
     29:   ctx = (ST_CAYLEY*)st->data;
     30:   nu = ctx->nu;

     32:   if (st->matmode == ST_MATMODE_INPLACE) nu = nu + st->sigma;

     34:   if (st->nmat>1) {
     35:     /* generalized eigenproblem: y = (A + tB)x */
     36:     PetscCall(MatMult(st->A[0],x,y));
     37:     PetscCall(MatMult(st->A[1],x,st->work[1]));
     38:     PetscCall(VecAXPY(y,nu,st->work[1]));
     39:   } else {
     40:     /* standard eigenproblem: y = (A + tI)x */
     41:     PetscCall(MatMult(st->A[0],x,y));
     42:     PetscCall(VecAXPY(y,nu,x));
     43:   }
     44:   PetscFunctionReturn(PETSC_SUCCESS);
     45: }

     47: static PetscErrorCode MatMultTranspose_Cayley(Mat B,Vec x,Vec y)
     48: {
     49:   ST             st;
     50:   ST_CAYLEY      *ctx;
     51:   PetscScalar    nu;

     53:   PetscFunctionBegin;
     54:   PetscCall(MatShellGetContext(B,&st));
     55:   ctx = (ST_CAYLEY*)st->data;
     56:   nu = ctx->nu;

     58:   if (st->matmode == ST_MATMODE_INPLACE) nu = nu + st->sigma;
     59:   nu = PetscConj(nu);

     61:   if (st->nmat>1) {
     62:     /* generalized eigenproblem: y = (A + tB)x */
     63:     PetscCall(MatMultTranspose(st->A[0],x,y));
     64:     PetscCall(MatMultTranspose(st->A[1],x,st->work[1]));
     65:     PetscCall(VecAXPY(y,nu,st->work[1]));
     66:   } else {
     67:     /* standard eigenproblem: y = (A + tI)x */
     68:     PetscCall(MatMultTranspose(st->A[0],x,y));
     69:     PetscCall(VecAXPY(y,nu,x));
     70:   }
     71:   PetscFunctionReturn(PETSC_SUCCESS);
     72: }

     74: static PetscErrorCode STGetBilinearForm_Cayley(ST st,Mat *B)
     75: {
     76:   PetscFunctionBegin;
     77:   PetscCall(STSetUp(st));
     78:   *B = st->T[0];
     79:   PetscCall(PetscObjectReference((PetscObject)*B));
     80:   PetscFunctionReturn(PETSC_SUCCESS);
     81: }

     83: static PetscErrorCode STBackTransform_Cayley(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
     84: {
     85:   ST_CAYLEY   *ctx = (ST_CAYLEY*)st->data;
     86:   PetscInt    j;
     87: #if !defined(PETSC_USE_COMPLEX)
     88:   PetscScalar t,i,r;
     89: #endif

     91:   PetscFunctionBegin;
     92: #if !defined(PETSC_USE_COMPLEX)
     93:   for (j=0;j<n;j++) {
     94:     if (eigi[j] == 0.0) eigr[j] = (ctx->nu + eigr[j] * st->sigma) / (eigr[j] - 1.0);
     95:     else {
     96:       r = eigr[j];
     97:       i = eigi[j];
     98:       r = st->sigma * (r * r + i * i - r) + ctx->nu * (r - 1);
     99:       i = - st->sigma * i - ctx->nu * i;
    100:       t = i * i + r * (r - 2.0) + 1.0;
    101:       eigr[j] = r / t;
    102:       eigi[j] = i / t;
    103:     }
    104:   }
    105: #else
    106:   for (j=0;j<n;j++) {
    107:     eigr[j] = (ctx->nu + eigr[j] * st->sigma) / (eigr[j] - 1.0);
    108:   }
    109: #endif
    110:   PetscFunctionReturn(PETSC_SUCCESS);
    111: }

    113: static PetscErrorCode STPostSolve_Cayley(ST st)
    114: {
    115:   PetscFunctionBegin;
    116:   if (st->matmode == ST_MATMODE_INPLACE) {
    117:     if (st->nmat>1) PetscCall(MatAXPY(st->A[0],st->sigma,st->A[1],st->str));
    118:     else PetscCall(MatShift(st->A[0],st->sigma));
    119:     st->Astate[0] = ((PetscObject)st->A[0])->state;
    120:     st->state   = ST_STATE_INITIAL;
    121:     st->opready = PETSC_FALSE;
    122:   }
    123:   PetscFunctionReturn(PETSC_SUCCESS);
    124: }

    126: /*
    127:    Operator (cayley):
    128:                Op                  P         M
    129:    if nmat=1:  (A-sI)^-1 (A+tI)    A-sI      A+tI
    130:    if nmat=2:  (A-sB)^-1 (A+tB)    A-sB      A+tI
    131: */
    132: static PetscErrorCode STComputeOperator_Cayley(ST st)
    133: {
    134:   PetscInt       n,m;
    135:   ST_CAYLEY      *ctx = (ST_CAYLEY*)st->data;

    137:   PetscFunctionBegin;
    138:   /* if the user did not set the shift, use the target value */
    139:   if (!st->sigma_set) st->sigma = st->defsigma;

    141:   if (!ctx->nu_set) ctx->nu = st->sigma;
    142:   PetscCheck(ctx->nu!=0.0 || st->sigma!=0.0,PetscObjectComm((PetscObject)st),PETSC_ERR_USER_INPUT,"Values of shift and antishift cannot be zero simultaneously");
    143:   PetscCheck(ctx->nu!=-st->sigma,PetscObjectComm((PetscObject)st),PETSC_ERR_USER_INPUT,"It is not allowed to set the antishift equal to minus the shift (the target)");

    145:   /* T[0] = A+nu*B */
    146:   if (st->matmode==ST_MATMODE_INPLACE) {
    147:     PetscCall(MatGetLocalSize(st->A[0],&n,&m));
    148:     PetscCall(MatCreateShell(PetscObjectComm((PetscObject)st),n,m,PETSC_DETERMINE,PETSC_DETERMINE,st,&st->T[0]));
    149:     PetscCall(MatShellSetOperation(st->T[0],MATOP_MULT,(PetscErrorCodeFn*)MatMult_Cayley));
    150:     PetscCall(MatShellSetOperation(st->T[0],MATOP_MULT_TRANSPOSE,(PetscErrorCodeFn*)MatMultTranspose_Cayley));
    151:   } else PetscCall(STMatMAXPY_Private(st,ctx->nu,0.0,0,NULL,PetscNot(st->state==ST_STATE_UPDATED),PETSC_FALSE,&st->T[0]));
    152:   st->M = st->T[0];

    154:   /* T[1] = A-sigma*B */
    155:   PetscCall(STMatMAXPY_Private(st,-st->sigma,0.0,0,NULL,PetscNot(st->state==ST_STATE_UPDATED),PETSC_FALSE,&st->T[1]));
    156:   PetscCall(PetscObjectReference((PetscObject)st->T[1]));
    157:   PetscCall(MatDestroy(&st->P));
    158:   st->P = st->T[1];
    159:   if (st->Psplit) {  /* build custom preconditioner from the split matrices */
    160:     PetscCall(STMatMAXPY_Private(st,-st->sigma,0.0,0,NULL,PETSC_TRUE,PETSC_TRUE,&st->Pmat));
    161:   }
    162:   PetscFunctionReturn(PETSC_SUCCESS);
    163: }

    165: static PetscErrorCode STSetUp_Cayley(ST st)
    166: {
    167:   PetscFunctionBegin;
    168:   PetscCheck(st->nmat<=2,PetscObjectComm((PetscObject)st),PETSC_ERR_SUP,"Cayley transform cannot be used in polynomial eigenproblems");
    169:   PetscCall(STSetWorkVecs(st,2));
    170:   PetscCall(KSPSetUp(st->ksp));
    171:   PetscFunctionReturn(PETSC_SUCCESS);
    172: }

    174: static PetscErrorCode STSetShift_Cayley(ST st,PetscScalar newshift)
    175: {
    176:   ST_CAYLEY      *ctx = (ST_CAYLEY*)st->data;

    178:   PetscFunctionBegin;
    179:   PetscCheck(newshift!=0.0 || (ctx->nu_set && ctx->nu!=0.0),PetscObjectComm((PetscObject)st),PETSC_ERR_USER_INPUT,"Values of shift and antishift cannot be zero simultaneously");
    180:   PetscCheck(ctx->nu!=-newshift,PetscObjectComm((PetscObject)st),PETSC_ERR_USER_INPUT,"It is not allowed to set the shift equal to minus the antishift");

    182:   if (!ctx->nu_set) {
    183:     if (st->matmode!=ST_MATMODE_INPLACE) PetscCall(STMatMAXPY_Private(st,newshift,ctx->nu,0,NULL,PETSC_FALSE,PETSC_FALSE,&st->T[0]));
    184:     ctx->nu = newshift;
    185:   }
    186:   PetscCall(STMatMAXPY_Private(st,-newshift,-st->sigma,0,NULL,PETSC_FALSE,PETSC_FALSE,&st->T[1]));
    187:   if (st->P!=st->T[1]) {
    188:     PetscCall(PetscObjectReference((PetscObject)st->T[1]));
    189:     PetscCall(MatDestroy(&st->P));
    190:     st->P = st->T[1];
    191:   }
    192:   if (st->Psplit) {  /* build custom preconditioner from the split matrices */
    193:     PetscCall(STMatMAXPY_Private(st,-newshift,-st->sigma,0,NULL,PETSC_FALSE,PETSC_TRUE,&st->Pmat));
    194:   }
    195:   PetscCall(ST_KSPSetOperators(st,st->P,st->Pmat?st->Pmat:st->P));
    196:   PetscCall(KSPSetUp(st->ksp));
    197:   PetscFunctionReturn(PETSC_SUCCESS);
    198: }

    200: static PetscErrorCode STSetFromOptions_Cayley(ST st,PetscOptionItems PetscOptionsObject)
    201: {
    202:   PetscScalar    nu;
    203:   PetscBool      flg;
    204:   ST_CAYLEY      *ctx = (ST_CAYLEY*)st->data;

    206:   PetscFunctionBegin;
    207:   PetscOptionsHeadBegin(PetscOptionsObject,"ST Cayley Options");

    209:     PetscCall(PetscOptionsScalar("-st_cayley_antishift","Value of the antishift","STCayleySetAntishift",ctx->nu,&nu,&flg));
    210:     if (flg) PetscCall(STCayleySetAntishift(st,nu));

    212:   PetscOptionsHeadEnd();
    213:   PetscFunctionReturn(PETSC_SUCCESS);
    214: }

    216: static PetscErrorCode STCayleySetAntishift_Cayley(ST st,PetscScalar newshift)
    217: {
    218:   ST_CAYLEY *ctx = (ST_CAYLEY*)st->data;

    220:   PetscFunctionBegin;
    221:   if (ctx->nu != newshift) {
    222:     STCheckNotSeized(st,1);
    223:     if (st->state && st->matmode!=ST_MATMODE_INPLACE) PetscCall(STMatMAXPY_Private(st,newshift,ctx->nu,0,NULL,PETSC_FALSE,PETSC_FALSE,&st->T[0]));
    224:     ctx->nu = newshift;
    225:   }
    226:   ctx->nu_set = PETSC_TRUE;
    227:   PetscFunctionReturn(PETSC_SUCCESS);
    228: }

    230: /*@
    231:    STCayleySetAntishift - Sets the value of the anti-shift for the Cayley
    232:    spectral transformation.

    234:    Logically Collective

    236:    Input Parameters:
    237: +  st  - the spectral transformation context
    238: -  nu  - the anti-shift

    240:    Options Database Key:
    241: .  -st_cayley_antishift nu - sets the value of the anti-shift

    243:    Level: intermediate

    245:    Note:
    246:    In the generalized Cayley transform, the operator can be expressed as
    247:    $(A - \sigma B)^{-1}(A + \nu B)$. This function sets the value of $\nu$.
    248:    Use `STSetShift()` for setting $\sigma$. The value $\nu=-\sigma$ is not allowed.

    250: .seealso: [](ch:st), `STCAYLEY`, `STSetShift()`, `STCayleyGetAntishift()`
    251: @*/
    252: PetscErrorCode STCayleySetAntishift(ST st,PetscScalar nu)
    253: {
    254:   PetscFunctionBegin;
    257:   PetscTryMethod(st,"STCayleySetAntishift_C",(ST,PetscScalar),(st,nu));
    258:   PetscFunctionReturn(PETSC_SUCCESS);
    259: }

    261: static PetscErrorCode STCayleyGetAntishift_Cayley(ST st,PetscScalar *nu)
    262: {
    263:   ST_CAYLEY *ctx = (ST_CAYLEY*)st->data;

    265:   PetscFunctionBegin;
    266:   *nu = ctx->nu;
    267:   PetscFunctionReturn(PETSC_SUCCESS);
    268: }

    270: /*@
    271:    STCayleyGetAntishift - Gets the value of the anti-shift used in the Cayley
    272:    spectral transformation.

    274:    Not Collective

    276:    Input Parameter:
    277: .  st  - the spectral transformation context

    279:    Output Parameter:
    280: .  nu  - the anti-shift

    282:    Level: intermediate

    284: .seealso: [](ch:st), `STCAYLEY`, `STGetShift()`, `STCayleySetAntishift()`
    285: @*/
    286: PetscErrorCode STCayleyGetAntishift(ST st,PetscScalar *nu)
    287: {
    288:   PetscFunctionBegin;
    290:   PetscAssertPointer(nu,2);
    291:   PetscUseMethod(st,"STCayleyGetAntishift_C",(ST,PetscScalar*),(st,nu));
    292:   PetscFunctionReturn(PETSC_SUCCESS);
    293: }

    295: static PetscErrorCode STView_Cayley(ST st,PetscViewer viewer)
    296: {
    297:   char           str[50];
    298:   ST_CAYLEY      *ctx = (ST_CAYLEY*)st->data;
    299:   PetscBool      isascii;

    301:   PetscFunctionBegin;
    302:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
    303:   if (isascii) {
    304:     PetscCall(SlepcSNPrintfScalar(str,sizeof(str),ctx->nu,PETSC_FALSE));
    305:     PetscCall(PetscViewerASCIIPrintf(viewer,"  antishift: %s\n",str));
    306:   }
    307:   PetscFunctionReturn(PETSC_SUCCESS);
    308: }

    310: static PetscErrorCode STDestroy_Cayley(ST st)
    311: {
    312:   PetscFunctionBegin;
    313:   PetscCall(PetscFree(st->data));
    314:   PetscCall(PetscObjectComposeFunction((PetscObject)st,"STCayleySetAntishift_C",NULL));
    315:   PetscCall(PetscObjectComposeFunction((PetscObject)st,"STCayleyGetAntishift_C",NULL));
    316:   PetscFunctionReturn(PETSC_SUCCESS);
    317: }

    319: /*MC
    320:    STCAYLEY - STCAYLEY = "cayley" - The generalized Cayley transform.

    322:    Level: beginner

    324:    Notes:
    325:    The operator in this `ST` is $(A-\sigma B)^{-1}(A+\nu B)$, where $A$ and $B$
    326:    are the matrices set with `STSetMatrices()`, the shift $\sigma$ is given in
    327:    `STSetShift()`, and $\nu$ is given in `STCayleySetAntishift()`.

    329:    This transformation is mathematically equivalent to `STSINVERT`.

    331: .seealso: [](ch:st), `ST`, `STType`, `STSetType()`, `STSetMatrices()`, `STSetShift()`, `STCayleySetAntishift()`
    332: M*/

    334: SLEPC_EXTERN PetscErrorCode STCreate_Cayley(ST st)
    335: {
    336:   ST_CAYLEY      *ctx;

    338:   PetscFunctionBegin;
    339:   PetscCall(PetscNew(&ctx));
    340:   st->data = (void*)ctx;

    342:   st->usesksp = PETSC_TRUE;

    344:   st->ops->apply           = STApply_Generic;
    345:   st->ops->applytrans      = STApplyTranspose_Generic;
    346:   st->ops->applyhermtrans  = STApplyHermitianTranspose_Generic;
    347:   st->ops->backtransform   = STBackTransform_Cayley;
    348:   st->ops->setshift        = STSetShift_Cayley;
    349:   st->ops->getbilinearform = STGetBilinearForm_Cayley;
    350:   st->ops->setup           = STSetUp_Cayley;
    351:   st->ops->computeoperator = STComputeOperator_Cayley;
    352:   st->ops->setfromoptions  = STSetFromOptions_Cayley;
    353:   st->ops->postsolve       = STPostSolve_Cayley;
    354:   st->ops->destroy         = STDestroy_Cayley;
    355:   st->ops->view            = STView_Cayley;
    356:   st->ops->checknullspace  = STCheckNullSpace_Default;
    357:   st->ops->setdefaultksp   = STSetDefaultKSP_Default;

    359:   PetscCall(PetscObjectComposeFunction((PetscObject)st,"STCayleySetAntishift_C",STCayleySetAntishift_Cayley));
    360:   PetscCall(PetscObjectComposeFunction((PetscObject)st,"STCayleyGetAntishift_C",STCayleyGetAntishift_Cayley));
    361:   PetscFunctionReturn(PETSC_SUCCESS);
    362: }