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