1 /* 2 Provides an interface to the Tufo-Fischer parallel direct solver 3 */ 4 5 #include "src/ksp/pc/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 #undef __FUNCT__ 17 #define __FUNCT__ "PCDestroy_TFS" 18 PetscErrorCode PCDestroy_TFS(PC pc) 19 { 20 PC_TFS *tfs = (PC_TFS*)pc->data; 21 PetscErrorCode ierr; 22 23 PetscFunctionBegin; 24 /* free the XXT datastructures */ 25 if (tfs->xxt) { 26 ierr = XXT_free(tfs->xxt);CHKERRQ(ierr); 27 } 28 if (tfs->xyt) { 29 ierr = XYT_free(tfs->xyt);CHKERRQ(ierr); 30 } 31 if (tfs->b) { 32 ierr = VecDestroy(tfs->b);CHKERRQ(ierr); 33 } 34 if (tfs->xd) { 35 ierr = VecDestroy(tfs->xd);CHKERRQ(ierr); 36 } 37 if (tfs->xo) { 38 ierr = VecDestroy(tfs->xo);CHKERRQ(ierr); 39 } 40 ierr = PetscFree(tfs);CHKERRQ(ierr); 41 PetscFunctionReturn(0); 42 } 43 44 #undef __FUNCT__ 45 #define __FUNCT__ "PCApply_TFS_XXT" 46 static PetscErrorCode PCApply_TFS_XXT(PC pc,Vec x,Vec y) 47 { 48 PC_TFS *tfs = (PC_TFS*)pc->data; 49 PetscScalar *xx,*yy; 50 PetscErrorCode ierr; 51 52 PetscFunctionBegin; 53 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 54 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 55 ierr = XXT_solve(tfs->xxt,yy,xx);CHKERRQ(ierr); 56 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 57 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 58 PetscFunctionReturn(0); 59 } 60 61 #undef __FUNCT__ 62 #define __FUNCT__ "PCApply_TFS_XYT" 63 static PetscErrorCode PCApply_TFS_XYT(PC pc,Vec x,Vec y) 64 { 65 PC_TFS *tfs = (PC_TFS*)pc->data; 66 PetscScalar *xx,*yy; 67 PetscErrorCode ierr; 68 69 PetscFunctionBegin; 70 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 71 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 72 ierr = XYT_solve(tfs->xyt,yy,xx);CHKERRQ(ierr); 73 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 74 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 75 PetscFunctionReturn(0); 76 } 77 78 #undef __FUNCT__ 79 #define __FUNCT__ "LocalMult_TFS" 80 static PetscErrorCode LocalMult_TFS(PC pc,PetscScalar *xin,PetscScalar *xout) 81 { 82 PC_TFS *tfs = (PC_TFS*)pc->data; 83 Mat A = pc->pmat; 84 Mat_MPIAIJ *a = (Mat_MPIAIJ*)A->data; 85 PetscErrorCode ierr; 86 87 PetscFunctionBegin; 88 ierr = VecPlaceArray(tfs->b,xout);CHKERRQ(ierr); 89 ierr = VecPlaceArray(tfs->xd,xin);CHKERRQ(ierr); 90 ierr = VecPlaceArray(tfs->xo,xin+tfs->nd);CHKERRQ(ierr); 91 ierr = MatMult(a->A,tfs->xd,tfs->b);CHKERRQ(ierr); 92 ierr = MatMultAdd(a->B,tfs->xo,tfs->b,tfs->b);CHKERRQ(ierr); 93 PetscFunctionReturn(0); 94 } 95 96 #undef __FUNCT__ 97 #define __FUNCT__ "PCSetUp_TFS" 98 static PetscErrorCode PCSetUp_TFS(PC pc) 99 { 100 PC_TFS *tfs = (PC_TFS*)pc->data; 101 Mat A = pc->pmat; 102 Mat_MPIAIJ *a = (Mat_MPIAIJ*)A->data; 103 PetscErrorCode ierr; 104 PetscInt *localtoglobal,ncol,i; 105 PetscTruth ismpiaij; 106 107 /* 108 PetscTruth issymmetric; 109 Petsc Real tol = 0.0; 110 */ 111 112 PetscFunctionBegin; 113 if (A->N != A->M) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 114 ierr = PetscTypeCompare((PetscObject)pc->pmat,MATMPIAIJ,&ismpiaij);CHKERRQ(ierr); 115 if (!ismpiaij) { 116 SETERRQ(PETSC_ERR_SUP,"Currently only supports MPIAIJ matrices"); 117 } 118 119 /* generate the local to global mapping */ 120 ncol = a->A->n + a->B->n; 121 ierr = PetscMalloc((ncol)*sizeof(int),&localtoglobal);CHKERRQ(ierr); 122 for (i=0; i<a->A->n; i++) { 123 localtoglobal[i] = a->cstart + i + 1; 124 } 125 for (i=0; i<a->B->n; i++) { 126 localtoglobal[i+a->A->n] = a->garray[i] + 1; 127 } 128 /* generate the vectors needed for the local solves */ 129 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,a->A->m,PETSC_NULL,&tfs->b);CHKERRQ(ierr); 130 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,a->A->n,PETSC_NULL,&tfs->xd);CHKERRQ(ierr); 131 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,a->B->n,PETSC_NULL,&tfs->xo);CHKERRQ(ierr); 132 tfs->nd = a->A->n; 133 134 135 /* ierr = MatIsSymmetric(A,tol,&issymmetric); */ 136 /* if (issymmetric) { */ 137 ierr = PetscBarrier((PetscObject)pc);CHKERRQ(ierr); 138 if (A->symmetric) { 139 tfs->xxt = XXT_new(); 140 ierr = XXT_factor(tfs->xxt,localtoglobal,A->m,ncol,(void*)LocalMult_TFS,pc);CHKERRQ(ierr); 141 pc->ops->apply = PCApply_TFS_XXT; 142 } else { 143 tfs->xyt = XYT_new(); 144 ierr = XYT_factor(tfs->xyt,localtoglobal,A->m,ncol,(void*)LocalMult_TFS,pc);CHKERRQ(ierr); 145 pc->ops->apply = PCApply_TFS_XYT; 146 } 147 148 ierr = PetscFree(localtoglobal);CHKERRQ(ierr); 149 PetscFunctionReturn(0); 150 } 151 152 #undef __FUNCT__ 153 #define __FUNCT__ "PCSetFromOptions_TFS" 154 static PetscErrorCode PCSetFromOptions_TFS(PC pc) 155 { 156 PetscFunctionBegin; 157 PetscFunctionReturn(0); 158 } 159 #undef __FUNCT__ 160 #define __FUNCT__ "PCView_TFS" 161 static PetscErrorCode PCView_TFS(PC pc,PetscViewer viewer) 162 { 163 PetscFunctionBegin; 164 PetscFunctionReturn(0); 165 } 166 167 EXTERN_C_BEGIN 168 #undef __FUNCT__ 169 #define __FUNCT__ "PCCreate_TFS" 170 PetscErrorCode PCCreate_TFS(PC pc) 171 { 172 PetscErrorCode ierr; 173 PC_TFS *tfs; 174 175 PetscFunctionBegin; 176 ierr = PetscNew(PC_TFS,&tfs);CHKERRQ(ierr); 177 ierr = PetscLogObjectMemory(pc,sizeof(PC_TFS));CHKERRQ(ierr); 178 179 tfs->xxt = 0; 180 tfs->xyt = 0; 181 tfs->b = 0; 182 tfs->xd = 0; 183 tfs->xo = 0; 184 tfs->nd = 0; 185 186 pc->ops->apply = 0; 187 pc->ops->applytranspose = 0; 188 pc->ops->setup = PCSetUp_TFS; 189 pc->ops->destroy = PCDestroy_TFS; 190 pc->ops->setfromoptions = PCSetFromOptions_TFS; 191 pc->ops->view = PCView_TFS; 192 pc->ops->applyrichardson = 0; 193 pc->ops->applysymmetricleft = 0; 194 pc->ops->applysymmetricright = 0; 195 pc->data = (void*)tfs; 196 PetscFunctionReturn(0); 197 } 198 EXTERN_C_END 199 200