147c6ae99SBarry Smith 2af0996ceSBarry Smith #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/ 307475bc1SBarry Smith #include <petscmat.h> 447c6ae99SBarry Smith 5e727c939SJed Brown extern PetscErrorCode DMCreateColoring_DA_1d_MPIAIJ(DM,ISColoringType,ISColoring*); 6e727c939SJed Brown extern PetscErrorCode DMCreateColoring_DA_2d_MPIAIJ(DM,ISColoringType,ISColoring*); 7e727c939SJed Brown extern PetscErrorCode DMCreateColoring_DA_2d_5pt_MPIAIJ(DM,ISColoringType,ISColoring*); 8e727c939SJed Brown extern PetscErrorCode DMCreateColoring_DA_3d_MPIAIJ(DM,ISColoringType,ISColoring*); 947c6ae99SBarry Smith 1047c6ae99SBarry Smith /* 1147c6ae99SBarry Smith For ghost i that may be negative or greater than the upper bound this 1247c6ae99SBarry Smith maps it into the 0:m-1 range using periodicity 1347c6ae99SBarry Smith */ 1447c6ae99SBarry Smith #define SetInRange(i,m) ((i < 0) ? m+i : ((i >= m) ? i-m : i)) 1547c6ae99SBarry Smith 16ce308e1dSBarry Smith static PetscErrorCode DMDASetBlockFills_Private(const PetscInt *dfill,PetscInt w,PetscInt **rfill) 1747c6ae99SBarry Smith { 1847c6ae99SBarry Smith PetscErrorCode ierr; 1947c6ae99SBarry Smith PetscInt i,j,nz,*fill; 2047c6ae99SBarry Smith 2147c6ae99SBarry Smith PetscFunctionBegin; 2247c6ae99SBarry Smith if (!dfill) PetscFunctionReturn(0); 2347c6ae99SBarry Smith 2447c6ae99SBarry Smith /* count number nonzeros */ 2547c6ae99SBarry Smith nz = 0; 2647c6ae99SBarry Smith for (i=0; i<w; i++) { 2747c6ae99SBarry Smith for (j=0; j<w; j++) { 2847c6ae99SBarry Smith if (dfill[w*i+j]) nz++; 2947c6ae99SBarry Smith } 3047c6ae99SBarry Smith } 31854ce69bSBarry Smith ierr = PetscMalloc1(nz + w + 1,&fill);CHKERRQ(ierr); 3247c6ae99SBarry Smith /* construct modified CSR storage of nonzero structure */ 33ce308e1dSBarry Smith /* fill[0 -- w] marks starts of each row of column indices (and end of last row) 34ce308e1dSBarry Smith so fill[1] - fill[0] gives number of nonzeros in first row etc */ 3547c6ae99SBarry Smith nz = w + 1; 3647c6ae99SBarry Smith for (i=0; i<w; i++) { 3747c6ae99SBarry Smith fill[i] = nz; 3847c6ae99SBarry Smith for (j=0; j<w; j++) { 3947c6ae99SBarry Smith if (dfill[w*i+j]) { 4047c6ae99SBarry Smith fill[nz] = j; 4147c6ae99SBarry Smith nz++; 4247c6ae99SBarry Smith } 4347c6ae99SBarry Smith } 4447c6ae99SBarry Smith } 4547c6ae99SBarry Smith fill[w] = nz; 4647c6ae99SBarry Smith 4747c6ae99SBarry Smith *rfill = fill; 4847c6ae99SBarry Smith PetscFunctionReturn(0); 4947c6ae99SBarry Smith } 5047c6ae99SBarry Smith 5147c6ae99SBarry Smith /*@ 52aa219208SBarry Smith DMDASetBlockFills - Sets the fill pattern in each block for a multi-component problem 53950540a4SJed Brown of the matrix returned by DMCreateMatrix(). 5447c6ae99SBarry Smith 55aa219208SBarry Smith Logically Collective on DMDA 5647c6ae99SBarry Smith 5747c6ae99SBarry Smith Input Parameter: 5847c6ae99SBarry Smith + da - the distributed array 590298fd71SBarry Smith . dfill - the fill pattern in the diagonal block (may be NULL, means use dense block) 6047c6ae99SBarry Smith - ofill - the fill pattern in the off-diagonal blocks 6147c6ae99SBarry Smith 6247c6ae99SBarry Smith 6347c6ae99SBarry Smith Level: developer 6447c6ae99SBarry Smith 6547c6ae99SBarry Smith Notes: This only makes sense when you are doing multicomponent problems but using the 6647c6ae99SBarry Smith MPIAIJ matrix format 6747c6ae99SBarry Smith 6847c6ae99SBarry Smith The format for dfill and ofill is a 2 dimensional dof by dof matrix with 1 entries 6947c6ae99SBarry Smith representing coupling and 0 entries for missing coupling. For example 7047c6ae99SBarry Smith $ dfill[9] = {1, 0, 0, 7147c6ae99SBarry Smith $ 1, 1, 0, 7247c6ae99SBarry Smith $ 0, 1, 1} 7347c6ae99SBarry Smith means that row 0 is coupled with only itself in the diagonal block, row 1 is coupled with 7447c6ae99SBarry Smith itself and row 0 (in the diagonal block) and row 2 is coupled with itself and row 1 (in the 7547c6ae99SBarry Smith diagonal block). 7647c6ae99SBarry Smith 77aa219208SBarry Smith DMDASetGetMatrix() allows you to provide general code for those more complicated nonzero patterns then 7847c6ae99SBarry Smith can be represented in the dfill, ofill format 7947c6ae99SBarry Smith 8047c6ae99SBarry Smith Contributed by Glenn Hammond 8147c6ae99SBarry Smith 828ddb5d8bSBarry Smith .seealso DMCreateMatrix(), DMDASetGetMatrix(), DMSetMatrixPreallocateOnly() 8347c6ae99SBarry Smith 8447c6ae99SBarry Smith @*/ 85ce308e1dSBarry Smith PetscErrorCode DMDASetBlockFills(DM da,const PetscInt *dfill,const PetscInt *ofill) 8647c6ae99SBarry Smith { 8747c6ae99SBarry Smith DM_DA *dd = (DM_DA*)da->data; 8847c6ae99SBarry Smith PetscErrorCode ierr; 89ae4f298aSBarry Smith PetscInt i,k,cnt = 1; 9047c6ae99SBarry Smith 9147c6ae99SBarry Smith PetscFunctionBegin; 92aa219208SBarry Smith ierr = DMDASetBlockFills_Private(dfill,dd->w,&dd->dfill);CHKERRQ(ierr); 93aa219208SBarry Smith ierr = DMDASetBlockFills_Private(ofill,dd->w,&dd->ofill);CHKERRQ(ierr); 94ae4f298aSBarry Smith 95ae4f298aSBarry Smith /* ofillcount tracks the columns of ofill that have any nonzero in thems; the value in each location is the number of 96ae4f298aSBarry Smith columns to the left with any nonzeros in them plus 1 */ 971795a4d1SJed Brown ierr = PetscCalloc1(dd->w,&dd->ofillcols);CHKERRQ(ierr); 98ae4f298aSBarry Smith for (i=0; i<dd->w; i++) { 99ae4f298aSBarry Smith for (k=dd->ofill[i]; k<dd->ofill[i+1]; k++) dd->ofillcols[dd->ofill[k]] = 1; 100ae4f298aSBarry Smith } 101ae4f298aSBarry Smith for (i=0; i<dd->w; i++) { 102ae4f298aSBarry Smith if (dd->ofillcols[i]) { 103ae4f298aSBarry Smith dd->ofillcols[i] = cnt++; 104ae4f298aSBarry Smith } 105ae4f298aSBarry Smith } 10647c6ae99SBarry Smith PetscFunctionReturn(0); 10747c6ae99SBarry Smith } 10847c6ae99SBarry Smith 10947c6ae99SBarry Smith 110b412c318SBarry Smith PetscErrorCode DMCreateColoring_DA(DM da,ISColoringType ctype,ISColoring *coloring) 11147c6ae99SBarry Smith { 11247c6ae99SBarry Smith PetscErrorCode ierr; 11347c6ae99SBarry Smith PetscInt dim,m,n,p,nc; 114bff4a2f0SMatthew G. Knepley DMBoundaryType bx,by,bz; 11547c6ae99SBarry Smith MPI_Comm comm; 11647c6ae99SBarry Smith PetscMPIInt size; 11747c6ae99SBarry Smith PetscBool isBAIJ; 11847c6ae99SBarry Smith DM_DA *dd = (DM_DA*)da->data; 11947c6ae99SBarry Smith 12047c6ae99SBarry Smith PetscFunctionBegin; 12147c6ae99SBarry Smith /* 12247c6ae99SBarry Smith m 12347c6ae99SBarry Smith ------------------------------------------------------ 12447c6ae99SBarry Smith | | 12547c6ae99SBarry Smith | | 12647c6ae99SBarry Smith | ---------------------- | 12747c6ae99SBarry Smith | | | | 12847c6ae99SBarry Smith n | yn | | | 12947c6ae99SBarry Smith | | | | 13047c6ae99SBarry Smith | .--------------------- | 13147c6ae99SBarry Smith | (xs,ys) xn | 13247c6ae99SBarry Smith | . | 13347c6ae99SBarry Smith | (gxs,gys) | 13447c6ae99SBarry Smith | | 13547c6ae99SBarry Smith ----------------------------------------------------- 13647c6ae99SBarry Smith */ 13747c6ae99SBarry Smith 13847c6ae99SBarry Smith /* 13947c6ae99SBarry Smith nc - number of components per grid point 14047c6ae99SBarry Smith col - number of colors needed in one direction for single component problem 14147c6ae99SBarry Smith 14247c6ae99SBarry Smith */ 1431321219cSEthan Coon ierr = DMDAGetInfo(da,&dim,0,0,0,&m,&n,&p,&nc,0,&bx,&by,&bz,0);CHKERRQ(ierr); 14447c6ae99SBarry Smith 14547c6ae99SBarry Smith ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 14647c6ae99SBarry Smith ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1475bdb020cSBarry Smith if (ctype == IS_COLORING_LOCAL) { 14847c6ae99SBarry Smith if (size == 1) { 14947c6ae99SBarry Smith ctype = IS_COLORING_GLOBAL; 15047c6ae99SBarry Smith } else if (dim > 1) { 151bff4a2f0SMatthew G. Knepley if ((m==1 && bx == DM_BOUNDARY_PERIODIC) || (n==1 && by == DM_BOUNDARY_PERIODIC) || (p==1 && bz == DM_BOUNDARY_PERIODIC)) { 1525bdb020cSBarry Smith SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"IS_COLORING_LOCAL cannot be used for periodic boundary condition having both ends of the domain on the same process"); 15347c6ae99SBarry Smith } 15447c6ae99SBarry Smith } 15547c6ae99SBarry Smith } 15647c6ae99SBarry Smith 157aa219208SBarry Smith /* Tell the DMDA it has 1 degree of freedom per grid point so that the coloring for BAIJ 15847c6ae99SBarry Smith matrices is for the blocks, not the individual matrix elements */ 159b412c318SBarry Smith ierr = PetscStrcmp(da->mattype,MATBAIJ,&isBAIJ);CHKERRQ(ierr); 160b412c318SBarry Smith if (!isBAIJ) {ierr = PetscStrcmp(da->mattype,MATMPIBAIJ,&isBAIJ);CHKERRQ(ierr);} 161b412c318SBarry Smith if (!isBAIJ) {ierr = PetscStrcmp(da->mattype,MATSEQBAIJ,&isBAIJ);CHKERRQ(ierr);} 16247c6ae99SBarry Smith if (isBAIJ) { 16347c6ae99SBarry Smith dd->w = 1; 16447c6ae99SBarry Smith dd->xs = dd->xs/nc; 16547c6ae99SBarry Smith dd->xe = dd->xe/nc; 16647c6ae99SBarry Smith dd->Xs = dd->Xs/nc; 16747c6ae99SBarry Smith dd->Xe = dd->Xe/nc; 16847c6ae99SBarry Smith } 16947c6ae99SBarry Smith 17047c6ae99SBarry Smith /* 171aa219208SBarry Smith We do not provide a getcoloring function in the DMDA operations because 172aa219208SBarry Smith the basic DMDA does not know about matrices. We think of DMDA as being more 17347c6ae99SBarry Smith more low-level then matrices. 17447c6ae99SBarry Smith */ 17547c6ae99SBarry Smith if (dim == 1) { 176e727c939SJed Brown ierr = DMCreateColoring_DA_1d_MPIAIJ(da,ctype,coloring);CHKERRQ(ierr); 17747c6ae99SBarry Smith } else if (dim == 2) { 178e727c939SJed Brown ierr = DMCreateColoring_DA_2d_MPIAIJ(da,ctype,coloring);CHKERRQ(ierr); 17947c6ae99SBarry Smith } else if (dim == 3) { 180e727c939SJed Brown ierr = DMCreateColoring_DA_3d_MPIAIJ(da,ctype,coloring);CHKERRQ(ierr); 181ce94432eSBarry Smith } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Not done for %D dimension, send us mail petsc-maint@mcs.anl.gov for code",dim); 18247c6ae99SBarry Smith if (isBAIJ) { 18347c6ae99SBarry Smith dd->w = nc; 18447c6ae99SBarry Smith dd->xs = dd->xs*nc; 18547c6ae99SBarry Smith dd->xe = dd->xe*nc; 18647c6ae99SBarry Smith dd->Xs = dd->Xs*nc; 18747c6ae99SBarry Smith dd->Xe = dd->Xe*nc; 18847c6ae99SBarry Smith } 18947c6ae99SBarry Smith PetscFunctionReturn(0); 19047c6ae99SBarry Smith } 19147c6ae99SBarry Smith 19247c6ae99SBarry Smith /* ---------------------------------------------------------------------------------*/ 19347c6ae99SBarry Smith 194e727c939SJed Brown PetscErrorCode DMCreateColoring_DA_2d_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring) 19547c6ae99SBarry Smith { 19647c6ae99SBarry Smith PetscErrorCode ierr; 19747c6ae99SBarry Smith PetscInt xs,ys,nx,ny,i,j,ii,gxs,gys,gnx,gny,m,n,M,N,dim,s,k,nc,col; 19847c6ae99SBarry Smith PetscInt ncolors; 19947c6ae99SBarry Smith MPI_Comm comm; 200bff4a2f0SMatthew G. Knepley DMBoundaryType bx,by; 201aa219208SBarry Smith DMDAStencilType st; 20247c6ae99SBarry Smith ISColoringValue *colors; 20347c6ae99SBarry Smith DM_DA *dd = (DM_DA*)da->data; 20447c6ae99SBarry Smith 20547c6ae99SBarry Smith PetscFunctionBegin; 20647c6ae99SBarry Smith /* 20747c6ae99SBarry Smith nc - number of components per grid point 20847c6ae99SBarry Smith col - number of colors needed in one direction for single component problem 20947c6ae99SBarry Smith 21047c6ae99SBarry Smith */ 2111321219cSEthan Coon ierr = DMDAGetInfo(da,&dim,&m,&n,0,&M,&N,0,&nc,&s,&bx,&by,0,&st);CHKERRQ(ierr); 21247c6ae99SBarry Smith col = 2*s + 1; 213aa219208SBarry Smith ierr = DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);CHKERRQ(ierr); 214aa219208SBarry Smith ierr = DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);CHKERRQ(ierr); 21547c6ae99SBarry Smith ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 21647c6ae99SBarry Smith 21747c6ae99SBarry Smith /* special case as taught to us by Paul Hovland */ 218aa219208SBarry Smith if (st == DMDA_STENCIL_STAR && s == 1) { 219e727c939SJed Brown ierr = DMCreateColoring_DA_2d_5pt_MPIAIJ(da,ctype,coloring);CHKERRQ(ierr); 22047c6ae99SBarry Smith } else { 22147c6ae99SBarry Smith 222bff4a2f0SMatthew G. Knepley if (bx == DM_BOUNDARY_PERIODIC && (m % col)) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X (%d) is divisible\n\ 22347c6ae99SBarry Smith by 2*stencil_width + 1 (%d)\n", m, col); 224bff4a2f0SMatthew G. Knepley if (by == DM_BOUNDARY_PERIODIC && (n % col)) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y (%d) is divisible\n\ 22547c6ae99SBarry Smith by 2*stencil_width + 1 (%d)\n", n, col); 22647c6ae99SBarry Smith if (ctype == IS_COLORING_GLOBAL) { 22747c6ae99SBarry Smith if (!dd->localcoloring) { 228785e854fSJed Brown ierr = PetscMalloc1(nc*nx*ny,&colors);CHKERRQ(ierr); 22947c6ae99SBarry Smith ii = 0; 23047c6ae99SBarry Smith for (j=ys; j<ys+ny; j++) { 23147c6ae99SBarry Smith for (i=xs; i<xs+nx; i++) { 23247c6ae99SBarry Smith for (k=0; k<nc; k++) { 23347c6ae99SBarry Smith colors[ii++] = k + nc*((i % col) + col*(j % col)); 23447c6ae99SBarry Smith } 23547c6ae99SBarry Smith } 23647c6ae99SBarry Smith } 23747c6ae99SBarry Smith ncolors = nc + nc*(col-1 + col*(col-1)); 238aaf3ff59SMatthew G. Knepley ierr = ISColoringCreate(comm,ncolors,nc*nx*ny,colors,PETSC_OWN_POINTER,&dd->localcoloring);CHKERRQ(ierr); 23947c6ae99SBarry Smith } 24047c6ae99SBarry Smith *coloring = dd->localcoloring; 2415bdb020cSBarry Smith } else if (ctype == IS_COLORING_LOCAL) { 24247c6ae99SBarry Smith if (!dd->ghostedcoloring) { 243785e854fSJed Brown ierr = PetscMalloc1(nc*gnx*gny,&colors);CHKERRQ(ierr); 24447c6ae99SBarry Smith ii = 0; 24547c6ae99SBarry Smith for (j=gys; j<gys+gny; j++) { 24647c6ae99SBarry Smith for (i=gxs; i<gxs+gnx; i++) { 24747c6ae99SBarry Smith for (k=0; k<nc; k++) { 24847c6ae99SBarry Smith /* the complicated stuff is to handle periodic boundaries */ 24947c6ae99SBarry Smith colors[ii++] = k + nc*((SetInRange(i,m) % col) + col*(SetInRange(j,n) % col)); 25047c6ae99SBarry Smith } 25147c6ae99SBarry Smith } 25247c6ae99SBarry Smith } 25347c6ae99SBarry Smith ncolors = nc + nc*(col - 1 + col*(col-1)); 254aaf3ff59SMatthew G. Knepley ierr = ISColoringCreate(comm,ncolors,nc*gnx*gny,colors,PETSC_OWN_POINTER,&dd->ghostedcoloring);CHKERRQ(ierr); 25547c6ae99SBarry Smith /* PetscIntView(ncolors,(PetscInt*)colors,0); */ 25647c6ae99SBarry Smith 2575bdb020cSBarry Smith ierr = ISColoringSetType(dd->ghostedcoloring,IS_COLORING_LOCAL);CHKERRQ(ierr); 25847c6ae99SBarry Smith } 25947c6ae99SBarry Smith *coloring = dd->ghostedcoloring; 260ce94432eSBarry Smith } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype); 26147c6ae99SBarry Smith } 26247c6ae99SBarry Smith ierr = ISColoringReference(*coloring);CHKERRQ(ierr); 26347c6ae99SBarry Smith PetscFunctionReturn(0); 26447c6ae99SBarry Smith } 26547c6ae99SBarry Smith 26647c6ae99SBarry Smith /* ---------------------------------------------------------------------------------*/ 26747c6ae99SBarry Smith 268e727c939SJed Brown PetscErrorCode DMCreateColoring_DA_3d_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring) 26947c6ae99SBarry Smith { 27047c6ae99SBarry Smith PetscErrorCode ierr; 27147c6ae99SBarry Smith PetscInt xs,ys,nx,ny,i,j,gxs,gys,gnx,gny,m,n,p,dim,s,k,nc,col,zs,gzs,ii,l,nz,gnz,M,N,P; 27247c6ae99SBarry Smith PetscInt ncolors; 27347c6ae99SBarry Smith MPI_Comm comm; 274bff4a2f0SMatthew G. Knepley DMBoundaryType bx,by,bz; 275aa219208SBarry Smith DMDAStencilType st; 27647c6ae99SBarry Smith ISColoringValue *colors; 27747c6ae99SBarry Smith DM_DA *dd = (DM_DA*)da->data; 27847c6ae99SBarry Smith 27947c6ae99SBarry Smith PetscFunctionBegin; 28047c6ae99SBarry Smith /* 28147c6ae99SBarry Smith nc - number of components per grid point 28247c6ae99SBarry Smith col - number of colors needed in one direction for single component problem 28347c6ae99SBarry Smith 28447c6ae99SBarry Smith */ 2851321219cSEthan Coon ierr = DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&bx,&by,&bz,&st);CHKERRQ(ierr); 28647c6ae99SBarry Smith col = 2*s + 1; 287bff4a2f0SMatthew G. Knepley if (bx == DM_BOUNDARY_PERIODIC && (m % col)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisible\n\ 28847c6ae99SBarry Smith by 2*stencil_width + 1\n"); 289bff4a2f0SMatthew G. Knepley if (by == DM_BOUNDARY_PERIODIC && (n % col)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisible\n\ 29047c6ae99SBarry Smith by 2*stencil_width + 1\n"); 291bff4a2f0SMatthew G. Knepley if (bz == DM_BOUNDARY_PERIODIC && (p % col)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Z is divisible\n\ 29247c6ae99SBarry Smith by 2*stencil_width + 1\n"); 29347c6ae99SBarry Smith 294aa219208SBarry Smith ierr = DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);CHKERRQ(ierr); 295aa219208SBarry Smith ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);CHKERRQ(ierr); 29647c6ae99SBarry Smith ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 29747c6ae99SBarry Smith 29847c6ae99SBarry Smith /* create the coloring */ 29947c6ae99SBarry Smith if (ctype == IS_COLORING_GLOBAL) { 30047c6ae99SBarry Smith if (!dd->localcoloring) { 301785e854fSJed Brown ierr = PetscMalloc1(nc*nx*ny*nz,&colors);CHKERRQ(ierr); 30247c6ae99SBarry Smith ii = 0; 30347c6ae99SBarry Smith for (k=zs; k<zs+nz; k++) { 30447c6ae99SBarry Smith for (j=ys; j<ys+ny; j++) { 30547c6ae99SBarry Smith for (i=xs; i<xs+nx; i++) { 30647c6ae99SBarry Smith for (l=0; l<nc; l++) { 30747c6ae99SBarry Smith colors[ii++] = l + nc*((i % col) + col*(j % col) + col*col*(k % col)); 30847c6ae99SBarry Smith } 30947c6ae99SBarry Smith } 31047c6ae99SBarry Smith } 31147c6ae99SBarry Smith } 31247c6ae99SBarry Smith ncolors = nc + nc*(col-1 + col*(col-1)+ col*col*(col-1)); 313aaf3ff59SMatthew G. Knepley ierr = ISColoringCreate(comm,ncolors,nc*nx*ny*nz,colors,PETSC_OWN_POINTER,&dd->localcoloring);CHKERRQ(ierr); 31447c6ae99SBarry Smith } 31547c6ae99SBarry Smith *coloring = dd->localcoloring; 3165bdb020cSBarry Smith } else if (ctype == IS_COLORING_LOCAL) { 31747c6ae99SBarry Smith if (!dd->ghostedcoloring) { 318785e854fSJed Brown ierr = PetscMalloc1(nc*gnx*gny*gnz,&colors);CHKERRQ(ierr); 31947c6ae99SBarry Smith ii = 0; 32047c6ae99SBarry Smith for (k=gzs; k<gzs+gnz; k++) { 32147c6ae99SBarry Smith for (j=gys; j<gys+gny; j++) { 32247c6ae99SBarry Smith for (i=gxs; i<gxs+gnx; i++) { 32347c6ae99SBarry Smith for (l=0; l<nc; l++) { 32447c6ae99SBarry Smith /* the complicated stuff is to handle periodic boundaries */ 32547c6ae99SBarry Smith colors[ii++] = l + nc*((SetInRange(i,m) % col) + col*(SetInRange(j,n) % col) + col*col*(SetInRange(k,p) % col)); 32647c6ae99SBarry Smith } 32747c6ae99SBarry Smith } 32847c6ae99SBarry Smith } 32947c6ae99SBarry Smith } 33047c6ae99SBarry Smith ncolors = nc + nc*(col-1 + col*(col-1)+ col*col*(col-1)); 331aaf3ff59SMatthew G. Knepley ierr = ISColoringCreate(comm,ncolors,nc*gnx*gny*gnz,colors,PETSC_OWN_POINTER,&dd->ghostedcoloring);CHKERRQ(ierr); 3325bdb020cSBarry Smith ierr = ISColoringSetType(dd->ghostedcoloring,IS_COLORING_LOCAL);CHKERRQ(ierr); 33347c6ae99SBarry Smith } 33447c6ae99SBarry Smith *coloring = dd->ghostedcoloring; 335ce94432eSBarry Smith } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype); 33647c6ae99SBarry Smith ierr = ISColoringReference(*coloring);CHKERRQ(ierr); 33747c6ae99SBarry Smith PetscFunctionReturn(0); 33847c6ae99SBarry Smith } 33947c6ae99SBarry Smith 34047c6ae99SBarry Smith /* ---------------------------------------------------------------------------------*/ 34147c6ae99SBarry Smith 342e727c939SJed Brown PetscErrorCode DMCreateColoring_DA_1d_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring) 34347c6ae99SBarry Smith { 34447c6ae99SBarry Smith PetscErrorCode ierr; 34547c6ae99SBarry Smith PetscInt xs,nx,i,i1,gxs,gnx,l,m,M,dim,s,nc,col; 34647c6ae99SBarry Smith PetscInt ncolors; 34747c6ae99SBarry Smith MPI_Comm comm; 348bff4a2f0SMatthew G. Knepley DMBoundaryType bx; 34947c6ae99SBarry Smith ISColoringValue *colors; 35047c6ae99SBarry Smith DM_DA *dd = (DM_DA*)da->data; 35147c6ae99SBarry Smith 35247c6ae99SBarry Smith PetscFunctionBegin; 35347c6ae99SBarry Smith /* 35447c6ae99SBarry Smith nc - number of components per grid point 35547c6ae99SBarry Smith col - number of colors needed in one direction for single component problem 35647c6ae99SBarry Smith 35747c6ae99SBarry Smith */ 3581321219cSEthan Coon ierr = DMDAGetInfo(da,&dim,&m,0,0,&M,0,0,&nc,&s,&bx,0,0,0);CHKERRQ(ierr); 35947c6ae99SBarry Smith col = 2*s + 1; 36047c6ae99SBarry Smith 361bff4a2f0SMatthew G. Knepley if (bx == DM_BOUNDARY_PERIODIC && (m % col)) SETERRQ2(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points %d is divisible\n\ 36231e6f798SBarry Smith by 2*stencil_width + 1 %d\n",(int)m,(int)col); 36347c6ae99SBarry Smith 364aa219208SBarry Smith ierr = DMDAGetCorners(da,&xs,0,0,&nx,0,0);CHKERRQ(ierr); 365aa219208SBarry Smith ierr = DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);CHKERRQ(ierr); 36647c6ae99SBarry Smith ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 36747c6ae99SBarry Smith 36847c6ae99SBarry Smith /* create the coloring */ 36947c6ae99SBarry Smith if (ctype == IS_COLORING_GLOBAL) { 37047c6ae99SBarry Smith if (!dd->localcoloring) { 371785e854fSJed Brown ierr = PetscMalloc1(nc*nx,&colors);CHKERRQ(ierr); 372ae4f298aSBarry Smith if (dd->ofillcols) { 373ae4f298aSBarry Smith PetscInt tc = 0; 374ae4f298aSBarry Smith for (i=0; i<nc; i++) tc += (PetscInt) (dd->ofillcols[i] > 0); 375ae4f298aSBarry Smith i1 = 0; 376ae4f298aSBarry Smith for (i=xs; i<xs+nx; i++) { 377ae4f298aSBarry Smith for (l=0; l<nc; l++) { 378ae4f298aSBarry Smith if (dd->ofillcols[l] && (i % col)) { 379ae4f298aSBarry Smith colors[i1++] = nc - 1 + tc*((i % col) - 1) + dd->ofillcols[l]; 380ae4f298aSBarry Smith } else { 381ae4f298aSBarry Smith colors[i1++] = l; 382ae4f298aSBarry Smith } 383ae4f298aSBarry Smith } 384ae4f298aSBarry Smith } 385ae4f298aSBarry Smith ncolors = nc + 2*s*tc; 386ae4f298aSBarry Smith } else { 38747c6ae99SBarry Smith i1 = 0; 38847c6ae99SBarry Smith for (i=xs; i<xs+nx; i++) { 38947c6ae99SBarry Smith for (l=0; l<nc; l++) { 39047c6ae99SBarry Smith colors[i1++] = l + nc*(i % col); 39147c6ae99SBarry Smith } 39247c6ae99SBarry Smith } 39347c6ae99SBarry Smith ncolors = nc + nc*(col-1); 394ae4f298aSBarry Smith } 395aaf3ff59SMatthew G. Knepley ierr = ISColoringCreate(comm,ncolors,nc*nx,colors,PETSC_OWN_POINTER,&dd->localcoloring);CHKERRQ(ierr); 39647c6ae99SBarry Smith } 39747c6ae99SBarry Smith *coloring = dd->localcoloring; 3985bdb020cSBarry Smith } else if (ctype == IS_COLORING_LOCAL) { 39947c6ae99SBarry Smith if (!dd->ghostedcoloring) { 400785e854fSJed Brown ierr = PetscMalloc1(nc*gnx,&colors);CHKERRQ(ierr); 40147c6ae99SBarry Smith i1 = 0; 40247c6ae99SBarry Smith for (i=gxs; i<gxs+gnx; i++) { 40347c6ae99SBarry Smith for (l=0; l<nc; l++) { 40447c6ae99SBarry Smith /* the complicated stuff is to handle periodic boundaries */ 40547c6ae99SBarry Smith colors[i1++] = l + nc*(SetInRange(i,m) % col); 40647c6ae99SBarry Smith } 40747c6ae99SBarry Smith } 40847c6ae99SBarry Smith ncolors = nc + nc*(col-1); 409aaf3ff59SMatthew G. Knepley ierr = ISColoringCreate(comm,ncolors,nc*gnx,colors,PETSC_OWN_POINTER,&dd->ghostedcoloring);CHKERRQ(ierr); 4105bdb020cSBarry Smith ierr = ISColoringSetType(dd->ghostedcoloring,IS_COLORING_LOCAL);CHKERRQ(ierr); 41147c6ae99SBarry Smith } 41247c6ae99SBarry Smith *coloring = dd->ghostedcoloring; 413ce94432eSBarry Smith } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype); 41447c6ae99SBarry Smith ierr = ISColoringReference(*coloring);CHKERRQ(ierr); 41547c6ae99SBarry Smith PetscFunctionReturn(0); 41647c6ae99SBarry Smith } 41747c6ae99SBarry Smith 418e727c939SJed Brown PetscErrorCode DMCreateColoring_DA_2d_5pt_MPIAIJ(DM da,ISColoringType ctype,ISColoring *coloring) 41947c6ae99SBarry Smith { 42047c6ae99SBarry Smith PetscErrorCode ierr; 42147c6ae99SBarry Smith PetscInt xs,ys,nx,ny,i,j,ii,gxs,gys,gnx,gny,m,n,dim,s,k,nc; 42247c6ae99SBarry Smith PetscInt ncolors; 42347c6ae99SBarry Smith MPI_Comm comm; 424bff4a2f0SMatthew G. Knepley DMBoundaryType bx,by; 42547c6ae99SBarry Smith ISColoringValue *colors; 42647c6ae99SBarry Smith DM_DA *dd = (DM_DA*)da->data; 42747c6ae99SBarry Smith 42847c6ae99SBarry Smith PetscFunctionBegin; 42947c6ae99SBarry Smith /* 43047c6ae99SBarry Smith nc - number of components per grid point 43147c6ae99SBarry Smith col - number of colors needed in one direction for single component problem 43247c6ae99SBarry Smith 43347c6ae99SBarry Smith */ 4341321219cSEthan Coon ierr = DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,0);CHKERRQ(ierr); 435aa219208SBarry Smith ierr = DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);CHKERRQ(ierr); 436aa219208SBarry Smith ierr = DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);CHKERRQ(ierr); 43747c6ae99SBarry Smith ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 43847c6ae99SBarry Smith 439bff4a2f0SMatthew G. Knepley if (bx == DM_BOUNDARY_PERIODIC && (m % 5)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisible by 5\n"); 440bff4a2f0SMatthew G. Knepley if (by == DM_BOUNDARY_PERIODIC && (n % 5)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisible by 5\n"); 44147c6ae99SBarry Smith 44247c6ae99SBarry Smith /* create the coloring */ 44347c6ae99SBarry Smith if (ctype == IS_COLORING_GLOBAL) { 44447c6ae99SBarry Smith if (!dd->localcoloring) { 445785e854fSJed Brown ierr = PetscMalloc1(nc*nx*ny,&colors);CHKERRQ(ierr); 44647c6ae99SBarry Smith ii = 0; 44747c6ae99SBarry Smith for (j=ys; j<ys+ny; j++) { 44847c6ae99SBarry Smith for (i=xs; i<xs+nx; i++) { 44947c6ae99SBarry Smith for (k=0; k<nc; k++) { 45047c6ae99SBarry Smith colors[ii++] = k + nc*((3*j+i) % 5); 45147c6ae99SBarry Smith } 45247c6ae99SBarry Smith } 45347c6ae99SBarry Smith } 45447c6ae99SBarry Smith ncolors = 5*nc; 455aaf3ff59SMatthew G. Knepley ierr = ISColoringCreate(comm,ncolors,nc*nx*ny,colors,PETSC_OWN_POINTER,&dd->localcoloring);CHKERRQ(ierr); 45647c6ae99SBarry Smith } 45747c6ae99SBarry Smith *coloring = dd->localcoloring; 4585bdb020cSBarry Smith } else if (ctype == IS_COLORING_LOCAL) { 45947c6ae99SBarry Smith if (!dd->ghostedcoloring) { 460785e854fSJed Brown ierr = PetscMalloc1(nc*gnx*gny,&colors);CHKERRQ(ierr); 46147c6ae99SBarry Smith ii = 0; 46247c6ae99SBarry Smith for (j=gys; j<gys+gny; j++) { 46347c6ae99SBarry Smith for (i=gxs; i<gxs+gnx; i++) { 46447c6ae99SBarry Smith for (k=0; k<nc; k++) { 46547c6ae99SBarry Smith colors[ii++] = k + nc*((3*SetInRange(j,n) + SetInRange(i,m)) % 5); 46647c6ae99SBarry Smith } 46747c6ae99SBarry Smith } 46847c6ae99SBarry Smith } 46947c6ae99SBarry Smith ncolors = 5*nc; 470aaf3ff59SMatthew G. Knepley ierr = ISColoringCreate(comm,ncolors,nc*gnx*gny,colors,PETSC_OWN_POINTER,&dd->ghostedcoloring);CHKERRQ(ierr); 4715bdb020cSBarry Smith ierr = ISColoringSetType(dd->ghostedcoloring,IS_COLORING_LOCAL);CHKERRQ(ierr); 47247c6ae99SBarry Smith } 47347c6ae99SBarry Smith *coloring = dd->ghostedcoloring; 474ce94432eSBarry Smith } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype); 47547c6ae99SBarry Smith PetscFunctionReturn(0); 47647c6ae99SBarry Smith } 47747c6ae99SBarry Smith 47847c6ae99SBarry Smith /* =========================================================================== */ 479950540a4SJed Brown extern PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ(DM,Mat); 480ce308e1dSBarry Smith extern PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ_Fill(DM,Mat); 481950540a4SJed Brown extern PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ(DM,Mat); 482950540a4SJed Brown extern PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ_Fill(DM,Mat); 483950540a4SJed Brown extern PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ(DM,Mat); 484950540a4SJed Brown extern PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ_Fill(DM,Mat); 485950540a4SJed Brown extern PetscErrorCode DMCreateMatrix_DA_2d_MPIBAIJ(DM,Mat); 486950540a4SJed Brown extern PetscErrorCode DMCreateMatrix_DA_3d_MPIBAIJ(DM,Mat); 487950540a4SJed Brown extern PetscErrorCode DMCreateMatrix_DA_2d_MPISBAIJ(DM,Mat); 488950540a4SJed Brown extern PetscErrorCode DMCreateMatrix_DA_3d_MPISBAIJ(DM,Mat); 48947c6ae99SBarry Smith 4908bbdbebaSMatthew G Knepley /*@C 491c688c046SMatthew G Knepley MatSetupDM - Sets the DMDA that is to be used by the HYPRE_StructMatrix PETSc matrix 49247c6ae99SBarry Smith 49347c6ae99SBarry Smith Logically Collective on Mat 49447c6ae99SBarry Smith 49547c6ae99SBarry Smith Input Parameters: 49647c6ae99SBarry Smith + mat - the matrix 49747c6ae99SBarry Smith - da - the da 49847c6ae99SBarry Smith 49947c6ae99SBarry Smith Level: intermediate 50047c6ae99SBarry Smith 50147c6ae99SBarry Smith @*/ 502c688c046SMatthew G Knepley PetscErrorCode MatSetupDM(Mat mat,DM da) 50347c6ae99SBarry Smith { 50447c6ae99SBarry Smith PetscErrorCode ierr; 50547c6ae99SBarry Smith 50647c6ae99SBarry Smith PetscFunctionBegin; 50747c6ae99SBarry Smith PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 50847c6ae99SBarry Smith PetscValidHeaderSpecific(da,DM_CLASSID,1); 509c688c046SMatthew G Knepley ierr = PetscTryMethod(mat,"MatSetupDM_C",(Mat,DM),(mat,da));CHKERRQ(ierr); 51047c6ae99SBarry Smith PetscFunctionReturn(0); 51147c6ae99SBarry Smith } 51247c6ae99SBarry Smith 5137087cfbeSBarry Smith PetscErrorCode MatView_MPI_DA(Mat A,PetscViewer viewer) 51447c6ae99SBarry Smith { 5159a42bb27SBarry Smith DM da; 51647c6ae99SBarry Smith PetscErrorCode ierr; 51747c6ae99SBarry Smith const char *prefix; 51847c6ae99SBarry Smith Mat Anatural; 51947c6ae99SBarry Smith AO ao; 52047c6ae99SBarry Smith PetscInt rstart,rend,*petsc,i; 52147c6ae99SBarry Smith IS is; 52247c6ae99SBarry Smith MPI_Comm comm; 52374388724SJed Brown PetscViewerFormat format; 52447c6ae99SBarry Smith 52547c6ae99SBarry Smith PetscFunctionBegin; 52674388724SJed Brown /* Check whether we are just printing info, in which case MatView() already viewed everything we wanted to view */ 52774388724SJed Brown ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 52874388724SJed Brown if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(0); 52974388724SJed Brown 53047c6ae99SBarry Smith ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 531c688c046SMatthew G Knepley ierr = MatGetDM(A, &da);CHKERRQ(ierr); 532ce94432eSBarry Smith if (!da) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix not generated from a DMDA"); 53347c6ae99SBarry Smith 534aa219208SBarry Smith ierr = DMDAGetAO(da,&ao);CHKERRQ(ierr); 53547c6ae99SBarry Smith ierr = MatGetOwnershipRange(A,&rstart,&rend);CHKERRQ(ierr); 536854ce69bSBarry Smith ierr = PetscMalloc1(rend-rstart,&petsc);CHKERRQ(ierr); 53747c6ae99SBarry Smith for (i=rstart; i<rend; i++) petsc[i-rstart] = i; 53847c6ae99SBarry Smith ierr = AOApplicationToPetsc(ao,rend-rstart,petsc);CHKERRQ(ierr); 53947c6ae99SBarry Smith ierr = ISCreateGeneral(comm,rend-rstart,petsc,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 54047c6ae99SBarry Smith 54147c6ae99SBarry Smith /* call viewer on natural ordering */ 542*7dae84e0SHong Zhang ierr = MatCreateSubMatrix(A,is,is,MAT_INITIAL_MATRIX,&Anatural);CHKERRQ(ierr); 543fcfd50ebSBarry Smith ierr = ISDestroy(&is);CHKERRQ(ierr); 54447c6ae99SBarry Smith ierr = PetscObjectGetOptionsPrefix((PetscObject)A,&prefix);CHKERRQ(ierr); 54547c6ae99SBarry Smith ierr = PetscObjectSetOptionsPrefix((PetscObject)Anatural,prefix);CHKERRQ(ierr); 54647c6ae99SBarry Smith ierr = PetscObjectSetName((PetscObject)Anatural,((PetscObject)A)->name);CHKERRQ(ierr); 547539c167fSBarry Smith ierr = MatView(Anatural,viewer);CHKERRQ(ierr); 548fcfd50ebSBarry Smith ierr = MatDestroy(&Anatural);CHKERRQ(ierr); 54947c6ae99SBarry Smith PetscFunctionReturn(0); 55047c6ae99SBarry Smith } 55147c6ae99SBarry Smith 5527087cfbeSBarry Smith PetscErrorCode MatLoad_MPI_DA(Mat A,PetscViewer viewer) 55347c6ae99SBarry Smith { 5549a42bb27SBarry Smith DM da; 55547c6ae99SBarry Smith PetscErrorCode ierr; 55647c6ae99SBarry Smith Mat Anatural,Aapp; 55747c6ae99SBarry Smith AO ao; 558539c167fSBarry Smith PetscInt rstart,rend,*app,i,m,n,M,N; 55947c6ae99SBarry Smith IS is; 56047c6ae99SBarry Smith MPI_Comm comm; 56147c6ae99SBarry Smith 56247c6ae99SBarry Smith PetscFunctionBegin; 56347c6ae99SBarry Smith ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr); 564c688c046SMatthew G Knepley ierr = MatGetDM(A, &da);CHKERRQ(ierr); 565ce94432eSBarry Smith if (!da) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix not generated from a DMDA"); 56647c6ae99SBarry Smith 56747c6ae99SBarry Smith /* Load the matrix in natural ordering */ 568ce94432eSBarry Smith ierr = MatCreate(PetscObjectComm((PetscObject)A),&Anatural);CHKERRQ(ierr); 56947c6ae99SBarry Smith ierr = MatSetType(Anatural,((PetscObject)A)->type_name);CHKERRQ(ierr); 570539c167fSBarry Smith ierr = MatGetSize(A,&M,&N);CHKERRQ(ierr); 571539c167fSBarry Smith ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 572539c167fSBarry Smith ierr = MatSetSizes(Anatural,m,n,M,N);CHKERRQ(ierr); 57347c6ae99SBarry Smith ierr = MatLoad(Anatural,viewer);CHKERRQ(ierr); 57447c6ae99SBarry Smith 57547c6ae99SBarry Smith /* Map natural ordering to application ordering and create IS */ 576aa219208SBarry Smith ierr = DMDAGetAO(da,&ao);CHKERRQ(ierr); 57747c6ae99SBarry Smith ierr = MatGetOwnershipRange(Anatural,&rstart,&rend);CHKERRQ(ierr); 578854ce69bSBarry Smith ierr = PetscMalloc1(rend-rstart,&app);CHKERRQ(ierr); 57947c6ae99SBarry Smith for (i=rstart; i<rend; i++) app[i-rstart] = i; 58047c6ae99SBarry Smith ierr = AOPetscToApplication(ao,rend-rstart,app);CHKERRQ(ierr); 58147c6ae99SBarry Smith ierr = ISCreateGeneral(comm,rend-rstart,app,PETSC_OWN_POINTER,&is);CHKERRQ(ierr); 58247c6ae99SBarry Smith 58347c6ae99SBarry Smith /* Do permutation and replace header */ 584*7dae84e0SHong Zhang ierr = MatCreateSubMatrix(Anatural,is,is,MAT_INITIAL_MATRIX,&Aapp);CHKERRQ(ierr); 58528be2f97SBarry Smith ierr = MatHeaderReplace(A,&Aapp);CHKERRQ(ierr); 586fcfd50ebSBarry Smith ierr = ISDestroy(&is);CHKERRQ(ierr); 587fcfd50ebSBarry Smith ierr = MatDestroy(&Anatural);CHKERRQ(ierr); 58847c6ae99SBarry Smith PetscFunctionReturn(0); 58947c6ae99SBarry Smith } 59047c6ae99SBarry Smith 591b412c318SBarry Smith PetscErrorCode DMCreateMatrix_DA(DM da, Mat *J) 59247c6ae99SBarry Smith { 59347c6ae99SBarry Smith PetscErrorCode ierr; 59447c6ae99SBarry Smith PetscInt dim,dof,nx,ny,nz,dims[3],starts[3],M,N,P; 59547c6ae99SBarry Smith Mat A; 59647c6ae99SBarry Smith MPI_Comm comm; 59719fd82e9SBarry Smith MatType Atype; 59837d0c07bSMatthew G Knepley PetscSection section, sectionGlobal; 5990298fd71SBarry Smith void (*aij)(void)=NULL,(*baij)(void)=NULL,(*sbaij)(void)=NULL; 600b412c318SBarry Smith MatType mtype; 60147c6ae99SBarry Smith PetscMPIInt size; 60247c6ae99SBarry Smith DM_DA *dd = (DM_DA*)da->data; 60347c6ae99SBarry Smith 60447c6ae99SBarry Smith PetscFunctionBegin; 605607a6623SBarry Smith ierr = MatInitializePackage();CHKERRQ(ierr); 606b412c318SBarry Smith mtype = da->mattype; 60747c6ae99SBarry Smith 60837d0c07bSMatthew G Knepley ierr = DMGetDefaultSection(da, §ion);CHKERRQ(ierr); 60937d0c07bSMatthew G Knepley if (section) { 61037d0c07bSMatthew G Knepley PetscInt bs = -1; 61137d0c07bSMatthew G Knepley PetscInt localSize; 61237d0c07bSMatthew G Knepley PetscBool isShell, isBlock, isSeqBlock, isMPIBlock, isSymBlock, isSymSeqBlock, isSymMPIBlock, isSymmetric; 61337d0c07bSMatthew G Knepley 61437d0c07bSMatthew G Knepley ierr = DMGetDefaultGlobalSection(da, §ionGlobal);CHKERRQ(ierr); 61537d0c07bSMatthew G Knepley ierr = PetscSectionGetConstrainedStorageSize(sectionGlobal, &localSize);CHKERRQ(ierr); 61682f516ccSBarry Smith ierr = MatCreate(PetscObjectComm((PetscObject)da), J);CHKERRQ(ierr); 61737d0c07bSMatthew G Knepley ierr = MatSetSizes(*J, localSize, localSize, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 61837d0c07bSMatthew G Knepley ierr = MatSetType(*J, mtype);CHKERRQ(ierr); 61937d0c07bSMatthew G Knepley ierr = MatSetFromOptions(*J);CHKERRQ(ierr); 62037d0c07bSMatthew G Knepley ierr = PetscStrcmp(mtype, MATSHELL, &isShell);CHKERRQ(ierr); 62137d0c07bSMatthew G Knepley ierr = PetscStrcmp(mtype, MATBAIJ, &isBlock);CHKERRQ(ierr); 62237d0c07bSMatthew G Knepley ierr = PetscStrcmp(mtype, MATSEQBAIJ, &isSeqBlock);CHKERRQ(ierr); 62337d0c07bSMatthew G Knepley ierr = PetscStrcmp(mtype, MATMPIBAIJ, &isMPIBlock);CHKERRQ(ierr); 62437d0c07bSMatthew G Knepley ierr = PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);CHKERRQ(ierr); 62537d0c07bSMatthew G Knepley ierr = PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);CHKERRQ(ierr); 62637d0c07bSMatthew G Knepley ierr = PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);CHKERRQ(ierr); 62737d0c07bSMatthew G Knepley /* Check for symmetric storage */ 62837d0c07bSMatthew G Knepley isSymmetric = (PetscBool) (isSymBlock || isSymSeqBlock || isSymMPIBlock); 62937d0c07bSMatthew G Knepley if (isSymmetric) { 63037d0c07bSMatthew G Knepley ierr = MatSetOption(*J, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);CHKERRQ(ierr); 63137d0c07bSMatthew G Knepley } 63237d0c07bSMatthew G Knepley if (!isShell) { 63337d0c07bSMatthew G Knepley PetscInt *dnz, *onz, *dnzu, *onzu, bsLocal; 63437d0c07bSMatthew G Knepley 63537d0c07bSMatthew G Knepley if (bs < 0) { 63637d0c07bSMatthew G Knepley if (isBlock || isSeqBlock || isMPIBlock || isSymBlock || isSymSeqBlock || isSymMPIBlock) { 63737d0c07bSMatthew G Knepley PetscInt pStart, pEnd, p, dof; 63837d0c07bSMatthew G Knepley 63937d0c07bSMatthew G Knepley ierr = PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);CHKERRQ(ierr); 64037d0c07bSMatthew G Knepley for (p = pStart; p < pEnd; ++p) { 64137d0c07bSMatthew G Knepley ierr = PetscSectionGetDof(sectionGlobal, p, &dof);CHKERRQ(ierr); 64237d0c07bSMatthew G Knepley if (dof) { 64337d0c07bSMatthew G Knepley bs = dof; 64437d0c07bSMatthew G Knepley break; 64537d0c07bSMatthew G Knepley } 64637d0c07bSMatthew G Knepley } 64737d0c07bSMatthew G Knepley } else { 64837d0c07bSMatthew G Knepley bs = 1; 64937d0c07bSMatthew G Knepley } 65037d0c07bSMatthew G Knepley /* Must have same blocksize on all procs (some might have no points) */ 65137d0c07bSMatthew G Knepley bsLocal = bs; 652b2566f29SBarry Smith ierr = MPIU_Allreduce(&bsLocal, &bs, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)da));CHKERRQ(ierr); 65337d0c07bSMatthew G Knepley } 6541795a4d1SJed Brown ierr = PetscCalloc4(localSize/bs, &dnz, localSize/bs, &onz, localSize/bs, &dnzu, localSize/bs, &onzu);CHKERRQ(ierr); 655552f7358SJed Brown /* ierr = DMPlexPreallocateOperator(dm, bs, section, sectionGlobal, dnz, onz, dnzu, onzu, *J, fillMatrix);CHKERRQ(ierr); */ 65637d0c07bSMatthew G Knepley ierr = PetscFree4(dnz, onz, dnzu, onzu);CHKERRQ(ierr); 65737d0c07bSMatthew G Knepley } 65837d0c07bSMatthew G Knepley } 65947c6ae99SBarry Smith /* 66047c6ae99SBarry Smith m 66147c6ae99SBarry Smith ------------------------------------------------------ 66247c6ae99SBarry Smith | | 66347c6ae99SBarry Smith | | 66447c6ae99SBarry Smith | ---------------------- | 66547c6ae99SBarry Smith | | | | 66647c6ae99SBarry Smith n | ny | | | 66747c6ae99SBarry Smith | | | | 66847c6ae99SBarry Smith | .--------------------- | 66947c6ae99SBarry Smith | (xs,ys) nx | 67047c6ae99SBarry Smith | . | 67147c6ae99SBarry Smith | (gxs,gys) | 67247c6ae99SBarry Smith | | 67347c6ae99SBarry Smith ----------------------------------------------------- 67447c6ae99SBarry Smith */ 67547c6ae99SBarry Smith 67647c6ae99SBarry Smith /* 67747c6ae99SBarry Smith nc - number of components per grid point 67847c6ae99SBarry Smith col - number of colors needed in one direction for single component problem 67947c6ae99SBarry Smith 68047c6ae99SBarry Smith */ 681e30e807fSPeter Brune M = dd->M; 682e30e807fSPeter Brune N = dd->N; 683e30e807fSPeter Brune P = dd->P; 684c73cfb54SMatthew G. Knepley dim = da->dim; 685e30e807fSPeter Brune dof = dd->w; 686e30e807fSPeter Brune /* ierr = DMDAGetInfo(da,&dim,&M,&N,&P,0,0,0,&dof,0,0,0,0,0);CHKERRQ(ierr); */ 687aa219208SBarry Smith ierr = DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);CHKERRQ(ierr); 68847c6ae99SBarry Smith ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 68947c6ae99SBarry Smith ierr = MatCreate(comm,&A);CHKERRQ(ierr); 69047c6ae99SBarry Smith ierr = MatSetSizes(A,dof*nx*ny*nz,dof*nx*ny*nz,dof*M*N*P,dof*M*N*P);CHKERRQ(ierr); 691b412c318SBarry Smith ierr = MatSetType(A,mtype);CHKERRQ(ierr); 69295ee5b0eSBarry Smith ierr = MatSetDM(A,da);CHKERRQ(ierr); 69347c6ae99SBarry Smith ierr = MatSetFromOptions(A);CHKERRQ(ierr); 69447c6ae99SBarry Smith ierr = MatGetType(A,&Atype);CHKERRQ(ierr); 69547c6ae99SBarry Smith /* 696aa219208SBarry Smith We do not provide a getmatrix function in the DMDA operations because 697aa219208SBarry Smith the basic DMDA does not know about matrices. We think of DMDA as being more 69847c6ae99SBarry Smith more low-level than matrices. This is kind of cheating but, cause sometimes 699aa219208SBarry Smith we think of DMDA has higher level than matrices. 70047c6ae99SBarry Smith 70147c6ae99SBarry Smith We could switch based on Atype (or mtype), but we do not since the 70247c6ae99SBarry Smith specialized setting routines depend only the particular preallocation 70347c6ae99SBarry Smith details of the matrix, not the type itself. 70447c6ae99SBarry Smith */ 70547c6ae99SBarry Smith ierr = PetscObjectQueryFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",&aij);CHKERRQ(ierr); 70647c6ae99SBarry Smith if (!aij) { 70747c6ae99SBarry Smith ierr = PetscObjectQueryFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",&aij);CHKERRQ(ierr); 70847c6ae99SBarry Smith } 70947c6ae99SBarry Smith if (!aij) { 71047c6ae99SBarry Smith ierr = PetscObjectQueryFunction((PetscObject)A,"MatMPIBAIJSetPreallocation_C",&baij);CHKERRQ(ierr); 71147c6ae99SBarry Smith if (!baij) { 71247c6ae99SBarry Smith ierr = PetscObjectQueryFunction((PetscObject)A,"MatSeqBAIJSetPreallocation_C",&baij);CHKERRQ(ierr); 71347c6ae99SBarry Smith } 71447c6ae99SBarry Smith if (!baij) { 71547c6ae99SBarry Smith ierr = PetscObjectQueryFunction((PetscObject)A,"MatMPISBAIJSetPreallocation_C",&sbaij);CHKERRQ(ierr); 71647c6ae99SBarry Smith if (!sbaij) { 71747c6ae99SBarry Smith ierr = PetscObjectQueryFunction((PetscObject)A,"MatSeqSBAIJSetPreallocation_C",&sbaij);CHKERRQ(ierr); 71847c6ae99SBarry Smith } 71947c6ae99SBarry Smith } 72047c6ae99SBarry Smith } 72147c6ae99SBarry Smith if (aij) { 72247c6ae99SBarry Smith if (dim == 1) { 723ce308e1dSBarry Smith if (dd->ofill) { 724ce308e1dSBarry Smith ierr = DMCreateMatrix_DA_1d_MPIAIJ_Fill(da,A);CHKERRQ(ierr); 725ce308e1dSBarry Smith } else { 726950540a4SJed Brown ierr = DMCreateMatrix_DA_1d_MPIAIJ(da,A);CHKERRQ(ierr); 727ce308e1dSBarry Smith } 72847c6ae99SBarry Smith } else if (dim == 2) { 72947c6ae99SBarry Smith if (dd->ofill) { 730950540a4SJed Brown ierr = DMCreateMatrix_DA_2d_MPIAIJ_Fill(da,A);CHKERRQ(ierr); 73147c6ae99SBarry Smith } else { 732950540a4SJed Brown ierr = DMCreateMatrix_DA_2d_MPIAIJ(da,A);CHKERRQ(ierr); 73347c6ae99SBarry Smith } 73447c6ae99SBarry Smith } else if (dim == 3) { 73547c6ae99SBarry Smith if (dd->ofill) { 736950540a4SJed Brown ierr = DMCreateMatrix_DA_3d_MPIAIJ_Fill(da,A);CHKERRQ(ierr); 73747c6ae99SBarry Smith } else { 738950540a4SJed Brown ierr = DMCreateMatrix_DA_3d_MPIAIJ(da,A);CHKERRQ(ierr); 73947c6ae99SBarry Smith } 74047c6ae99SBarry Smith } 74147c6ae99SBarry Smith } else if (baij) { 74247c6ae99SBarry Smith if (dim == 2) { 743950540a4SJed Brown ierr = DMCreateMatrix_DA_2d_MPIBAIJ(da,A);CHKERRQ(ierr); 74447c6ae99SBarry Smith } else if (dim == 3) { 745950540a4SJed Brown ierr = DMCreateMatrix_DA_3d_MPIBAIJ(da,A);CHKERRQ(ierr); 746ce94432eSBarry Smith } else SETERRQ3(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Not implemented for %D dimension and Matrix Type: %s in %D dimension! Send mail to petsc-maint@mcs.anl.gov for code",dim,Atype,dim); 74747c6ae99SBarry Smith } else if (sbaij) { 74847c6ae99SBarry Smith if (dim == 2) { 749950540a4SJed Brown ierr = DMCreateMatrix_DA_2d_MPISBAIJ(da,A);CHKERRQ(ierr); 75047c6ae99SBarry Smith } else if (dim == 3) { 751950540a4SJed Brown ierr = DMCreateMatrix_DA_3d_MPISBAIJ(da,A);CHKERRQ(ierr); 752ce94432eSBarry Smith } else SETERRQ3(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Not implemented for %D dimension and Matrix Type: %s in %D dimension! Send mail to petsc-maint@mcs.anl.gov for code",dim,Atype,dim); 753869776cdSLisandro Dalcin } else { 75445b6f7e9SBarry Smith ISLocalToGlobalMapping ltog; 755869776cdSLisandro Dalcin ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 7562949035bSJed Brown ierr = MatSetUp(A);CHKERRQ(ierr); 757869776cdSLisandro Dalcin ierr = MatSetLocalToGlobalMapping(A,ltog,ltog);CHKERRQ(ierr); 75847c6ae99SBarry Smith } 759aa219208SBarry Smith ierr = DMDAGetGhostCorners(da,&starts[0],&starts[1],&starts[2],&dims[0],&dims[1],&dims[2]);CHKERRQ(ierr); 76047c6ae99SBarry Smith ierr = MatSetStencil(A,dim,dims,starts,dof);CHKERRQ(ierr); 761c688c046SMatthew G Knepley ierr = MatSetDM(A,da);CHKERRQ(ierr); 76247c6ae99SBarry Smith ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 76347c6ae99SBarry Smith if (size > 1) { 76447c6ae99SBarry Smith /* change viewer to display matrix in natural ordering */ 76547c6ae99SBarry Smith ierr = MatShellSetOperation(A, MATOP_VIEW, (void (*)(void))MatView_MPI_DA);CHKERRQ(ierr); 76647c6ae99SBarry Smith ierr = MatShellSetOperation(A, MATOP_LOAD, (void (*)(void))MatLoad_MPI_DA);CHKERRQ(ierr); 76747c6ae99SBarry Smith } 76847c6ae99SBarry Smith *J = A; 76947c6ae99SBarry Smith PetscFunctionReturn(0); 77047c6ae99SBarry Smith } 77147c6ae99SBarry Smith 77247c6ae99SBarry Smith /* ---------------------------------------------------------------------------------*/ 773950540a4SJed Brown PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ(DM da,Mat J) 77447c6ae99SBarry Smith { 77547c6ae99SBarry Smith PetscErrorCode ierr; 776c1154cd5SBarry Smith PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny,m,n,dim,s,*cols = NULL,k,nc,*rows = NULL,col,cnt,l,p,M,N; 77747c6ae99SBarry Smith PetscInt lstart,lend,pstart,pend,*dnz,*onz; 77847c6ae99SBarry Smith MPI_Comm comm; 77947c6ae99SBarry Smith PetscScalar *values; 780bff4a2f0SMatthew G. Knepley DMBoundaryType bx,by; 78145b6f7e9SBarry Smith ISLocalToGlobalMapping ltog; 782aa219208SBarry Smith DMDAStencilType st; 783c1154cd5SBarry Smith PetscBool removedups = PETSC_FALSE; 78447c6ae99SBarry Smith 78547c6ae99SBarry Smith PetscFunctionBegin; 78647c6ae99SBarry Smith /* 78747c6ae99SBarry Smith nc - number of components per grid point 78847c6ae99SBarry Smith col - number of colors needed in one direction for single component problem 78947c6ae99SBarry Smith 79047c6ae99SBarry Smith */ 791c1154cd5SBarry Smith ierr = DMDAGetInfo(da,&dim,&m,&n,&M,&N,0,0,&nc,&s,&bx,&by,0,&st);CHKERRQ(ierr); 79247c6ae99SBarry Smith col = 2*s + 1; 793c1154cd5SBarry Smith /* 794c1154cd5SBarry Smith With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times 795c1154cd5SBarry Smith because of "wrapping" around the end of the domain hitting an entry already counted in the other direction. 796c1154cd5SBarry Smith */ 797c1154cd5SBarry Smith if (M == 1 && 2*s >= m) removedups = PETSC_TRUE; 798c1154cd5SBarry Smith if (N == 1 && 2*s >= n) removedups = PETSC_TRUE; 799aa219208SBarry Smith ierr = DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);CHKERRQ(ierr); 800aa219208SBarry Smith ierr = DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);CHKERRQ(ierr); 80147c6ae99SBarry Smith ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 80247c6ae99SBarry Smith 803dcca6d9dSJed Brown ierr = PetscMalloc2(nc,&rows,col*col*nc*nc,&cols);CHKERRQ(ierr); 8041411c6eeSJed Brown ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 80547c6ae99SBarry Smith 80606ca8cadSBarry Smith ierr = MatSetBlockSize(J,nc);CHKERRQ(ierr); 80747c6ae99SBarry Smith /* determine the matrix preallocation information */ 80847c6ae99SBarry Smith ierr = MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);CHKERRQ(ierr); 80947c6ae99SBarry Smith for (i=xs; i<xs+nx; i++) { 81047c6ae99SBarry Smith 811bff4a2f0SMatthew G. Knepley pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 812bff4a2f0SMatthew G. Knepley pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 81347c6ae99SBarry Smith 81447c6ae99SBarry Smith for (j=ys; j<ys+ny; j++) { 81547c6ae99SBarry Smith slot = i - gxs + gnx*(j - gys); 81647c6ae99SBarry Smith 817bff4a2f0SMatthew G. Knepley lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 818bff4a2f0SMatthew G. Knepley lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 81947c6ae99SBarry Smith 82047c6ae99SBarry Smith cnt = 0; 82147c6ae99SBarry Smith for (k=0; k<nc; k++) { 82247c6ae99SBarry Smith for (l=lstart; l<lend+1; l++) { 82347c6ae99SBarry Smith for (p=pstart; p<pend+1; p++) { 824aa219208SBarry Smith if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */ 82547c6ae99SBarry Smith cols[cnt++] = k + nc*(slot + gnx*l + p); 82647c6ae99SBarry Smith } 82747c6ae99SBarry Smith } 82847c6ae99SBarry Smith } 82947c6ae99SBarry Smith rows[k] = k + nc*(slot); 83047c6ae99SBarry Smith } 831c1154cd5SBarry Smith if (removedups) { 832c1154cd5SBarry Smith ierr = MatPreallocateSetLocalRemoveDups(ltog,nc,rows,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 833c1154cd5SBarry Smith } else { 834784ac674SJed Brown ierr = MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 83547c6ae99SBarry Smith } 83647c6ae99SBarry Smith } 837c1154cd5SBarry Smith } 838f73d5cc4SBarry Smith ierr = MatSetBlockSize(J,nc);CHKERRQ(ierr); 83947c6ae99SBarry Smith ierr = MatSeqAIJSetPreallocation(J,0,dnz);CHKERRQ(ierr); 84047c6ae99SBarry Smith ierr = MatMPIAIJSetPreallocation(J,0,dnz,0,onz);CHKERRQ(ierr); 84147c6ae99SBarry Smith ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 84247c6ae99SBarry Smith 843784ac674SJed Brown ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 84447c6ae99SBarry Smith 84547c6ae99SBarry Smith /* 84647c6ae99SBarry Smith For each node in the grid: we get the neighbors in the local (on processor ordering 84747c6ae99SBarry Smith that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 84847c6ae99SBarry Smith PETSc ordering. 84947c6ae99SBarry Smith */ 850fcfd50ebSBarry Smith if (!da->prealloc_only) { 8511795a4d1SJed Brown ierr = PetscCalloc1(col*col*nc*nc,&values);CHKERRQ(ierr); 85247c6ae99SBarry Smith for (i=xs; i<xs+nx; i++) { 85347c6ae99SBarry Smith 854bff4a2f0SMatthew G. Knepley pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 855bff4a2f0SMatthew G. Knepley pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 85647c6ae99SBarry Smith 85747c6ae99SBarry Smith for (j=ys; j<ys+ny; j++) { 85847c6ae99SBarry Smith slot = i - gxs + gnx*(j - gys); 85947c6ae99SBarry Smith 860bff4a2f0SMatthew G. Knepley lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 861bff4a2f0SMatthew G. Knepley lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 86247c6ae99SBarry Smith 86347c6ae99SBarry Smith cnt = 0; 86447c6ae99SBarry Smith for (k=0; k<nc; k++) { 86547c6ae99SBarry Smith for (l=lstart; l<lend+1; l++) { 86647c6ae99SBarry Smith for (p=pstart; p<pend+1; p++) { 867aa219208SBarry Smith if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */ 86847c6ae99SBarry Smith cols[cnt++] = k + nc*(slot + gnx*l + p); 86947c6ae99SBarry Smith } 87047c6ae99SBarry Smith } 87147c6ae99SBarry Smith } 87247c6ae99SBarry Smith rows[k] = k + nc*(slot); 87347c6ae99SBarry Smith } 87447c6ae99SBarry Smith ierr = MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 87547c6ae99SBarry Smith } 87647c6ae99SBarry Smith } 87747c6ae99SBarry Smith ierr = PetscFree(values);CHKERRQ(ierr); 87847c6ae99SBarry Smith ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 87947c6ae99SBarry Smith ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 880189e4007SBarry Smith ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 88147c6ae99SBarry Smith } 88247c6ae99SBarry Smith ierr = PetscFree2(rows,cols);CHKERRQ(ierr); 88347c6ae99SBarry Smith PetscFunctionReturn(0); 88447c6ae99SBarry Smith } 88547c6ae99SBarry Smith 886950540a4SJed Brown PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ_Fill(DM da,Mat J) 88747c6ae99SBarry Smith { 88847c6ae99SBarry Smith PetscErrorCode ierr; 88947c6ae99SBarry Smith PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 890c1154cd5SBarry Smith PetscInt m,n,dim,s,*cols,k,nc,row,col,cnt,maxcnt = 0,l,p,M,N; 89147c6ae99SBarry Smith PetscInt lstart,lend,pstart,pend,*dnz,*onz; 89247c6ae99SBarry Smith DM_DA *dd = (DM_DA*)da->data; 89347c6ae99SBarry Smith PetscInt ifill_col,*ofill = dd->ofill, *dfill = dd->dfill; 89447c6ae99SBarry Smith MPI_Comm comm; 89547c6ae99SBarry Smith PetscScalar *values; 896bff4a2f0SMatthew G. Knepley DMBoundaryType bx,by; 89745b6f7e9SBarry Smith ISLocalToGlobalMapping ltog; 898aa219208SBarry Smith DMDAStencilType st; 899c1154cd5SBarry Smith PetscBool removedups = PETSC_FALSE; 90047c6ae99SBarry Smith 90147c6ae99SBarry Smith PetscFunctionBegin; 90247c6ae99SBarry Smith /* 90347c6ae99SBarry Smith nc - number of components per grid point 90447c6ae99SBarry Smith col - number of colors needed in one direction for single component problem 90547c6ae99SBarry Smith 90647c6ae99SBarry Smith */ 907c1154cd5SBarry Smith ierr = DMDAGetInfo(da,&dim,&m,&n,&M,&N,0,0,&nc,&s,&bx,&by,0,&st);CHKERRQ(ierr); 90847c6ae99SBarry Smith col = 2*s + 1; 909c1154cd5SBarry Smith /* 910c1154cd5SBarry Smith With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times 911c1154cd5SBarry Smith because of "wrapping" around the end of the domain hitting an entry already counted in the other direction. 912c1154cd5SBarry Smith */ 913c1154cd5SBarry Smith if (M == 1 && 2*s >= m) removedups = PETSC_TRUE; 914c1154cd5SBarry Smith if (N == 1 && 2*s >= n) removedups = PETSC_TRUE; 915aa219208SBarry Smith ierr = DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);CHKERRQ(ierr); 916aa219208SBarry Smith ierr = DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);CHKERRQ(ierr); 91747c6ae99SBarry Smith ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 91847c6ae99SBarry Smith 9194b26d1cfSBarry Smith ierr = PetscMalloc1(col*col*nc,&cols);CHKERRQ(ierr); 9201411c6eeSJed Brown ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 92147c6ae99SBarry Smith 92206ca8cadSBarry Smith ierr = MatSetBlockSize(J,nc);CHKERRQ(ierr); 92347c6ae99SBarry Smith /* determine the matrix preallocation information */ 92447c6ae99SBarry Smith ierr = MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);CHKERRQ(ierr); 92547c6ae99SBarry Smith for (i=xs; i<xs+nx; i++) { 92647c6ae99SBarry Smith 927bff4a2f0SMatthew G. Knepley pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 928bff4a2f0SMatthew G. Knepley pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 92947c6ae99SBarry Smith 93047c6ae99SBarry Smith for (j=ys; j<ys+ny; j++) { 93147c6ae99SBarry Smith slot = i - gxs + gnx*(j - gys); 93247c6ae99SBarry Smith 933bff4a2f0SMatthew G. Knepley lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 934bff4a2f0SMatthew G. Knepley lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 93547c6ae99SBarry Smith 93647c6ae99SBarry Smith for (k=0; k<nc; k++) { 93747c6ae99SBarry Smith cnt = 0; 93847c6ae99SBarry Smith for (l=lstart; l<lend+1; l++) { 93947c6ae99SBarry Smith for (p=pstart; p<pend+1; p++) { 94047c6ae99SBarry Smith if (l || p) { 941aa219208SBarry Smith if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star */ 9428865f1eaSKarl Rupp for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++) cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p); 94347c6ae99SBarry Smith } 94447c6ae99SBarry Smith } else { 94547c6ae99SBarry Smith if (dfill) { 9468865f1eaSKarl Rupp for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++) cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p); 94747c6ae99SBarry Smith } else { 9488865f1eaSKarl Rupp for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + gnx*l + p); 94947c6ae99SBarry Smith } 95047c6ae99SBarry Smith } 95147c6ae99SBarry Smith } 95247c6ae99SBarry Smith } 95347c6ae99SBarry Smith row = k + nc*(slot); 954c0ab637bSBarry Smith maxcnt = PetscMax(maxcnt,cnt); 955c1154cd5SBarry Smith if (removedups) { 956c1154cd5SBarry Smith ierr = MatPreallocateSetLocalRemoveDups(ltog,1,&row,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 957c1154cd5SBarry Smith } else { 958784ac674SJed Brown ierr = MatPreallocateSetLocal(ltog,1,&row,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 95947c6ae99SBarry Smith } 96047c6ae99SBarry Smith } 96147c6ae99SBarry Smith } 962c1154cd5SBarry Smith } 96347c6ae99SBarry Smith ierr = MatSeqAIJSetPreallocation(J,0,dnz);CHKERRQ(ierr); 96447c6ae99SBarry Smith ierr = MatMPIAIJSetPreallocation(J,0,dnz,0,onz);CHKERRQ(ierr); 96547c6ae99SBarry Smith ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 966784ac674SJed Brown ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 96747c6ae99SBarry Smith 96847c6ae99SBarry Smith /* 96947c6ae99SBarry Smith For each node in the grid: we get the neighbors in the local (on processor ordering 97047c6ae99SBarry Smith that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 97147c6ae99SBarry Smith PETSc ordering. 97247c6ae99SBarry Smith */ 973fcfd50ebSBarry Smith if (!da->prealloc_only) { 974c0ab637bSBarry Smith ierr = PetscCalloc1(maxcnt,&values);CHKERRQ(ierr); 97547c6ae99SBarry Smith for (i=xs; i<xs+nx; i++) { 97647c6ae99SBarry Smith 977bff4a2f0SMatthew G. Knepley pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 978bff4a2f0SMatthew G. Knepley pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 97947c6ae99SBarry Smith 98047c6ae99SBarry Smith for (j=ys; j<ys+ny; j++) { 98147c6ae99SBarry Smith slot = i - gxs + gnx*(j - gys); 98247c6ae99SBarry Smith 983bff4a2f0SMatthew G. Knepley lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 984bff4a2f0SMatthew G. Knepley lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 98547c6ae99SBarry Smith 98647c6ae99SBarry Smith for (k=0; k<nc; k++) { 98747c6ae99SBarry Smith cnt = 0; 98847c6ae99SBarry Smith for (l=lstart; l<lend+1; l++) { 98947c6ae99SBarry Smith for (p=pstart; p<pend+1; p++) { 99047c6ae99SBarry Smith if (l || p) { 991aa219208SBarry Smith if ((st == DMDA_STENCIL_BOX) || (!l || !p)) { /* entries on star */ 9928865f1eaSKarl Rupp for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++) cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p); 99347c6ae99SBarry Smith } 99447c6ae99SBarry Smith } else { 99547c6ae99SBarry Smith if (dfill) { 9968865f1eaSKarl Rupp for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++) cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p); 99747c6ae99SBarry Smith } else { 9988865f1eaSKarl Rupp for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + gnx*l + p); 99947c6ae99SBarry Smith } 100047c6ae99SBarry Smith } 100147c6ae99SBarry Smith } 100247c6ae99SBarry Smith } 100347c6ae99SBarry Smith row = k + nc*(slot); 100447c6ae99SBarry Smith ierr = MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 100547c6ae99SBarry Smith } 100647c6ae99SBarry Smith } 100747c6ae99SBarry Smith } 100847c6ae99SBarry Smith ierr = PetscFree(values);CHKERRQ(ierr); 100947c6ae99SBarry Smith ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 101047c6ae99SBarry Smith ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1011189e4007SBarry Smith ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 101247c6ae99SBarry Smith } 101347c6ae99SBarry Smith ierr = PetscFree(cols);CHKERRQ(ierr); 101447c6ae99SBarry Smith PetscFunctionReturn(0); 101547c6ae99SBarry Smith } 101647c6ae99SBarry Smith 101747c6ae99SBarry Smith /* ---------------------------------------------------------------------------------*/ 101847c6ae99SBarry Smith 1019950540a4SJed Brown PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ(DM da,Mat J) 102047c6ae99SBarry Smith { 102147c6ae99SBarry Smith PetscErrorCode ierr; 102247c6ae99SBarry Smith PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 10230298fd71SBarry Smith PetscInt m,n,dim,s,*cols = NULL,k,nc,*rows = NULL,col,cnt,l,p,*dnz = NULL,*onz = NULL; 1024c1154cd5SBarry Smith PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk,M,N,P; 102547c6ae99SBarry Smith MPI_Comm comm; 102647c6ae99SBarry Smith PetscScalar *values; 1027bff4a2f0SMatthew G. Knepley DMBoundaryType bx,by,bz; 102845b6f7e9SBarry Smith ISLocalToGlobalMapping ltog; 1029aa219208SBarry Smith DMDAStencilType st; 1030c1154cd5SBarry Smith PetscBool removedups = PETSC_FALSE; 103147c6ae99SBarry Smith 103247c6ae99SBarry Smith PetscFunctionBegin; 103347c6ae99SBarry Smith /* 103447c6ae99SBarry Smith nc - number of components per grid point 103547c6ae99SBarry Smith col - number of colors needed in one direction for single component problem 103647c6ae99SBarry Smith 103747c6ae99SBarry Smith */ 1038c1154cd5SBarry Smith ierr = DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&bx,&by,&bz,&st);CHKERRQ(ierr); 103947c6ae99SBarry Smith col = 2*s + 1; 104047c6ae99SBarry Smith 1041c1154cd5SBarry Smith /* 1042c1154cd5SBarry Smith With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times 1043c1154cd5SBarry Smith because of "wrapping" around the end of the domain hitting an entry already counted in the other direction. 1044c1154cd5SBarry Smith */ 1045c1154cd5SBarry Smith if (M == 1 && 2*s >= m) removedups = PETSC_TRUE; 1046c1154cd5SBarry Smith if (N == 1 && 2*s >= n) removedups = PETSC_TRUE; 1047c1154cd5SBarry Smith if (P == 1 && 2*s >= p) removedups = PETSC_TRUE; 1048c1154cd5SBarry Smith 1049aa219208SBarry Smith ierr = DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);CHKERRQ(ierr); 1050aa219208SBarry Smith ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);CHKERRQ(ierr); 105147c6ae99SBarry Smith ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 105247c6ae99SBarry Smith 1053dcca6d9dSJed Brown ierr = PetscMalloc2(nc,&rows,col*col*col*nc*nc,&cols);CHKERRQ(ierr); 10541411c6eeSJed Brown ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 105547c6ae99SBarry Smith 105606ca8cadSBarry Smith ierr = MatSetBlockSize(J,nc);CHKERRQ(ierr); 105747c6ae99SBarry Smith /* determine the matrix preallocation information */ 105847c6ae99SBarry Smith ierr = MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);CHKERRQ(ierr); 105947c6ae99SBarry Smith for (i=xs; i<xs+nx; i++) { 1060bff4a2f0SMatthew G. Knepley istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1061bff4a2f0SMatthew G. Knepley iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 106247c6ae99SBarry Smith for (j=ys; j<ys+ny; j++) { 1063bff4a2f0SMatthew G. Knepley jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1064bff4a2f0SMatthew G. Knepley jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 106547c6ae99SBarry Smith for (k=zs; k<zs+nz; k++) { 1066bff4a2f0SMatthew G. Knepley kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1067bff4a2f0SMatthew G. Knepley kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 106847c6ae99SBarry Smith 106947c6ae99SBarry Smith slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 107047c6ae99SBarry Smith 107147c6ae99SBarry Smith cnt = 0; 107247c6ae99SBarry Smith for (l=0; l<nc; l++) { 107347c6ae99SBarry Smith for (ii=istart; ii<iend+1; ii++) { 107447c6ae99SBarry Smith for (jj=jstart; jj<jend+1; jj++) { 107547c6ae99SBarry Smith for (kk=kstart; kk<kend+1; kk++) { 1076aa219208SBarry Smith if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 107747c6ae99SBarry Smith cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk); 107847c6ae99SBarry Smith } 107947c6ae99SBarry Smith } 108047c6ae99SBarry Smith } 108147c6ae99SBarry Smith } 108247c6ae99SBarry Smith rows[l] = l + nc*(slot); 108347c6ae99SBarry Smith } 1084c1154cd5SBarry Smith if (removedups) { 1085c1154cd5SBarry Smith ierr = MatPreallocateSetLocalRemoveDups(ltog,nc,rows,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 1086c1154cd5SBarry Smith } else { 1087784ac674SJed Brown ierr = MatPreallocateSetLocal(ltog,nc,rows,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 108847c6ae99SBarry Smith } 108947c6ae99SBarry Smith } 109047c6ae99SBarry Smith } 1091c1154cd5SBarry Smith } 1092f73d5cc4SBarry Smith ierr = MatSetBlockSize(J,nc);CHKERRQ(ierr); 109347c6ae99SBarry Smith ierr = MatSeqAIJSetPreallocation(J,0,dnz);CHKERRQ(ierr); 109447c6ae99SBarry Smith ierr = MatMPIAIJSetPreallocation(J,0,dnz,0,onz);CHKERRQ(ierr); 109547c6ae99SBarry Smith ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1096784ac674SJed Brown ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 109747c6ae99SBarry Smith 109847c6ae99SBarry Smith /* 109947c6ae99SBarry Smith For each node in the grid: we get the neighbors in the local (on processor ordering 110047c6ae99SBarry Smith that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 110147c6ae99SBarry Smith PETSc ordering. 110247c6ae99SBarry Smith */ 1103fcfd50ebSBarry Smith if (!da->prealloc_only) { 11041795a4d1SJed Brown ierr = PetscCalloc1(col*col*col*nc*nc*nc,&values);CHKERRQ(ierr); 110547c6ae99SBarry Smith for (i=xs; i<xs+nx; i++) { 1106bff4a2f0SMatthew G. Knepley istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1107bff4a2f0SMatthew G. Knepley iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 110847c6ae99SBarry Smith for (j=ys; j<ys+ny; j++) { 1109bff4a2f0SMatthew G. Knepley jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1110bff4a2f0SMatthew G. Knepley jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 111147c6ae99SBarry Smith for (k=zs; k<zs+nz; k++) { 1112bff4a2f0SMatthew G. Knepley kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1113bff4a2f0SMatthew G. Knepley kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 111447c6ae99SBarry Smith 111547c6ae99SBarry Smith slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 111647c6ae99SBarry Smith 111747c6ae99SBarry Smith cnt = 0; 111847c6ae99SBarry Smith for (l=0; l<nc; l++) { 111947c6ae99SBarry Smith for (ii=istart; ii<iend+1; ii++) { 112047c6ae99SBarry Smith for (jj=jstart; jj<jend+1; jj++) { 112147c6ae99SBarry Smith for (kk=kstart; kk<kend+1; kk++) { 1122aa219208SBarry Smith if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 112347c6ae99SBarry Smith cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk); 112447c6ae99SBarry Smith } 112547c6ae99SBarry Smith } 112647c6ae99SBarry Smith } 112747c6ae99SBarry Smith } 112847c6ae99SBarry Smith rows[l] = l + nc*(slot); 112947c6ae99SBarry Smith } 113047c6ae99SBarry Smith ierr = MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 113147c6ae99SBarry Smith } 113247c6ae99SBarry Smith } 113347c6ae99SBarry Smith } 113447c6ae99SBarry Smith ierr = PetscFree(values);CHKERRQ(ierr); 113547c6ae99SBarry Smith ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 113647c6ae99SBarry Smith ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1137189e4007SBarry Smith ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 113847c6ae99SBarry Smith } 113947c6ae99SBarry Smith ierr = PetscFree2(rows,cols);CHKERRQ(ierr); 114047c6ae99SBarry Smith PetscFunctionReturn(0); 114147c6ae99SBarry Smith } 114247c6ae99SBarry Smith 114347c6ae99SBarry Smith /* ---------------------------------------------------------------------------------*/ 114447c6ae99SBarry Smith 1145ce308e1dSBarry Smith PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ_Fill(DM da,Mat J) 1146ce308e1dSBarry Smith { 1147ce308e1dSBarry Smith PetscErrorCode ierr; 1148ce308e1dSBarry Smith DM_DA *dd = (DM_DA*)da->data; 1149ce308e1dSBarry Smith PetscInt xs,nx,i,j,gxs,gnx,row,k,l; 11508d4c968fSBarry Smith PetscInt m,dim,s,*cols = NULL,nc,cnt,maxcnt = 0,*ocols; 11510acb5bebSBarry Smith PetscInt *ofill = dd->ofill,*dfill = dd->dfill; 1152ce308e1dSBarry Smith PetscScalar *values; 1153bff4a2f0SMatthew G. Knepley DMBoundaryType bx; 115445b6f7e9SBarry Smith ISLocalToGlobalMapping ltog; 1155ce308e1dSBarry Smith PetscMPIInt rank,size; 1156ce308e1dSBarry Smith 1157ce308e1dSBarry Smith PetscFunctionBegin; 1158bff4a2f0SMatthew G. Knepley if (dd->bx == DM_BOUNDARY_PERIODIC) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"With fill provided not implemented with periodic boundary conditions"); 1159ce94432eSBarry Smith ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)da),&rank);CHKERRQ(ierr); 1160ce94432eSBarry Smith ierr = MPI_Comm_size(PetscObjectComm((PetscObject)da),&size);CHKERRQ(ierr); 1161ce308e1dSBarry Smith 1162ce308e1dSBarry Smith /* 1163ce308e1dSBarry Smith nc - number of components per grid point 1164ce308e1dSBarry Smith 1165ce308e1dSBarry Smith */ 1166ce308e1dSBarry Smith ierr = DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);CHKERRQ(ierr); 1167ce308e1dSBarry Smith ierr = DMDAGetCorners(da,&xs,0,0,&nx,0,0);CHKERRQ(ierr); 1168ce308e1dSBarry Smith ierr = DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);CHKERRQ(ierr); 1169ce308e1dSBarry Smith 1170ce308e1dSBarry Smith ierr = MatSetBlockSize(J,nc);CHKERRQ(ierr); 11711795a4d1SJed Brown ierr = PetscCalloc2(nx*nc,&cols,nx*nc,&ocols);CHKERRQ(ierr); 1172ce308e1dSBarry Smith 1173ce308e1dSBarry Smith /* 1174ce308e1dSBarry Smith note should be smaller for first and last process with no periodic 1175ce308e1dSBarry Smith does not handle dfill 1176ce308e1dSBarry Smith */ 1177ce308e1dSBarry Smith cnt = 0; 1178ce308e1dSBarry Smith /* coupling with process to the left */ 1179ce308e1dSBarry Smith for (i=0; i<s; i++) { 1180ce308e1dSBarry Smith for (j=0; j<nc; j++) { 1181ce308e1dSBarry Smith ocols[cnt] = ((!rank) ? 0 : (s - i)*(ofill[j+1] - ofill[j])); 11820acb5bebSBarry Smith cols[cnt] = dfill[j+1] - dfill[j] + (s + i)*(ofill[j+1] - ofill[j]); 1183c0ab637bSBarry Smith maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]); 1184ce308e1dSBarry Smith cnt++; 1185ce308e1dSBarry Smith } 1186ce308e1dSBarry Smith } 1187ce308e1dSBarry Smith for (i=s; i<nx-s; i++) { 1188ce308e1dSBarry Smith for (j=0; j<nc; j++) { 11890acb5bebSBarry Smith cols[cnt] = dfill[j+1] - dfill[j] + 2*s*(ofill[j+1] - ofill[j]); 1190c0ab637bSBarry Smith maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]); 1191ce308e1dSBarry Smith cnt++; 1192ce308e1dSBarry Smith } 1193ce308e1dSBarry Smith } 1194ce308e1dSBarry Smith /* coupling with process to the right */ 1195ce308e1dSBarry Smith for (i=nx-s; i<nx; i++) { 1196ce308e1dSBarry Smith for (j=0; j<nc; j++) { 1197ce308e1dSBarry Smith ocols[cnt] = ((rank == (size-1)) ? 0 : (i - nx + s + 1)*(ofill[j+1] - ofill[j])); 11980acb5bebSBarry Smith cols[cnt] = dfill[j+1] - dfill[j] + (s + nx - i - 1)*(ofill[j+1] - ofill[j]); 1199c0ab637bSBarry Smith maxcnt = PetscMax(maxcnt,ocols[cnt]+cols[cnt]); 1200ce308e1dSBarry Smith cnt++; 1201ce308e1dSBarry Smith } 1202ce308e1dSBarry Smith } 1203ce308e1dSBarry Smith 1204ce308e1dSBarry Smith ierr = MatSeqAIJSetPreallocation(J,0,cols);CHKERRQ(ierr); 1205ce308e1dSBarry Smith ierr = MatMPIAIJSetPreallocation(J,0,cols,0,ocols);CHKERRQ(ierr); 1206ce308e1dSBarry Smith ierr = PetscFree2(cols,ocols);CHKERRQ(ierr); 1207ce308e1dSBarry Smith 1208ce308e1dSBarry Smith ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1209ce308e1dSBarry Smith ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 1210ce308e1dSBarry Smith 1211ce308e1dSBarry Smith /* 1212ce308e1dSBarry Smith For each node in the grid: we get the neighbors in the local (on processor ordering 1213ce308e1dSBarry Smith that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 1214ce308e1dSBarry Smith PETSc ordering. 1215ce308e1dSBarry Smith */ 1216ce308e1dSBarry Smith if (!da->prealloc_only) { 1217c0ab637bSBarry Smith ierr = PetscCalloc2(maxcnt,&values,maxcnt,&cols);CHKERRQ(ierr); 1218ce308e1dSBarry Smith 1219ce308e1dSBarry Smith row = xs*nc; 1220ce308e1dSBarry Smith /* coupling with process to the left */ 1221ce308e1dSBarry Smith for (i=xs; i<xs+s; i++) { 1222ce308e1dSBarry Smith for (j=0; j<nc; j++) { 1223ce308e1dSBarry Smith cnt = 0; 1224ce308e1dSBarry Smith if (rank) { 1225ce308e1dSBarry Smith for (l=0; l<s; l++) { 1226ce308e1dSBarry Smith for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k]; 1227ce308e1dSBarry Smith } 1228ce308e1dSBarry Smith } 12290acb5bebSBarry Smith if (dfill) { 12300acb5bebSBarry Smith for (k=dfill[j]; k<dfill[j+1]; k++) { 12310acb5bebSBarry Smith cols[cnt++] = i*nc + dfill[k]; 12320acb5bebSBarry Smith } 12330acb5bebSBarry Smith } else { 1234ce308e1dSBarry Smith for (k=0; k<nc; k++) { 1235ce308e1dSBarry Smith cols[cnt++] = i*nc + k; 1236ce308e1dSBarry Smith } 12370acb5bebSBarry Smith } 1238ce308e1dSBarry Smith for (l=0; l<s; l++) { 1239ce308e1dSBarry Smith for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k]; 1240ce308e1dSBarry Smith } 1241ce308e1dSBarry Smith ierr = MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1242ce308e1dSBarry Smith row++; 1243ce308e1dSBarry Smith } 1244ce308e1dSBarry Smith } 1245ce308e1dSBarry Smith for (i=xs+s; i<xs+nx-s; i++) { 1246ce308e1dSBarry Smith for (j=0; j<nc; j++) { 1247ce308e1dSBarry Smith cnt = 0; 1248ce308e1dSBarry Smith for (l=0; l<s; l++) { 1249ce308e1dSBarry Smith for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k]; 1250ce308e1dSBarry Smith } 12510acb5bebSBarry Smith if (dfill) { 12520acb5bebSBarry Smith for (k=dfill[j]; k<dfill[j+1]; k++) { 12530acb5bebSBarry Smith cols[cnt++] = i*nc + dfill[k]; 12540acb5bebSBarry Smith } 12550acb5bebSBarry Smith } else { 1256ce308e1dSBarry Smith for (k=0; k<nc; k++) { 1257ce308e1dSBarry Smith cols[cnt++] = i*nc + k; 1258ce308e1dSBarry Smith } 12590acb5bebSBarry Smith } 1260ce308e1dSBarry Smith for (l=0; l<s; l++) { 1261ce308e1dSBarry Smith for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k]; 1262ce308e1dSBarry Smith } 1263ce308e1dSBarry Smith ierr = MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1264ce308e1dSBarry Smith row++; 1265ce308e1dSBarry Smith } 1266ce308e1dSBarry Smith } 1267ce308e1dSBarry Smith /* coupling with process to the right */ 1268ce308e1dSBarry Smith for (i=xs+nx-s; i<xs+nx; i++) { 1269ce308e1dSBarry Smith for (j=0; j<nc; j++) { 1270ce308e1dSBarry Smith cnt = 0; 1271ce308e1dSBarry Smith for (l=0; l<s; l++) { 1272ce308e1dSBarry Smith for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i - s + l)*nc + ofill[k]; 1273ce308e1dSBarry Smith } 12740acb5bebSBarry Smith if (dfill) { 12750acb5bebSBarry Smith for (k=dfill[j]; k<dfill[j+1]; k++) { 12760acb5bebSBarry Smith cols[cnt++] = i*nc + dfill[k]; 12770acb5bebSBarry Smith } 12780acb5bebSBarry Smith } else { 1279ce308e1dSBarry Smith for (k=0; k<nc; k++) { 1280ce308e1dSBarry Smith cols[cnt++] = i*nc + k; 1281ce308e1dSBarry Smith } 12820acb5bebSBarry Smith } 1283ce308e1dSBarry Smith if (rank < size-1) { 1284ce308e1dSBarry Smith for (l=0; l<s; l++) { 1285ce308e1dSBarry Smith for (k=ofill[j]; k<ofill[j+1]; k++) cols[cnt++] = (i + s - l)*nc + ofill[k]; 1286ce308e1dSBarry Smith } 1287ce308e1dSBarry Smith } 1288ce308e1dSBarry Smith ierr = MatSetValues(J,1,&row,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 1289ce308e1dSBarry Smith row++; 1290ce308e1dSBarry Smith } 1291ce308e1dSBarry Smith } 1292c0ab637bSBarry Smith ierr = PetscFree2(values,cols);CHKERRQ(ierr); 1293ce308e1dSBarry Smith ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1294ce308e1dSBarry Smith ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1295189e4007SBarry Smith ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1296ce308e1dSBarry Smith } 1297ce308e1dSBarry Smith PetscFunctionReturn(0); 1298ce308e1dSBarry Smith } 1299ce308e1dSBarry Smith 1300ce308e1dSBarry Smith /* ---------------------------------------------------------------------------------*/ 1301ce308e1dSBarry Smith 1302950540a4SJed Brown PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ(DM da,Mat J) 130347c6ae99SBarry Smith { 130447c6ae99SBarry Smith PetscErrorCode ierr; 130547c6ae99SBarry Smith PetscInt xs,nx,i,i1,slot,gxs,gnx; 13060298fd71SBarry Smith PetscInt m,dim,s,*cols = NULL,nc,*rows = NULL,col,cnt,l; 130747c6ae99SBarry Smith PetscInt istart,iend; 130847c6ae99SBarry Smith PetscScalar *values; 1309bff4a2f0SMatthew G. Knepley DMBoundaryType bx; 131045b6f7e9SBarry Smith ISLocalToGlobalMapping ltog; 131147c6ae99SBarry Smith 131247c6ae99SBarry Smith PetscFunctionBegin; 131347c6ae99SBarry Smith /* 131447c6ae99SBarry Smith nc - number of components per grid point 131547c6ae99SBarry Smith col - number of colors needed in one direction for single component problem 131647c6ae99SBarry Smith 131747c6ae99SBarry Smith */ 13181321219cSEthan Coon ierr = DMDAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&bx,0,0,0);CHKERRQ(ierr); 131947c6ae99SBarry Smith col = 2*s + 1; 132047c6ae99SBarry Smith 1321aa219208SBarry Smith ierr = DMDAGetCorners(da,&xs,0,0,&nx,0,0);CHKERRQ(ierr); 1322aa219208SBarry Smith ierr = DMDAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);CHKERRQ(ierr); 132347c6ae99SBarry Smith 1324f73d5cc4SBarry Smith ierr = MatSetBlockSize(J,nc);CHKERRQ(ierr); 132547c6ae99SBarry Smith ierr = MatSeqAIJSetPreallocation(J,col*nc,0);CHKERRQ(ierr); 132647c6ae99SBarry Smith ierr = MatMPIAIJSetPreallocation(J,col*nc,0,col*nc,0);CHKERRQ(ierr); 132747c6ae99SBarry Smith 13281411c6eeSJed Brown ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 1329784ac674SJed Brown ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 133047c6ae99SBarry Smith 133147c6ae99SBarry Smith /* 133247c6ae99SBarry Smith For each node in the grid: we get the neighbors in the local (on processor ordering 133347c6ae99SBarry Smith that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 133447c6ae99SBarry Smith PETSc ordering. 133547c6ae99SBarry Smith */ 1336fcfd50ebSBarry Smith if (!da->prealloc_only) { 1337dcca6d9dSJed Brown ierr = PetscMalloc2(nc,&rows,col*nc*nc,&cols);CHKERRQ(ierr); 13381795a4d1SJed Brown ierr = PetscCalloc1(col*nc*nc,&values);CHKERRQ(ierr); 133947c6ae99SBarry Smith for (i=xs; i<xs+nx; i++) { 134047c6ae99SBarry Smith istart = PetscMax(-s,gxs - i); 134147c6ae99SBarry Smith iend = PetscMin(s,gxs + gnx - i - 1); 134247c6ae99SBarry Smith slot = i - gxs; 134347c6ae99SBarry Smith 134447c6ae99SBarry Smith cnt = 0; 134547c6ae99SBarry Smith for (l=0; l<nc; l++) { 134647c6ae99SBarry Smith for (i1=istart; i1<iend+1; i1++) { 134747c6ae99SBarry Smith cols[cnt++] = l + nc*(slot + i1); 134847c6ae99SBarry Smith } 134947c6ae99SBarry Smith rows[l] = l + nc*(slot); 135047c6ae99SBarry Smith } 135147c6ae99SBarry Smith ierr = MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 135247c6ae99SBarry Smith } 135347c6ae99SBarry Smith ierr = PetscFree(values);CHKERRQ(ierr); 135447c6ae99SBarry Smith ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 135547c6ae99SBarry Smith ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1356189e4007SBarry Smith ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 135747c6ae99SBarry Smith ierr = PetscFree2(rows,cols);CHKERRQ(ierr); 1358ce308e1dSBarry Smith } 135947c6ae99SBarry Smith PetscFunctionReturn(0); 136047c6ae99SBarry Smith } 136147c6ae99SBarry Smith 1362950540a4SJed Brown PetscErrorCode DMCreateMatrix_DA_2d_MPIBAIJ(DM da,Mat J) 136347c6ae99SBarry Smith { 136447c6ae99SBarry Smith PetscErrorCode ierr; 136547c6ae99SBarry Smith PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 136647c6ae99SBarry Smith PetscInt m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz; 136747c6ae99SBarry Smith PetscInt istart,iend,jstart,jend,ii,jj; 136847c6ae99SBarry Smith MPI_Comm comm; 136947c6ae99SBarry Smith PetscScalar *values; 1370bff4a2f0SMatthew G. Knepley DMBoundaryType bx,by; 1371aa219208SBarry Smith DMDAStencilType st; 137245b6f7e9SBarry Smith ISLocalToGlobalMapping ltog; 137347c6ae99SBarry Smith 137447c6ae99SBarry Smith PetscFunctionBegin; 137547c6ae99SBarry Smith /* 137647c6ae99SBarry Smith nc - number of components per grid point 137747c6ae99SBarry Smith col - number of colors needed in one direction for single component problem 137847c6ae99SBarry Smith */ 13791321219cSEthan Coon ierr = DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);CHKERRQ(ierr); 138047c6ae99SBarry Smith col = 2*s + 1; 138147c6ae99SBarry Smith 1382aa219208SBarry Smith ierr = DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);CHKERRQ(ierr); 1383aa219208SBarry Smith ierr = DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);CHKERRQ(ierr); 138447c6ae99SBarry Smith ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 138547c6ae99SBarry Smith 1386785e854fSJed Brown ierr = PetscMalloc1(col*col*nc*nc,&cols);CHKERRQ(ierr); 138747c6ae99SBarry Smith 13881411c6eeSJed Brown ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 138947c6ae99SBarry Smith 139047c6ae99SBarry Smith /* determine the matrix preallocation information */ 139147c6ae99SBarry Smith ierr = MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);CHKERRQ(ierr); 139247c6ae99SBarry Smith for (i=xs; i<xs+nx; i++) { 1393bff4a2f0SMatthew G. Knepley istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1394bff4a2f0SMatthew G. Knepley iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 139547c6ae99SBarry Smith for (j=ys; j<ys+ny; j++) { 1396bff4a2f0SMatthew G. Knepley jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1397bff4a2f0SMatthew G. Knepley jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 139847c6ae99SBarry Smith slot = i - gxs + gnx*(j - gys); 139947c6ae99SBarry Smith 140047c6ae99SBarry Smith /* Find block columns in block row */ 140147c6ae99SBarry Smith cnt = 0; 140247c6ae99SBarry Smith for (ii=istart; ii<iend+1; ii++) { 140347c6ae99SBarry Smith for (jj=jstart; jj<jend+1; jj++) { 1404aa219208SBarry Smith if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */ 140547c6ae99SBarry Smith cols[cnt++] = slot + ii + gnx*jj; 140647c6ae99SBarry Smith } 140747c6ae99SBarry Smith } 140847c6ae99SBarry Smith } 1409d6e23781SBarry Smith ierr = MatPreallocateSetLocalBlock(ltog,1,&slot,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 141047c6ae99SBarry Smith } 141147c6ae99SBarry Smith } 141247c6ae99SBarry Smith ierr = MatSeqBAIJSetPreallocation(J,nc,0,dnz);CHKERRQ(ierr); 141347c6ae99SBarry Smith ierr = MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);CHKERRQ(ierr); 141447c6ae99SBarry Smith ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 141547c6ae99SBarry Smith 1416784ac674SJed Brown ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 141747c6ae99SBarry Smith 141847c6ae99SBarry Smith /* 141947c6ae99SBarry Smith For each node in the grid: we get the neighbors in the local (on processor ordering 142047c6ae99SBarry Smith that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 142147c6ae99SBarry Smith PETSc ordering. 142247c6ae99SBarry Smith */ 1423fcfd50ebSBarry Smith if (!da->prealloc_only) { 14241795a4d1SJed Brown ierr = PetscCalloc1(col*col*nc*nc,&values);CHKERRQ(ierr); 142547c6ae99SBarry Smith for (i=xs; i<xs+nx; i++) { 1426bff4a2f0SMatthew G. Knepley istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1427bff4a2f0SMatthew G. Knepley iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 142847c6ae99SBarry Smith for (j=ys; j<ys+ny; j++) { 1429bff4a2f0SMatthew G. Knepley jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1430bff4a2f0SMatthew G. Knepley jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 143147c6ae99SBarry Smith slot = i - gxs + gnx*(j - gys); 143247c6ae99SBarry Smith cnt = 0; 143347c6ae99SBarry Smith for (ii=istart; ii<iend+1; ii++) { 143447c6ae99SBarry Smith for (jj=jstart; jj<jend+1; jj++) { 1435aa219208SBarry Smith if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */ 143647c6ae99SBarry Smith cols[cnt++] = slot + ii + gnx*jj; 143747c6ae99SBarry Smith } 143847c6ae99SBarry Smith } 143947c6ae99SBarry Smith } 144047c6ae99SBarry Smith ierr = MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 144147c6ae99SBarry Smith } 144247c6ae99SBarry Smith } 144347c6ae99SBarry Smith ierr = PetscFree(values);CHKERRQ(ierr); 144447c6ae99SBarry Smith ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 144547c6ae99SBarry Smith ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1446189e4007SBarry Smith ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 144747c6ae99SBarry Smith } 144847c6ae99SBarry Smith ierr = PetscFree(cols);CHKERRQ(ierr); 144947c6ae99SBarry Smith PetscFunctionReturn(0); 145047c6ae99SBarry Smith } 145147c6ae99SBarry Smith 1452950540a4SJed Brown PetscErrorCode DMCreateMatrix_DA_3d_MPIBAIJ(DM da,Mat J) 145347c6ae99SBarry Smith { 145447c6ae99SBarry Smith PetscErrorCode ierr; 145547c6ae99SBarry Smith PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 145647c6ae99SBarry Smith PetscInt m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz; 145747c6ae99SBarry Smith PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk; 145847c6ae99SBarry Smith MPI_Comm comm; 145947c6ae99SBarry Smith PetscScalar *values; 1460bff4a2f0SMatthew G. Knepley DMBoundaryType bx,by,bz; 1461aa219208SBarry Smith DMDAStencilType st; 146245b6f7e9SBarry Smith ISLocalToGlobalMapping ltog; 146347c6ae99SBarry Smith 146447c6ae99SBarry Smith PetscFunctionBegin; 146547c6ae99SBarry Smith /* 146647c6ae99SBarry Smith nc - number of components per grid point 146747c6ae99SBarry Smith col - number of colors needed in one direction for single component problem 146847c6ae99SBarry Smith 146947c6ae99SBarry Smith */ 14701321219cSEthan Coon ierr = DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);CHKERRQ(ierr); 147147c6ae99SBarry Smith col = 2*s + 1; 147247c6ae99SBarry Smith 1473aa219208SBarry Smith ierr = DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);CHKERRQ(ierr); 1474aa219208SBarry Smith ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);CHKERRQ(ierr); 147547c6ae99SBarry Smith ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 147647c6ae99SBarry Smith 1477785e854fSJed Brown ierr = PetscMalloc1(col*col*col,&cols);CHKERRQ(ierr); 147847c6ae99SBarry Smith 14791411c6eeSJed Brown ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 148047c6ae99SBarry Smith 148147c6ae99SBarry Smith /* determine the matrix preallocation information */ 148247c6ae99SBarry Smith ierr = MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);CHKERRQ(ierr); 148347c6ae99SBarry Smith for (i=xs; i<xs+nx; i++) { 1484bff4a2f0SMatthew G. Knepley istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1485bff4a2f0SMatthew G. Knepley iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 148647c6ae99SBarry Smith for (j=ys; j<ys+ny; j++) { 1487bff4a2f0SMatthew G. Knepley jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1488bff4a2f0SMatthew G. Knepley jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 148947c6ae99SBarry Smith for (k=zs; k<zs+nz; k++) { 1490bff4a2f0SMatthew G. Knepley kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1491bff4a2f0SMatthew G. Knepley kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 149247c6ae99SBarry Smith 149347c6ae99SBarry Smith slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 149447c6ae99SBarry Smith 149547c6ae99SBarry Smith /* Find block columns in block row */ 149647c6ae99SBarry Smith cnt = 0; 149747c6ae99SBarry Smith for (ii=istart; ii<iend+1; ii++) { 149847c6ae99SBarry Smith for (jj=jstart; jj<jend+1; jj++) { 149947c6ae99SBarry Smith for (kk=kstart; kk<kend+1; kk++) { 1500aa219208SBarry Smith if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 150147c6ae99SBarry Smith cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk; 150247c6ae99SBarry Smith } 150347c6ae99SBarry Smith } 150447c6ae99SBarry Smith } 150547c6ae99SBarry Smith } 1506d6e23781SBarry Smith ierr = MatPreallocateSetLocalBlock(ltog,1,&slot,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 150747c6ae99SBarry Smith } 150847c6ae99SBarry Smith } 150947c6ae99SBarry Smith } 151047c6ae99SBarry Smith ierr = MatSeqBAIJSetPreallocation(J,nc,0,dnz);CHKERRQ(ierr); 151147c6ae99SBarry Smith ierr = MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);CHKERRQ(ierr); 151247c6ae99SBarry Smith ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 151347c6ae99SBarry Smith 1514784ac674SJed Brown ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 151547c6ae99SBarry Smith 151647c6ae99SBarry Smith /* 151747c6ae99SBarry Smith For each node in the grid: we get the neighbors in the local (on processor ordering 151847c6ae99SBarry Smith that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 151947c6ae99SBarry Smith PETSc ordering. 152047c6ae99SBarry Smith */ 1521fcfd50ebSBarry Smith if (!da->prealloc_only) { 15221795a4d1SJed Brown ierr = PetscCalloc1(col*col*col*nc*nc,&values);CHKERRQ(ierr); 152347c6ae99SBarry Smith for (i=xs; i<xs+nx; i++) { 1524bff4a2f0SMatthew G. Knepley istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1525bff4a2f0SMatthew G. Knepley iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 152647c6ae99SBarry Smith for (j=ys; j<ys+ny; j++) { 1527bff4a2f0SMatthew G. Knepley jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1528bff4a2f0SMatthew G. Knepley jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 152947c6ae99SBarry Smith for (k=zs; k<zs+nz; k++) { 1530bff4a2f0SMatthew G. Knepley kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1531bff4a2f0SMatthew G. Knepley kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 153247c6ae99SBarry Smith 153347c6ae99SBarry Smith slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 153447c6ae99SBarry Smith 153547c6ae99SBarry Smith cnt = 0; 153647c6ae99SBarry Smith for (ii=istart; ii<iend+1; ii++) { 153747c6ae99SBarry Smith for (jj=jstart; jj<jend+1; jj++) { 153847c6ae99SBarry Smith for (kk=kstart; kk<kend+1; kk++) { 1539aa219208SBarry Smith if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 154047c6ae99SBarry Smith cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk; 154147c6ae99SBarry Smith } 154247c6ae99SBarry Smith } 154347c6ae99SBarry Smith } 154447c6ae99SBarry Smith } 154547c6ae99SBarry Smith ierr = MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 154647c6ae99SBarry Smith } 154747c6ae99SBarry Smith } 154847c6ae99SBarry Smith } 154947c6ae99SBarry Smith ierr = PetscFree(values);CHKERRQ(ierr); 155047c6ae99SBarry Smith ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 155147c6ae99SBarry Smith ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1552189e4007SBarry Smith ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 155347c6ae99SBarry Smith } 155447c6ae99SBarry Smith ierr = PetscFree(cols);CHKERRQ(ierr); 155547c6ae99SBarry Smith PetscFunctionReturn(0); 155647c6ae99SBarry Smith } 155747c6ae99SBarry Smith 155847c6ae99SBarry Smith /* 155947c6ae99SBarry Smith This helper is for of SBAIJ preallocation, to discard the lower-triangular values which are difficult to 156047c6ae99SBarry Smith identify in the local ordering with periodic domain. 156147c6ae99SBarry Smith */ 156247c6ae99SBarry Smith static PetscErrorCode L2GFilterUpperTriangular(ISLocalToGlobalMapping ltog,PetscInt *row,PetscInt *cnt,PetscInt col[]) 156347c6ae99SBarry Smith { 156447c6ae99SBarry Smith PetscErrorCode ierr; 156547c6ae99SBarry Smith PetscInt i,n; 156647c6ae99SBarry Smith 156747c6ae99SBarry Smith PetscFunctionBegin; 1568d6e23781SBarry Smith ierr = ISLocalToGlobalMappingApplyBlock(ltog,1,row,row);CHKERRQ(ierr); 1569d6e23781SBarry Smith ierr = ISLocalToGlobalMappingApplyBlock(ltog,*cnt,col,col);CHKERRQ(ierr); 157047c6ae99SBarry Smith for (i=0,n=0; i<*cnt; i++) { 157147c6ae99SBarry Smith if (col[i] >= *row) col[n++] = col[i]; 157247c6ae99SBarry Smith } 157347c6ae99SBarry Smith *cnt = n; 157447c6ae99SBarry Smith PetscFunctionReturn(0); 157547c6ae99SBarry Smith } 157647c6ae99SBarry Smith 1577950540a4SJed Brown PetscErrorCode DMCreateMatrix_DA_2d_MPISBAIJ(DM da,Mat J) 157847c6ae99SBarry Smith { 157947c6ae99SBarry Smith PetscErrorCode ierr; 158047c6ae99SBarry Smith PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 158147c6ae99SBarry Smith PetscInt m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz; 158247c6ae99SBarry Smith PetscInt istart,iend,jstart,jend,ii,jj; 158347c6ae99SBarry Smith MPI_Comm comm; 158447c6ae99SBarry Smith PetscScalar *values; 1585bff4a2f0SMatthew G. Knepley DMBoundaryType bx,by; 1586aa219208SBarry Smith DMDAStencilType st; 158745b6f7e9SBarry Smith ISLocalToGlobalMapping ltog; 158847c6ae99SBarry Smith 158947c6ae99SBarry Smith PetscFunctionBegin; 159047c6ae99SBarry Smith /* 159147c6ae99SBarry Smith nc - number of components per grid point 159247c6ae99SBarry Smith col - number of colors needed in one direction for single component problem 159347c6ae99SBarry Smith */ 15941321219cSEthan Coon ierr = DMDAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&bx,&by,0,&st);CHKERRQ(ierr); 159547c6ae99SBarry Smith col = 2*s + 1; 159647c6ae99SBarry Smith 1597aa219208SBarry Smith ierr = DMDAGetCorners(da,&xs,&ys,0,&nx,&ny,0);CHKERRQ(ierr); 1598aa219208SBarry Smith ierr = DMDAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);CHKERRQ(ierr); 159947c6ae99SBarry Smith ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 160047c6ae99SBarry Smith 1601785e854fSJed Brown ierr = PetscMalloc1(col*col*nc*nc,&cols);CHKERRQ(ierr); 160247c6ae99SBarry Smith 16031411c6eeSJed Brown ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 160447c6ae99SBarry Smith 160547c6ae99SBarry Smith /* determine the matrix preallocation information */ 1606eabe889fSLisandro Dalcin ierr = MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);CHKERRQ(ierr); 160747c6ae99SBarry Smith for (i=xs; i<xs+nx; i++) { 1608bff4a2f0SMatthew G. Knepley istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1609bff4a2f0SMatthew G. Knepley iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 161047c6ae99SBarry Smith for (j=ys; j<ys+ny; j++) { 1611bff4a2f0SMatthew G. Knepley jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1612bff4a2f0SMatthew G. Knepley jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 161347c6ae99SBarry Smith slot = i - gxs + gnx*(j - gys); 161447c6ae99SBarry Smith 161547c6ae99SBarry Smith /* Find block columns in block row */ 161647c6ae99SBarry Smith cnt = 0; 161747c6ae99SBarry Smith for (ii=istart; ii<iend+1; ii++) { 161847c6ae99SBarry Smith for (jj=jstart; jj<jend+1; jj++) { 1619aa219208SBarry Smith if (st == DMDA_STENCIL_BOX || !ii || !jj) { 162047c6ae99SBarry Smith cols[cnt++] = slot + ii + gnx*jj; 162147c6ae99SBarry Smith } 162247c6ae99SBarry Smith } 162347c6ae99SBarry Smith } 162445b6f7e9SBarry Smith ierr = L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);CHKERRQ(ierr); 1625d6e23781SBarry Smith ierr = MatPreallocateSymmetricSetBlock(slot,cnt,cols,dnz,onz);CHKERRQ(ierr); 162647c6ae99SBarry Smith } 162747c6ae99SBarry Smith } 162847c6ae99SBarry Smith ierr = MatSeqSBAIJSetPreallocation(J,nc,0,dnz);CHKERRQ(ierr); 162947c6ae99SBarry Smith ierr = MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);CHKERRQ(ierr); 163047c6ae99SBarry Smith ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 163147c6ae99SBarry Smith 1632784ac674SJed Brown ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 163347c6ae99SBarry Smith 163447c6ae99SBarry Smith /* 163547c6ae99SBarry Smith For each node in the grid: we get the neighbors in the local (on processor ordering 163647c6ae99SBarry Smith that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 163747c6ae99SBarry Smith PETSc ordering. 163847c6ae99SBarry Smith */ 1639fcfd50ebSBarry Smith if (!da->prealloc_only) { 16401795a4d1SJed Brown ierr = PetscCalloc1(col*col*nc*nc,&values);CHKERRQ(ierr); 164147c6ae99SBarry Smith for (i=xs; i<xs+nx; i++) { 1642bff4a2f0SMatthew G. Knepley istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1643bff4a2f0SMatthew G. Knepley iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 164447c6ae99SBarry Smith for (j=ys; j<ys+ny; j++) { 1645bff4a2f0SMatthew G. Knepley jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1646bff4a2f0SMatthew G. Knepley jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 164747c6ae99SBarry Smith slot = i - gxs + gnx*(j - gys); 164847c6ae99SBarry Smith 164947c6ae99SBarry Smith /* Find block columns in block row */ 165047c6ae99SBarry Smith cnt = 0; 165147c6ae99SBarry Smith for (ii=istart; ii<iend+1; ii++) { 165247c6ae99SBarry Smith for (jj=jstart; jj<jend+1; jj++) { 1653aa219208SBarry Smith if (st == DMDA_STENCIL_BOX || !ii || !jj) { 165447c6ae99SBarry Smith cols[cnt++] = slot + ii + gnx*jj; 165547c6ae99SBarry Smith } 165647c6ae99SBarry Smith } 165747c6ae99SBarry Smith } 165845b6f7e9SBarry Smith ierr = L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);CHKERRQ(ierr); 165947c6ae99SBarry Smith ierr = MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 166047c6ae99SBarry Smith } 166147c6ae99SBarry Smith } 166247c6ae99SBarry Smith ierr = PetscFree(values);CHKERRQ(ierr); 166347c6ae99SBarry Smith ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 166447c6ae99SBarry Smith ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1665189e4007SBarry Smith ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 166647c6ae99SBarry Smith } 166747c6ae99SBarry Smith ierr = PetscFree(cols);CHKERRQ(ierr); 166847c6ae99SBarry Smith PetscFunctionReturn(0); 166947c6ae99SBarry Smith } 167047c6ae99SBarry Smith 1671950540a4SJed Brown PetscErrorCode DMCreateMatrix_DA_3d_MPISBAIJ(DM da,Mat J) 167247c6ae99SBarry Smith { 167347c6ae99SBarry Smith PetscErrorCode ierr; 167447c6ae99SBarry Smith PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 167547c6ae99SBarry Smith PetscInt m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz; 167647c6ae99SBarry Smith PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk; 167747c6ae99SBarry Smith MPI_Comm comm; 167847c6ae99SBarry Smith PetscScalar *values; 1679bff4a2f0SMatthew G. Knepley DMBoundaryType bx,by,bz; 1680aa219208SBarry Smith DMDAStencilType st; 168145b6f7e9SBarry Smith ISLocalToGlobalMapping ltog; 168247c6ae99SBarry Smith 168347c6ae99SBarry Smith PetscFunctionBegin; 168447c6ae99SBarry Smith /* 168547c6ae99SBarry Smith nc - number of components per grid point 168647c6ae99SBarry Smith col - number of colors needed in one direction for single component problem 168747c6ae99SBarry Smith */ 16881321219cSEthan Coon ierr = DMDAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&bx,&by,&bz,&st);CHKERRQ(ierr); 168947c6ae99SBarry Smith col = 2*s + 1; 169047c6ae99SBarry Smith 1691aa219208SBarry Smith ierr = DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);CHKERRQ(ierr); 1692aa219208SBarry Smith ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);CHKERRQ(ierr); 169347c6ae99SBarry Smith ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 169447c6ae99SBarry Smith 169547c6ae99SBarry Smith /* create the matrix */ 1696785e854fSJed Brown ierr = PetscMalloc1(col*col*col,&cols);CHKERRQ(ierr); 169747c6ae99SBarry Smith 16981411c6eeSJed Brown ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 169947c6ae99SBarry Smith 170047c6ae99SBarry Smith /* determine the matrix preallocation information */ 1701eabe889fSLisandro Dalcin ierr = MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);CHKERRQ(ierr); 170247c6ae99SBarry Smith for (i=xs; i<xs+nx; i++) { 1703bff4a2f0SMatthew G. Knepley istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1704bff4a2f0SMatthew G. Knepley iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 170547c6ae99SBarry Smith for (j=ys; j<ys+ny; j++) { 1706bff4a2f0SMatthew G. Knepley jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1707bff4a2f0SMatthew G. Knepley jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 170847c6ae99SBarry Smith for (k=zs; k<zs+nz; k++) { 1709bff4a2f0SMatthew G. Knepley kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1710bff4a2f0SMatthew G. Knepley kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 171147c6ae99SBarry Smith 171247c6ae99SBarry Smith slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 171347c6ae99SBarry Smith 171447c6ae99SBarry Smith /* Find block columns in block row */ 171547c6ae99SBarry Smith cnt = 0; 171647c6ae99SBarry Smith for (ii=istart; ii<iend+1; ii++) { 171747c6ae99SBarry Smith for (jj=jstart; jj<jend+1; jj++) { 171847c6ae99SBarry Smith for (kk=kstart; kk<kend+1; kk++) { 1719aa219208SBarry Smith if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) { 172047c6ae99SBarry Smith cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk; 172147c6ae99SBarry Smith } 172247c6ae99SBarry Smith } 172347c6ae99SBarry Smith } 172447c6ae99SBarry Smith } 172545b6f7e9SBarry Smith ierr = L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);CHKERRQ(ierr); 1726d6e23781SBarry Smith ierr = MatPreallocateSymmetricSetBlock(slot,cnt,cols,dnz,onz);CHKERRQ(ierr); 172747c6ae99SBarry Smith } 172847c6ae99SBarry Smith } 172947c6ae99SBarry Smith } 173047c6ae99SBarry Smith ierr = MatSeqSBAIJSetPreallocation(J,nc,0,dnz);CHKERRQ(ierr); 173147c6ae99SBarry Smith ierr = MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);CHKERRQ(ierr); 173247c6ae99SBarry Smith ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 173347c6ae99SBarry Smith 1734784ac674SJed Brown ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 173547c6ae99SBarry Smith 173647c6ae99SBarry Smith /* 173747c6ae99SBarry Smith For each node in the grid: we get the neighbors in the local (on processor ordering 173847c6ae99SBarry Smith that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 173947c6ae99SBarry Smith PETSc ordering. 174047c6ae99SBarry Smith */ 1741fcfd50ebSBarry Smith if (!da->prealloc_only) { 17421795a4d1SJed Brown ierr = PetscCalloc1(col*col*col*nc*nc,&values);CHKERRQ(ierr); 174347c6ae99SBarry Smith for (i=xs; i<xs+nx; i++) { 1744bff4a2f0SMatthew G. Knepley istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1745bff4a2f0SMatthew G. Knepley iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 174647c6ae99SBarry Smith for (j=ys; j<ys+ny; j++) { 1747bff4a2f0SMatthew G. Knepley jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1748bff4a2f0SMatthew G. Knepley jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 174947c6ae99SBarry Smith for (k=zs; k<zs+nz; k++) { 1750bff4a2f0SMatthew G. Knepley kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1751bff4a2f0SMatthew G. Knepley kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 175247c6ae99SBarry Smith 175347c6ae99SBarry Smith slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 175447c6ae99SBarry Smith 175547c6ae99SBarry Smith cnt = 0; 175647c6ae99SBarry Smith for (ii=istart; ii<iend+1; ii++) { 175747c6ae99SBarry Smith for (jj=jstart; jj<jend+1; jj++) { 175847c6ae99SBarry Smith for (kk=kstart; kk<kend+1; kk++) { 1759aa219208SBarry Smith if ((st == DMDA_STENCIL_BOX) || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) { 176047c6ae99SBarry Smith cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk; 176147c6ae99SBarry Smith } 176247c6ae99SBarry Smith } 176347c6ae99SBarry Smith } 176447c6ae99SBarry Smith } 176545b6f7e9SBarry Smith ierr = L2GFilterUpperTriangular(ltog,&slot,&cnt,cols);CHKERRQ(ierr); 176647c6ae99SBarry Smith ierr = MatSetValuesBlocked(J,1,&slot,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 176747c6ae99SBarry Smith } 176847c6ae99SBarry Smith } 176947c6ae99SBarry Smith } 177047c6ae99SBarry Smith ierr = PetscFree(values);CHKERRQ(ierr); 177147c6ae99SBarry Smith ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 177247c6ae99SBarry Smith ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1773189e4007SBarry Smith ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 177447c6ae99SBarry Smith } 177547c6ae99SBarry Smith ierr = PetscFree(cols);CHKERRQ(ierr); 177647c6ae99SBarry Smith PetscFunctionReturn(0); 177747c6ae99SBarry Smith } 177847c6ae99SBarry Smith 177947c6ae99SBarry Smith /* ---------------------------------------------------------------------------------*/ 178047c6ae99SBarry Smith 1781950540a4SJed Brown PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ_Fill(DM da,Mat J) 178247c6ae99SBarry Smith { 178347c6ae99SBarry Smith PetscErrorCode ierr; 178447c6ae99SBarry Smith PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny; 1785c0ab637bSBarry Smith PetscInt m,n,dim,s,*cols,k,nc,row,col,cnt, maxcnt = 0,l,p,*dnz,*onz; 1786c1154cd5SBarry Smith PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk,M,N,P; 178747c6ae99SBarry Smith DM_DA *dd = (DM_DA*)da->data; 178847c6ae99SBarry Smith PetscInt ifill_col,*dfill = dd->dfill,*ofill = dd->ofill; 178947c6ae99SBarry Smith MPI_Comm comm; 179047c6ae99SBarry Smith PetscScalar *values; 1791bff4a2f0SMatthew G. Knepley DMBoundaryType bx,by,bz; 179245b6f7e9SBarry Smith ISLocalToGlobalMapping ltog; 1793aa219208SBarry Smith DMDAStencilType st; 1794c1154cd5SBarry Smith PetscBool removedups = PETSC_FALSE; 179547c6ae99SBarry Smith 179647c6ae99SBarry Smith PetscFunctionBegin; 179747c6ae99SBarry Smith /* 179847c6ae99SBarry Smith nc - number of components per grid point 179947c6ae99SBarry Smith col - number of colors needed in one direction for single component problem 180047c6ae99SBarry Smith 180147c6ae99SBarry Smith */ 1802c1154cd5SBarry Smith ierr = DMDAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&bx,&by,&bz,&st);CHKERRQ(ierr); 180347c6ae99SBarry Smith col = 2*s + 1; 1804bff4a2f0SMatthew G. Knepley if (bx == DM_BOUNDARY_PERIODIC && (m % col)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisible\n\ 180547c6ae99SBarry Smith by 2*stencil_width + 1\n"); 1806bff4a2f0SMatthew G. Knepley if (by == DM_BOUNDARY_PERIODIC && (n % col)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisible\n\ 180747c6ae99SBarry Smith by 2*stencil_width + 1\n"); 1808bff4a2f0SMatthew G. Knepley if (bz == DM_BOUNDARY_PERIODIC && (p % col)) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Z is divisible\n\ 180947c6ae99SBarry Smith by 2*stencil_width + 1\n"); 181047c6ae99SBarry Smith 1811c1154cd5SBarry Smith /* 1812c1154cd5SBarry Smith With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times 1813c1154cd5SBarry Smith because of "wrapping" around the end of the domain hitting an entry already counted in the other direction. 1814c1154cd5SBarry Smith */ 1815c1154cd5SBarry Smith if (M == 1 && 2*s >= m) removedups = PETSC_TRUE; 1816c1154cd5SBarry Smith if (N == 1 && 2*s >= n) removedups = PETSC_TRUE; 1817c1154cd5SBarry Smith if (P == 1 && 2*s >= p) removedups = PETSC_TRUE; 1818c1154cd5SBarry Smith 1819aa219208SBarry Smith ierr = DMDAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);CHKERRQ(ierr); 1820aa219208SBarry Smith ierr = DMDAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);CHKERRQ(ierr); 182147c6ae99SBarry Smith ierr = PetscObjectGetComm((PetscObject)da,&comm);CHKERRQ(ierr); 182247c6ae99SBarry Smith 1823785e854fSJed Brown ierr = PetscMalloc1(col*col*col*nc,&cols);CHKERRQ(ierr); 18241411c6eeSJed Brown ierr = DMGetLocalToGlobalMapping(da,<og);CHKERRQ(ierr); 182547c6ae99SBarry Smith 182647c6ae99SBarry Smith /* determine the matrix preallocation information */ 182747c6ae99SBarry Smith ierr = MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);CHKERRQ(ierr); 182847c6ae99SBarry Smith 182906ca8cadSBarry Smith ierr = MatSetBlockSize(J,nc);CHKERRQ(ierr); 183047c6ae99SBarry Smith for (i=xs; i<xs+nx; i++) { 1831bff4a2f0SMatthew G. Knepley istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1832bff4a2f0SMatthew G. Knepley iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 183347c6ae99SBarry Smith for (j=ys; j<ys+ny; j++) { 1834bff4a2f0SMatthew G. Knepley jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1835bff4a2f0SMatthew G. Knepley jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 183647c6ae99SBarry Smith for (k=zs; k<zs+nz; k++) { 1837bff4a2f0SMatthew G. Knepley kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1838bff4a2f0SMatthew G. Knepley kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 183947c6ae99SBarry Smith 184047c6ae99SBarry Smith slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 184147c6ae99SBarry Smith 184247c6ae99SBarry Smith for (l=0; l<nc; l++) { 184347c6ae99SBarry Smith cnt = 0; 184447c6ae99SBarry Smith for (ii=istart; ii<iend+1; ii++) { 184547c6ae99SBarry Smith for (jj=jstart; jj<jend+1; jj++) { 184647c6ae99SBarry Smith for (kk=kstart; kk<kend+1; kk++) { 184747c6ae99SBarry Smith if (ii || jj || kk) { 1848aa219208SBarry Smith if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 18498865f1eaSKarl Rupp for (ifill_col=ofill[l]; ifill_col<ofill[l+1]; ifill_col++) cols[cnt++] = ofill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk); 185047c6ae99SBarry Smith } 185147c6ae99SBarry Smith } else { 185247c6ae99SBarry Smith if (dfill) { 18538865f1eaSKarl Rupp for (ifill_col=dfill[l]; ifill_col<dfill[l+1]; ifill_col++) cols[cnt++] = dfill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk); 185447c6ae99SBarry Smith } else { 18558865f1eaSKarl Rupp for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk); 185647c6ae99SBarry Smith } 185747c6ae99SBarry Smith } 185847c6ae99SBarry Smith } 185947c6ae99SBarry Smith } 186047c6ae99SBarry Smith } 186147c6ae99SBarry Smith row = l + nc*(slot); 1862c0ab637bSBarry Smith maxcnt = PetscMax(maxcnt,cnt); 1863c1154cd5SBarry Smith if (removedups) { 1864c1154cd5SBarry Smith ierr = MatPreallocateSetLocalRemoveDups(ltog,1,&row,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 1865c1154cd5SBarry Smith } else { 1866784ac674SJed Brown ierr = MatPreallocateSetLocal(ltog,1,&row,ltog,cnt,cols,dnz,onz);CHKERRQ(ierr); 186747c6ae99SBarry Smith } 186847c6ae99SBarry Smith } 186947c6ae99SBarry Smith } 187047c6ae99SBarry Smith } 1871c1154cd5SBarry Smith } 187247c6ae99SBarry Smith ierr = MatSeqAIJSetPreallocation(J,0,dnz);CHKERRQ(ierr); 187347c6ae99SBarry Smith ierr = MatMPIAIJSetPreallocation(J,0,dnz,0,onz);CHKERRQ(ierr); 187447c6ae99SBarry Smith ierr = MatPreallocateFinalize(dnz,onz);CHKERRQ(ierr); 1875784ac674SJed Brown ierr = MatSetLocalToGlobalMapping(J,ltog,ltog);CHKERRQ(ierr); 187647c6ae99SBarry Smith 187747c6ae99SBarry Smith /* 187847c6ae99SBarry Smith For each node in the grid: we get the neighbors in the local (on processor ordering 187947c6ae99SBarry Smith that includes the ghost points) then MatSetValuesLocal() maps those indices to the global 188047c6ae99SBarry Smith PETSc ordering. 188147c6ae99SBarry Smith */ 1882fcfd50ebSBarry Smith if (!da->prealloc_only) { 1883c0ab637bSBarry Smith ierr = PetscCalloc1(maxcnt,&values);CHKERRQ(ierr); 188447c6ae99SBarry Smith for (i=xs; i<xs+nx; i++) { 1885bff4a2f0SMatthew G. Knepley istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-i)); 1886bff4a2f0SMatthew G. Knepley iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,m-i-1)); 188747c6ae99SBarry Smith for (j=ys; j<ys+ny; j++) { 1888bff4a2f0SMatthew G. Knepley jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-j)); 1889bff4a2f0SMatthew G. Knepley jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,n-j-1)); 189047c6ae99SBarry Smith for (k=zs; k<zs+nz; k++) { 1891bff4a2f0SMatthew G. Knepley kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s,-k)); 1892bff4a2f0SMatthew G. Knepley kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s,p-k-1)); 189347c6ae99SBarry Smith 189447c6ae99SBarry Smith slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs); 189547c6ae99SBarry Smith 189647c6ae99SBarry Smith for (l=0; l<nc; l++) { 189747c6ae99SBarry Smith cnt = 0; 189847c6ae99SBarry Smith for (ii=istart; ii<iend+1; ii++) { 189947c6ae99SBarry Smith for (jj=jstart; jj<jend+1; jj++) { 190047c6ae99SBarry Smith for (kk=kstart; kk<kend+1; kk++) { 190147c6ae99SBarry Smith if (ii || jj || kk) { 1902aa219208SBarry Smith if ((st == DMDA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/ 19038865f1eaSKarl Rupp for (ifill_col=ofill[l]; ifill_col<ofill[l+1]; ifill_col++) cols[cnt++] = ofill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk); 190447c6ae99SBarry Smith } 190547c6ae99SBarry Smith } else { 190647c6ae99SBarry Smith if (dfill) { 19078865f1eaSKarl Rupp for (ifill_col=dfill[l]; ifill_col<dfill[l+1]; ifill_col++) cols[cnt++] = dfill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk); 190847c6ae99SBarry Smith } else { 19098865f1eaSKarl Rupp for (ifill_col=0; ifill_col<nc; ifill_col++) cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk); 191047c6ae99SBarry Smith } 191147c6ae99SBarry Smith } 191247c6ae99SBarry Smith } 191347c6ae99SBarry Smith } 191447c6ae99SBarry Smith } 191547c6ae99SBarry Smith row = l + nc*(slot); 191647c6ae99SBarry Smith ierr = MatSetValuesLocal(J,1,&row,cnt,cols,values,INSERT_VALUES);CHKERRQ(ierr); 191747c6ae99SBarry Smith } 191847c6ae99SBarry Smith } 191947c6ae99SBarry Smith } 192047c6ae99SBarry Smith } 192147c6ae99SBarry Smith ierr = PetscFree(values);CHKERRQ(ierr); 192247c6ae99SBarry Smith ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 192347c6ae99SBarry Smith ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1924189e4007SBarry Smith ierr = MatSetOption(J,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 192547c6ae99SBarry Smith } 192647c6ae99SBarry Smith ierr = PetscFree(cols);CHKERRQ(ierr); 192747c6ae99SBarry Smith PetscFunctionReturn(0); 192847c6ae99SBarry Smith } 1929