Actual source code: fdda.c
1: #define PETSCDM_DLL
2:
3: #include src/dm/da/daimpl.h
4: #include petscmat.h
7: EXTERN PetscErrorCode DAGetColoring1d_MPIAIJ(DA,ISColoringType,ISColoring *);
8: EXTERN PetscErrorCode DAGetColoring2d_MPIAIJ(DA,ISColoringType,ISColoring *);
9: EXTERN PetscErrorCode DAGetColoring2d_5pt_MPIAIJ(DA,ISColoringType,ISColoring *);
10: EXTERN PetscErrorCode DAGetColoring3d_MPIAIJ(DA,ISColoringType,ISColoring *);
12: /*
13: For ghost i that may be negative or greater than the upper bound this
14: maps it into the 0:m-1 range using periodicity
15: */
16: #define SetInRange(i,m) ((i < 0) ? m+i:((i >= m) ? i-m:i))
20: static PetscErrorCode DASetBlockFills_Private(PetscInt *dfill,PetscInt w,PetscInt **rfill)
21: {
23: PetscInt i,j,nz,*fill;
26: /* count number nonzeros */
27: nz = 0;
28: for (i=0; i<w; i++) {
29: for (j=0; j<w; j++) {
30: if (dfill[w*i+j]) nz++;
31: }
32: }
33: PetscMalloc((nz + w + 1)*sizeof(PetscInt),&fill);
34: /* construct modified CSR storage of nonzero structure */
35: nz = w + 1;
36: for (i=0; i<w; i++) {
37: fill[i] = nz;
38: for (j=0; j<w; j++) {
39: if (dfill[w*i+j]) {
40: fill[nz] = j;
41: nz++;
42: }
43: }
44: }
45: fill[w] = nz;
46:
47: *rfill = fill;
48: return(0);
49: }
53: /*@C
54: DASetBlockFills - Sets the fill pattern in each block for a multi-component problem
55: of the matrix returned by DAGetMatrix().
57: Collective on DA
59: Input Parameter:
60: + da - the distributed array
61: . dfill - the fill pattern in the diagonal block (may be PETSC_NULL, means use dense block)
62: - ofill - the fill pattern in the off-diagonal blocks
65: Level: developer
67: Notes: This only makes sense when you are doing multicomponent problems but using the
68: MPIAIJ matrix format
70: The format for dfill and ofill is a 2 dimensional dof by dof matrix with 1 entries
71: representing coupling and 0 entries for missing coupling. For example
72: $ dfill[3][3] = {1, 0, 0,
73: $ 1, 1, 0,
74: $ 0, 1, 1}
75: means that row 0 is coupled with only itself in the diagonal block, row 1 is coupled with
76: itself and row 0 (in the diagonal block) and row 2 is coupled with itself and row 1 (in the
77: diagonal block.
79: DASetGetMatrix() allows you to provide general code for those more complicated nonzero patterns then
80: can be represented in the dfill, ofill format
82: Contributed by Glenn Hammond
84: .seealso DAGetMatrix(), DASetGetMatrix()
86: @*/
87: PetscErrorCode PETSCDM_DLLEXPORT DASetBlockFills(DA da,PetscInt *dfill,PetscInt *ofill)
88: {
92: if (dfill) {
93: DASetBlockFills_Private(dfill,da->w,&da->dfill);
94: }
95: DASetBlockFills_Private(ofill,da->w,&da->ofill);
96: return(0);
97: }
102: /*@C
103: DAGetColoring - Gets the coloring required for computing the Jacobian via
104: finite differences on a function defined using a stencil on the DA.
106: Collective on DA
108: Input Parameter:
109: + da - the distributed array
110: - ctype - IS_COLORING_LOCAL or IS_COLORING_GHOSTED
112: Output Parameters:
113: . coloring - matrix coloring for use in computing Jacobians (or PETSC_NULL if not needed)
115: Level: advanced
117: Notes: These compute the graph coloring of the graph of A^{T}A. The coloring used
118: for efficient (parallel or thread based) triangular solves etc is NOT yet
119: available.
122: .seealso ISColoringView(), ISColoringGetIS(), MatFDColoringCreate(), ISColoringType, ISColoring
124: @*/
125: PetscErrorCode PETSCDM_DLLEXPORT DAGetColoring(DA da,ISColoringType ctype,ISColoring *coloring)
126: {
128: PetscInt dim;
131: /*
132: m
133: ------------------------------------------------------
134: | |
135: | |
136: | ---------------------- |
137: | | | |
138: n | yn | | |
139: | | | |
140: | .--------------------- |
141: | (xs,ys) xn |
142: | . |
143: | (gxs,gys) |
144: | |
145: -----------------------------------------------------
146: */
148: /*
149: nc - number of components per grid point
150: col - number of colors needed in one direction for single component problem
151:
152: */
153: DAGetInfo(da,&dim,0,0,0,0,0,0,0,0,0,0);
155: /*
156: We do not provide a getcoloring function in the DA operations because
157: the basic DA does not know about matrices. We think of DA as being more
158: more low-level then matrices.
159: */
160: if (dim == 1) {
161: DAGetColoring1d_MPIAIJ(da,ctype,coloring);
162: } else if (dim == 2) {
163: DAGetColoring2d_MPIAIJ(da,ctype,coloring);
164: } else if (dim == 3) {
165: DAGetColoring3d_MPIAIJ(da,ctype,coloring);
166: } else {
167: SETERRQ1(PETSC_ERR_SUP,"Not done for %D dimension, send us mail petsc-maint@mcs.anl.gov for code",dim);
168: }
169: return(0);
170: }
172: /* ---------------------------------------------------------------------------------*/
176: PetscErrorCode DAGetColoring2d_MPIAIJ(DA da,ISColoringType ctype,ISColoring *coloring)
177: {
178: PetscErrorCode ierr;
179: PetscInt xs,ys,nx,ny,i,j,ii,gxs,gys,gnx,gny,m,n,M,N,dim,s,k,nc,col;
180: PetscMPIInt size;
181: MPI_Comm comm;
182: DAPeriodicType wrap;
183: DAStencilType st;
184: ISColoringValue *colors;
187: /*
188: nc - number of components per grid point
189: col - number of colors needed in one direction for single component problem
190:
191: */
192: DAGetInfo(da,&dim,&m,&n,0,&M,&N,0,&nc,&s,&wrap,&st);
193: col = 2*s + 1;
194: DAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
195: DAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
196: PetscObjectGetComm((PetscObject)da,&comm);
197: MPI_Comm_size(comm,&size);
199: /* special case as taught to us by Paul Hovland */
200: if (st == DA_STENCIL_STAR && s == 1) {
201: DAGetColoring2d_5pt_MPIAIJ(da,ctype,coloring);
202: } else {
204: if (DAXPeriodic(wrap) && (m % col)){
205: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisible\n\
206: by 2*stencil_width + 1\n");
207: }
208: if (DAYPeriodic(wrap) && (n % col)){
209: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisible\n\
210: by 2*stencil_width + 1\n");
211: }
212: if (ctype == IS_COLORING_LOCAL) {
213: if (!da->localcoloring) {
214: PetscMalloc(nc*nx*ny*sizeof(ISColoringValue),&colors);
215: ii = 0;
216: for (j=ys; j<ys+ny; j++) {
217: for (i=xs; i<xs+nx; i++) {
218: for (k=0; k<nc; k++) {
219: colors[ii++] = k + nc*((i % col) + col*(j % col));
220: }
221: }
222: }
223: ISColoringCreate(comm,nc*nx*ny,colors,&da->localcoloring);
224: }
225: *coloring = da->localcoloring;
226: } else if (ctype == IS_COLORING_GHOSTED) {
227: if (!da->ghostedcoloring) {
228: PetscMalloc(nc*gnx*gny*sizeof(ISColoringValue),&colors);
229: ii = 0;
230: for (j=gys; j<gys+gny; j++) {
231: for (i=gxs; i<gxs+gnx; i++) {
232: for (k=0; k<nc; k++) {
233: /* the complicated stuff is to handle periodic boundaries */
234: colors[ii++] = k + nc*((SetInRange(i,m) % col) + col*(SetInRange(j,n) % col));
235: }
236: }
237: }
238: ISColoringCreate(comm,nc*gnx*gny,colors,&da->ghostedcoloring);
239: ISColoringSetType(da->ghostedcoloring,IS_COLORING_GHOSTED);
240: }
241: *coloring = da->ghostedcoloring;
242: } else SETERRQ1(PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
243: }
244: ISColoringReference(*coloring);
245: return(0);
246: }
248: /* ---------------------------------------------------------------------------------*/
252: PetscErrorCode DAGetColoring3d_MPIAIJ(DA da,ISColoringType ctype,ISColoring *coloring)
253: {
254: PetscErrorCode ierr;
255: 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;
256: PetscMPIInt size;
257: MPI_Comm comm;
258: DAPeriodicType wrap;
259: DAStencilType st;
260: ISColoringValue *colors;
263: /*
264: nc - number of components per grid point
265: col - number of colors needed in one direction for single component problem
266:
267: */
268: DAGetInfo(da,&dim,&m,&n,&p,&M,&N,&P,&nc,&s,&wrap,&st);
269: col = 2*s + 1;
270: if (DAXPeriodic(wrap) && (m % col)){
271: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisible\n\
272: by 2*stencil_width + 1\n");
273: }
274: if (DAYPeriodic(wrap) && (n % col)){
275: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisible\n\
276: by 2*stencil_width + 1\n");
277: }
278: if (DAZPeriodic(wrap) && (p % col)){
279: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Z is divisible\n\
280: by 2*stencil_width + 1\n");
281: }
283: DAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
284: DAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
285: PetscObjectGetComm((PetscObject)da,&comm);
286: MPI_Comm_size(comm,&size);
288: /* create the coloring */
289: if (ctype == IS_COLORING_LOCAL) {
290: if (!da->localcoloring) {
291: PetscMalloc(nc*nx*ny*nz*sizeof(ISColoringValue),&colors);
292: ii = 0;
293: for (k=zs; k<zs+nz; k++) {
294: for (j=ys; j<ys+ny; j++) {
295: for (i=xs; i<xs+nx; i++) {
296: for (l=0; l<nc; l++) {
297: colors[ii++] = l + nc*((i % col) + col*(j % col) + col*col*(k % col));
298: }
299: }
300: }
301: }
302: ISColoringCreate(comm,nc*nx*ny*nz,colors,&da->localcoloring);
303: }
304: *coloring = da->localcoloring;
305: } else if (ctype == IS_COLORING_GHOSTED) {
306: if (!da->ghostedcoloring) {
307: PetscMalloc(nc*gnx*gny*gnz*sizeof(ISColoringValue),&colors);
308: ii = 0;
309: for (k=gzs; k<gzs+gnz; k++) {
310: for (j=gys; j<gys+gny; j++) {
311: for (i=gxs; i<gxs+gnx; i++) {
312: for (l=0; l<nc; l++) {
313: /* the complicated stuff is to handle periodic boundaries */
314: colors[ii++] = l + nc*((SetInRange(i,m) % col) + col*(SetInRange(j,n) % col) + col*col*(SetInRange(k,p) % col));
315: }
316: }
317: }
318: }
319: ISColoringCreate(comm,nc*gnx*gny*gnz,colors,&da->ghostedcoloring);
320: ISColoringSetType(da->ghostedcoloring,IS_COLORING_GHOSTED);
321: }
322: *coloring = da->ghostedcoloring;
323: } else SETERRQ1(PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
324: ISColoringReference(*coloring);
325: return(0);
326: }
328: /* ---------------------------------------------------------------------------------*/
332: PetscErrorCode DAGetColoring1d_MPIAIJ(DA da,ISColoringType ctype,ISColoring *coloring)
333: {
334: PetscErrorCode ierr;
335: PetscInt xs,nx,i,i1,gxs,gnx,l,m,M,dim,s,nc,col;
336: PetscMPIInt size;
337: MPI_Comm comm;
338: DAPeriodicType wrap;
339: ISColoringValue *colors;
342: /*
343: nc - number of components per grid point
344: col - number of colors needed in one direction for single component problem
345:
346: */
347: DAGetInfo(da,&dim,&m,0,0,&M,0,0,&nc,&s,&wrap,0);
348: col = 2*s + 1;
350: if (DAXPeriodic(wrap) && (m % col)) {
351: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points is divisible\n\
352: by 2*stencil_width + 1\n");
353: }
355: DAGetCorners(da,&xs,0,0,&nx,0,0);
356: DAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);
357: PetscObjectGetComm((PetscObject)da,&comm);
358: MPI_Comm_size(comm,&size);
360: /* create the coloring */
361: if (ctype == IS_COLORING_LOCAL) {
362: if (!da->localcoloring) {
363: PetscMalloc(nc*nx*sizeof(ISColoringValue),&colors);
364: i1 = 0;
365: for (i=xs; i<xs+nx; i++) {
366: for (l=0; l<nc; l++) {
367: colors[i1++] = l + nc*(i % col);
368: }
369: }
370: ISColoringCreate(comm,nc*nx,colors,&da->localcoloring);
371: }
372: *coloring = da->localcoloring;
373: } else if (ctype == IS_COLORING_GHOSTED) {
374: if (!da->ghostedcoloring) {
375: PetscMalloc(nc*gnx*sizeof(ISColoringValue),&colors);
376: i1 = 0;
377: for (i=gxs; i<gxs+gnx; i++) {
378: for (l=0; l<nc; l++) {
379: /* the complicated stuff is to handle periodic boundaries */
380: colors[i1++] = l + nc*(SetInRange(i,m) % col);
381: }
382: }
383: ISColoringCreate(comm,nc*gnx,colors,&da->ghostedcoloring);
384: ISColoringSetType(da->ghostedcoloring,IS_COLORING_GHOSTED);
385: }
386: *coloring = da->ghostedcoloring;
387: } else SETERRQ1(PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
388: ISColoringReference(*coloring);
389: return(0);
390: }
394: PetscErrorCode DAGetColoring2d_5pt_MPIAIJ(DA da,ISColoringType ctype,ISColoring *coloring)
395: {
396: PetscErrorCode ierr;
397: PetscInt xs,ys,nx,ny,i,j,ii,gxs,gys,gnx,gny,m,n,dim,s,k,nc;
398: MPI_Comm comm;
399: DAPeriodicType wrap;
400: ISColoringValue *colors;
403: /*
404: nc - number of components per grid point
405: col - number of colors needed in one direction for single component problem
406:
407: */
408: DAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&wrap,0);
409: DAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
410: DAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
411: PetscObjectGetComm((PetscObject)da,&comm);
413: if (DAXPeriodic(wrap) && (m % 5)){
414: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisible\n\
415: by 5\n");
416: }
417: if (DAYPeriodic(wrap) && (n % 5)){
418: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisible\n\
419: by 5\n");
420: }
422: /* create the coloring */
423: if (ctype == IS_COLORING_LOCAL) {
424: if (!da->localcoloring) {
425: PetscMalloc(nc*nx*ny*sizeof(ISColoringValue),&colors);
426: ii = 0;
427: for (j=ys; j<ys+ny; j++) {
428: for (i=xs; i<xs+nx; i++) {
429: for (k=0; k<nc; k++) {
430: colors[ii++] = k + nc*((3*j+i) % 5);
431: }
432: }
433: }
434: ISColoringCreate(comm,nc*nx*ny,colors,&da->localcoloring);
435: }
436: *coloring = da->localcoloring;
437: } else if (ctype == IS_COLORING_GHOSTED) {
438: if (!da->ghostedcoloring) {
439: PetscMalloc(nc*gnx*gny*sizeof(ISColoringValue),&colors);
440: ii = 0;
441: for (j=gys; j<gys+gny; j++) {
442: for (i=gxs; i<gxs+gnx; i++) {
443: for (k=0; k<nc; k++) {
444: colors[ii++] = k + nc*((3*SetInRange(j,n) + SetInRange(i,m)) % 5);
445: }
446: }
447: }
448: ISColoringCreate(comm,nc*gnx*gny,colors,&da->ghostedcoloring);
449: ISColoringSetType(da->ghostedcoloring,IS_COLORING_GHOSTED);
450: }
451: *coloring = da->ghostedcoloring;
452: } else SETERRQ1(PETSC_ERR_ARG_WRONG,"Unknown ISColoringType %d",(int)ctype);
453: return(0);
454: }
456: /* =========================================================================== */
457: EXTERN PetscErrorCode DAGetMatrix1d_MPIAIJ(DA,Mat);
458: EXTERN PetscErrorCode DAGetMatrix2d_MPIAIJ(DA,Mat);
459: EXTERN PetscErrorCode DAGetMatrix2d_MPIAIJ_Fill(DA,Mat);
460: EXTERN PetscErrorCode DAGetMatrix3d_MPIAIJ(DA,Mat);
461: EXTERN PetscErrorCode DAGetMatrix3d_MPIAIJ_Fill(DA,Mat);
462: EXTERN PetscErrorCode DAGetMatrix2d_MPIBAIJ(DA,Mat);
463: EXTERN PetscErrorCode DAGetMatrix3d_MPIBAIJ(DA,Mat);
464: EXTERN PetscErrorCode DAGetMatrix2d_MPISBAIJ(DA,Mat);
465: EXTERN PetscErrorCode DAGetMatrix3d_MPISBAIJ(DA,Mat);
469: /*@C
470: DAGetMatrix - Creates a matrix with the correct parallel layout and nonzero structure required for
471: computing the Jacobian on a function defined using the stencil set in the DA.
473: Collective on DA
475: Input Parameter:
476: + da - the distributed array
477: - mtype - Supported types are MATSEQAIJ, MATMPIAIJ, MATSEQBAIJ, MATMPIBAIJ, MATSEQSBAIJ, MATMPISBAIJ,
478: or any type which inherits from one of these (such as MATAIJ, MATLUSOL, etc.).
480: Output Parameters:
481: . J - matrix with the correct nonzero structure
482: (obviously without the correct Jacobian values)
484: Level: advanced
486: Notes: This properly preallocates the number of nonzeros in the sparse matrix so you
487: do not need to do it yourself.
489: .seealso ISColoringView(), ISColoringGetIS(), MatFDColoringCreate(), DASetBlockFills()
491: @*/
492: PetscErrorCode PETSCDM_DLLEXPORT DAGetMatrix(DA da, MatType mtype,Mat *J)
493: {
495: PetscInt dim,dof,nx,ny,nz,dims[3],starts[3];
496: Mat A;
497: MPI_Comm comm;
498: MatType Atype;
499: void (*aij)(void)=PETSC_NULL,(*baij)(void)=PETSC_NULL,(*sbaij)(void)=PETSC_NULL;
502: /*
503: m
504: ------------------------------------------------------
505: | |
506: | |
507: | ---------------------- |
508: | | | |
509: n | ny | | |
510: | | | |
511: | .--------------------- |
512: | (xs,ys) nx |
513: | . |
514: | (gxs,gys) |
515: | |
516: -----------------------------------------------------
517: */
519: /*
520: nc - number of components per grid point
521: col - number of colors needed in one direction for single component problem
522:
523: */
524: DAGetInfo(da,&dim,0,0,0,0,0,0,&dof,0,0,0);
525: DAGetCorners(da,0,0,0,&nx,&ny,&nz);
526: PetscObjectGetComm((PetscObject)da,&comm);
527: MatCreate(comm,&A);
528: MatSetSizes(A,dof*nx*ny*nz,dof*nx*ny*nz,PETSC_DECIDE,PETSC_DECIDE);
529: MatSetType(A,mtype);
530: MatSetFromOptions(A);
531: MatGetType(A,&Atype);
532: /*
533: We do not provide a getmatrix function in the DA operations because
534: the basic DA does not know about matrices. We think of DA as being more
535: more low-level than matrices. This is kind of cheating but, cause sometimes
536: we think of DA has higher level than matrices.
538: We could switch based on Atype (or mtype), but we do not since the
539: specialized setting routines depend only the particular preallocation
540: details of the matrix, not the type itself.
541: */
542: PetscObjectQueryFunction((PetscObject)A,"MatMPIAIJSetPreallocation_C",&aij);
543: if (!aij) {
544: PetscObjectQueryFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",&aij);
545: }
546: if (!aij) {
547: PetscObjectQueryFunction((PetscObject)A,"MatMPIBAIJSetPreallocation_C",&baij);
548: if (!baij) {
549: PetscObjectQueryFunction((PetscObject)A,"MatSeqBAIJSetPreallocation_C",&baij);
550: }
551: if (!baij){
552: PetscObjectQueryFunction((PetscObject)A,"MatMPISBAIJSetPreallocation_C",&sbaij);
553: if (!sbaij) {
554: PetscObjectQueryFunction((PetscObject)A,"MatSeqSBAIJSetPreallocation_C",&sbaij);
555: }
556: if (!sbaij) SETERRQ2(PETSC_ERR_SUP,"Not implemented for the matrix type: %s!\n" \
557: "Send mail to petsc-maint@mcs.anl.gov for code",dim,Atype);
558: }
559: }
560: if (aij) {
561: if (dim == 1) {
562: DAGetMatrix1d_MPIAIJ(da,A);
563: } else if (dim == 2) {
564: if (da->ofill) {
565: DAGetMatrix2d_MPIAIJ_Fill(da,A);
566: } else {
567: DAGetMatrix2d_MPIAIJ(da,A);
568: }
569: } else if (dim == 3) {
570: if (da->ofill) {
571: DAGetMatrix3d_MPIAIJ_Fill(da,A);
572: } else {
573: DAGetMatrix3d_MPIAIJ(da,A);
574: }
575: }
576: } else if (baij) {
577: if (dim == 2) {
578: DAGetMatrix2d_MPIBAIJ(da,A);
579: } else if (dim == 3) {
580: DAGetMatrix3d_MPIBAIJ(da,A);
581: } else {
582: SETERRQ2(PETSC_ERR_SUP,"Not implemented for %D dimension and Matrix Type: %s!\n" \
583: "Send mail to petsc-maint@mcs.anl.gov for code",dim,Atype);
584: }
585: } else if (sbaij) {
586: if (dim == 2) {
587: DAGetMatrix2d_MPISBAIJ(da,A);
588: } else if (dim == 3) {
589: DAGetMatrix3d_MPISBAIJ(da,A);
590: } else {
591: SETERRQ2(PETSC_ERR_SUP,"Not implemented for %D dimension and Matrix Type: %s!\n" \
592: "Send mail to petsc-maint@mcs.anl.gov for code",dim,Atype);
593: }
594: }
595: DAGetGhostCorners(da,&starts[0],&starts[1],&starts[2],&dims[0],&dims[1],&dims[2]);
596: MatSetStencil(A,dim,dims,starts,dof);
597: *J = A;
598: return(0);
599: }
601: /* ---------------------------------------------------------------------------------*/
604: PetscErrorCode DAGetMatrix2d_MPIAIJ(DA da,Mat J)
605: {
606: PetscErrorCode ierr;
607: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny,m,n,dim,s,*cols,k,nc,*rows,col,cnt,l,p;
608: PetscInt lstart,lend,pstart,pend,*dnz,*onz;
609: MPI_Comm comm;
610: PetscScalar *values;
611: DAPeriodicType wrap;
612: ISLocalToGlobalMapping ltog;
613: DAStencilType st;
616: /*
617: nc - number of components per grid point
618: col - number of colors needed in one direction for single component problem
619:
620: */
621: DAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&wrap,&st);
622: col = 2*s + 1;
623: DAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
624: DAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
625: PetscObjectGetComm((PetscObject)da,&comm);
627: PetscMalloc(col*col*nc*nc*sizeof(PetscScalar),&values);
628: PetscMemzero(values,col*col*nc*nc*sizeof(PetscScalar));
629: PetscMalloc(nc*sizeof(PetscInt),&rows);
630: PetscMalloc(col*col*nc*nc*sizeof(PetscInt),&cols);
631: DAGetISLocalToGlobalMapping(da,<og);
632:
633: /* determine the matrix preallocation information */
634: MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
635: for (i=xs; i<xs+nx; i++) {
637: pstart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
638: pend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
640: for (j=ys; j<ys+ny; j++) {
641: slot = i - gxs + gnx*(j - gys);
643: lstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
644: lend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
646: cnt = 0;
647: for (k=0; k<nc; k++) {
648: for (l=lstart; l<lend+1; l++) {
649: for (p=pstart; p<pend+1; p++) {
650: if ((st == DA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */
651: cols[cnt++] = k + nc*(slot + gnx*l + p);
652: }
653: }
654: }
655: rows[k] = k + nc*(slot);
656: }
657: MatPreallocateSetLocal(ltog,nc,rows,cnt,cols,dnz,onz);
658: }
659: }
660: MatSeqAIJSetPreallocation(J,0,dnz);
661: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
662: MatSetBlockSize(J,nc);
663: MatPreallocateFinalize(dnz,onz);
665: MatSetLocalToGlobalMapping(J,ltog);
667: /*
668: For each node in the grid: we get the neighbors in the local (on processor ordering
669: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
670: PETSc ordering.
671: */
672: for (i=xs; i<xs+nx; i++) {
673:
674: pstart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
675: pend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
676:
677: for (j=ys; j<ys+ny; j++) {
678: slot = i - gxs + gnx*(j - gys);
679:
680: lstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
681: lend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
683: cnt = 0;
684: for (k=0; k<nc; k++) {
685: for (l=lstart; l<lend+1; l++) {
686: for (p=pstart; p<pend+1; p++) {
687: if ((st == DA_STENCIL_BOX) || (!l || !p)) { /* entries on star have either l = 0 or p = 0 */
688: cols[cnt++] = k + nc*(slot + gnx*l + p);
689: }
690: }
691: }
692: rows[k] = k + nc*(slot);
693: }
694: MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
695: }
696: }
697: PetscFree(values);
698: PetscFree(rows);
699: PetscFree(cols);
700: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
701: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
702: return(0);
703: }
707: PetscErrorCode DAGetMatrix2d_MPIAIJ_Fill(DA da,Mat J)
708: {
709: PetscErrorCode ierr;
710: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
711: PetscInt m,n,dim,s,*cols,k,nc,*rows,col,cnt,l,p;
712: PetscInt lstart,lend,pstart,pend,*dnz,*onz;
713: PetscInt ifill_col,*ofill = da->ofill, *dfill = da->dfill;
714: MPI_Comm comm;
715: PetscScalar *values;
716: DAPeriodicType wrap;
717: ISLocalToGlobalMapping ltog;
718: DAStencilType st;
721: /*
722: nc - number of components per grid point
723: col - number of colors needed in one direction for single component problem
724:
725: */
726: DAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&wrap,&st);
727: col = 2*s + 1;
728: DAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
729: DAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
730: PetscObjectGetComm((PetscObject)da,&comm);
732: PetscMalloc(col*col*nc*nc*sizeof(PetscScalar),&values);
733: PetscMemzero(values,col*col*nc*nc*sizeof(PetscScalar));
734: PetscMalloc(nc*sizeof(PetscInt),&rows);
735: PetscMalloc(col*col*nc*nc*sizeof(PetscInt),&cols);
736: DAGetISLocalToGlobalMapping(da,<og);
737:
738: /* determine the matrix preallocation information */
739: MatPreallocateInitialize(comm,nc*nx*ny,nc*nx*ny,dnz,onz);
740: for (i=xs; i<xs+nx; i++) {
742: pstart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
743: pend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
745: for (j=ys; j<ys+ny; j++) {
746: slot = i - gxs + gnx*(j - gys);
748: lstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
749: lend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
751: for (k=0; k<nc; k++) {
752: cnt = 0;
753: for (l=lstart; l<lend+1; l++) {
754: for (p=pstart; p<pend+1; p++) {
755: if (l || p) {
756: if ((st == DA_STENCIL_BOX) || (!l || !p)) { /* entries on star */
757: for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++)
758: cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p);
759: }
760: } else {
761: if (dfill) {
762: for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++)
763: cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p);
764: } else {
765: for (ifill_col=0; ifill_col<nc; ifill_col++)
766: cols[cnt++] = ifill_col + nc*(slot + gnx*l + p);
767: }
768: }
769: }
770: }
771: rows[0] = k + nc*(slot);
772: MatPreallocateSetLocal(ltog,1,rows,cnt,cols,dnz,onz);
773: }
774: }
775: }
776: MatSeqAIJSetPreallocation(J,0,dnz);
777: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
778: MatPreallocateFinalize(dnz,onz);
779: MatSetLocalToGlobalMapping(J,ltog);
781: /*
782: For each node in the grid: we get the neighbors in the local (on processor ordering
783: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
784: PETSc ordering.
785: */
786: for (i=xs; i<xs+nx; i++) {
787:
788: pstart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
789: pend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
790:
791: for (j=ys; j<ys+ny; j++) {
792: slot = i - gxs + gnx*(j - gys);
793:
794: lstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
795: lend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
797: for (k=0; k<nc; k++) {
798: cnt = 0;
799: for (l=lstart; l<lend+1; l++) {
800: for (p=pstart; p<pend+1; p++) {
801: if (l || p) {
802: if ((st == DA_STENCIL_BOX) || (!l || !p)) { /* entries on star */
803: for (ifill_col=ofill[k]; ifill_col<ofill[k+1]; ifill_col++)
804: cols[cnt++] = ofill[ifill_col] + nc*(slot + gnx*l + p);
805: }
806: } else {
807: if (dfill) {
808: for (ifill_col=dfill[k]; ifill_col<dfill[k+1]; ifill_col++)
809: cols[cnt++] = dfill[ifill_col] + nc*(slot + gnx*l + p);
810: } else {
811: for (ifill_col=0; ifill_col<nc; ifill_col++)
812: cols[cnt++] = ifill_col + nc*(slot + gnx*l + p);
813: }
814: }
815: }
816: }
817: rows[0] = k + nc*(slot);
818: MatSetValuesLocal(J,1,rows,cnt,cols,values,INSERT_VALUES);
819: }
820: }
821: }
822: PetscFree(values);
823: PetscFree(rows);
824: PetscFree(cols);
825: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
826: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
827: return(0);
828: }
830: /* ---------------------------------------------------------------------------------*/
834: PetscErrorCode DAGetMatrix3d_MPIAIJ(DA da,Mat J)
835: {
836: PetscErrorCode ierr;
837: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
838: PetscInt m,n,dim,s,*cols,k,nc,*rows,col,cnt,l,p,*dnz,*onz;
839: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
840: MPI_Comm comm;
841: PetscScalar *values;
842: DAPeriodicType wrap;
843: ISLocalToGlobalMapping ltog;
844: DAStencilType st;
847: /*
848: nc - number of components per grid point
849: col - number of colors needed in one direction for single component problem
850:
851: */
852: DAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&wrap,&st);
853: col = 2*s + 1;
855: DAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
856: DAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
857: PetscObjectGetComm((PetscObject)da,&comm);
859: PetscMalloc(col*col*col*nc*nc*nc*sizeof(PetscScalar),&values);
860: PetscMemzero(values,col*col*col*nc*nc*nc*sizeof(PetscScalar));
861: PetscMalloc(nc*sizeof(PetscInt),&rows);
862: PetscMalloc(col*col*col*nc*sizeof(PetscInt),&cols);
863: DAGetISLocalToGlobalMapping(da,<og);
865: /* determine the matrix preallocation information */
866: MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);
867: for (i=xs; i<xs+nx; i++) {
868: istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
869: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
870: for (j=ys; j<ys+ny; j++) {
871: jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
872: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
873: for (k=zs; k<zs+nz; k++) {
874: kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
875: kend = DAZPeriodic(wrap) ? s : (PetscMin(s,p-k-1));
876:
877: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
878:
879: cnt = 0;
880: for (l=0; l<nc; l++) {
881: for (ii=istart; ii<iend+1; ii++) {
882: for (jj=jstart; jj<jend+1; jj++) {
883: for (kk=kstart; kk<kend+1; kk++) {
884: if ((st == DA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
885: cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
886: }
887: }
888: }
889: }
890: rows[l] = l + nc*(slot);
891: }
892: MatPreallocateSetLocal(ltog,nc,rows,cnt,cols,dnz,onz);
893: }
894: }
895: }
896: MatSeqAIJSetPreallocation(J,0,dnz);
897: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
898: MatPreallocateFinalize(dnz,onz);
899: MatSetBlockSize(J,nc);CHKERRQ(ierr)
900: MatSetLocalToGlobalMapping(J,ltog);
902: /*
903: For each node in the grid: we get the neighbors in the local (on processor ordering
904: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
905: PETSc ordering.
906: */
907: for (i=xs; i<xs+nx; i++) {
908: istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
909: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
910: for (j=ys; j<ys+ny; j++) {
911: jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
912: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
913: for (k=zs; k<zs+nz; k++) {
914: kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
915: kend = DAZPeriodic(wrap) ? s : (PetscMin(s,p-k-1));
916:
917: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
918:
919: cnt = 0;
920: for (l=0; l<nc; l++) {
921: for (ii=istart; ii<iend+1; ii++) {
922: for (jj=jstart; jj<jend+1; jj++) {
923: for (kk=kstart; kk<kend+1; kk++) {
924: if ((st == DA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
925: cols[cnt++] = l + nc*(slot + ii + gnx*jj + gnx*gny*kk);
926: }
927: }
928: }
929: }
930: rows[l] = l + nc*(slot);
931: }
932: MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
933: }
934: }
935: }
936: PetscFree(values);
937: PetscFree(rows);
938: PetscFree(cols);
939: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
940: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
941: return(0);
942: }
944: /* ---------------------------------------------------------------------------------*/
948: PetscErrorCode DAGetMatrix1d_MPIAIJ(DA da,Mat J)
949: {
950: PetscErrorCode ierr;
951: PetscInt xs,nx,i,i1,slot,gxs,gnx;
952: PetscInt m,dim,s,*cols,nc,*rows,col,cnt,l;
953: PetscInt istart,iend;
954: PetscScalar *values;
955: DAPeriodicType wrap;
956: ISLocalToGlobalMapping ltog;
959: /*
960: nc - number of components per grid point
961: col - number of colors needed in one direction for single component problem
962:
963: */
964: DAGetInfo(da,&dim,&m,0,0,0,0,0,&nc,&s,&wrap,0);
965: col = 2*s + 1;
967: DAGetCorners(da,&xs,0,0,&nx,0,0);
968: DAGetGhostCorners(da,&gxs,0,0,&gnx,0,0);
970: MatSeqAIJSetPreallocation(J,col*nc,0);
971: MatMPIAIJSetPreallocation(J,col*nc,0,0,0);
972: MatSetBlockSize(J,nc);
973: PetscMalloc(col*nc*nc*sizeof(PetscScalar),&values);
974: PetscMemzero(values,col*nc*nc*sizeof(PetscScalar));
975: PetscMalloc(nc*sizeof(PetscInt),&rows);
976: PetscMalloc(col*nc*sizeof(PetscInt),&cols);
977:
978: DAGetISLocalToGlobalMapping(da,<og);
979: MatSetLocalToGlobalMapping(J,ltog);
980:
981: /*
982: For each node in the grid: we get the neighbors in the local (on processor ordering
983: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
984: PETSc ordering.
985: */
986: for (i=xs; i<xs+nx; i++) {
987: istart = PetscMax(-s,gxs - i);
988: iend = PetscMin(s,gxs + gnx - i - 1);
989: slot = i - gxs;
990:
991: cnt = 0;
992: for (l=0; l<nc; l++) {
993: for (i1=istart; i1<iend+1; i1++) {
994: cols[cnt++] = l + nc*(slot + i1);
995: }
996: rows[l] = l + nc*(slot);
997: }
998: MatSetValuesLocal(J,nc,rows,cnt,cols,values,INSERT_VALUES);
999: }
1000: PetscFree(values);
1001: PetscFree(rows);
1002: PetscFree(cols);
1003: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1004: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1005: return(0);
1006: }
1010: PetscErrorCode DAGetMatrix2d_MPIBAIJ(DA da,Mat J)
1011: {
1012: PetscErrorCode ierr;
1013: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1014: PetscInt m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz;
1015: PetscInt istart,iend,jstart,jend,ii,jj;
1016: MPI_Comm comm;
1017: PetscScalar *values;
1018: DAPeriodicType wrap;
1019: DAStencilType st;
1020: ISLocalToGlobalMapping ltog,ltogb;
1023: /*
1024: nc - number of components per grid point
1025: col - number of colors needed in one direction for single component problem
1026: */
1027: DAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&wrap,&st);
1028: if (wrap != DA_NONPERIODIC) SETERRQ(PETSC_ERR_SUP,"Currently no support for periodic");
1029: col = 2*s + 1;
1031: DAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
1032: DAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
1033: PetscObjectGetComm((PetscObject)da,&comm);
1035: PetscMalloc(col*col*nc*nc*sizeof(PetscScalar),&values);
1036: PetscMemzero(values,col*col*nc*nc*sizeof(PetscScalar));
1037: PetscMalloc(col*col*nc*nc*sizeof(PetscInt),&cols);
1039: DAGetISLocalToGlobalMapping(da,<og);
1040: DAGetISLocalToGlobalMappingBlck(da,<ogb);
1042: /* determine the matrix preallocation information */
1043: MatPreallocateInitialize(comm,nx*ny,nx*ny,dnz,onz);
1044: for (i=xs; i<xs+nx; i++) {
1045: istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
1046: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
1047: for (j=ys; j<ys+ny; j++) {
1048: jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
1049: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
1050: slot = i - gxs + gnx*(j - gys);
1052: /* Find block columns in block row */
1053: cnt = 0;
1054: if (st == DA_STENCIL_BOX) { /* if using BOX stencil */
1055: for (ii=istart; ii<iend+1; ii++) {
1056: for (jj=jstart; jj<jend+1; jj++) {
1057: cols[cnt++] = slot + ii + gnx*jj;
1058: }
1059: }
1060: } else { /* Star stencil */
1061: cnt = 0;
1062: for (ii=istart; ii<iend+1; ii++) {
1063: if (ii) { /* jj must be zero */
1064: cols[cnt++] = slot + ii;
1065: } else {
1066: for (jj=jstart; jj<jend+1; jj++) { /* ii must be zero */
1067: cols[cnt++] = slot + gnx*jj;
1068: }
1069: }
1070: }
1071: }
1072: MatPreallocateSetLocal(ltogb,1,&slot,cnt,cols,dnz,onz);
1073: }
1074: }
1075: MatSeqBAIJSetPreallocation(J,nc,0,dnz);
1076: MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1077: MatPreallocateFinalize(dnz,onz);
1079: MatSetLocalToGlobalMapping(J,ltog);
1080: MatSetLocalToGlobalMappingBlock(J,ltogb);
1082: /*
1083: For each node in the grid: we get the neighbors in the local (on processor ordering
1084: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1085: PETSc ordering.
1086: */
1087: for (i=xs; i<xs+nx; i++) {
1088: istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
1089: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
1090: for (j=ys; j<ys+ny; j++) {
1091: jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
1092: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
1093: slot = i - gxs + gnx*(j - gys);
1094: cnt = 0;
1095: if (st == DA_STENCIL_BOX) { /* if using BOX stencil */
1096: for (ii=istart; ii<iend+1; ii++) {
1097: for (jj=jstart; jj<jend+1; jj++) {
1098: cols[cnt++] = slot + ii + gnx*jj;
1099: }
1100: }
1101: } else { /* Star stencil */
1102: cnt = 0;
1103: for (ii=istart; ii<iend+1; ii++) {
1104: if (ii) { /* jj must be zero */
1105: cols[cnt++] = slot + ii;
1106: } else {
1107: for (jj=jstart; jj<jend+1; jj++) {/* ii must be zero */
1108: cols[cnt++] = slot + gnx*jj;
1109: }
1110: }
1111: }
1112: }
1113: MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1114: }
1115: }
1116: PetscFree(values);
1117: PetscFree(cols);
1118: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1119: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1120: return(0);
1121: }
1125: PetscErrorCode DAGetMatrix3d_MPIBAIJ(DA da,Mat J)
1126: {
1127: PetscErrorCode ierr;
1128: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1129: PetscInt m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
1130: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
1131: MPI_Comm comm;
1132: PetscScalar *values;
1133: DAPeriodicType wrap;
1134: DAStencilType st;
1135: ISLocalToGlobalMapping ltog,ltogb;
1138: /*
1139: nc - number of components per grid point
1140: col - number of colors needed in one direction for single component problem
1141:
1142: */
1143: DAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&wrap,&st);
1144: if (wrap != DA_NONPERIODIC) SETERRQ(PETSC_ERR_SUP,"Currently no support for periodic");
1145: col = 2*s + 1;
1147: DAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1148: DAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1149: PetscObjectGetComm((PetscObject)da,&comm);
1151: PetscMalloc(col*col*col*nc*nc*sizeof(PetscScalar),&values);
1152: PetscMemzero(values,col*col*col*nc*nc*sizeof(PetscScalar));
1153: PetscMalloc(col*col*col*sizeof(PetscInt),&cols);
1155: DAGetISLocalToGlobalMapping(da,<og);
1156: DAGetISLocalToGlobalMappingBlck(da,<ogb);
1158: /* determine the matrix preallocation information */
1159: MatPreallocateInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
1160: for (i=xs; i<xs+nx; i++) {
1161: istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
1162: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
1163: for (j=ys; j<ys+ny; j++) {
1164: jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
1165: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
1166: for (k=zs; k<zs+nz; k++) {
1167: kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
1168: kend = DAZPeriodic(wrap) ? s : (PetscMin(s,p-k-1));
1170: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1172: /* Find block columns in block row */
1173: cnt = 0;
1174: if (st == DA_STENCIL_BOX) { /* if using BOX stencil */
1175: for (ii=istart; ii<iend+1; ii++) {
1176: for (jj=jstart; jj<jend+1; jj++) {
1177: for (kk=kstart; kk<kend+1; kk++) {
1178: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1179: }
1180: }
1181: }
1182: } else { /* Star stencil */
1183: cnt = 0;
1184: for (ii=istart; ii<iend+1; ii++) {
1185: if (ii) {
1186: /* jj and kk must be zero */
1187: /* cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk; */
1188: cols[cnt++] = slot + ii;
1189: } else {
1190: for (jj=jstart; jj<jend+1; jj++) {
1191: if (jj) {
1192: /* ii and kk must be zero */
1193: cols[cnt++] = slot + gnx*jj;
1194: } else {
1195: /* ii and jj must be zero */
1196: for (kk=kstart; kk<kend+1; kk++) {
1197: cols[cnt++] = slot + gnx*gny*kk;
1198: }
1199: }
1200: }
1201: }
1202: }
1203: }
1204: MatPreallocateSetLocal(ltogb,1,&slot,cnt,cols,dnz,onz);
1205: }
1206: }
1207: }
1208: MatSeqBAIJSetPreallocation(J,nc,0,dnz);
1209: MatMPIBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1210: MatPreallocateFinalize(dnz,onz);
1212: MatSetLocalToGlobalMapping(J,ltog);
1213: MatSetLocalToGlobalMappingBlock(J,ltogb);
1215: /*
1216: For each node in the grid: we get the neighbors in the local (on processor ordering
1217: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1218: PETSc ordering.
1219: */
1221: for (i=xs; i<xs+nx; i++) {
1222: istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
1223: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
1224: for (j=ys; j<ys+ny; j++) {
1225: jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
1226: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
1227: for (k=zs; k<zs+nz; k++) {
1228: kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
1229: kend = DAZPeriodic(wrap) ? s : (PetscMin(s,p-k-1));
1230:
1231: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1232:
1233: cnt = 0;
1234: if (st == DA_STENCIL_BOX) { /* if using BOX stencil */
1235: for (ii=istart; ii<iend+1; ii++) {
1236: for (jj=jstart; jj<jend+1; jj++) {
1237: for (kk=kstart; kk<kend+1; kk++) {
1238: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1239: }
1240: }
1241: }
1242: } else { /* Star stencil */
1243: cnt = 0;
1244: for (ii=istart; ii<iend+1; ii++) {
1245: if (ii) {
1246: /* jj and kk must be zero */
1247: cols[cnt++] = slot + ii;
1248: } else {
1249: for (jj=jstart; jj<jend+1; jj++) {
1250: if (jj) {
1251: /* ii and kk must be zero */
1252: cols[cnt++] = slot + gnx*jj;
1253: } else {
1254: /* ii and jj must be zero */
1255: for (kk=kstart; kk<kend+1; kk++) {
1256: cols[cnt++] = slot + gnx*gny*kk;
1257: }
1258: }
1259: }
1260: }
1261: }
1262: }
1263: MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1264: }
1265: }
1266: }
1267: PetscFree(values);
1268: PetscFree(cols);
1269: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1270: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1271: return(0);
1272: }
1275: PetscErrorCode DAGetMatrix2d_MPISBAIJ(DA da,Mat J)
1276: {
1277: PetscErrorCode ierr;
1278: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1279: PetscInt m,n,dim,s,*cols,nc,col,cnt,*dnz,*onz;
1280: PetscInt iend,jend,ii,jj;
1281: MPI_Comm comm;
1282: PetscScalar *values;
1283: DAPeriodicType wrap;
1284: DAStencilType st;
1285: ISLocalToGlobalMapping ltog,ltogb;
1288: /*
1289: nc - number of components per grid point
1290: col - number of colors needed in one direction for single component problem
1291: */
1292: DAGetInfo(da,&dim,&m,&n,0,0,0,0,&nc,&s,&wrap,&st);
1293: if (wrap != DA_NONPERIODIC) SETERRQ(PETSC_ERR_SUP,"Currently no support for periodic");
1294: col = 2*s + 1;
1296: DAGetCorners(da,&xs,&ys,0,&nx,&ny,0);
1297: DAGetGhostCorners(da,&gxs,&gys,0,&gnx,&gny,0);
1298: PetscObjectGetComm((PetscObject)da,&comm);
1300: PetscMalloc(col*col*nc*nc*sizeof(PetscScalar),&values);
1301: PetscMemzero(values,col*col*nc*nc*sizeof(PetscScalar));
1302: PetscMalloc(col*col*nc*nc*sizeof(PetscInt),&cols);
1304: DAGetISLocalToGlobalMapping(da,<og);
1305: DAGetISLocalToGlobalMappingBlck(da,<ogb);
1307: /* determine the matrix preallocation information */
1308: MatPreallocateSymmetricInitialize(comm,nx*ny,nx*ny,dnz,onz);
1309: for (i=xs; i<xs+nx; i++) {
1310: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
1311: for (j=ys; j<ys+ny; j++) {
1312: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
1313: slot = i - gxs + gnx*(j - gys);
1315: /* Find block columns in block row */
1316: cnt = 0;
1317: if (st == DA_STENCIL_BOX) { /* if using BOX stencil */
1318: for (ii=0; ii<iend+1; ii++) {
1319: for (jj=0; jj<jend+1; jj++) {
1320: cols[cnt++] = slot + ii + gnx*jj;
1321: }
1322: }
1323: } else { /* Star stencil */
1324: cnt = 0;
1325: for (ii=0; ii<iend+1; ii++) {
1326: if (ii) { /* jj must be zero */
1327: cols[cnt++] = slot + ii;
1328: } else {
1329: for (jj=0; jj<jend+1; jj++) { /* ii must be zero */
1330: cols[cnt++] = slot + gnx*jj;
1331: }
1332: }
1333: }
1334: }
1335: MatPreallocateSymmetricSetLocal(ltogb,1,&slot,cnt,cols,dnz,onz);
1336: }
1337: }
1338: MatSeqSBAIJSetPreallocation(J,nc,0,dnz);
1339: MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1340: MatPreallocateFinalize(dnz,onz);
1342: MatSetLocalToGlobalMapping(J,ltog);
1343: MatSetLocalToGlobalMappingBlock(J,ltogb);
1345: /*
1346: For each node in the grid: we get the neighbors in the local (on processor ordering
1347: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1348: PETSc ordering.
1349: */
1350: for (i=xs; i<xs+nx; i++) {
1351: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
1352: for (j=ys; j<ys+ny; j++) {
1353: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
1354: slot = i - gxs + gnx*(j - gys);
1355: cnt = 0;
1356: if (st == DA_STENCIL_BOX) { /* if using BOX stencil */
1357: for (ii=0; ii<iend+1; ii++) {
1358: for (jj=0; jj<jend+1; jj++) {
1359: cols[cnt++] = slot + ii + gnx*jj;
1360: }
1361: }
1362: } else { /* Star stencil */
1363: cnt = 0;
1364: for (ii=0; ii<iend+1; ii++) {
1365: if (ii) { /* jj must be zero */
1366: cols[cnt++] = slot + ii;
1367: } else {
1368: for (jj=0; jj<jend+1; jj++) {/* ii must be zero */
1369: cols[cnt++] = slot + gnx*jj;
1370: }
1371: }
1372: }
1373: }
1374: MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1375: }
1376: }
1377: PetscFree(values);
1378: PetscFree(cols);
1379: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1380: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1381: return(0);
1382: }
1386: PetscErrorCode DAGetMatrix3d_MPISBAIJ(DA da,Mat J)
1387: {
1388: PetscErrorCode ierr;
1389: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1390: PetscInt m,n,dim,s,*cols,k,nc,col,cnt,p,*dnz,*onz;
1391: PetscInt iend,jend,kend,zs,nz,gzs,gnz,ii,jj,kk;
1392: MPI_Comm comm;
1393: PetscScalar *values;
1394: DAPeriodicType wrap;
1395: DAStencilType st;
1396: ISLocalToGlobalMapping ltog,ltogb;
1399: /*
1400: nc - number of components per grid point
1401: col - number of colors needed in one direction for single component problem
1402: */
1403: DAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&wrap,&st);
1404: if (wrap != DA_NONPERIODIC) SETERRQ(PETSC_ERR_SUP,"Currently no support for periodic");
1405: col = 2*s + 1;
1407: DAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1408: DAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1409: PetscObjectGetComm((PetscObject)da,&comm);
1411: /* create the matrix */
1412: PetscMalloc(col*col*col*nc*nc*sizeof(PetscScalar),&values);
1413: PetscMemzero(values,col*col*col*nc*nc*sizeof(PetscScalar));
1414: PetscMalloc(col*col*col*sizeof(PetscInt),&cols);
1416: DAGetISLocalToGlobalMapping(da,<og);
1417: DAGetISLocalToGlobalMappingBlck(da,<ogb);
1419: /* determine the matrix preallocation information */
1420: MatPreallocateSymmetricInitialize(comm,nx*ny*nz,nx*ny*nz,dnz,onz);
1421: for (i=xs; i<xs+nx; i++) {
1422: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
1423: for (j=ys; j<ys+ny; j++) {
1424: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
1425: for (k=zs; k<zs+nz; k++) {
1426: kend = DAZPeriodic(wrap) ? s : (PetscMin(s,p-k-1));
1428: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1430: /* Find block columns in block row */
1431: cnt = 0;
1432: if (st == DA_STENCIL_BOX) { /* if using BOX stencil */
1433: for (ii=0; ii<iend+1; ii++) {
1434: for (jj=0; jj<jend+1; jj++) {
1435: for (kk=0; kk<kend+1; kk++) {
1436: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1437: }
1438: }
1439: }
1440: } else { /* Star stencil */
1441: cnt = 0;
1442: for (ii=0; ii<iend+1; ii++) {
1443: if (ii) {
1444: /* jj and kk must be zero */
1445: cols[cnt++] = slot + ii;
1446: } else {
1447: for (jj=0; jj<jend+1; jj++) {
1448: if (jj) {
1449: /* ii and kk must be zero */
1450: cols[cnt++] = slot + gnx*jj;
1451: } else {
1452: /* ii and jj must be zero */
1453: for (kk=0; kk<kend+1; kk++) {
1454: cols[cnt++] = slot + gnx*gny*kk;
1455: }
1456: }
1457: }
1458: }
1459: }
1460: }
1461: MatPreallocateSymmetricSetLocal(ltogb,1,&slot,cnt,cols,dnz,onz);
1462: }
1463: }
1464: }
1465: MatSeqSBAIJSetPreallocation(J,nc,0,dnz);
1466: MatMPISBAIJSetPreallocation(J,nc,0,dnz,0,onz);
1467: MatPreallocateFinalize(dnz,onz);
1469: MatSetLocalToGlobalMapping(J,ltog);
1470: MatSetLocalToGlobalMappingBlock(J,ltogb);
1472: /*
1473: For each node in the grid: we get the neighbors in the local (on processor ordering
1474: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1475: PETSc ordering.
1476: */
1478: for (i=xs; i<xs+nx; i++) {
1479: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
1480: for (j=ys; j<ys+ny; j++) {
1481: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
1482: for (k=zs; k<zs+nz; k++) {
1483: kend = DAZPeriodic(wrap) ? s : (PetscMin(s,p-k-1));
1484:
1485: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1486:
1487: cnt = 0;
1488: if (st == DA_STENCIL_BOX) { /* if using BOX stencil */
1489: for (ii=0; ii<iend+1; ii++) {
1490: for (jj=0; jj<jend+1; jj++) {
1491: for (kk=0; kk<kend+1; kk++) {
1492: cols[cnt++] = slot + ii + gnx*jj + gnx*gny*kk;
1493: }
1494: }
1495: }
1496: } else { /* Star stencil */
1497: cnt = 0;
1498: for (ii=0; ii<iend+1; ii++) {
1499: if (ii) {
1500: /* jj and kk must be zero */
1501: cols[cnt++] = slot + ii;
1502: } else {
1503: for (jj=0; jj<jend+1; jj++) {
1504: if (jj) {
1505: /* ii and kk must be zero */
1506: cols[cnt++] = slot + gnx*jj;
1507: } else {
1508: /* ii and jj must be zero */
1509: for (kk=0; kk<kend+1; kk++) {
1510: cols[cnt++] = slot + gnx*gny*kk;
1511: }
1512: }
1513: }
1514: }
1515: }
1516: }
1517: MatSetValuesBlockedLocal(J,1,&slot,cnt,cols,values,INSERT_VALUES);
1518: }
1519: }
1520: }
1521: PetscFree(values);
1522: PetscFree(cols);
1523: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1524: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1525: return(0);
1526: }
1528: /* ---------------------------------------------------------------------------------*/
1532: PetscErrorCode DAGetMatrix3d_MPIAIJ_Fill(DA da,Mat J)
1533: {
1534: PetscErrorCode ierr;
1535: PetscInt xs,ys,nx,ny,i,j,slot,gxs,gys,gnx,gny;
1536: PetscInt m,n,dim,s,*cols,k,nc,*rows,col,cnt,l,p,*dnz,*onz;
1537: PetscInt istart,iend,jstart,jend,kstart,kend,zs,nz,gzs,gnz,ii,jj,kk;
1538: PetscInt ifill_col,*dfill = da->dfill,*ofill = da->ofill;
1539: MPI_Comm comm;
1540: PetscScalar *values;
1541: DAPeriodicType wrap;
1542: ISLocalToGlobalMapping ltog;
1543: DAStencilType st;
1546: /*
1547: nc - number of components per grid point
1548: col - number of colors needed in one direction for single component problem
1549:
1550: */
1551: DAGetInfo(da,&dim,&m,&n,&p,0,0,0,&nc,&s,&wrap,&st);
1552: col = 2*s + 1;
1553: if (DAXPeriodic(wrap) && (m % col)){
1554: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in X is divisible\n\
1555: by 2*stencil_width + 1\n");
1556: }
1557: if (DAYPeriodic(wrap) && (n % col)){
1558: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Y is divisible\n\
1559: by 2*stencil_width + 1\n");
1560: }
1561: if (DAZPeriodic(wrap) && (p % col)){
1562: SETERRQ(PETSC_ERR_SUP,"For coloring efficiency ensure number of grid points in Z is divisible\n\
1563: by 2*stencil_width + 1\n");
1564: }
1566: DAGetCorners(da,&xs,&ys,&zs,&nx,&ny,&nz);
1567: DAGetGhostCorners(da,&gxs,&gys,&gzs,&gnx,&gny,&gnz);
1568: PetscObjectGetComm((PetscObject)da,&comm);
1570: PetscMalloc(col*col*col*nc*nc*nc*sizeof(PetscScalar),&values);
1571: PetscMemzero(values,col*col*col*nc*nc*nc*sizeof(PetscScalar));
1572: PetscMalloc(nc*sizeof(PetscInt),&rows);
1573: PetscMalloc(col*col*col*nc*sizeof(PetscInt),&cols);
1574: DAGetISLocalToGlobalMapping(da,<og);
1576: /* determine the matrix preallocation information */
1577: MatPreallocateInitialize(comm,nc*nx*ny*nz,nc*nx*ny*nz,dnz,onz);
1580: for (i=xs; i<xs+nx; i++) {
1581: istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
1582: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
1583: for (j=ys; j<ys+ny; j++) {
1584: jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
1585: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
1586: for (k=zs; k<zs+nz; k++) {
1587: kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
1588: kend = DAZPeriodic(wrap) ? s : (PetscMin(s,p-k-1));
1589:
1590: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1591:
1592: for (l=0; l<nc; l++) {
1593: cnt = 0;
1594: for (ii=istart; ii<iend+1; ii++) {
1595: for (jj=jstart; jj<jend+1; jj++) {
1596: for (kk=kstart; kk<kend+1; kk++) {
1597: if (ii || jj || kk) {
1598: if ((st == DA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1599: for (ifill_col=ofill[l]; ifill_col<ofill[l+1]; ifill_col++)
1600: cols[cnt++] = ofill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1601: }
1602: } else {
1603: if (dfill) {
1604: for (ifill_col=dfill[l]; ifill_col<dfill[l+1]; ifill_col++)
1605: cols[cnt++] = dfill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1606: } else {
1607: for (ifill_col=0; ifill_col<nc; ifill_col++)
1608: cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1609: }
1610: }
1611: }
1612: }
1613: }
1614: rows[0] = l + nc*(slot);
1615: MatPreallocateSetLocal(ltog,1,rows,cnt,cols,dnz,onz);
1616: }
1617: }
1618: }
1619: }
1620: MatSeqAIJSetPreallocation(J,0,dnz);
1621: MatMPIAIJSetPreallocation(J,0,dnz,0,onz);
1622: MatPreallocateFinalize(dnz,onz);
1623: MatSetLocalToGlobalMapping(J,ltog);
1625: /*
1626: For each node in the grid: we get the neighbors in the local (on processor ordering
1627: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1628: PETSc ordering.
1629: */
1630: for (i=xs; i<xs+nx; i++) {
1631: istart = DAXPeriodic(wrap) ? -s : (PetscMax(-s,-i));
1632: iend = DAXPeriodic(wrap) ? s : (PetscMin(s,m-i-1));
1633: for (j=ys; j<ys+ny; j++) {
1634: jstart = DAYPeriodic(wrap) ? -s : (PetscMax(-s,-j));
1635: jend = DAYPeriodic(wrap) ? s : (PetscMin(s,n-j-1));
1636: for (k=zs; k<zs+nz; k++) {
1637: kstart = DAZPeriodic(wrap) ? -s : (PetscMax(-s,-k));
1638: kend = DAZPeriodic(wrap) ? s : (PetscMin(s,p-k-1));
1639:
1640: slot = i - gxs + gnx*(j - gys) + gnx*gny*(k - gzs);
1641:
1642: for (l=0; l<nc; l++) {
1643: cnt = 0;
1644: for (ii=istart; ii<iend+1; ii++) {
1645: for (jj=jstart; jj<jend+1; jj++) {
1646: for (kk=kstart; kk<kend+1; kk++) {
1647: if (ii || jj || kk) {
1648: if ((st == DA_STENCIL_BOX) || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) {/* entries on star*/
1649: for (ifill_col=ofill[l]; ifill_col<ofill[l+1]; ifill_col++)
1650: cols[cnt++] = ofill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1651: }
1652: } else {
1653: if (dfill) {
1654: for (ifill_col=dfill[l]; ifill_col<dfill[l+1]; ifill_col++)
1655: cols[cnt++] = dfill[ifill_col] + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1656: } else {
1657: for (ifill_col=0; ifill_col<nc; ifill_col++)
1658: cols[cnt++] = ifill_col + nc*(slot + ii + gnx*jj + gnx*gny*kk);
1659: }
1660: }
1661: }
1662: }
1663: }
1664: rows[0] = l + nc*(slot);
1665: MatSetValuesLocal(J,1,rows,cnt,cols,values,INSERT_VALUES);
1666: }
1667: }
1668: }
1669: }
1670: PetscFree(values);
1671: PetscFree(rows);
1672: PetscFree(cols);
1673: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
1674: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
1675: return(0);
1676: }