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 119 /* ierr = MatIsSymmetric(A,tol,&issymmetric); */ 120 /* if (issymmetric) { */ 121 ierr = PetscBarrier((PetscObject)pc);CHKERRQ(ierr); 122 if (A->symmetric) { 123 tfs->xxt = XXT_new(); 124 ierr = XXT_factor(tfs->xxt,localtoglobal,A->rmap->n,ncol,(PetscErrorCode (*)(void*,PetscScalar*,PetscScalar*))PCTFSLocalMult_TFS,pc);CHKERRQ(ierr); 125 pc->ops->apply = PCApply_TFS_XXT; 126 } else { 127 tfs->xyt = XYT_new(); 128 ierr = XYT_factor(tfs->xyt,localtoglobal,A->rmap->n,ncol,(PetscErrorCode (*)(void*,PetscScalar*,PetscScalar*))PCTFSLocalMult_TFS,pc);CHKERRQ(ierr); 129 pc->ops->apply = PCApply_TFS_XYT; 130 } 131 132 ierr = PetscFree(localtoglobal);CHKERRQ(ierr); 133 PetscFunctionReturn(0); 134 } 135 136 static PetscErrorCode PCSetFromOptions_TFS(PetscOptionItems *PetscOptionsObject,PC pc) 137 { 138 PetscFunctionBegin; 139 PetscFunctionReturn(0); 140 } 141 static PetscErrorCode PCView_TFS(PC pc,PetscViewer viewer) 142 { 143 PetscFunctionBegin; 144 PetscFunctionReturn(0); 145 } 146 147 /*MC 148 PCTFS - A parallel direct solver intended for problems with very few unknowns (like the 149 coarse grid in multigrid). 150 151 Implemented by Henry M. Tufo III and Paul Fischer 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 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC 161 M*/ 162 PETSC_EXTERN PetscErrorCode PCCreate_TFS(PC pc) 163 { 164 PetscErrorCode ierr; 165 PC_TFS *tfs; 166 PetscMPIInt cmp; 167 168 PetscFunctionBegin; 169 ierr = MPI_Comm_compare(PETSC_COMM_WORLD,PetscObjectComm((PetscObject)pc),&cmp);CHKERRQ(ierr); 170 if (cmp != MPI_IDENT && cmp != MPI_CONGRUENT) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"TFS only works with PETSC_COMM_WORLD objects"); 171 ierr = PetscNewLog(pc,&tfs);CHKERRQ(ierr); 172 173 tfs->xxt = 0; 174 tfs->xyt = 0; 175 tfs->b = 0; 176 tfs->xd = 0; 177 tfs->xo = 0; 178 tfs->nd = 0; 179 180 pc->ops->apply = 0; 181 pc->ops->applytranspose = 0; 182 pc->ops->setup = PCSetUp_TFS; 183 pc->ops->destroy = PCDestroy_TFS; 184 pc->ops->setfromoptions = PCSetFromOptions_TFS; 185 pc->ops->view = PCView_TFS; 186 pc->ops->applyrichardson = 0; 187 pc->ops->applysymmetricleft = 0; 188 pc->ops->applysymmetricright = 0; 189 pc->data = (void*)tfs; 190 PetscFunctionReturn(0); 191 } 192 193