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