1 2 /* 3 Creates hypre ijmatrix from PETSc matrix 4 */ 5 #include <petscsys.h> 6 #include <petsc/private/petschypre.h> 7 #include <petsc/private/matimpl.h> 8 #include <petscdmda.h> /*I "petscdmda.h" I*/ 9 #include <../src/dm/impls/da/hypre/mhyp.h> 10 11 /* -----------------------------------------------------------------------------------------------------------------*/ 12 13 /*MC 14 MATHYPRESTRUCT - MATHYPRESTRUCT = "hyprestruct" - A matrix type to be used for parallel sparse matrices 15 based on the hypre HYPRE_StructMatrix. 16 17 Level: intermediate 18 19 Notes: 20 Unlike the more general support for blocks in hypre this allows only one block per process and requires the block 21 be defined by a DMDA. 22 23 The matrix needs a DMDA associated with it by either a call to MatSetDM() or if the matrix is obtained from DMCreateMatrix() 24 25 .seealso: MatCreate(), PCPFMG, MatSetDM(), DMCreateMatrix() 26 M*/ 27 28 PetscErrorCode MatSetValuesLocal_HYPREStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv) 29 { 30 HYPRE_Int index[3],entries[7]; 31 PetscInt i,j,stencil,row; 32 HYPRE_Complex *values = (HYPRE_Complex*)y; 33 Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data; 34 35 PetscFunctionBegin; 36 if (PetscUnlikely(ncol > 7)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"ncol %D > 7 too large",ncol); 37 for (i=0; i<nrow; i++) { 38 for (j=0; j<ncol; j++) { 39 stencil = icol[j] - irow[i]; 40 if (!stencil) { 41 entries[j] = 3; 42 } else if (stencil == -1) { 43 entries[j] = 2; 44 } else if (stencil == 1) { 45 entries[j] = 4; 46 } else if (stencil == -ex->gnx) { 47 entries[j] = 1; 48 } else if (stencil == ex->gnx) { 49 entries[j] = 5; 50 } else if (stencil == -ex->gnxgny) { 51 entries[j] = 0; 52 } else if (stencil == ex->gnxgny) { 53 entries[j] = 6; 54 } else SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row %D local column %D have bad stencil %D",irow[i],icol[j],stencil); 55 } 56 row = ex->gindices[irow[i]] - ex->rstart; 57 index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx)); 58 index[1] = (HYPRE_Int)(ex->ys + ((row/ex->nx) % ex->ny)); 59 index[2] = (HYPRE_Int)(ex->zs + (row/(ex->nxny))); 60 if (addv == ADD_VALUES) { 61 PetscStackCallStandard(HYPRE_StructMatrixAddToValues,(ex->hmat,index,ncol,entries,values)); 62 } else { 63 PetscStackCallStandard(HYPRE_StructMatrixSetValues,(ex->hmat,index,ncol,entries,values)); 64 } 65 values += ncol; 66 } 67 PetscFunctionReturn(0); 68 } 69 70 PetscErrorCode MatZeroRowsLocal_HYPREStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscScalar d,Vec x,Vec b) 71 { 72 PetscErrorCode ierr; 73 HYPRE_Int index[3],entries[7] = {0,1,2,3,4,5,6}; 74 PetscInt row,i; 75 HYPRE_Complex values[7]; 76 Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data; 77 78 PetscFunctionBegin; 79 if (x && b) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"No support"); 80 ierr = PetscArrayzero(values,7);CHKERRQ(ierr); 81 ierr = PetscHYPREScalarCast(d,&values[3]);CHKERRQ(ierr); 82 for (i=0; i<nrow; i++) { 83 row = ex->gindices[irow[i]] - ex->rstart; 84 index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx)); 85 index[1] = (HYPRE_Int)(ex->ys + ((row/ex->nx) % ex->ny)); 86 index[2] = (HYPRE_Int)(ex->zs + (row/(ex->nxny))); 87 PetscStackCallStandard(HYPRE_StructMatrixSetValues,(ex->hmat,index,7,entries,values)); 88 } 89 PetscStackCallStandard(HYPRE_StructMatrixAssemble,(ex->hmat)); 90 PetscFunctionReturn(0); 91 } 92 93 PetscErrorCode MatZeroEntries_HYPREStruct_3d(Mat mat) 94 { 95 HYPRE_Int indices[7] = {0,1,2,3,4,5,6}; 96 Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data; 97 98 PetscFunctionBegin; 99 /* hypre has no public interface to do this */ 100 PetscStackCallStandard(hypre_StructMatrixClearBoxValues,(ex->hmat,&ex->hbox,7,indices,0,1)); 101 PetscStackCallStandard(HYPRE_StructMatrixAssemble,(ex->hmat)); 102 PetscFunctionReturn(0); 103 } 104 105 static PetscErrorCode MatSetUp_HYPREStruct(Mat mat) 106 { 107 PetscErrorCode ierr; 108 Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data; 109 HYPRE_Int sw[6]; 110 HYPRE_Int hlower[3],hupper[3]; 111 PetscInt dim,dof,psw,nx,ny,nz,ilower[3],iupper[3],ssize,i; 112 DMBoundaryType px,py,pz; 113 DMDAStencilType st; 114 ISLocalToGlobalMapping ltog; 115 DM da; 116 117 PetscFunctionBegin; 118 ierr = MatGetDM(mat,(DM*)&da);CHKERRQ(ierr); 119 ex->da = da; 120 ierr = PetscObjectReference((PetscObject)da);CHKERRQ(ierr); 121 122 ierr = DMDAGetInfo(ex->da,&dim,0,0,0,0,0,0,&dof,&psw,&px,&py,&pz,&st);CHKERRQ(ierr); 123 ierr = DMDAGetCorners(ex->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);CHKERRQ(ierr); 124 125 /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */ 126 iupper[0] += ilower[0] - 1; 127 iupper[1] += ilower[1] - 1; 128 iupper[2] += ilower[2] - 1; 129 hlower[0] = (HYPRE_Int)ilower[0]; 130 hlower[1] = (HYPRE_Int)ilower[1]; 131 hlower[2] = (HYPRE_Int)ilower[2]; 132 hupper[0] = (HYPRE_Int)iupper[0]; 133 hupper[1] = (HYPRE_Int)iupper[1]; 134 hupper[2] = (HYPRE_Int)iupper[2]; 135 136 /* the hypre_Box is used to zero out the matrix entries in MatZeroValues() */ 137 ex->hbox.imin[0] = hlower[0]; 138 ex->hbox.imin[1] = hlower[1]; 139 ex->hbox.imin[2] = hlower[2]; 140 ex->hbox.imax[0] = hupper[0]; 141 ex->hbox.imax[1] = hupper[1]; 142 ex->hbox.imax[2] = hupper[2]; 143 144 /* create the hypre grid object and set its information */ 145 if (dof > 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Currently only support for scalar problems"); 146 if (px || py || pz) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add periodic support by calling HYPRE_StructGridSetPeriodic()"); 147 PetscStackCallStandard(HYPRE_StructGridCreate,(ex->hcomm,dim,&ex->hgrid)); 148 PetscStackCallStandard(HYPRE_StructGridSetExtents,(ex->hgrid,hlower,hupper)); 149 PetscStackCallStandard(HYPRE_StructGridAssemble,(ex->hgrid)); 150 151 sw[5] = sw[4] = sw[3] = sw[2] = sw[1] = sw[0] = psw; 152 PetscStackCallStandard(HYPRE_StructGridSetNumGhost,(ex->hgrid,sw)); 153 154 /* create the hypre stencil object and set its information */ 155 if (sw[0] > 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add support for wider stencils"); 156 if (st == DMDA_STENCIL_BOX) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add support for box stencils"); 157 if (dim == 1) { 158 HYPRE_Int offsets[3][1] = {{-1},{0},{1}}; 159 ssize = 3; 160 PetscStackCallStandard(HYPRE_StructStencilCreate,(dim,ssize,&ex->hstencil)); 161 for (i=0; i<ssize; i++) { 162 PetscStackCallStandard(HYPRE_StructStencilSetElement,(ex->hstencil,i,offsets[i])); 163 } 164 } else if (dim == 2) { 165 HYPRE_Int offsets[5][2] = {{0,-1},{-1,0},{0,0},{1,0},{0,1}}; 166 ssize = 5; 167 PetscStackCallStandard(HYPRE_StructStencilCreate,(dim,ssize,&ex->hstencil)); 168 for (i=0; i<ssize; i++) { 169 PetscStackCallStandard(HYPRE_StructStencilSetElement,(ex->hstencil,i,offsets[i])); 170 } 171 } else if (dim == 3) { 172 HYPRE_Int offsets[7][3] = {{0,0,-1},{0,-1,0},{-1,0,0},{0,0,0},{1,0,0},{0,1,0},{0,0,1}}; 173 ssize = 7; 174 PetscStackCallStandard(HYPRE_StructStencilCreate,(dim,ssize,&ex->hstencil)); 175 for (i=0; i<ssize; i++) { 176 PetscStackCallStandard(HYPRE_StructStencilSetElement,(ex->hstencil,i,offsets[i])); 177 } 178 } 179 180 /* create the HYPRE vector for rhs and solution */ 181 PetscStackCallStandard(HYPRE_StructVectorCreate,(ex->hcomm,ex->hgrid,&ex->hb)); 182 PetscStackCallStandard(HYPRE_StructVectorCreate,(ex->hcomm,ex->hgrid,&ex->hx)); 183 PetscStackCallStandard(HYPRE_StructVectorInitialize,(ex->hb)); 184 PetscStackCallStandard(HYPRE_StructVectorInitialize,(ex->hx)); 185 PetscStackCallStandard(HYPRE_StructVectorAssemble,(ex->hb)); 186 PetscStackCallStandard(HYPRE_StructVectorAssemble,(ex->hx)); 187 188 /* create the hypre matrix object and set its information */ 189 PetscStackCallStandard(HYPRE_StructMatrixCreate,(ex->hcomm,ex->hgrid,ex->hstencil,&ex->hmat)); 190 PetscStackCallStandard(HYPRE_StructGridDestroy,(ex->hgrid)); 191 PetscStackCallStandard(HYPRE_StructStencilDestroy,(ex->hstencil)); 192 if (ex->needsinitialization) { 193 PetscStackCallStandard(HYPRE_StructMatrixInitialize,(ex->hmat)); 194 ex->needsinitialization = PETSC_FALSE; 195 } 196 197 /* set the global and local sizes of the matrix */ 198 ierr = DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);CHKERRQ(ierr); 199 ierr = MatSetSizes(mat,dof*nx*ny*nz,dof*nx*ny*nz,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 200 ierr = PetscLayoutSetBlockSize(mat->rmap,1);CHKERRQ(ierr); 201 ierr = PetscLayoutSetBlockSize(mat->cmap,1);CHKERRQ(ierr); 202 ierr = PetscLayoutSetUp(mat->rmap);CHKERRQ(ierr); 203 ierr = PetscLayoutSetUp(mat->cmap);CHKERRQ(ierr); 204 mat->preallocated = PETSC_TRUE; 205 206 if (dim == 3) { 207 mat->ops->setvalueslocal = MatSetValuesLocal_HYPREStruct_3d; 208 mat->ops->zerorowslocal = MatZeroRowsLocal_HYPREStruct_3d; 209 mat->ops->zeroentries = MatZeroEntries_HYPREStruct_3d; 210 211 /* ierr = MatZeroEntries_HYPREStruct_3d(mat);CHKERRQ(ierr); */ 212 } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Only support for 3d DMDA currently"); 213 214 /* get values that will be used repeatedly in MatSetValuesLocal() and MatZeroRowsLocal() repeatedly */ 215 ierr = MatGetOwnershipRange(mat,&ex->rstart,NULL);CHKERRQ(ierr); 216 ierr = DMGetLocalToGlobalMapping(ex->da,<og);CHKERRQ(ierr); 217 ierr = ISLocalToGlobalMappingGetIndices(ltog, (const PetscInt **) &ex->gindices);CHKERRQ(ierr); 218 ierr = DMDAGetGhostCorners(ex->da,0,0,0,&ex->gnx,&ex->gnxgny,0);CHKERRQ(ierr); 219 ex->gnxgny *= ex->gnx; 220 ierr = DMDAGetCorners(ex->da,&ex->xs,&ex->ys,&ex->zs,&ex->nx,&ex->ny,0);CHKERRQ(ierr); 221 ex->nxny = ex->nx*ex->ny; 222 PetscFunctionReturn(0); 223 } 224 225 PetscErrorCode MatMult_HYPREStruct(Mat A,Vec x,Vec y) 226 { 227 PetscErrorCode ierr; 228 const PetscScalar *xx; 229 PetscScalar *yy; 230 PetscInt ilower[3],iupper[3]; 231 HYPRE_Int hlower[3],hupper[3]; 232 Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(A->data); 233 234 PetscFunctionBegin; 235 ierr = DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);CHKERRQ(ierr); 236 /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */ 237 iupper[0] += ilower[0] - 1; 238 iupper[1] += ilower[1] - 1; 239 iupper[2] += ilower[2] - 1; 240 hlower[0] = (HYPRE_Int)ilower[0]; 241 hlower[1] = (HYPRE_Int)ilower[1]; 242 hlower[2] = (HYPRE_Int)ilower[2]; 243 hupper[0] = (HYPRE_Int)iupper[0]; 244 hupper[1] = (HYPRE_Int)iupper[1]; 245 hupper[2] = (HYPRE_Int)iupper[2]; 246 247 /* copy x values over to hypre */ 248 PetscStackCallStandard(HYPRE_StructVectorSetConstantValues,(mx->hb,0.0)); 249 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 250 PetscStackCallStandard(HYPRE_StructVectorSetBoxValues,(mx->hb,hlower,hupper,(HYPRE_Complex*)xx)); 251 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 252 PetscStackCallStandard(HYPRE_StructVectorAssemble,(mx->hb)); 253 PetscStackCallStandard(HYPRE_StructMatrixMatvec,(1.0,mx->hmat,mx->hb,0.0,mx->hx)); 254 255 /* copy solution values back to PETSc */ 256 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 257 PetscStackCallStandard(HYPRE_StructVectorGetBoxValues,(mx->hx,hlower,hupper,(HYPRE_Complex*)yy)); 258 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 259 PetscFunctionReturn(0); 260 } 261 262 PetscErrorCode MatAssemblyEnd_HYPREStruct(Mat mat,MatAssemblyType mode) 263 { 264 Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data; 265 266 PetscFunctionBegin; 267 PetscStackCallStandard(HYPRE_StructMatrixAssemble,(ex->hmat)); 268 /* PetscStackCallStandard(HYPRE_StructMatrixPrint,("dummy",ex->hmat,0)); */ 269 PetscFunctionReturn(0); 270 } 271 272 PetscErrorCode MatZeroEntries_HYPREStruct(Mat mat) 273 { 274 PetscFunctionBegin; 275 /* before the DMDA is set to the matrix the zero doesn't need to do anything */ 276 PetscFunctionReturn(0); 277 } 278 279 PetscErrorCode MatDestroy_HYPREStruct(Mat mat) 280 { 281 Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data; 282 PetscErrorCode ierr; 283 284 PetscFunctionBegin; 285 PetscStackCallStandard(HYPRE_StructMatrixDestroy,(ex->hmat)); 286 PetscStackCallStandard(HYPRE_StructVectorDestroy,(ex->hx)); 287 PetscStackCallStandard(HYPRE_StructVectorDestroy,(ex->hb)); 288 ierr = PetscObjectDereference((PetscObject)ex->da);CHKERRQ(ierr); 289 ierr = MPI_Comm_free(&(ex->hcomm));CHKERRQ(ierr); 290 ierr = PetscFree(ex);CHKERRQ(ierr); 291 PetscFunctionReturn(0); 292 } 293 294 PETSC_EXTERN PetscErrorCode MatCreate_HYPREStruct(Mat B) 295 { 296 Mat_HYPREStruct *ex; 297 PetscErrorCode ierr; 298 299 PetscFunctionBegin; 300 ierr = PetscNewLog(B,&ex);CHKERRQ(ierr); 301 B->data = (void*)ex; 302 B->rmap->bs = 1; 303 B->assembled = PETSC_FALSE; 304 305 B->insertmode = NOT_SET_VALUES; 306 307 B->ops->assemblyend = MatAssemblyEnd_HYPREStruct; 308 B->ops->mult = MatMult_HYPREStruct; 309 B->ops->zeroentries = MatZeroEntries_HYPREStruct; 310 B->ops->destroy = MatDestroy_HYPREStruct; 311 B->ops->setup = MatSetUp_HYPREStruct; 312 313 ex->needsinitialization = PETSC_TRUE; 314 315 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)B),&(ex->hcomm));CHKERRQ(ierr); 316 ierr = PetscObjectChangeTypeName((PetscObject)B,MATHYPRESTRUCT);CHKERRQ(ierr); 317 PetscFunctionReturn(0); 318 } 319 320 /*MC 321 MATHYPRESSTRUCT - MATHYPRESSTRUCT = "hypresstruct" - A matrix type to be used for parallel sparse matrices 322 based on the hypre HYPRE_SStructMatrix. 323 324 325 Level: intermediate 326 327 Notes: 328 Unlike hypre's general semi-struct object consisting of a collection of structured-grid objects and unstructured 329 grid objects, we restrict the semi-struct objects to consist of only structured-grid components. 330 331 Unlike the more general support for parts and blocks in hypre this allows only one part, and one block per process and requires the block 332 be defined by a DMDA. 333 334 The matrix needs a DMDA associated with it by either a call to MatSetDM() or if the matrix is obtained from DMCreateMatrix() 335 336 M*/ 337 338 PetscErrorCode MatSetValuesLocal_HYPRESStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv) 339 { 340 PetscErrorCode ierr; 341 HYPRE_Int index[3],*entries; 342 PetscInt i,j,stencil; 343 HYPRE_Complex *values = (HYPRE_Complex*)y; 344 Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data; 345 346 PetscInt part = 0; /* Petsc sstruct interface only allows 1 part */ 347 PetscInt ordering; 348 PetscInt grid_rank, to_grid_rank; 349 PetscInt var_type, to_var_type; 350 PetscInt to_var_entry = 0; 351 PetscInt nvars= ex->nvars; 352 PetscInt row; 353 354 PetscFunctionBegin; 355 ierr = PetscMalloc1(7*nvars,&entries);CHKERRQ(ierr); 356 ordering = ex->dofs_order; /* ordering= 0 nodal ordering 357 1 variable ordering */ 358 /* stencil entries are ordered by variables: var0_stencil0, var0_stencil1, ..., var0_stencil6, var1_stencil0, var1_stencil1, ... */ 359 if (!ordering) { 360 for (i=0; i<nrow; i++) { 361 grid_rank= irow[i]/nvars; 362 var_type = (irow[i] % nvars); 363 364 for (j=0; j<ncol; j++) { 365 to_grid_rank= icol[j]/nvars; 366 to_var_type = (icol[j] % nvars); 367 368 to_var_entry= to_var_entry*7; 369 entries[j] = to_var_entry; 370 371 stencil = to_grid_rank-grid_rank; 372 if (!stencil) { 373 entries[j] += 3; 374 } else if (stencil == -1) { 375 entries[j] += 2; 376 } else if (stencil == 1) { 377 entries[j] += 4; 378 } else if (stencil == -ex->gnx) { 379 entries[j] += 1; 380 } else if (stencil == ex->gnx) { 381 entries[j] += 5; 382 } else if (stencil == -ex->gnxgny) { 383 entries[j] += 0; 384 } else if (stencil == ex->gnxgny) { 385 entries[j] += 6; 386 } else SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row %D local column %D have bad stencil %D",irow[i],icol[j],stencil); 387 } 388 389 row = ex->gindices[grid_rank] - ex->rstart; 390 index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx)); 391 index[1] = (HYPRE_Int)(ex->ys + ((row/ex->nx) % ex->ny)); 392 index[2] = (HYPRE_Int)(ex->zs + (row/(ex->nxny))); 393 394 if (addv == ADD_VALUES) { 395 PetscStackCallStandard(HYPRE_SStructMatrixAddToValues,(ex->ss_mat,part,index,var_type,ncol,entries,values)); 396 } else { 397 PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,index,var_type,ncol,entries,values)); 398 } 399 values += ncol; 400 } 401 } else { 402 for (i=0; i<nrow; i++) { 403 var_type = irow[i]/(ex->gnxgnygnz); 404 grid_rank= irow[i] - var_type*(ex->gnxgnygnz); 405 406 for (j=0; j<ncol; j++) { 407 to_var_type = icol[j]/(ex->gnxgnygnz); 408 to_grid_rank= icol[j] - to_var_type*(ex->gnxgnygnz); 409 410 to_var_entry= to_var_entry*7; 411 entries[j] = to_var_entry; 412 413 stencil = to_grid_rank-grid_rank; 414 if (!stencil) { 415 entries[j] += 3; 416 } else if (stencil == -1) { 417 entries[j] += 2; 418 } else if (stencil == 1) { 419 entries[j] += 4; 420 } else if (stencil == -ex->gnx) { 421 entries[j] += 1; 422 } else if (stencil == ex->gnx) { 423 entries[j] += 5; 424 } else if (stencil == -ex->gnxgny) { 425 entries[j] += 0; 426 } else if (stencil == ex->gnxgny) { 427 entries[j] += 6; 428 } else SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row %D local column %D have bad stencil %D",irow[i],icol[j],stencil); 429 } 430 431 row = ex->gindices[grid_rank] - ex->rstart; 432 index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx)); 433 index[1] = (HYPRE_Int)(ex->ys + ((row/ex->nx) % ex->ny)); 434 index[2] = (HYPRE_Int)(ex->zs + (row/(ex->nxny))); 435 436 if (addv == ADD_VALUES) { 437 PetscStackCallStandard(HYPRE_SStructMatrixAddToValues,(ex->ss_mat,part,index,var_type,ncol,entries,values)); 438 } else { 439 PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,index,var_type,ncol,entries,values)); 440 } 441 values += ncol; 442 } 443 } 444 ierr = PetscFree(entries);CHKERRQ(ierr); 445 PetscFunctionReturn(0); 446 } 447 448 PetscErrorCode MatZeroRowsLocal_HYPRESStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscScalar d,Vec x,Vec b) 449 { 450 PetscErrorCode ierr; 451 HYPRE_Int index[3],*entries; 452 PetscInt i; 453 HYPRE_Complex **values; 454 Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data; 455 456 PetscInt part = 0; /* Petsc sstruct interface only allows 1 part */ 457 PetscInt ordering = ex->dofs_order; 458 PetscInt grid_rank; 459 PetscInt var_type; 460 PetscInt nvars= ex->nvars; 461 PetscInt row; 462 463 PetscFunctionBegin; 464 if (x && b) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"No support"); 465 ierr = PetscMalloc1(7*nvars,&entries);CHKERRQ(ierr); 466 467 ierr = PetscMalloc1(nvars,&values);CHKERRQ(ierr); 468 ierr = PetscMalloc1(7*nvars*nvars,&values[0]);CHKERRQ(ierr); 469 for (i=1; i<nvars; i++) { 470 values[i] = values[i-1] + nvars*7; 471 } 472 473 for (i=0; i< nvars; i++) { 474 ierr = PetscArrayzero(values[i],nvars*7*sizeof(HYPRE_Complex));CHKERRQ(ierr); 475 ierr = PetscHYPREScalarCast(d,values[i]+3);CHKERRQ(ierr); 476 } 477 478 for (i=0; i< nvars*7; i++) entries[i] = i; 479 480 if (!ordering) { 481 for (i=0; i<nrow; i++) { 482 grid_rank = irow[i] / nvars; 483 var_type = (irow[i] % nvars); 484 485 row = ex->gindices[grid_rank] - ex->rstart; 486 index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx)); 487 index[1] = (HYPRE_Int)(ex->ys + ((row/ex->nx) % ex->ny)); 488 index[2] = (HYPRE_Int)(ex->zs + (row/(ex->nxny))); 489 PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,index,var_type,7*nvars,entries,values[var_type])); 490 } 491 } else { 492 for (i=0; i<nrow; i++) { 493 var_type = irow[i]/(ex->gnxgnygnz); 494 grid_rank= irow[i] - var_type*(ex->gnxgnygnz); 495 496 row = ex->gindices[grid_rank] - ex->rstart; 497 index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx)); 498 index[1] = (HYPRE_Int)(ex->ys + ((row/ex->nx) % ex->ny)); 499 index[2] = (HYPRE_Int)(ex->zs + (row/(ex->nxny))); 500 PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,index,var_type,7*nvars,entries,values[var_type])); 501 } 502 } 503 PetscStackCallStandard(HYPRE_SStructMatrixAssemble,(ex->ss_mat)); 504 ierr = PetscFree(values[0]);CHKERRQ(ierr); 505 ierr = PetscFree(values);CHKERRQ(ierr); 506 ierr = PetscFree(entries);CHKERRQ(ierr); 507 PetscFunctionReturn(0); 508 } 509 510 PetscErrorCode MatZeroEntries_HYPRESStruct_3d(Mat mat) 511 { 512 PetscErrorCode ierr; 513 Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data; 514 PetscInt nvars= ex->nvars; 515 PetscInt size; 516 PetscInt part= 0; /* only one part */ 517 518 PetscFunctionBegin; 519 size = ((ex->hbox.imax[0])-(ex->hbox.imin[0])+1)*((ex->hbox.imax[1])-(ex->hbox.imin[1])+1)*((ex->hbox.imax[2])-(ex->hbox.imin[2])+1); 520 { 521 HYPRE_Int i,*entries,iupper[3],ilower[3]; 522 HYPRE_Complex *values; 523 524 525 for (i= 0; i< 3; i++) { 526 ilower[i] = ex->hbox.imin[i]; 527 iupper[i] = ex->hbox.imax[i]; 528 } 529 530 ierr = PetscMalloc2(nvars*7,&entries,nvars*7*size,&values);CHKERRQ(ierr); 531 for (i= 0; i< nvars*7; i++) entries[i] = i; 532 ierr = PetscArrayzero(values,nvars*7*size);CHKERRQ(ierr); 533 534 for (i= 0; i< nvars; i++) { 535 PetscStackCallStandard(HYPRE_SStructMatrixSetBoxValues,(ex->ss_mat,part,ilower,iupper,i,nvars*7,entries,values)); 536 } 537 ierr = PetscFree2(entries,values);CHKERRQ(ierr); 538 } 539 PetscStackCallStandard(HYPRE_SStructMatrixAssemble,(ex->ss_mat)); 540 PetscFunctionReturn(0); 541 } 542 543 544 static PetscErrorCode MatSetUp_HYPRESStruct(Mat mat) 545 { 546 PetscErrorCode ierr; 547 Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data; 548 PetscInt dim,dof,sw[3],nx,ny,nz; 549 PetscInt ilower[3],iupper[3],ssize,i; 550 DMBoundaryType px,py,pz; 551 DMDAStencilType st; 552 PetscInt nparts= 1; /* assuming only one part */ 553 PetscInt part = 0; 554 ISLocalToGlobalMapping ltog; 555 DM da; 556 557 PetscFunctionBegin; 558 ierr = MatGetDM(mat,(DM*)&da);CHKERRQ(ierr); 559 ex->da = da; 560 ierr = PetscObjectReference((PetscObject)da);CHKERRQ(ierr); 561 562 ierr = DMDAGetInfo(ex->da,&dim,0,0,0,0,0,0,&dof,&sw[0],&px,&py,&pz,&st);CHKERRQ(ierr); 563 ierr = DMDAGetCorners(ex->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);CHKERRQ(ierr); 564 iupper[0] += ilower[0] - 1; 565 iupper[1] += ilower[1] - 1; 566 iupper[2] += ilower[2] - 1; 567 /* the hypre_Box is used to zero out the matrix entries in MatZeroValues() */ 568 ex->hbox.imin[0] = ilower[0]; 569 ex->hbox.imin[1] = ilower[1]; 570 ex->hbox.imin[2] = ilower[2]; 571 ex->hbox.imax[0] = iupper[0]; 572 ex->hbox.imax[1] = iupper[1]; 573 ex->hbox.imax[2] = iupper[2]; 574 575 ex->dofs_order = 0; 576 577 /* assuming that the same number of dofs on each gridpoint. Also assume all cell-centred based */ 578 ex->nvars= dof; 579 580 /* create the hypre grid object and set its information */ 581 if (px || py || pz) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add periodic support by calling HYPRE_SStructGridSetPeriodic()"); 582 PetscStackCallStandard(HYPRE_SStructGridCreate,(ex->hcomm,dim,nparts,&ex->ss_grid)); 583 PetscStackCallStandard(HYPRE_SStructGridSetExtents,(ex->ss_grid,part,ex->hbox.imin,ex->hbox.imax)); 584 { 585 HYPRE_SStructVariable *vartypes; 586 ierr = PetscMalloc1(ex->nvars,&vartypes);CHKERRQ(ierr); 587 for (i= 0; i< ex->nvars; i++) vartypes[i]= HYPRE_SSTRUCT_VARIABLE_CELL; 588 PetscStackCallStandard(HYPRE_SStructGridSetVariables,(ex->ss_grid, part, ex->nvars,vartypes)); 589 ierr = PetscFree(vartypes);CHKERRQ(ierr); 590 } 591 PetscStackCallStandard(HYPRE_SStructGridAssemble,(ex->ss_grid)); 592 593 sw[1] = sw[0]; 594 sw[2] = sw[1]; 595 /* PetscStackCallStandard(HYPRE_SStructGridSetNumGhost,(ex->ss_grid,sw)); */ 596 597 /* create the hypre stencil object and set its information */ 598 if (sw[0] > 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add support for wider stencils"); 599 if (st == DMDA_STENCIL_BOX) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add support for box stencils"); 600 601 if (dim == 1) { 602 HYPRE_Int offsets[3][1] = {{-1},{0},{1}}; 603 PetscInt j, cnt; 604 605 ssize = 3*(ex->nvars); 606 PetscStackCallStandard(HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil)); 607 cnt= 0; 608 for (i = 0; i < (ex->nvars); i++) { 609 for (j = 0; j < 3; j++) { 610 PetscStackCallStandard(HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, offsets[j], i)); 611 cnt++; 612 } 613 } 614 } else if (dim == 2) { 615 HYPRE_Int offsets[5][2] = {{0,-1},{-1,0},{0,0},{1,0},{0,1}}; 616 PetscInt j, cnt; 617 618 ssize = 5*(ex->nvars); 619 PetscStackCallStandard(HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil)); 620 cnt= 0; 621 for (i= 0; i< (ex->nvars); i++) { 622 for (j= 0; j< 5; j++) { 623 PetscStackCallStandard(HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, offsets[j], i)); 624 cnt++; 625 } 626 } 627 } else if (dim == 3) { 628 HYPRE_Int offsets[7][3] = {{0,0,-1},{0,-1,0},{-1,0,0},{0,0,0},{1,0,0},{0,1,0},{0,0,1}}; 629 PetscInt j, cnt; 630 631 ssize = 7*(ex->nvars); 632 PetscStackCallStandard(HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil)); 633 cnt= 0; 634 for (i= 0; i< (ex->nvars); i++) { 635 for (j= 0; j< 7; j++) { 636 PetscStackCallStandard(HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, offsets[j], i)); 637 cnt++; 638 } 639 } 640 } 641 642 /* create the HYPRE graph */ 643 PetscStackCallStandard(HYPRE_SStructGraphCreate,(ex->hcomm, ex->ss_grid, &(ex->ss_graph))); 644 645 /* set the stencil graph. Note that each variable has the same graph. This means that each 646 variable couples to all the other variable and with the same stencil pattern. */ 647 for (i= 0; i< (ex->nvars); i++) { 648 PetscStackCallStandard(HYPRE_SStructGraphSetStencil,(ex->ss_graph,part,i,ex->ss_stencil)); 649 } 650 PetscStackCallStandard(HYPRE_SStructGraphAssemble,(ex->ss_graph)); 651 652 /* create the HYPRE sstruct vectors for rhs and solution */ 653 PetscStackCallStandard(HYPRE_SStructVectorCreate,(ex->hcomm,ex->ss_grid,&ex->ss_b)); 654 PetscStackCallStandard(HYPRE_SStructVectorCreate,(ex->hcomm,ex->ss_grid,&ex->ss_x)); 655 PetscStackCallStandard(HYPRE_SStructVectorInitialize,(ex->ss_b)); 656 PetscStackCallStandard(HYPRE_SStructVectorInitialize,(ex->ss_x)); 657 PetscStackCallStandard(HYPRE_SStructVectorAssemble,(ex->ss_b)); 658 PetscStackCallStandard(HYPRE_SStructVectorAssemble,(ex->ss_x)); 659 660 /* create the hypre matrix object and set its information */ 661 PetscStackCallStandard(HYPRE_SStructMatrixCreate,(ex->hcomm,ex->ss_graph,&ex->ss_mat)); 662 PetscStackCallStandard(HYPRE_SStructGridDestroy,(ex->ss_grid)); 663 PetscStackCallStandard(HYPRE_SStructStencilDestroy,(ex->ss_stencil)); 664 if (ex->needsinitialization) { 665 PetscStackCallStandard(HYPRE_SStructMatrixInitialize,(ex->ss_mat)); 666 ex->needsinitialization = PETSC_FALSE; 667 } 668 669 /* set the global and local sizes of the matrix */ 670 ierr = DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);CHKERRQ(ierr); 671 ierr = MatSetSizes(mat,dof*nx*ny*nz,dof*nx*ny*nz,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 672 ierr = PetscLayoutSetBlockSize(mat->rmap,dof);CHKERRQ(ierr); 673 ierr = PetscLayoutSetBlockSize(mat->cmap,dof);CHKERRQ(ierr); 674 ierr = PetscLayoutSetUp(mat->rmap);CHKERRQ(ierr); 675 ierr = PetscLayoutSetUp(mat->cmap);CHKERRQ(ierr); 676 mat->preallocated = PETSC_TRUE; 677 678 if (dim == 3) { 679 mat->ops->setvalueslocal = MatSetValuesLocal_HYPRESStruct_3d; 680 mat->ops->zerorowslocal = MatZeroRowsLocal_HYPRESStruct_3d; 681 mat->ops->zeroentries = MatZeroEntries_HYPRESStruct_3d; 682 683 /* ierr = MatZeroEntries_HYPRESStruct_3d(mat);CHKERRQ(ierr); */ 684 } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Only support for 3d DMDA currently"); 685 686 /* get values that will be used repeatedly in MatSetValuesLocal() and MatZeroRowsLocal() repeatedly */ 687 ierr = MatGetOwnershipRange(mat,&ex->rstart,NULL);CHKERRQ(ierr); 688 ierr = DMGetLocalToGlobalMapping(ex->da,<og);CHKERRQ(ierr); 689 ierr = ISLocalToGlobalMappingGetIndices(ltog, (const PetscInt **) &ex->gindices);CHKERRQ(ierr); 690 ierr = DMDAGetGhostCorners(ex->da,0,0,0,&ex->gnx,&ex->gnxgny,&ex->gnxgnygnz);CHKERRQ(ierr); 691 692 ex->gnxgny *= ex->gnx; 693 ex->gnxgnygnz *= ex->gnxgny; 694 695 ierr = DMDAGetCorners(ex->da,&ex->xs,&ex->ys,&ex->zs,&ex->nx,&ex->ny,&ex->nz);CHKERRQ(ierr); 696 697 ex->nxny = ex->nx*ex->ny; 698 ex->nxnynz = ex->nz*ex->nxny; 699 PetscFunctionReturn(0); 700 } 701 702 PetscErrorCode MatMult_HYPRESStruct(Mat A,Vec x,Vec y) 703 { 704 PetscErrorCode ierr; 705 const PetscScalar *xx; 706 PetscScalar *yy; 707 HYPRE_Int hlower[3],hupper[3]; 708 PetscInt ilower[3],iupper[3]; 709 Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(A->data); 710 PetscInt ordering= mx->dofs_order; 711 PetscInt nvars = mx->nvars; 712 PetscInt part = 0; 713 PetscInt size; 714 PetscInt i; 715 716 PetscFunctionBegin; 717 ierr = DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);CHKERRQ(ierr); 718 719 /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */ 720 iupper[0] += ilower[0] - 1; 721 iupper[1] += ilower[1] - 1; 722 iupper[2] += ilower[2] - 1; 723 hlower[0] = (HYPRE_Int)ilower[0]; 724 hlower[1] = (HYPRE_Int)ilower[1]; 725 hlower[2] = (HYPRE_Int)ilower[2]; 726 hupper[0] = (HYPRE_Int)iupper[0]; 727 hupper[1] = (HYPRE_Int)iupper[1]; 728 hupper[2] = (HYPRE_Int)iupper[2]; 729 730 size= 1; 731 for (i = 0; i < 3; i++) size *= (iupper[i]-ilower[i]+1); 732 733 /* copy x values over to hypre for variable ordering */ 734 if (ordering) { 735 PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0)); 736 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 737 for (i= 0; i< nvars; i++) { 738 PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,hlower,hupper,i,(HYPRE_Complex*)(xx+(size*i)))); 739 } 740 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 741 PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b)); 742 PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x)); 743 744 /* copy solution values back to PETSc */ 745 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 746 for (i= 0; i< nvars; i++) { 747 PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,hlower,hupper,i,(HYPRE_Complex*)(yy+(size*i)))); 748 } 749 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 750 } else { /* nodal ordering must be mapped to variable ordering for sys_pfmg */ 751 PetscScalar *z; 752 PetscInt j, k; 753 754 ierr = PetscMalloc1(nvars*size,&z);CHKERRQ(ierr); 755 PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0)); 756 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 757 758 /* transform nodal to hypre's variable ordering for sys_pfmg */ 759 for (i= 0; i< size; i++) { 760 k= i*nvars; 761 for (j= 0; j< nvars; j++) z[j*size+i]= xx[k+j]; 762 } 763 for (i= 0; i< nvars; i++) { 764 PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,hlower,hupper,i,(HYPRE_Complex*)(z+(size*i)))); 765 } 766 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 767 PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b)); 768 PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x)); 769 770 /* copy solution values back to PETSc */ 771 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 772 for (i= 0; i< nvars; i++) { 773 PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,hlower,hupper,i,(HYPRE_Complex*)(z+(size*i)))); 774 } 775 /* transform hypre's variable ordering for sys_pfmg to nodal ordering */ 776 for (i= 0; i< size; i++) { 777 k= i*nvars; 778 for (j= 0; j< nvars; j++) yy[k+j]= z[j*size+i]; 779 } 780 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 781 ierr = PetscFree(z);CHKERRQ(ierr); 782 } 783 PetscFunctionReturn(0); 784 } 785 786 PetscErrorCode MatAssemblyEnd_HYPRESStruct(Mat mat,MatAssemblyType mode) 787 { 788 Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data; 789 790 PetscFunctionBegin; 791 PetscStackCallStandard(HYPRE_SStructMatrixAssemble,(ex->ss_mat)); 792 PetscFunctionReturn(0); 793 } 794 795 PetscErrorCode MatZeroEntries_HYPRESStruct(Mat mat) 796 { 797 PetscFunctionBegin; 798 /* before the DMDA is set to the matrix the zero doesn't need to do anything */ 799 PetscFunctionReturn(0); 800 } 801 802 PetscErrorCode MatDestroy_HYPRESStruct(Mat mat) 803 { 804 Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data; 805 PetscErrorCode ierr; 806 ISLocalToGlobalMapping ltog; 807 808 PetscFunctionBegin; 809 ierr = DMGetLocalToGlobalMapping(ex->da,<og);CHKERRQ(ierr); 810 ierr = ISLocalToGlobalMappingRestoreIndices(ltog, (const PetscInt **) &ex->gindices);CHKERRQ(ierr); 811 PetscStackCallStandard(HYPRE_SStructGraphDestroy,(ex->ss_graph)); 812 PetscStackCallStandard(HYPRE_SStructMatrixDestroy,(ex->ss_mat)); 813 PetscStackCallStandard(HYPRE_SStructVectorDestroy,(ex->ss_x)); 814 PetscStackCallStandard(HYPRE_SStructVectorDestroy,(ex->ss_b)); 815 ierr = PetscObjectDereference((PetscObject)ex->da);CHKERRQ(ierr); 816 ierr = MPI_Comm_free(&(ex->hcomm));CHKERRQ(ierr); 817 ierr = PetscFree(ex);CHKERRQ(ierr); 818 PetscFunctionReturn(0); 819 } 820 821 PETSC_EXTERN PetscErrorCode MatCreate_HYPRESStruct(Mat B) 822 { 823 Mat_HYPRESStruct *ex; 824 PetscErrorCode ierr; 825 826 PetscFunctionBegin; 827 ierr = PetscNewLog(B,&ex);CHKERRQ(ierr); 828 B->data = (void*)ex; 829 B->rmap->bs = 1; 830 B->assembled = PETSC_FALSE; 831 832 B->insertmode = NOT_SET_VALUES; 833 834 B->ops->assemblyend = MatAssemblyEnd_HYPRESStruct; 835 B->ops->mult = MatMult_HYPRESStruct; 836 B->ops->zeroentries = MatZeroEntries_HYPRESStruct; 837 B->ops->destroy = MatDestroy_HYPRESStruct; 838 B->ops->setup = MatSetUp_HYPRESStruct; 839 840 ex->needsinitialization = PETSC_TRUE; 841 842 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)B),&(ex->hcomm));CHKERRQ(ierr); 843 ierr = PetscObjectChangeTypeName((PetscObject)B,MATHYPRESSTRUCT);CHKERRQ(ierr); 844 PetscFunctionReturn(0); 845 } 846 847