Actual source code: matstruct.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: */

     11: #include <slepc/private/slepcimpl.h>

     13: /*@
     14:    MatCreateBSE - Create a matrix that can be used to define a structured eigenvalue
     15:    problem of type BSE (Bethe-Salpeter Equation).

     17:    Collective

     19:    Input Parameters:
     20: +  R - matrix for the diagonal block (resonant)
     21: -  C - matrix for the off-diagonal block (coupling)

     23:    Output Parameter:
     24: .  H  - the resulting matrix

     26:    Notes:
     27:    The resulting matrix has the block form $$H = \begin{bmatrix} R & C \\ -C^* & -R^T \end{bmatrix},$$
     28:    where $R$ is assumed
     29:    to be (complex) Hermitian and $C$ complex symmetric. Note that this function does
     30:    not check these properties, so if the matrices provided by the user do not satisfy
     31:    them, then the solver will not behave as expected.

     33:    The obtained matrix can be used as an input matrix to `EPS` eigensolvers via
     34:    `EPSSetOperators()` for the case that the problem type is `EPS_BSE`. Note that the user
     35:    cannot just build a matrix with the required structure, it must be done via this
     36:    function.

     38:    In the current implementation, `H` is a `MATNEST` matrix, where `R` and `C` form the top
     39:    block row, while the bottom block row is composed of matrices of type
     40:    `MATTRANSPOSEVIRTUAL` and `MATHERMITIANTRANSPOSEVIRTUAL` scaled by -1.

     42:    Level: intermediate

     44: .seealso: [](sec:structured), `MatCreateNest()`, `EPSSetOperators()`, `EPSSetProblemType()`, `MatCreateHamiltonian()`, `MatCreateLREP()`
     45: @*/
     46: PetscErrorCode MatCreateBSE(Mat R,Mat C,Mat *H)
     47: {
     48:   PetscInt       Mr,Mc,Nr,Nc,mr,mc,nr,nc;
     49:   Mat            block[4] = { R, C, NULL, NULL };
     50:   SlepcMatStruct mctx;

     52:   PetscFunctionBegin;
     55:   PetscCheckSameTypeAndComm(R,1,C,2);
     56:   PetscAssertPointer(H,3);

     58:   /* check sizes */
     59:   PetscCall(MatGetSize(R,&Mr,&Nr));
     60:   PetscCall(MatGetLocalSize(R,&mr,&nr));
     61:   PetscCall(MatGetSize(C,&Mc,&Nc));
     62:   PetscCall(MatGetLocalSize(C,&mc,&nc));
     63:   PetscCheck(Mc==Mr && mc==mr,PetscObjectComm((PetscObject)R),PETSC_ERR_ARG_INCOMP,"Incompatible row dimensions");
     64:   PetscCheck(Nc==Nr && nc==nr,PetscObjectComm((PetscObject)R),PETSC_ERR_ARG_INCOMP,"Incompatible column dimensions");

     66:   /* bottom block row */
     67:   PetscCall(MatCreateHermitianTranspose(C,&block[2]));
     68:   PetscCall(MatScale(block[2],-1.0));
     69:   PetscCall(MatCreateTranspose(R,&block[3]));
     70:   PetscCall(MatScale(block[3],-1.0));

     72:   /* create nest matrix and compose context */
     73:   PetscCall(MatCreateNest(PetscObjectComm((PetscObject)R),2,NULL,2,NULL,block,H));
     74:   PetscCall(MatDestroy(&block[2]));
     75:   PetscCall(MatDestroy(&block[3]));

     77:   PetscCall(PetscNew(&mctx));
     78:   mctx->cookie = SLEPC_MAT_STRUCT_BSE;
     79:   PetscCall(PetscObjectContainerCompose((PetscObject)*H,"SlepcMatStruct",mctx,PetscCtxDestroyDefault));
     80:   PetscFunctionReturn(PETSC_SUCCESS);
     81: }

     83: /*@
     84:    MatCreateHamiltonian - Create a matrix that can be used to define a structured
     85:    eigenvalue problem of Hamiltonian type.

     87:    Collective

     89:    Input Parameters:
     90: +  A - matrix for the (0,0) block
     91: .  B - matrix for the (0,1) block, must be real symmetric or complex Hermitian
     92: -  C - matrix for the (1,0) block, must be real symmetric or complex Hermitian

     94:    Output Parameter:
     95: .  H  - the resulting matrix

     97:    Notes:
     98:    The resulting matrix has the block form $$H = \begin{bmatrix} A & B \\ C & -A^* \end{bmatrix},$$
     99:    where $B$ and $C$ are
    100:    assumed to be symmetric in the real case or Hermitian in the complex case. Note
    101:    that this function does not check this property, so if the matrices provided by
    102:    the user do not satisfy it, then the solver will not behave as expected.

    104:    The obtained matrix can be used as an input matrix to `EPS` eigensolvers via
    105:    `EPSSetOperators()` for the case that the problem type is `EPS_HAMILT`. Note that the
    106:    user cannot just build a matrix with the required structure, it must be done via
    107:    this function.

    109:    In the current implementation, `H` is a `MATNEST` matrix, where the (1,1) block is
    110:    a matrix of type `MATHERMITIANTRANSPOSEVIRTUAL` obtained from `A` and scaled by -1.

    112:    Level: intermediate

    114: .seealso: [](sec:structured), `MatCreateNest()`, `EPSSetOperators()`, `EPSSetProblemType()`, `MatCreateBSE()`, `MatCreateLREP()`
    115: @*/
    116: PetscErrorCode MatCreateHamiltonian(Mat A,Mat B,Mat C,Mat *H)
    117: {
    118:   PetscInt       Ma,Mb,Mc,Na,Nb,Nc,ma,mb,mc,na,nb,nc;
    119:   Mat            block[4] = { A, B, C, NULL };
    120:   SlepcMatStruct mctx;

    122:   PetscFunctionBegin;
    126:   PetscCheckSameTypeAndComm(A,1,B,2);
    127:   PetscCheckSameTypeAndComm(A,1,C,3);
    128:   PetscAssertPointer(H,4);

    130:   /* check sizes */
    131:   PetscCall(MatGetSize(A,&Ma,&Na));
    132:   PetscCall(MatGetLocalSize(A,&ma,&na));
    133:   PetscCall(MatGetSize(B,&Mb,&Nb));
    134:   PetscCall(MatGetLocalSize(B,&mb,&nb));
    135:   PetscCall(MatGetSize(C,&Mc,&Nc));
    136:   PetscCall(MatGetLocalSize(C,&mc,&nc));
    137:   PetscCheck(Mb==Ma && mb==ma,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Incompatible row dimensions");
    138:   PetscCheck(Nc==Na && nc==na,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Incompatible column dimensions");

    140:   /* (1,1) block */
    141:   PetscCall(MatCreateHermitianTranspose(A,&block[3]));
    142:   PetscCall(MatScale(block[3],-1.0));

    144:   /* create nest matrix and compose context */
    145:   PetscCall(MatCreateNest(PetscObjectComm((PetscObject)A),2,NULL,2,NULL,block,H));
    146:   PetscCall(MatDestroy(&block[3]));

    148:   PetscCall(PetscNew(&mctx));
    149:   mctx->cookie = SLEPC_MAT_STRUCT_HAMILT;
    150:   PetscCall(PetscObjectContainerCompose((PetscObject)*H,"SlepcMatStruct",mctx,PetscCtxDestroyDefault));
    151:   PetscFunctionReturn(PETSC_SUCCESS);
    152: }

    154: /*@
    155:    MatCreateLREP - Create a matrix that can be used to define a structured Linear
    156:    Response eigenvalue problem.

    158:    Collective

    160:    Input Parameters:
    161: +  AK  - matrix for the diagonal block ($A$) or the top block ($K$)
    162: .  BM  - matrix for the off-diagonal block ($B$) or the bottom block ($M$)
    163: -  red - whether the reduced form should be built

    165:    Output Parameter:
    166: .  H  - the resulting matrix

    168:    Notes:
    169:    The resulting matrix has the block form $$H = \begin{bmatrix} A & B \\ -B & -A \end{bmatrix},$$
    170:    where both $A$ and $B$ are assumed to be real symmetric. An alternative form,
    171:    called reduced form (`red`=`PETSC_TRUE`) is $$H = \begin{bmatrix} 0 & K \\ M & 0 \end{bmatrix},$$
    172:    where both $K$ and $M$ are also assumed to be real symmetric. Note that this function
    173:    does not check these properties, so if the matrices provided by the user do not
    174:    satisfy them, then the solver will not behave as expected.

    176:    The reduced form can be obtained from the original formulation by setting $K=A-B$
    177:    and $M=A+B$. However, $K$ and $M$ need not have such relation, in which case the
    178:    computed eigenvalues will be $\pm\lambda_i$, where $\lambda_i^2$ are the eigenvalues
    179:    of the _product eigenvalue problem_ $KM$ (or $MK$).

    181:    The obtained matrix can be used as an input matrix to `EPS` eigensolvers via
    182:    `EPSSetOperators()` for the case that the problem type is `EPS_LREP`. Note that the
    183:    user cannot just build a matrix with the required structure, it must be done via
    184:    this function.

    186:    In the current implementation, `H` is a `MATNEST` matrix, where `A` and `B` form the
    187:    top block row, while the bottom block row is composed of shell matrices scaled by -1.
    188:    Similarly, in the reduced formulation, the `MATNEST` has two zero blocks.

    190:    Internally, the eigensolver will work with the reduced form, and it will carry out
    191:    the conversion transparently if the first form is provided by the user.

    193:    Level: intermediate

    195: .seealso: [](sec:structured), `MatCreateNest()`, `EPSSetOperators()`, `EPSSetProblemType()`, `MatCreateBSE()`, `MatCreateHamiltonian()`
    196: @*/
    197: PetscErrorCode MatCreateLREP(Mat AK,Mat BM,PetscBool red,Mat *H)
    198: {
    199:   PetscInt       Ma,Mb,Na,Nb,ma,mb,na,nb;
    200:   SlepcMatStruct mbtx;

    202:   PetscFunctionBegin;
    205:   PetscCheckSameTypeAndComm(AK,1,BM,2);
    207:   PetscAssertPointer(H,4);

    209:   /* check sizes */
    210:   PetscCall(MatGetSize(AK,&Ma,&Na));
    211:   PetscCall(MatGetLocalSize(AK,&ma,&na));
    212:   PetscCall(MatGetSize(BM,&Mb,&Nb));
    213:   PetscCall(MatGetLocalSize(BM,&mb,&nb));
    214:   PetscCheck(Mb==Ma && mb==ma,PetscObjectComm((PetscObject)AK),PETSC_ERR_ARG_INCOMP,"Incompatible row dimensions");
    215:   PetscCheck(Nb==Na && nb==na,PetscObjectComm((PetscObject)AK),PETSC_ERR_ARG_INCOMP,"Incompatible column dimensions");

    217:   if (red) {  /* reduced form */

    219:     Mat block[4] = { NULL, AK, BM, NULL };

    221:     /* create nest matrix */
    222:     PetscCall(MatCreateNest(PetscObjectComm((PetscObject)AK),2,NULL,2,NULL,block,H));

    224:   } else {  /* original LREP form */

    226:     Mat block[4] = { AK, BM, NULL, NULL };
    227:     const PetscScalar scal[] = { -1.0 };

    229:     /* bottom block row, using MATCOMPOSITE */
    230:     PetscCall(MatCreate(PetscObjectComm((PetscObject)AK),&block[2]));
    231:     PetscCall(MatSetSizes(block[2],mb,nb,Mb,Nb));
    232:     PetscCall(MatSetType(block[2],MATCOMPOSITE));
    233:     PetscCall(MatCompositeAddMat(block[2],BM));
    234:     PetscCall(MatAssemblyBegin(block[2],MAT_FINAL_ASSEMBLY));
    235:     PetscCall(MatAssemblyEnd(block[2],MAT_FINAL_ASSEMBLY));
    236:     PetscCall(MatCompositeSetScalings(block[2],scal));

    238:     PetscCall(MatCreate(PetscObjectComm((PetscObject)AK),&block[3]));
    239:     PetscCall(MatSetSizes(block[3],ma,na,Ma,Na));
    240:     PetscCall(MatSetType(block[3],MATCOMPOSITE));
    241:     PetscCall(MatCompositeAddMat(block[3],AK));
    242:     PetscCall(MatAssemblyBegin(block[3],MAT_FINAL_ASSEMBLY));
    243:     PetscCall(MatAssemblyEnd(block[3],MAT_FINAL_ASSEMBLY));
    244:     PetscCall(MatCompositeSetScalings(block[3],scal));

    246:     /* create nest matrix */
    247:     PetscCall(MatCreateNest(PetscObjectComm((PetscObject)AK),2,NULL,2,NULL,block,H));
    248:     PetscCall(MatDestroy(&block[2]));
    249:     PetscCall(MatDestroy(&block[3]));

    251:   }

    253:   /* compose context */
    254:   PetscCall(PetscNew(&mbtx));
    255:   mbtx->cookie = SLEPC_MAT_STRUCT_LREP;
    256:   PetscCall(PetscObjectContainerCompose((PetscObject)*H,"SlepcMatStruct",mbtx,PetscCtxDestroyDefault));
    257:   PetscFunctionReturn(PETSC_SUCCESS);
    258: }