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