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