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 Implemented by Henry M. Tufo III and Paul Fischer originally for Nek5000 and called XXT or XYT 136 137 Level: beginner 138 139 Notes: 140 Only implemented for the MPIAIJ matrices 141 142 Only works on a solver object that lives on all of PETSC_COMM_WORLD! 143 144 Only works for real numbers (is not built if PetscScalar is complex) 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(PetscNewLog(pc, &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