1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: mpiaijpc.c,v 1.36 1998/04/03 23:15:06 bsmith Exp bsmith $"; 3 #endif 4 /* 5 Defines a block Jacobi preconditioner for the SeqAIJ/MPIAIJ format. 6 Handles special case of single block per processor. 7 This file knows about storage formats for SeqMPI/MPIAIJ matrices. 8 The general case is handled in aijpc.c 9 */ 10 #include "src/mat/impls/aij/mpi/mpiaij.h" 11 #include "src/pc/pcimpl.h" 12 #include "src/vec/vecimpl.h" 13 #include "src/pc/impls/bjacobi/bjacobi.h" 14 #include "sles.h" 15 16 typedef struct { 17 Vec x, y; 18 } PC_BJacobi_MPIAIJ; 19 20 #undef __FUNC__ 21 #define __FUNC__ "PCDestroy_BJacobi_MPIAIJ" 22 int PCDestroy_BJacobi_MPIAIJ(PC pc) 23 { 24 PC_BJacobi *jac = (PC_BJacobi *) pc->data; 25 PC_BJacobi_MPIAIJ *bjac = (PC_BJacobi_MPIAIJ *) jac->data; 26 int ierr; 27 28 PetscFunctionBegin; 29 ierr = SLESDestroy(jac->sles[0]); CHKERRQ(ierr); 30 PetscFree(jac->sles); 31 ierr = VecDestroy(bjac->x); CHKERRQ(ierr); 32 ierr = VecDestroy(bjac->y); CHKERRQ(ierr); 33 if (jac->l_lens) PetscFree(jac->l_lens); 34 if (jac->g_lens) PetscFree(jac->g_lens); 35 PetscFree(bjac); PetscFree(jac); 36 PetscFunctionReturn(0); 37 } 38 39 40 #undef __FUNC__ 41 #define __FUNC__ "PCSetUpOnBlocks_BJacobi_MPIAIJ" 42 int PCSetUpOnBlocks_BJacobi_MPIAIJ(PC pc) 43 { 44 int ierr; 45 PC_BJacobi *jac = (PC_BJacobi *) pc->data; 46 PC_BJacobi_MPIAIJ *bjac = (PC_BJacobi_MPIAIJ *) jac->data; 47 48 PetscFunctionBegin; 49 ierr = SLESSetUp(jac->sles[0],bjac->x,bjac->y); CHKERRQ(ierr); 50 PetscFunctionReturn(0); 51 } 52 53 #undef __FUNC__ 54 #define __FUNC__ "PCApply_BJacobi_MPIAIJ" 55 int PCApply_BJacobi_MPIAIJ(PC pc,Vec x, Vec y) 56 { 57 int ierr,its; 58 PC_BJacobi *jac = (PC_BJacobi *) pc->data; 59 PC_BJacobi_MPIAIJ *bjac = (PC_BJacobi_MPIAIJ *) jac->data; 60 Scalar *x_array,*x_true_array, *y_array,*y_true_array; 61 62 PetscFunctionBegin; 63 /* 64 The VecPlaceArray() is to avoid having to copy the 65 y vector into the bjac->x vector. The reason for 66 the bjac->x vector is that we need a sequential vector 67 for the sequential solve. 68 */ 69 ierr = VecGetArray(x,&x_array); CHKERRQ(ierr); 70 ierr = VecGetArray(y,&y_array); CHKERRQ(ierr); 71 ierr = VecGetArray(bjac->x,&x_true_array); CHKERRQ(ierr); 72 ierr = VecGetArray(bjac->y,&y_true_array); CHKERRQ(ierr); 73 ierr = VecPlaceArray(bjac->x,x_array); CHKERRQ(ierr); 74 ierr = VecPlaceArray(bjac->y,y_array); CHKERRQ(ierr); 75 ierr = SLESSolve(jac->sles[0],bjac->x,bjac->y,&its); CHKERRQ(ierr); 76 ierr = VecPlaceArray(bjac->x,x_true_array); CHKERRQ(ierr); 77 ierr = VecPlaceArray(bjac->y,y_true_array);CHKERRQ(ierr); 78 PetscFunctionReturn(0); 79 } 80 81 #undef __FUNC__ 82 #define __FUNC__ "PCApplyTrans_BJacobi_MPIAIJ" 83 int PCApplyTrans_BJacobi_MPIAIJ(PC pc,Vec x, Vec y) 84 { 85 int ierr,its; 86 PC_BJacobi *jac = (PC_BJacobi *) pc->data; 87 PC_BJacobi_MPIAIJ *bjac = (PC_BJacobi_MPIAIJ *) jac->data; 88 Scalar *x_array,*x_true_array, *y_array,*y_true_array; 89 90 PetscFunctionBegin; 91 /* 92 The VecPlaceArray() is to avoid having to copy the 93 y vector into the bjac->x vector. The reason for 94 the bjac->x vector is that we need a sequential vector 95 for the sequential solve. 96 */ 97 ierr = VecGetArray(x,&x_array); CHKERRQ(ierr); 98 ierr = VecGetArray(y,&y_array); CHKERRQ(ierr); 99 ierr = VecGetArray(bjac->x,&x_true_array); CHKERRQ(ierr); 100 ierr = VecGetArray(bjac->y,&y_true_array); CHKERRQ(ierr); 101 ierr = VecPlaceArray(bjac->x,x_array); CHKERRQ(ierr); 102 ierr = VecPlaceArray(bjac->y,y_array); CHKERRQ(ierr); 103 ierr = SLESSolveTrans(jac->sles[0],bjac->x,bjac->y,&its); CHKERRQ(ierr); 104 ierr = VecPlaceArray(bjac->x,x_true_array); CHKERRQ(ierr); 105 ierr = VecPlaceArray(bjac->y,y_true_array);CHKERRQ(ierr); 106 PetscFunctionReturn(0); 107 } 108 109 #undef __FUNC__ 110 #define __FUNC__ "PCSetUp_BJacobi_MPIAIJ" 111 int PCSetUp_BJacobi_MPIAIJ(PC pc) 112 { 113 PC_BJacobi *jac = (PC_BJacobi *) pc->data; 114 Mat mat = pc->mat, pmat = pc->pmat; 115 Mat_MPIAIJ *pmatin = (Mat_MPIAIJ *) pmat->data; 116 Mat_MPIAIJ *matin = 0; 117 int ierr, m; 118 SLES sles; 119 Vec x,y; 120 PC_BJacobi_MPIAIJ *bjac; 121 KSP subksp; 122 PC subpc; 123 MatType type; 124 125 PetscFunctionBegin; 126 if (jac->use_true_local) { 127 ierr = MatGetType(pc->mat,&type,PETSC_NULL);CHKERRQ(ierr); 128 if (type != MATMPIAIJ) SETERRQ(PETSC_ERR_ARG_INCOMP,0,"Incompatible matrix type."); 129 matin = (Mat_MPIAIJ *) mat->data; 130 } 131 132 /* set default direct solver with no Krylov method */ 133 if (!pc->setupcalled) { 134 char *prefix; 135 ierr = SLESCreate(PETSC_COMM_SELF,&sles); CHKERRQ(ierr); 136 PLogObjectParent(pc,sles); 137 ierr = SLESGetKSP(sles,&subksp); CHKERRQ(ierr); 138 ierr = KSPSetType(subksp,KSPPREONLY); CHKERRQ(ierr); 139 ierr = SLESGetPC(sles,&subpc); CHKERRQ(ierr); 140 ierr = PCSetType(subpc,PCILU); CHKERRQ(ierr); 141 ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr); 142 ierr = SLESSetOptionsPrefix(sles,prefix); CHKERRQ(ierr); 143 ierr = SLESAppendOptionsPrefix(sles,"sub_"); CHKERRQ(ierr); 144 ierr = SLESSetFromOptions(sles); CHKERRQ(ierr); 145 /* 146 This is not so good. The only reason we need to generate this vector 147 is so KSP may generate seq vectors for the local solves 148 */ 149 ierr = MatGetSize(pmatin->A,&m,&m); CHKERRQ(ierr); 150 ierr = VecCreateSeq(PETSC_COMM_SELF,m,&x); CHKERRQ(ierr); 151 ierr = VecCreateSeq(PETSC_COMM_SELF,m,&y); CHKERRQ(ierr); 152 PLogObjectParent(pc,x); 153 PLogObjectParent(pc,y); 154 155 pc->destroy = PCDestroy_BJacobi_MPIAIJ; 156 pc->apply = PCApply_BJacobi_MPIAIJ; 157 pc->applytrans = PCApplyTrans_BJacobi_MPIAIJ; 158 pc->setuponblocks = PCSetUpOnBlocks_BJacobi_MPIAIJ; 159 160 bjac = (PC_BJacobi_MPIAIJ *) PetscMalloc(sizeof(PC_BJacobi_MPIAIJ));CHKPTRQ(bjac); 161 PLogObjectMemory(pc,sizeof(PC_BJacobi_MPIAIJ)); 162 bjac->x = x; 163 bjac->y = y; 164 165 jac->sles = (SLES*) PetscMalloc( sizeof(SLES) ); CHKPTRQ(jac->sles); 166 jac->sles[0] = sles; 167 jac->data = (void *) bjac; 168 } else { 169 sles = jac->sles[0]; 170 bjac = (PC_BJacobi_MPIAIJ *)jac->data; 171 } 172 if (jac->use_true_local) { 173 ierr = SLESSetOperators(sles,matin->A,pmatin->A,pc->flag); CHKERRQ(ierr); 174 } else { 175 ierr = SLESSetOperators(sles,pmatin->A,pmatin->A,pc->flag); CHKERRQ(ierr); 176 } 177 PetscFunctionReturn(0); 178 } 179 180 #undef __FUNC__ 181 #define __FUNC__ "PCSetUp_BJacobi_SeqAIJ" 182 int PCSetUp_BJacobi_SeqAIJ(PC pc) 183 { 184 PC_BJacobi *jac = (PC_BJacobi *) pc->data; 185 Mat mat = pc->mat, pmat = pc->pmat; 186 int ierr, m; 187 SLES sles; 188 Vec x,y; 189 PC_BJacobi_MPIAIJ *bjac; 190 KSP subksp; 191 PC subpc; 192 MatType type; 193 194 PetscFunctionBegin; 195 if (jac->use_true_local) { 196 MatGetType(mat,&type,PETSC_NULL); 197 if (type != MATSEQAIJ) SETERRQ(PETSC_ERR_ARG_INCOMP,0,"Incompatible matrix type."); 198 } 199 200 /* set default direct solver with no Krylov method */ 201 if (!pc->setupcalled) { 202 char *prefix; 203 ierr = SLESCreate(PETSC_COMM_SELF,&sles); CHKERRQ(ierr); 204 PLogObjectParent(pc,sles); 205 ierr = SLESGetKSP(sles,&subksp); CHKERRQ(ierr); 206 ierr = KSPSetType(subksp,KSPPREONLY); CHKERRQ(ierr); 207 ierr = SLESGetPC(sles,&subpc); CHKERRQ(ierr); 208 ierr = PCSetType(subpc,PCILU); CHKERRQ(ierr); 209 ierr = PCGetOptionsPrefix(pc,&prefix);CHKERRQ(ierr); 210 ierr = SLESSetOptionsPrefix(sles,prefix); CHKERRQ(ierr); 211 ierr = SLESAppendOptionsPrefix(sles,"sub_"); CHKERRQ(ierr); 212 ierr = SLESSetFromOptions(sles); CHKERRQ(ierr); 213 /* 214 This is not so good. The only reason we need to generate this vector 215 is so KSP may generate seq vectors for the local solves 216 */ 217 ierr = MatGetSize(pmat,&m,&m); CHKERRQ(ierr); 218 ierr = VecCreateSeq(PETSC_COMM_SELF,m,&x); CHKERRQ(ierr); 219 ierr = VecCreateSeq(PETSC_COMM_SELF,m,&y); CHKERRQ(ierr); 220 PLogObjectParent(pc,x); 221 PLogObjectParent(pc,y); 222 223 pc->destroy = PCDestroy_BJacobi_MPIAIJ; 224 pc->apply = PCApply_BJacobi_MPIAIJ; 225 pc->applytrans = PCApplyTrans_BJacobi_MPIAIJ; 226 pc->setuponblocks = PCSetUpOnBlocks_BJacobi_MPIAIJ; 227 228 bjac = (PC_BJacobi_MPIAIJ *) PetscMalloc(sizeof(PC_BJacobi_MPIAIJ));CHKPTRQ(bjac); 229 PLogObjectMemory(pc,sizeof(PC_BJacobi_MPIAIJ)); 230 bjac->x = x; 231 bjac->y = y; 232 233 jac->sles = (SLES*) PetscMalloc( sizeof(SLES) ); CHKPTRQ(jac->sles); 234 jac->sles[0] = sles; 235 jac->data = (void *) bjac; 236 } else { 237 sles = jac->sles[0]; 238 bjac = (PC_BJacobi_MPIAIJ *)jac->data; 239 } 240 if (jac->use_true_local) { 241 ierr = SLESSetOperators(sles,mat,pmat,pc->flag); CHKERRQ(ierr); 242 } else { 243 ierr = SLESSetOperators(sles,pmat,pmat,pc->flag); CHKERRQ(ierr); 244 } 245 PetscFunctionReturn(0); 246 } 247 248