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 ierr = VecResetArray(tfs->b);CHKERRQ(ierr); 95 ierr = VecResetArray(tfs->xd);CHKERRQ(ierr); 96 ierr = VecResetArray(tfs->xo);CHKERRQ(ierr); 97 PetscFunctionReturn(0); 98 } 99 100 #undef __FUNCT__ 101 #define __FUNCT__ "PCSetUp_TFS" 102 static PetscErrorCode PCSetUp_TFS(PC pc) 103 { 104 PC_TFS *tfs = (PC_TFS*)pc->data; 105 Mat A = pc->pmat; 106 Mat_MPIAIJ *a = (Mat_MPIAIJ*)A->data; 107 PetscErrorCode ierr; 108 PetscInt *localtoglobal,ncol,i; 109 PetscTruth ismpiaij; 110 111 /* 112 PetscTruth issymmetric; 113 Petsc Real tol = 0.0; 114 */ 115 116 PetscFunctionBegin; 117 if (A->cmap.N != A->rmap.N) SETERRQ(PETSC_ERR_ARG_SIZ,"matrix must be square"); 118 ierr = PetscTypeCompare((PetscObject)pc->pmat,MATMPIAIJ,&ismpiaij);CHKERRQ(ierr); 119 if (!ismpiaij) { 120 SETERRQ(PETSC_ERR_SUP,"Currently only supports MPIAIJ matrices"); 121 } 122 123 /* generate the local to global mapping */ 124 ncol = a->A->cmap.n + a->B->cmap.n; 125 ierr = PetscMalloc((ncol)*sizeof(int),&localtoglobal);CHKERRQ(ierr); 126 for (i=0; i<a->A->cmap.n; i++) { 127 localtoglobal[i] = A->cmap.rstart + i + 1; 128 } 129 for (i=0; i<a->B->cmap.n; i++) { 130 localtoglobal[i+a->A->cmap.n] = a->garray[i] + 1; 131 } 132 /* generate the vectors needed for the local solves */ 133 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,a->A->rmap.n,PETSC_NULL,&tfs->b);CHKERRQ(ierr); 134 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,a->A->cmap.n,PETSC_NULL,&tfs->xd);CHKERRQ(ierr); 135 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,a->B->cmap.n,PETSC_NULL,&tfs->xo);CHKERRQ(ierr); 136 tfs->nd = a->A->cmap.n; 137 138 139 /* ierr = MatIsSymmetric(A,tol,&issymmetric); */ 140 /* if (issymmetric) { */ 141 ierr = PetscBarrier((PetscObject)pc);CHKERRQ(ierr); 142 if (A->symmetric) { 143 tfs->xxt = XXT_new(); 144 ierr = XXT_factor(tfs->xxt,localtoglobal,A->rmap.n,ncol,(void*)LocalMult_TFS,pc);CHKERRQ(ierr); 145 pc->ops->apply = PCApply_TFS_XXT; 146 } else { 147 tfs->xyt = XYT_new(); 148 ierr = XYT_factor(tfs->xyt,localtoglobal,A->rmap.n,ncol,(void*)LocalMult_TFS,pc);CHKERRQ(ierr); 149 pc->ops->apply = PCApply_TFS_XYT; 150 } 151 152 ierr = PetscFree(localtoglobal);CHKERRQ(ierr); 153 PetscFunctionReturn(0); 154 } 155 156 #undef __FUNCT__ 157 #define __FUNCT__ "PCSetFromOptions_TFS" 158 static PetscErrorCode PCSetFromOptions_TFS(PC pc) 159 { 160 PetscFunctionBegin; 161 PetscFunctionReturn(0); 162 } 163 #undef __FUNCT__ 164 #define __FUNCT__ "PCView_TFS" 165 static PetscErrorCode PCView_TFS(PC pc,PetscViewer viewer) 166 { 167 PetscFunctionBegin; 168 PetscFunctionReturn(0); 169 } 170 171 EXTERN_C_BEGIN 172 #undef __FUNCT__ 173 #define __FUNCT__ "PCCreate_TFS" 174 PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_TFS(PC pc) 175 { 176 PetscErrorCode ierr; 177 PC_TFS *tfs; 178 179 PetscFunctionBegin; 180 ierr = PetscNewLog(pc,PC_TFS,&tfs);CHKERRQ(ierr); 181 182 tfs->xxt = 0; 183 tfs->xyt = 0; 184 tfs->b = 0; 185 tfs->xd = 0; 186 tfs->xo = 0; 187 tfs->nd = 0; 188 189 pc->ops->apply = 0; 190 pc->ops->applytranspose = 0; 191 pc->ops->setup = PCSetUp_TFS; 192 pc->ops->destroy = PCDestroy_TFS; 193 pc->ops->setfromoptions = PCSetFromOptions_TFS; 194 pc->ops->view = PCView_TFS; 195 pc->ops->applyrichardson = 0; 196 pc->ops->applysymmetricleft = 0; 197 pc->ops->applysymmetricright = 0; 198 pc->data = (void*)tfs; 199 PetscFunctionReturn(0); 200 } 201 EXTERN_C_END 202 203