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()`
 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 comples 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()`
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: }