Actual source code: tfs.c

      1: /*
      2:         Provides an interface to the Tufo-Fischer parallel direct solver
      3: */

      5: #include <petsc/private/pcimpl.h>
      6: #include <../src/mat/impls/aij/mpi/mpiaij.h>
      7: #include <../src/ksp/pc/impls/tfs/tfs.h>

      9: typedef struct {
     10:   xxt_ADT  xxt;
     11:   xyt_ADT  xyt;
     12:   Vec      b, xd, xo;
     13:   PetscInt nd;
     14: } PC_TFS;

     16: static PetscErrorCode PCDestroy_TFS(PC pc)
     17: {
     18:   PC_TFS *tfs = (PC_TFS *)pc->data;

     20:   PetscFunctionBegin;
     21:   /* free the XXT datastructures */
     22:   if (tfs->xxt) PetscCall(XXT_free(tfs->xxt));
     23:   if (tfs->xyt) PetscCall(XYT_free(tfs->xyt));
     24:   PetscCall(VecDestroy(&tfs->b));
     25:   PetscCall(VecDestroy(&tfs->xd));
     26:   PetscCall(VecDestroy(&tfs->xo));
     27:   PetscCall(PetscFree(pc->data));
     28:   PetscFunctionReturn(PETSC_SUCCESS);
     29: }

     31: static PetscErrorCode PCApply_TFS_XXT(PC pc, Vec x, Vec y)
     32: {
     33:   PC_TFS            *tfs = (PC_TFS *)pc->data;
     34:   PetscScalar       *yy;
     35:   const PetscScalar *xx;

     37:   PetscFunctionBegin;
     38:   PetscCall(VecGetArrayRead(x, &xx));
     39:   PetscCall(VecGetArray(y, &yy));
     40:   PetscCall(XXT_solve(tfs->xxt, yy, (PetscScalar *)xx));
     41:   PetscCall(VecRestoreArrayRead(x, &xx));
     42:   PetscCall(VecRestoreArray(y, &yy));
     43:   PetscFunctionReturn(PETSC_SUCCESS);
     44: }

     46: static PetscErrorCode PCApply_TFS_XYT(PC pc, Vec x, Vec y)
     47: {
     48:   PC_TFS            *tfs = (PC_TFS *)pc->data;
     49:   PetscScalar       *yy;
     50:   const PetscScalar *xx;

     52:   PetscFunctionBegin;
     53:   PetscCall(VecGetArrayRead(x, &xx));
     54:   PetscCall(VecGetArray(y, &yy));
     55:   PetscCall(XYT_solve(tfs->xyt, yy, (PetscScalar *)xx));
     56:   PetscCall(VecRestoreArrayRead(x, &xx));
     57:   PetscCall(VecRestoreArray(y, &yy));
     58:   PetscFunctionReturn(PETSC_SUCCESS);
     59: }

     61: static PetscErrorCode PCTFSLocalMult_TFS(PC pc, PetscScalar *xin, PetscScalar *xout)
     62: {
     63:   PC_TFS     *tfs = (PC_TFS *)pc->data;
     64:   Mat         A   = pc->pmat;
     65:   Mat_MPIAIJ *a   = (Mat_MPIAIJ *)A->data;

     67:   PetscFunctionBegin;
     68:   PetscCall(VecPlaceArray(tfs->b, xout));
     69:   PetscCall(VecPlaceArray(tfs->xd, xin));
     70:   PetscCall(VecPlaceArray(tfs->xo, xin + tfs->nd));
     71:   PetscCall(MatMult(a->A, tfs->xd, tfs->b));
     72:   PetscCall(MatMultAdd(a->B, tfs->xo, tfs->b, tfs->b));
     73:   PetscCall(VecResetArray(tfs->b));
     74:   PetscCall(VecResetArray(tfs->xd));
     75:   PetscCall(VecResetArray(tfs->xo));
     76:   PetscFunctionReturn(PETSC_SUCCESS);
     77: }

     79: static PetscErrorCode PCSetUp_TFS(PC pc)
     80: {
     81:   PC_TFS     *tfs = (PC_TFS *)pc->data;
     82:   Mat         A   = pc->pmat;
     83:   Mat_MPIAIJ *a   = (Mat_MPIAIJ *)A->data;
     84:   PetscInt   *localtoglobal, ncol, i;
     85:   PetscBool   ismpiaij;

     87:   PetscFunctionBegin;
     88:   PetscCheck(A->cmap->N == A->rmap->N, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_SIZ, "matrix must be square");
     89:   PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATMPIAIJ, &ismpiaij));
     90:   PetscCheck(ismpiaij, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Currently only supports MPIAIJ matrices");

     92:   /* generate the local to global mapping */
     93:   ncol = a->A->cmap->n + a->B->cmap->n;
     94:   PetscCall(PetscMalloc1(ncol, &localtoglobal));
     95:   for (i = 0; i < a->A->cmap->n; i++) localtoglobal[i] = A->cmap->rstart + i + 1;
     96:   for (i = 0; i < a->B->cmap->n; i++) localtoglobal[i + a->A->cmap->n] = a->garray[i] + 1;

     98:   /* generate the vectors needed for the local solves */
     99:   PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, a->A->rmap->n, NULL, &tfs->b));
    100:   PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, a->A->cmap->n, NULL, &tfs->xd));
    101:   PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, a->B->cmap->n, NULL, &tfs->xo));
    102:   tfs->nd = a->A->cmap->n;

    104:   PetscCall(PetscBarrier((PetscObject)pc));
    105:   if (A->symmetric == PETSC_BOOL3_TRUE) {
    106:     tfs->xxt = XXT_new();
    107:     PetscCall(XXT_factor(tfs->xxt, localtoglobal, A->rmap->n, ncol, (PetscErrorCode (*)(void *, PetscScalar *, PetscScalar *))PCTFSLocalMult_TFS, pc));
    108:     pc->ops->apply = PCApply_TFS_XXT;
    109:   } else {
    110:     tfs->xyt = XYT_new();
    111:     PetscCall(XYT_factor(tfs->xyt, localtoglobal, A->rmap->n, ncol, (PetscErrorCode (*)(void *, PetscScalar *, PetscScalar *))PCTFSLocalMult_TFS, pc));
    112:     pc->ops->apply = PCApply_TFS_XYT;
    113:   }

    115:   PetscCall(PetscFree(localtoglobal));
    116:   PetscFunctionReturn(PETSC_SUCCESS);
    117: }

    119: /*MC
    120:      PCTFS - A parallel direct solver intended for problems with very few unknowns (like the
    121:          coarse grid in multigrid). Performs a Cholesky or LU factorization of a matrix defined by
    122:          its local matrix-vector product.

    124:    Level: beginner

    126:    Notes:
    127:     Only implemented for the `MATMPIAIJ` matrices

    129:     Only works on a solver object that lives on all of `PETSC_COMM_WORLD`!

    131:     Only works for real numbers (is not built if `PetscScalar` is complex)

    133:    Implemented by  Henry M. Tufo III and Paul Fischer originally for Nek5000 and called XXT or XYT

    135: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`
    136: M*/
    137: PETSC_EXTERN PetscErrorCode PCCreate_TFS(PC pc)
    138: {
    139:   PC_TFS     *tfs;
    140:   PetscMPIInt cmp;

    142:   PetscFunctionBegin;
    143:   PetscCallMPI(MPI_Comm_compare(PETSC_COMM_WORLD, PetscObjectComm((PetscObject)pc), &cmp));
    144:   PetscCheck(cmp == MPI_IDENT || cmp == MPI_CONGRUENT, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "TFS only works with PETSC_COMM_WORLD objects");
    145:   PetscCall(PetscNew(&tfs));

    147:   tfs->xxt = NULL;
    148:   tfs->xyt = NULL;
    149:   tfs->b   = NULL;
    150:   tfs->xd  = NULL;
    151:   tfs->xo  = NULL;
    152:   tfs->nd  = 0;

    154:   pc->ops->apply               = NULL;
    155:   pc->ops->applytranspose      = NULL;
    156:   pc->ops->setup               = PCSetUp_TFS;
    157:   pc->ops->destroy             = PCDestroy_TFS;
    158:   pc->ops->applyrichardson     = NULL;
    159:   pc->ops->applysymmetricleft  = NULL;
    160:   pc->ops->applysymmetricright = NULL;
    161:   pc->data                     = (void *)tfs;
    162:   PetscFunctionReturn(PETSC_SUCCESS);
    163: }