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             :    Implements the ST class for preconditioned eigenvalue methods
      12             : */
      13             : 
      14             : #include <slepc/private/stimpl.h>          /*I "slepcst.h" I*/
      15             : 
      16             : typedef struct {
      17             :   PetscBool ksphasmat;  /* the KSP must have the same matrix as PC */
      18             : } ST_PRECOND;
      19             : 
      20         294 : static PetscErrorCode STSetDefaultKSP_Precond(ST st)
      21             : {
      22         294 :   PC             pc;
      23         294 :   PCType         pctype;
      24         294 :   PetscBool      t0,t1;
      25             : 
      26         294 :   PetscFunctionBegin;
      27         294 :   PetscCall(KSPGetPC(st->ksp,&pc));
      28         294 :   PetscCall(PCGetType(pc,&pctype));
      29         294 :   if (!pctype && st->A && st->A[0]) {
      30         101 :     if (st->matmode == ST_MATMODE_SHELL) PetscCall(PCSetType(pc,PCJACOBI));
      31             :     else {
      32          99 :       PetscCall(MatHasOperation(st->A[0],MATOP_DUPLICATE,&t0));
      33          99 :       if (st->nmat>1) PetscCall(MatHasOperation(st->A[0],MATOP_AXPY,&t1));
      34          74 :       else t1 = PETSC_TRUE;
      35          99 :       PetscCall(PCSetType(pc,(t0 && t1)?PCBJACOBI:PCNONE));
      36             :     }
      37             :   }
      38         294 :   PetscCall(KSPSetErrorIfNotConverged(st->ksp,PETSC_FALSE));
      39         294 :   PetscFunctionReturn(PETSC_SUCCESS);
      40             : }
      41             : 
      42         204 : static PetscErrorCode STPostSolve_Precond(ST st)
      43             : {
      44         204 :   PetscFunctionBegin;
      45         204 :   if (st->matmode == ST_MATMODE_INPLACE && !(st->Pmat || (PetscAbsScalar(st->sigma)>=PETSC_MAX_REAL && st->nmat>1))) {
      46           1 :     if (st->nmat>1) PetscCall(MatAXPY(st->A[0],st->sigma,st->A[1],st->str));
      47           1 :     else PetscCall(MatShift(st->A[0],st->sigma));
      48           1 :     st->Astate[0] = ((PetscObject)st->A[0])->state;
      49           1 :     st->state   = ST_STATE_INITIAL;
      50           1 :     st->opready = PETSC_FALSE;
      51             :   }
      52         204 :   PetscFunctionReturn(PETSC_SUCCESS);
      53             : }
      54             : 
      55             : /*
      56             :    Operator (precond):
      57             :                Op        P         M
      58             :    if nmat=1:  ---       A-sI      NULL
      59             :    if nmat=2:  ---       A-sB      NULL
      60             : */
      61         179 : static PetscErrorCode STComputeOperator_Precond(ST st)
      62             : {
      63         179 :   PetscFunctionBegin;
      64             :   /* if the user did not set the shift, use the target value */
      65         179 :   if (!st->sigma_set) st->sigma = st->defsigma;
      66         179 :   st->M = NULL;
      67             : 
      68             :   /* build custom preconditioner from the split matrices */
      69         179 :   if (st->Psplit) {
      70           3 :     if (!(PetscAbsScalar(st->sigma) < PETSC_MAX_REAL) && st->nmat>1) {
      71           0 :       PetscCall(PetscObjectReference((PetscObject)st->Psplit[0]));
      72           0 :       PetscCall(MatDestroy(&st->Pmat));
      73           0 :       st->Pmat = st->Psplit[0];
      74           3 :     } else if (PetscAbsScalar(st->sigma)<PETSC_MAX_REAL) PetscCall(STMatMAXPY_Private(st,-st->sigma,0.0,0,NULL,PETSC_TRUE,PETSC_TRUE,&st->Pmat));
      75             :   }
      76             : 
      77             :   /* P = A-sigma*B */
      78         179 :   if (st->Pmat) {
      79          21 :     PetscCall(PetscObjectReference((PetscObject)st->Pmat));
      80          21 :     PetscCall(MatDestroy(&st->P));
      81          21 :     st->P = st->Pmat;
      82             :   } else {
      83         158 :     PetscCall(PetscObjectReference((PetscObject)st->A[1]));
      84         158 :     PetscCall(MatDestroy(&st->T[0]));
      85         158 :     st->T[0] = st->A[1];
      86         158 :     if (!(PetscAbsScalar(st->sigma) < PETSC_MAX_REAL) && st->nmat>1) {
      87           8 :       PetscCall(PetscObjectReference((PetscObject)st->T[0]));
      88           8 :       PetscCall(MatDestroy(&st->P));
      89           8 :       st->P = st->T[0];
      90         150 :     } else if (PetscAbsScalar(st->sigma)<PETSC_MAX_REAL) {
      91         113 :       PetscCall(STMatMAXPY_Private(st,-st->sigma,0.0,0,NULL,PetscNot(st->state==ST_STATE_UPDATED),PETSC_FALSE,&st->T[1]));
      92         113 :       PetscCall(PetscObjectReference((PetscObject)st->T[1]));
      93         113 :       PetscCall(MatDestroy(&st->P));
      94         113 :       st->P = st->T[1];
      95             :     }
      96             :   }
      97         179 :   PetscFunctionReturn(PETSC_SUCCESS);
      98             : }
      99             : 
     100         190 : static PetscErrorCode STSetUp_Precond(ST st)
     101             : {
     102         190 :   ST_PRECOND     *ctx = (ST_PRECOND*)st->data;
     103             : 
     104         190 :   PetscFunctionBegin;
     105         190 :   if (st->P) {
     106         223 :     PetscCall(ST_KSPSetOperators(st,ctx->ksphasmat?st->P:NULL,st->P));
     107             :     /* NOTE: we do not call KSPSetUp() here because some eigensolvers such as JD require a lazy setup */
     108             :   }
     109         190 :   PetscFunctionReturn(PETSC_SUCCESS);
     110             : }
     111             : 
     112           3 : static PetscErrorCode STSetShift_Precond(ST st,PetscScalar newshift)
     113             : {
     114           3 :   ST_PRECOND     *ctx = (ST_PRECOND*)st->data;
     115             : 
     116           3 :   PetscFunctionBegin;
     117           3 :   if (st->Psplit) { /* update custom preconditioner from the split matrices */
     118           0 :     if (PetscAbsScalar(st->sigma)<PETSC_MAX_REAL || st->nmat==1) PetscCall(STMatMAXPY_Private(st,-st->sigma,0.0,0,NULL,PETSC_FALSE,PETSC_TRUE,&st->Pmat));
     119             :   }
     120           3 :   if (st->transform && !st->Pmat) {
     121           3 :     PetscCall(STMatMAXPY_Private(st,-newshift,-st->sigma,0,NULL,PETSC_FALSE,PETSC_FALSE,&st->T[1]));
     122           3 :     if (st->P!=st->T[1]) {
     123           1 :       PetscCall(PetscObjectReference((PetscObject)st->T[1]));
     124           1 :       PetscCall(MatDestroy(&st->P));
     125           1 :       st->P = st->T[1];
     126             :     }
     127             :   }
     128           6 :   if (st->P) PetscCall(ST_KSPSetOperators(st,ctx->ksphasmat?st->P:NULL,st->P));
     129           3 :   PetscFunctionReturn(PETSC_SUCCESS);
     130             : }
     131             : 
     132         138 : static PetscErrorCode STPrecondSetKSPHasMat_Precond(ST st,PetscBool ksphasmat)
     133             : {
     134         138 :   ST_PRECOND *ctx = (ST_PRECOND*)st->data;
     135             : 
     136         138 :   PetscFunctionBegin;
     137         138 :   if (ctx->ksphasmat != ksphasmat) {
     138          82 :     ctx->ksphasmat = ksphasmat;
     139          82 :     st->state      = ST_STATE_INITIAL;
     140             :   }
     141         138 :   PetscFunctionReturn(PETSC_SUCCESS);
     142             : }
     143             : 
     144             : /*@
     145             :    STPrecondSetKSPHasMat - Sets a flag indicating that during STSetUp the coefficient
     146             :    matrix of the KSP linear solver (A) must be set to be the same matrix as the
     147             :    preconditioner (P).
     148             : 
     149             :    Collective
     150             : 
     151             :    Input Parameters:
     152             : +  st - the spectral transformation context
     153             : -  ksphasmat - the flag
     154             : 
     155             :    Notes:
     156             :    Often, the preconditioner matrix is used only in the PC object, but
     157             :    in some solvers this matrix must be provided also as the A-matrix in
     158             :    the KSP object.
     159             : 
     160             :    Level: developer
     161             : 
     162             : .seealso: STPrecondGetKSPHasMat(), STSetShift()
     163             : @*/
     164         140 : PetscErrorCode STPrecondSetKSPHasMat(ST st,PetscBool ksphasmat)
     165             : {
     166         140 :   PetscFunctionBegin;
     167         140 :   PetscValidHeaderSpecific(st,ST_CLASSID,1);
     168         420 :   PetscValidLogicalCollectiveBool(st,ksphasmat,2);
     169         140 :   PetscTryMethod(st,"STPrecondSetKSPHasMat_C",(ST,PetscBool),(st,ksphasmat));
     170         140 :   PetscFunctionReturn(PETSC_SUCCESS);
     171             : }
     172             : 
     173           2 : static PetscErrorCode STPrecondGetKSPHasMat_Precond(ST st,PetscBool *ksphasmat)
     174             : {
     175           2 :   ST_PRECOND *ctx = (ST_PRECOND*)st->data;
     176             : 
     177           2 :   PetscFunctionBegin;
     178           2 :   *ksphasmat = ctx->ksphasmat;
     179           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     180             : }
     181             : 
     182             : /*@
     183             :    STPrecondGetKSPHasMat - Returns the flag indicating if the coefficient
     184             :    matrix of the KSP linear system (A) is set to be the same matrix as the
     185             :    preconditioner (P).
     186             : 
     187             :    Not Collective
     188             : 
     189             :    Input Parameter:
     190             : .  st - the spectral transformation context
     191             : 
     192             :    Output Parameter:
     193             : .  ksphasmat - the flag
     194             : 
     195             :    Level: developer
     196             : 
     197             : .seealso: STPrecondSetKSPHasMat(), STSetShift()
     198             : @*/
     199           2 : PetscErrorCode STPrecondGetKSPHasMat(ST st,PetscBool *ksphasmat)
     200             : {
     201           2 :   PetscFunctionBegin;
     202           2 :   PetscValidHeaderSpecific(st,ST_CLASSID,1);
     203           2 :   PetscAssertPointer(ksphasmat,2);
     204           2 :   PetscUseMethod(st,"STPrecondGetKSPHasMat_C",(ST,PetscBool*),(st,ksphasmat));
     205           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     206             : }
     207             : 
     208         158 : static PetscErrorCode STDestroy_Precond(ST st)
     209             : {
     210         158 :   PetscFunctionBegin;
     211         158 :   PetscCall(PetscFree(st->data));
     212         158 :   PetscCall(PetscObjectComposeFunction((PetscObject)st,"STPrecondGetKSPHasMat_C",NULL));
     213         158 :   PetscCall(PetscObjectComposeFunction((PetscObject)st,"STPrecondSetKSPHasMat_C",NULL));
     214         158 :   PetscFunctionReturn(PETSC_SUCCESS);
     215             : }
     216             : 
     217         158 : SLEPC_EXTERN PetscErrorCode STCreate_Precond(ST st)
     218             : {
     219         158 :   ST_PRECOND     *ctx;
     220             : 
     221         158 :   PetscFunctionBegin;
     222         158 :   PetscCall(PetscNew(&ctx));
     223         158 :   st->data = (void*)ctx;
     224             : 
     225         158 :   st->usesksp = PETSC_TRUE;
     226             : 
     227         158 :   st->ops->apply           = STApply_Generic;
     228         158 :   st->ops->applymat        = STApplyMat_Generic;
     229         158 :   st->ops->applytrans      = STApplyTranspose_Generic;
     230         158 :   st->ops->applyhermtrans  = STApplyHermitianTranspose_Generic;
     231         158 :   st->ops->setshift        = STSetShift_Precond;
     232         158 :   st->ops->getbilinearform = STGetBilinearForm_Default;
     233         158 :   st->ops->setup           = STSetUp_Precond;
     234         158 :   st->ops->computeoperator = STComputeOperator_Precond;
     235         158 :   st->ops->postsolve       = STPostSolve_Precond;
     236         158 :   st->ops->destroy         = STDestroy_Precond;
     237         158 :   st->ops->setdefaultksp   = STSetDefaultKSP_Precond;
     238             : 
     239         158 :   PetscCall(PetscObjectComposeFunction((PetscObject)st,"STPrecondGetKSPHasMat_C",STPrecondGetKSPHasMat_Precond));
     240         158 :   PetscCall(PetscObjectComposeFunction((PetscObject)st,"STPrecondSetKSPHasMat_C",STPrecondSetKSPHasMat_Precond));
     241         158 :   PetscFunctionReturn(PETSC_SUCCESS);
     242             : }