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