1 /* 2 Provides an interface to the Tufo-Fischer parallel direct solver 3 */ 4 5 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/ 6 #include <../src/mat/impls/aij/mpi/mpiaij.h> 7 #include <../src/ksp/pc/impls/tfs/tfs.h> 8 9 typedef struct { 10 xxt_ADT xxt; 11 xyt_ADT xyt; 12 Vec b, xd, xo; 13 PetscInt nd; 14 } PC_TFS; 15 16 PetscErrorCode PCDestroy_TFS(PC pc) { 17 PC_TFS *tfs = (PC_TFS *)pc->data; 18 19 PetscFunctionBegin; 20 /* free the XXT datastructures */ 21 if (tfs->xxt) PetscCall(XXT_free(tfs->xxt)); 22 if (tfs->xyt) PetscCall(XYT_free(tfs->xyt)); 23 PetscCall(VecDestroy(&tfs->b)); 24 PetscCall(VecDestroy(&tfs->xd)); 25 PetscCall(VecDestroy(&tfs->xo)); 26 PetscCall(PetscFree(pc->data)); 27 PetscFunctionReturn(0); 28 } 29 30 static PetscErrorCode PCApply_TFS_XXT(PC pc, Vec x, Vec y) { 31 PC_TFS *tfs = (PC_TFS *)pc->data; 32 PetscScalar *yy; 33 const PetscScalar *xx; 34 35 PetscFunctionBegin; 36 PetscCall(VecGetArrayRead(x, &xx)); 37 PetscCall(VecGetArray(y, &yy)); 38 PetscCall(XXT_solve(tfs->xxt, yy, (PetscScalar *)xx)); 39 PetscCall(VecRestoreArrayRead(x, &xx)); 40 PetscCall(VecRestoreArray(y, &yy)); 41 PetscFunctionReturn(0); 42 } 43 44 static PetscErrorCode PCApply_TFS_XYT(PC pc, Vec x, Vec y) { 45 PC_TFS *tfs = (PC_TFS *)pc->data; 46 PetscScalar *yy; 47 const PetscScalar *xx; 48 49 PetscFunctionBegin; 50 PetscCall(VecGetArrayRead(x, &xx)); 51 PetscCall(VecGetArray(y, &yy)); 52 PetscCall(XYT_solve(tfs->xyt, yy, (PetscScalar *)xx)); 53 PetscCall(VecRestoreArrayRead(x, &xx)); 54 PetscCall(VecRestoreArray(y, &yy)); 55 PetscFunctionReturn(0); 56 } 57 58 static PetscErrorCode PCTFSLocalMult_TFS(PC pc, PetscScalar *xin, PetscScalar *xout) { 59 PC_TFS *tfs = (PC_TFS *)pc->data; 60 Mat A = pc->pmat; 61 Mat_MPIAIJ *a = (Mat_MPIAIJ *)A->data; 62 63 PetscFunctionBegin; 64 PetscCall(VecPlaceArray(tfs->b, xout)); 65 PetscCall(VecPlaceArray(tfs->xd, xin)); 66 PetscCall(VecPlaceArray(tfs->xo, xin + tfs->nd)); 67 PetscCall(MatMult(a->A, tfs->xd, tfs->b)); 68 PetscCall(MatMultAdd(a->B, tfs->xo, tfs->b, tfs->b)); 69 PetscCall(VecResetArray(tfs->b)); 70 PetscCall(VecResetArray(tfs->xd)); 71 PetscCall(VecResetArray(tfs->xo)); 72 PetscFunctionReturn(0); 73 } 74 75 static PetscErrorCode PCSetUp_TFS(PC pc) { 76 PC_TFS *tfs = (PC_TFS *)pc->data; 77 Mat A = pc->pmat; 78 Mat_MPIAIJ *a = (Mat_MPIAIJ *)A->data; 79 PetscInt *localtoglobal, ncol, i; 80 PetscBool ismpiaij; 81 82 /* 83 PetscBool issymmetric; 84 Petsc Real tol = 0.0; 85 */ 86 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"); 91 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; 97 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; 103 104 /* ierr = MatIsSymmetric(A,tol,&issymmetric); */ 105 /* if (issymmetric) { */ 106 PetscCall(PetscBarrier((PetscObject)pc)); 107 if (A->symmetric == PETSC_BOOL3_TRUE) { 108 tfs->xxt = XXT_new(); 109 PetscCall(XXT_factor(tfs->xxt, localtoglobal, A->rmap->n, ncol, (PetscErrorCode(*)(void *, PetscScalar *, PetscScalar *))PCTFSLocalMult_TFS, pc)); 110 pc->ops->apply = PCApply_TFS_XXT; 111 } else { 112 tfs->xyt = XYT_new(); 113 PetscCall(XYT_factor(tfs->xyt, localtoglobal, A->rmap->n, ncol, (PetscErrorCode(*)(void *, PetscScalar *, PetscScalar *))PCTFSLocalMult_TFS, pc)); 114 pc->ops->apply = PCApply_TFS_XYT; 115 } 116 117 PetscCall(PetscFree(localtoglobal)); 118 PetscFunctionReturn(0); 119 } 120 121 static PetscErrorCode PCSetFromOptions_TFS(PC pc, PetscOptionItems *PetscOptionsObject) { 122 PetscFunctionBegin; 123 PetscFunctionReturn(0); 124 } 125 static PetscErrorCode PCView_TFS(PC pc, PetscViewer viewer) { 126 PetscFunctionBegin; 127 PetscFunctionReturn(0); 128 } 129 130 /*MC 131 PCTFS - A parallel direct solver intended for problems with very few unknowns (like the 132 coarse grid in multigrid). Performs a Cholesky or LU factorization of a matrix defined by 133 its local matrix vector product. 134 135 Level: beginner 136 137 Notes: 138 Only implemented for the `MATMPIAIJ` matrices 139 140 Only works on a solver object that lives on all of `PETSC_COMM_WORLD`! 141 142 Only works for real numbers (is not built if `PetscScalar` is complex) 143 144 Implemented by Henry M. Tufo III and Paul Fischer originally for Nek5000 and called XXT or XYT 145 146 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC` 147 M*/ 148 PETSC_EXTERN PetscErrorCode PCCreate_TFS(PC pc) { 149 PC_TFS *tfs; 150 PetscMPIInt cmp; 151 152 PetscFunctionBegin; 153 PetscCallMPI(MPI_Comm_compare(PETSC_COMM_WORLD, PetscObjectComm((PetscObject)pc), &cmp)); 154 PetscCheck(cmp == MPI_IDENT || cmp == MPI_CONGRUENT, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "TFS only works with PETSC_COMM_WORLD objects"); 155 PetscCall(PetscNew(&tfs)); 156 157 tfs->xxt = NULL; 158 tfs->xyt = NULL; 159 tfs->b = NULL; 160 tfs->xd = NULL; 161 tfs->xo = NULL; 162 tfs->nd = 0; 163 164 pc->ops->apply = NULL; 165 pc->ops->applytranspose = NULL; 166 pc->ops->setup = PCSetUp_TFS; 167 pc->ops->destroy = PCDestroy_TFS; 168 pc->ops->setfromoptions = PCSetFromOptions_TFS; 169 pc->ops->view = PCView_TFS; 170 pc->ops->applyrichardson = NULL; 171 pc->ops->applysymmetricleft = NULL; 172 pc->ops->applysymmetricright = NULL; 173 pc->data = (void *)tfs; 174 PetscFunctionReturn(0); 175 } 176