Actual source code: factorschur.c

      1: #include <petsc/private/matimpl.h>
      2: #include <../src/mat/impls/dense/seq/dense.h>

      4: PETSC_INTERN PetscErrorCode MatFactorSetUpInPlaceSchur_Private(Mat F)
      5: {
      6:   Mat              St, S = F->schur;
      7:   MatFactorInfo    info;
      8:   PetscErrorCode   ierr;

     11:   MatSetUnfactored(S);
     12:   MatGetFactor(S,S->solvertype ? S->solvertype : MATSOLVERPETSC,F->factortype,&St);
     13:   if (St->factortype == MAT_FACTOR_CHOLESKY) { /* LDL^t regarded as Cholesky */
     14:     MatCholeskyFactorSymbolic(St,S,NULL,&info);
     15:   } else {
     16:     MatLUFactorSymbolic(St,S,NULL,NULL,&info);
     17:   }
     18:   S->ops->solve             = St->ops->solve;
     19:   S->ops->matsolve          = St->ops->matsolve;
     20:   S->ops->solvetranspose    = St->ops->solvetranspose;
     21:   S->ops->matsolvetranspose = St->ops->matsolvetranspose;
     22:   S->ops->solveadd          = St->ops->solveadd;
     23:   S->ops->solvetransposeadd = St->ops->solvetransposeadd;

     25:   MatDestroy(&St);
     26:   return(0);
     27: }

     29: PETSC_INTERN PetscErrorCode MatFactorUpdateSchurStatus_Private(Mat F)
     30: {
     31:   Mat            S = F->schur;

     35:   switch(F->schur_status) {
     36:   case MAT_FACTOR_SCHUR_UNFACTORED:
     37:   case MAT_FACTOR_SCHUR_INVERTED:
     38:     if (S) {
     39:       S->ops->solve             = NULL;
     40:       S->ops->matsolve          = NULL;
     41:       S->ops->solvetranspose    = NULL;
     42:       S->ops->matsolvetranspose = NULL;
     43:       S->ops->solveadd          = NULL;
     44:       S->ops->solvetransposeadd = NULL;
     45:       S->factortype             = MAT_FACTOR_NONE;
     46:       PetscFree(S->solvertype);
     47:     }
     48:     break;
     49:   case MAT_FACTOR_SCHUR_FACTORED:
     50:     break;
     51:   default:
     52:     SETERRQ1(PetscObjectComm((PetscObject)F),PETSC_ERR_SUP,"Unhandled MatFactorSchurStatus %D",F->schur_status);
     53:   }
     54:   return(0);
     55: }

     57: /* Schur status updated in the interface */
     58: PETSC_INTERN PetscErrorCode MatFactorFactorizeSchurComplement_Private(Mat F)
     59: {
     60:   MatFactorInfo  info;

     64:   PetscLogEventBegin(MAT_FactorFactS,F,0,0,0);
     65:   if (F->factortype == MAT_FACTOR_CHOLESKY) { /* LDL^t regarded as Cholesky */
     66:     MatCholeskyFactor(F->schur,NULL,&info);
     67:   } else {
     68:     MatLUFactor(F->schur,NULL,NULL,&info);
     69:   }
     70:   PetscLogEventEnd(MAT_FactorFactS,F,0,0,0);
     71:   return(0);
     72: }

     74: /* Schur status updated in the interface */
     75: PETSC_INTERN PetscErrorCode MatFactorInvertSchurComplement_Private(Mat F)
     76: {
     77:   Mat S = F->schur;

     80:   if (S) {
     81:     PetscMPIInt    size;
     82:     PetscBool      isdense,isdensecuda;

     85:     MPI_Comm_size(PetscObjectComm((PetscObject)S),&size);
     86:     if (size > 1) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_SUP,"Not yet implemented");
     87:     PetscObjectTypeCompare((PetscObject)S,MATSEQDENSE,&isdense);
     88:     PetscObjectTypeCompare((PetscObject)S,MATSEQDENSECUDA,&isdensecuda);
     89:     PetscLogEventBegin(MAT_FactorInvS,F,0,0,0);
     90:     if (isdense) {
     91:       MatSeqDenseInvertFactors_Private(S);
     92: #if defined(PETSC_HAVE_CUDA)
     93:     } else if (isdensecuda) {
     94:       MatSeqDenseCUDAInvertFactors_Private(S);
     95: #endif
     96:     } else SETERRQ1(PetscObjectComm((PetscObject)S),PETSC_ERR_SUP,"Not implemented for type %s",((PetscObject)S)->type_name);
     97:     PetscLogEventEnd(MAT_FactorInvS,F,0,0,0);
     98:   }
     99:   return(0);
    100: }