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