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