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