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