Line data    Source code

       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
       5             : 
       6             :    This file is part of SLEPc.
       7             :    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
       8             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       9             : */
      10             : 
      11             : #include <slepc/private/slepcimpl.h>            /*I "slepcsys.h" I*/
      12             : 
      13             : /*@
      14             :    MatCreateBSE - Create a matrix that can be used to define a structured eigenvalue
      15             :    problem of type BSE (Bethe-Salpeter Equation).
      16             : 
      17             :    Collective
      18             : 
      19             :    Input Parameters:
      20             : +  R - matrix for the diagonal block (resonant)
      21             : -  C - matrix for the off-diagonal block (coupling)
      22             : 
      23             :    Output Parameter:
      24             : .  H  - the resulting matrix
      25             : 
      26             :    Notes:
      27             :    The resulting matrix has the block form H = [ R C; -C^H -R^T ], where R is assumed
      28             :    to be (complex) Hermitian and C complex symmetric. Note that this function does
      29             :    not check these properties, so if the matrices provided by the user do not satisfy
      30             :    them, then the solver will not behave as expected.
      31             : 
      32             :    The obtained matrix can be used as an input matrix to EPS eigensolvers via
      33             :    EPSSetOperators() for the case that the problem type is EPS_BSE. Note that the user
      34             :    cannot just build a matrix with the required structure, it must be done via this
      35             :    function.
      36             : 
      37             :    In the current implementation, H is a MATNEST matrix, where R and C form the top
      38             :    block row, while the bottom block row is composed of matrices of type
      39             :    MATTRANSPOSEVIRTUAL and MATHERMITIANTRANSPOSEVIRTUAL scaled by -1.
      40             : 
      41             :    Level: intermediate
      42             : 
      43             : .seealso: MatCreateNest(), EPSSetOperators(), EPSSetProblemType()
      44             : @*/
      45          20 : PetscErrorCode MatCreateBSE(Mat R,Mat C,Mat *H)
      46             : {
      47          20 :   PetscInt       Mr,Mc,Nr,Nc,mr,mc,nr,nc;
      48          20 :   Mat            block[4] = { R, C, NULL, NULL };
      49          20 :   SlepcMatStruct mctx;
      50             : 
      51          20 :   PetscFunctionBegin;
      52          20 :   PetscValidHeaderSpecific(R,MAT_CLASSID,1);
      53          20 :   PetscValidHeaderSpecific(C,MAT_CLASSID,2);
      54          20 :   PetscCheckSameTypeAndComm(R,1,C,2);
      55          20 :   PetscAssertPointer(H,3);
      56             : 
      57             :   /* check sizes */
      58          20 :   PetscCall(MatGetSize(R,&Mr,&Nr));
      59          20 :   PetscCall(MatGetLocalSize(R,&mr,&nr));
      60          20 :   PetscCall(MatGetSize(C,&Mc,&Nc));
      61          20 :   PetscCall(MatGetLocalSize(C,&mc,&nc));
      62          20 :   PetscCheck(Mc==Mr && mc==mr,PetscObjectComm((PetscObject)R),PETSC_ERR_ARG_INCOMP,"Incompatible row dimensions");
      63          20 :   PetscCheck(Nc==Nr && nc==nr,PetscObjectComm((PetscObject)R),PETSC_ERR_ARG_INCOMP,"Incompatible column dimensions");
      64             : 
      65             :   /* bottom block row */
      66          20 :   PetscCall(MatCreateHermitianTranspose(C,&block[2]));
      67          20 :   PetscCall(MatScale(block[2],-1.0));
      68          20 :   PetscCall(MatCreateTranspose(R,&block[3]));
      69          20 :   PetscCall(MatScale(block[3],-1.0));
      70             : 
      71             :   /* create nest matrix and compose context */
      72          20 :   PetscCall(MatCreateNest(PetscObjectComm((PetscObject)R),2,NULL,2,NULL,block,H));
      73          20 :   PetscCall(MatDestroy(&block[2]));
      74          20 :   PetscCall(MatDestroy(&block[3]));
      75             : 
      76          20 :   PetscCall(PetscNew(&mctx));
      77          20 :   mctx->cookie = SLEPC_MAT_STRUCT_BSE;
      78          20 :   PetscCall(PetscObjectContainerCompose((PetscObject)*H,"SlepcMatStruct",mctx,PetscCtxDestroyDefault));
      79          20 :   PetscFunctionReturn(PETSC_SUCCESS);
      80             : }