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