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 { 18 PC_TFS *tfs = (PC_TFS*)pc->data; 19 20 PetscFunctionBegin; 21 /* free the XXT datastructures */ 22 if (tfs->xxt) { 23 PetscCall(XXT_free(tfs->xxt)); 24 } 25 if (tfs->xyt) { 26 PetscCall(XYT_free(tfs->xyt)); 27 } 28 PetscCall(VecDestroy(&tfs->b)); 29 PetscCall(VecDestroy(&tfs->xd)); 30 PetscCall(VecDestroy(&tfs->xo)); 31 PetscCall(PetscFree(pc->data)); 32 PetscFunctionReturn(0); 33 } 34 35 static PetscErrorCode PCApply_TFS_XXT(PC pc,Vec x,Vec y) 36 { 37 PC_TFS *tfs = (PC_TFS*)pc->data; 38 PetscScalar *yy; 39 const PetscScalar *xx; 40 41 PetscFunctionBegin; 42 PetscCall(VecGetArrayRead(x,&xx)); 43 PetscCall(VecGetArray(y,&yy)); 44 PetscCall(XXT_solve(tfs->xxt,yy,(PetscScalar*)xx)); 45 PetscCall(VecRestoreArrayRead(x,&xx)); 46 PetscCall(VecRestoreArray(y,&yy)); 47 PetscFunctionReturn(0); 48 } 49 50 static PetscErrorCode PCApply_TFS_XYT(PC pc,Vec x,Vec y) 51 { 52 PC_TFS *tfs = (PC_TFS*)pc->data; 53 PetscScalar *yy; 54 const PetscScalar *xx; 55 56 PetscFunctionBegin; 57 PetscCall(VecGetArrayRead(x,&xx)); 58 PetscCall(VecGetArray(y,&yy)); 59 PetscCall(XYT_solve(tfs->xyt,yy,(PetscScalar*)xx)); 60 PetscCall(VecRestoreArrayRead(x,&xx)); 61 PetscCall(VecRestoreArray(y,&yy)); 62 PetscFunctionReturn(0); 63 } 64 65 static PetscErrorCode PCTFSLocalMult_TFS(PC pc,PetscScalar *xin,PetscScalar *xout) 66 { 67 PC_TFS *tfs = (PC_TFS*)pc->data; 68 Mat A = pc->pmat; 69 Mat_MPIAIJ *a = (Mat_MPIAIJ*)A->data; 70 71 PetscFunctionBegin; 72 PetscCall(VecPlaceArray(tfs->b,xout)); 73 PetscCall(VecPlaceArray(tfs->xd,xin)); 74 PetscCall(VecPlaceArray(tfs->xo,xin+tfs->nd)); 75 PetscCall(MatMult(a->A,tfs->xd,tfs->b)); 76 PetscCall(MatMultAdd(a->B,tfs->xo,tfs->b,tfs->b)); 77 PetscCall(VecResetArray(tfs->b)); 78 PetscCall(VecResetArray(tfs->xd)); 79 PetscCall(VecResetArray(tfs->xo)); 80 PetscFunctionReturn(0); 81 } 82 83 static PetscErrorCode PCSetUp_TFS(PC pc) 84 { 85 PC_TFS *tfs = (PC_TFS*)pc->data; 86 Mat A = pc->pmat; 87 Mat_MPIAIJ *a = (Mat_MPIAIJ*)A->data; 88 PetscInt *localtoglobal,ncol,i; 89 PetscBool ismpiaij; 90 91 /* 92 PetscBool issymmetric; 93 Petsc Real tol = 0.0; 94 */ 95 96 PetscFunctionBegin; 97 PetscCheck(A->cmap->N == A->rmap->N,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_SIZ,"matrix must be square"); 98 PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat,MATMPIAIJ,&ismpiaij)); 99 PetscCheck(ismpiaij,PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Currently only supports MPIAIJ matrices"); 100 101 /* generate the local to global mapping */ 102 ncol = a->A->cmap->n + a->B->cmap->n; 103 PetscCall(PetscMalloc1(ncol,&localtoglobal)); 104 for (i=0; i<a->A->cmap->n; i++) localtoglobal[i] = A->cmap->rstart + i + 1; 105 for (i=0; i<a->B->cmap->n; i++) localtoglobal[i+a->A->cmap->n] = a->garray[i] + 1; 106 107 /* generate the vectors needed for the local solves */ 108 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF,1,a->A->rmap->n,NULL,&tfs->b)); 109 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF,1,a->A->cmap->n,NULL,&tfs->xd)); 110 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF,1,a->B->cmap->n,NULL,&tfs->xo)); 111 tfs->nd = a->A->cmap->n; 112 113 /* ierr = MatIsSymmetric(A,tol,&issymmetric); */ 114 /* if (issymmetric) { */ 115 PetscCall(PetscBarrier((PetscObject)pc)); 116 if (A->symmetric) { 117 tfs->xxt = XXT_new(); 118 PetscCall(XXT_factor(tfs->xxt,localtoglobal,A->rmap->n,ncol,(PetscErrorCode (*)(void*,PetscScalar*,PetscScalar*))PCTFSLocalMult_TFS,pc)); 119 pc->ops->apply = PCApply_TFS_XXT; 120 } else { 121 tfs->xyt = XYT_new(); 122 PetscCall(XYT_factor(tfs->xyt,localtoglobal,A->rmap->n,ncol,(PetscErrorCode (*)(void*,PetscScalar*,PetscScalar*))PCTFSLocalMult_TFS,pc)); 123 pc->ops->apply = PCApply_TFS_XYT; 124 } 125 126 PetscCall(PetscFree(localtoglobal)); 127 PetscFunctionReturn(0); 128 } 129 130 static PetscErrorCode PCSetFromOptions_TFS(PetscOptionItems *PetscOptionsObject,PC pc) 131 { 132 PetscFunctionBegin; 133 PetscFunctionReturn(0); 134 } 135 static PetscErrorCode PCView_TFS(PC pc,PetscViewer viewer) 136 { 137 PetscFunctionBegin; 138 PetscFunctionReturn(0); 139 } 140 141 /*MC 142 PCTFS - A parallel direct solver intended for problems with very few unknowns (like the 143 coarse grid in multigrid). Performs a Cholesky or LU factorization of a matrix defined by 144 its local matrix vector product. 145 146 Implemented by Henry M. Tufo III and Paul Fischer originally for Nek5000 and called XXT or XYT 147 148 Level: beginner 149 150 Notes: 151 Only implemented for the MPIAIJ matrices 152 153 Only works on a solver object that lives on all of PETSC_COMM_WORLD! 154 155 Only works for real numbers (is not built if PetscScalar is complex) 156 157 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC 158 M*/ 159 PETSC_EXTERN PetscErrorCode PCCreate_TFS(PC pc) 160 { 161 PC_TFS *tfs; 162 PetscMPIInt cmp; 163 164 PetscFunctionBegin; 165 PetscCallMPI(MPI_Comm_compare(PETSC_COMM_WORLD,PetscObjectComm((PetscObject)pc),&cmp)); 166 PetscCheckFalse(cmp != MPI_IDENT && cmp != MPI_CONGRUENT,PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"TFS only works with PETSC_COMM_WORLD objects"); 167 PetscCall(PetscNewLog(pc,&tfs)); 168 169 tfs->xxt = NULL; 170 tfs->xyt = NULL; 171 tfs->b = NULL; 172 tfs->xd = NULL; 173 tfs->xo = NULL; 174 tfs->nd = 0; 175 176 pc->ops->apply = NULL; 177 pc->ops->applytranspose = NULL; 178 pc->ops->setup = PCSetUp_TFS; 179 pc->ops->destroy = PCDestroy_TFS; 180 pc->ops->setfromoptions = PCSetFromOptions_TFS; 181 pc->ops->view = PCView_TFS; 182 pc->ops->applyrichardson = NULL; 183 pc->ops->applysymmetricleft = NULL; 184 pc->ops->applysymmetricright = NULL; 185 pc->data = (void*)tfs; 186 PetscFunctionReturn(0); 187 } 188