Actual source code: fdmpiaij.c
petsc-3.9.4 2018-09-11
1: #include <../src/mat/impls/sell/mpi/mpisell.h>
2: #include <../src/mat/impls/aij/mpi/mpiaij.h>
3: #include <../src/mat/impls/baij/mpi/mpibaij.h>
4: #include <petsc/private/isimpl.h>
6: PetscErrorCode MatFDColoringApply_BAIJ(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
7: {
8: PetscErrorCode (*f)(void*,Vec,Vec,void*)=(PetscErrorCode (*)(void*,Vec,Vec,void*))coloring->f;
9: PetscErrorCode ierr;
10: PetscInt k,cstart,cend,l,row,col,nz,spidx,i,j;
11: PetscScalar dx=0.0,*w3_array,*dy_i,*dy=coloring->dy;
12: PetscScalar *vscale_array;
13: const PetscScalar *xx;
14: PetscReal epsilon=coloring->error_rel,umin=coloring->umin,unorm;
15: Vec w1=coloring->w1,w2=coloring->w2,w3,vscale=coloring->vscale;
16: void *fctx=coloring->fctx;
17: PetscInt ctype=coloring->ctype,nxloc,nrows_k;
18: PetscScalar *valaddr;
19: MatEntry *Jentry=coloring->matentry;
20: MatEntry2 *Jentry2=coloring->matentry2;
21: const PetscInt ncolors=coloring->ncolors,*ncolumns=coloring->ncolumns,*nrows=coloring->nrows;
22: PetscInt bs=J->rmap->bs;
25: /* (1) Set w1 = F(x1) */
26: if (!coloring->fset) {
27: PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);
28: (*f)(sctx,x1,w1,fctx);
29: PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);
30: } else {
31: coloring->fset = PETSC_FALSE;
32: }
34: /* (2) Compute vscale = 1./dx - the local scale factors, including ghost points */
35: VecGetLocalSize(x1,&nxloc);
36: if (coloring->htype[0] == 'w') {
37: /* vscale = dx is a constant scalar */
38: VecNorm(x1,NORM_2,&unorm);
39: dx = 1.0/(PetscSqrtReal(1.0 + unorm)*epsilon);
40: } else {
41: VecGetArrayRead(x1,&xx);
42: VecGetArray(vscale,&vscale_array);
43: for (col=0; col<nxloc; col++) {
44: dx = xx[col];
45: if (PetscAbsScalar(dx) < umin) {
46: if (PetscRealPart(dx) >= 0.0) dx = umin;
47: else if (PetscRealPart(dx) < 0.0 ) dx = -umin;
48: }
49: dx *= epsilon;
50: vscale_array[col] = 1.0/dx;
51: }
52: VecRestoreArrayRead(x1,&xx);
53: VecRestoreArray(vscale,&vscale_array);
54: }
55: if (ctype == IS_COLORING_GLOBAL && coloring->htype[0] == 'd') {
56: VecGhostUpdateBegin(vscale,INSERT_VALUES,SCATTER_FORWARD);
57: VecGhostUpdateEnd(vscale,INSERT_VALUES,SCATTER_FORWARD);
58: }
60: /* (3) Loop over each color */
61: if (!coloring->w3) {
62: VecDuplicate(x1,&coloring->w3);
63: PetscLogObjectParent((PetscObject)coloring,(PetscObject)coloring->w3);
64: }
65: w3 = coloring->w3;
67: VecGetOwnershipRange(x1,&cstart,&cend); /* used by ghosted vscale */
68: if (vscale) {
69: VecGetArray(vscale,&vscale_array);
70: }
71: nz = 0;
72: for (k=0; k<ncolors; k++) {
73: coloring->currentcolor = k;
75: /*
76: (3-1) Loop over each column associated with color
77: adding the perturbation to the vector w3 = x1 + dx.
78: */
79: VecCopy(x1,w3);
80: dy_i = dy;
81: for (i=0; i<bs; i++) { /* Loop over a block of columns */
82: VecGetArray(w3,&w3_array);
83: if (ctype == IS_COLORING_GLOBAL) w3_array -= cstart; /* shift pointer so global index can be used */
84: if (coloring->htype[0] == 'w') {
85: for (l=0; l<ncolumns[k]; l++) {
86: col = i + bs*coloring->columns[k][l]; /* local column (in global index!) of the matrix we are probing for */
87: w3_array[col] += 1.0/dx;
88: if (i) w3_array[col-1] -= 1.0/dx; /* resume original w3[col-1] */
89: }
90: } else { /* htype == 'ds' */
91: vscale_array -= cstart; /* shift pointer so global index can be used */
92: for (l=0; l<ncolumns[k]; l++) {
93: col = i + bs*coloring->columns[k][l]; /* local column (in global index!) of the matrix we are probing for */
94: w3_array[col] += 1.0/vscale_array[col];
95: if (i) w3_array[col-1] -= 1.0/vscale_array[col-1]; /* resume original w3[col-1] */
96: }
97: vscale_array += cstart;
98: }
99: if (ctype == IS_COLORING_GLOBAL) w3_array += cstart;
100: VecRestoreArray(w3,&w3_array);
102: /*
103: (3-2) Evaluate function at w3 = x1 + dx (here dx is a vector of perturbations)
104: w2 = F(x1 + dx) - F(x1)
105: */
106: PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);
107: VecPlaceArray(w2,dy_i); /* place w2 to the array dy_i */
108: (*f)(sctx,w3,w2,fctx);
109: PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);
110: VecAXPY(w2,-1.0,w1);
111: VecResetArray(w2);
112: dy_i += nxloc; /* points to dy+i*nxloc */
113: }
115: /*
116: (3-3) Loop over rows of vector, putting results into Jacobian matrix
117: */
118: nrows_k = nrows[k];
119: if (coloring->htype[0] == 'w') {
120: for (l=0; l<nrows_k; l++) {
121: row = bs*Jentry2[nz].row; /* local row index */
122: valaddr = Jentry2[nz++].valaddr;
123: spidx = 0;
124: dy_i = dy;
125: for (i=0; i<bs; i++) { /* column of the block */
126: for (j=0; j<bs; j++) { /* row of the block */
127: valaddr[spidx++] = dy_i[row+j]*dx;
128: }
129: dy_i += nxloc; /* points to dy+i*nxloc */
130: }
131: }
132: } else { /* htype == 'ds' */
133: for (l=0; l<nrows_k; l++) {
134: row = bs*Jentry[nz].row; /* local row index */
135: col = bs*Jentry[nz].col; /* local column index */
136: valaddr = Jentry[nz++].valaddr;
137: spidx = 0;
138: dy_i = dy;
139: for (i=0; i<bs; i++) { /* column of the block */
140: for (j=0; j<bs; j++) { /* row of the block */
141: valaddr[spidx++] = dy_i[row+j]*vscale_array[col+i];
142: }
143: dy_i += nxloc; /* points to dy+i*nxloc */
144: }
145: }
146: }
147: }
148: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
149: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
150: if (vscale) {
151: VecRestoreArray(vscale,&vscale_array);
152: }
154: coloring->currentcolor = -1;
155: return(0);
156: }
158: /* this is declared PETSC_EXTERN because it is used by MatFDColoringUseDM() which is in the DM library */
159: PetscErrorCode MatFDColoringApply_AIJ(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
160: {
161: PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void*))coloring->f;
162: PetscErrorCode ierr;
163: PetscInt k,cstart,cend,l,row,col,nz;
164: PetscScalar dx=0.0,*y,*w3_array;
165: const PetscScalar *xx;
166: PetscScalar *vscale_array;
167: PetscReal epsilon=coloring->error_rel,umin=coloring->umin,unorm;
168: Vec w1=coloring->w1,w2=coloring->w2,w3,vscale=coloring->vscale;
169: void *fctx=coloring->fctx;
170: ISColoringType ctype=coloring->ctype;
171: PetscInt nxloc,nrows_k;
172: MatEntry *Jentry=coloring->matentry;
173: MatEntry2 *Jentry2=coloring->matentry2;
174: const PetscInt ncolors=coloring->ncolors,*ncolumns=coloring->ncolumns,*nrows=coloring->nrows;
177: if ((ctype == IS_COLORING_LOCAL) && (J->ops->fdcoloringapply == MatFDColoringApply_AIJ)) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_SUP,"Must call MatColoringUseDM() with IS_COLORING_LOCAL");
178: /* (1) Set w1 = F(x1) */
179: if (!coloring->fset) {
180: PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);
181: (*f)(sctx,x1,w1,fctx);
182: PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);
183: } else {
184: coloring->fset = PETSC_FALSE;
185: }
187: /* (2) Compute vscale = 1./dx - the local scale factors, including ghost points */
188: if (coloring->htype[0] == 'w') {
189: /* vscale = 1./dx is a constant scalar */
190: VecNorm(x1,NORM_2,&unorm);
191: dx = 1.0/(PetscSqrtReal(1.0 + unorm)*epsilon);
192: } else {
193: VecGetLocalSize(x1,&nxloc);
194: VecGetArrayRead(x1,&xx);
195: VecGetArray(vscale,&vscale_array);
196: for (col=0; col<nxloc; col++) {
197: dx = xx[col];
198: if (PetscAbsScalar(dx) < umin) {
199: if (PetscRealPart(dx) >= 0.0) dx = umin;
200: else if (PetscRealPart(dx) < 0.0 ) dx = -umin;
201: }
202: dx *= epsilon;
203: vscale_array[col] = 1.0/dx;
204: }
205: VecRestoreArrayRead(x1,&xx);
206: VecRestoreArray(vscale,&vscale_array);
207: }
208: if (ctype == IS_COLORING_GLOBAL && coloring->htype[0] == 'd') {
209: VecGhostUpdateBegin(vscale,INSERT_VALUES,SCATTER_FORWARD);
210: VecGhostUpdateEnd(vscale,INSERT_VALUES,SCATTER_FORWARD);
211: }
213: /* (3) Loop over each color */
214: if (!coloring->w3) {
215: VecDuplicate(x1,&coloring->w3);
216: PetscLogObjectParent((PetscObject)coloring,(PetscObject)coloring->w3);
217: }
218: w3 = coloring->w3;
220: VecGetOwnershipRange(x1,&cstart,&cend); /* used by ghosted vscale */
221: if (vscale) {
222: VecGetArray(vscale,&vscale_array);
223: }
224: nz = 0;
226: if (coloring->bcols > 1) { /* use blocked insertion of Jentry */
227: PetscInt i,m=J->rmap->n,nbcols,bcols=coloring->bcols;
228: PetscScalar *dy=coloring->dy,*dy_k;
230: nbcols = 0;
231: for (k=0; k<ncolors; k+=bcols) {
233: /*
234: (3-1) Loop over each column associated with color
235: adding the perturbation to the vector w3 = x1 + dx.
236: */
238: dy_k = dy;
239: if (k + bcols > ncolors) bcols = ncolors - k;
240: for (i=0; i<bcols; i++) {
241: coloring->currentcolor = k+i;
243: VecCopy(x1,w3);
244: VecGetArray(w3,&w3_array);
245: if (ctype == IS_COLORING_GLOBAL) w3_array -= cstart; /* shift pointer so global index can be used */
246: if (coloring->htype[0] == 'w') {
247: for (l=0; l<ncolumns[k+i]; l++) {
248: col = coloring->columns[k+i][l]; /* local column (in global index!) of the matrix we are probing for */
249: w3_array[col] += 1.0/dx;
250: }
251: } else { /* htype == 'ds' */
252: vscale_array -= cstart; /* shift pointer so global index can be used */
253: for (l=0; l<ncolumns[k+i]; l++) {
254: col = coloring->columns[k+i][l]; /* local column (in global index!) of the matrix we are probing for */
255: w3_array[col] += 1.0/vscale_array[col];
256: }
257: vscale_array += cstart;
258: }
259: if (ctype == IS_COLORING_GLOBAL) w3_array += cstart;
260: VecRestoreArray(w3,&w3_array);
262: /*
263: (3-2) Evaluate function at w3 = x1 + dx (here dx is a vector of perturbations)
264: w2 = F(x1 + dx) - F(x1)
265: */
266: PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);
267: VecPlaceArray(w2,dy_k); /* place w2 to the array dy_i */
268: (*f)(sctx,w3,w2,fctx);
269: PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);
270: VecAXPY(w2,-1.0,w1);
271: VecResetArray(w2);
272: dy_k += m; /* points to dy+i*nxloc */
273: }
275: /*
276: (3-3) Loop over block rows of vector, putting results into Jacobian matrix
277: */
278: nrows_k = nrows[nbcols++];
279: VecGetArray(w2,&y);
281: if (coloring->htype[0] == 'w') {
282: for (l=0; l<nrows_k; l++) {
283: row = Jentry2[nz].row; /* local row index */
284: *(Jentry2[nz++].valaddr) = dy[row]*dx;
285: }
286: } else { /* htype == 'ds' */
287: for (l=0; l<nrows_k; l++) {
288: row = Jentry[nz].row; /* local row index */
289: *(Jentry[nz].valaddr) = dy[row]*vscale_array[Jentry[nz].col];
290: nz++;
291: }
292: }
293: VecRestoreArray(w2,&y);
294: }
295: } else { /* bcols == 1 */
296: for (k=0; k<ncolors; k++) {
297: coloring->currentcolor = k;
299: /*
300: (3-1) Loop over each column associated with color
301: adding the perturbation to the vector w3 = x1 + dx.
302: */
303: VecCopy(x1,w3);
304: VecGetArray(w3,&w3_array);
305: if (ctype == IS_COLORING_GLOBAL) w3_array -= cstart; /* shift pointer so global index can be used */
306: if (coloring->htype[0] == 'w') {
307: for (l=0; l<ncolumns[k]; l++) {
308: col = coloring->columns[k][l]; /* local column (in global index!) of the matrix we are probing for */
309: w3_array[col] += 1.0/dx;
310: }
311: } else { /* htype == 'ds' */
312: vscale_array -= cstart; /* shift pointer so global index can be used */
313: for (l=0; l<ncolumns[k]; l++) {
314: col = coloring->columns[k][l]; /* local column (in global index!) of the matrix we are probing for */
315: w3_array[col] += 1.0/vscale_array[col];
316: }
317: vscale_array += cstart;
318: }
319: if (ctype == IS_COLORING_GLOBAL) w3_array += cstart;
320: VecRestoreArray(w3,&w3_array);
322: /*
323: (3-2) Evaluate function at w3 = x1 + dx (here dx is a vector of perturbations)
324: w2 = F(x1 + dx) - F(x1)
325: */
326: PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);
327: (*f)(sctx,w3,w2,fctx);
328: PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);
329: VecAXPY(w2,-1.0,w1);
331: /*
332: (3-3) Loop over rows of vector, putting results into Jacobian matrix
333: */
334: nrows_k = nrows[k];
335: VecGetArray(w2,&y);
336: if (coloring->htype[0] == 'w') {
337: for (l=0; l<nrows_k; l++) {
338: row = Jentry2[nz].row; /* local row index */
339: *(Jentry2[nz++].valaddr) = y[row]*dx;
340: }
341: } else { /* htype == 'ds' */
342: for (l=0; l<nrows_k; l++) {
343: row = Jentry[nz].row; /* local row index */
344: *(Jentry[nz].valaddr) = y[row]*vscale_array[Jentry[nz].col];
345: nz++;
346: }
347: }
348: VecRestoreArray(w2,&y);
349: }
350: }
352: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
353: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
354: if (vscale) {
355: VecRestoreArray(vscale,&vscale_array);
356: }
357: coloring->currentcolor = -1;
358: return(0);
359: }
361: PetscErrorCode MatFDColoringSetUp_MPIXAIJ(Mat mat,ISColoring iscoloring,MatFDColoring c)
362: {
363: PetscErrorCode ierr;
364: PetscMPIInt size,*ncolsonproc,*disp,nn;
365: PetscInt i,n,nrows,nrows_i,j,k,m,ncols,col,*rowhit,cstart,cend,colb;
366: const PetscInt *is,*A_ci,*A_cj,*B_ci,*B_cj,*row=NULL,*ltog=NULL;
367: PetscInt nis=iscoloring->n,nctot,*cols;
368: IS *isa;
369: ISLocalToGlobalMapping map=mat->cmap->mapping;
370: PetscInt ctype=c->ctype,*spidxA,*spidxB,nz,bs,bs2,spidx;
371: Mat A,B;
372: PetscScalar *A_val,*B_val,**valaddrhit;
373: MatEntry *Jentry;
374: MatEntry2 *Jentry2;
375: PetscBool isBAIJ,isSELL;
376: PetscInt bcols=c->bcols;
377: #if defined(PETSC_USE_CTABLE)
378: PetscTable colmap=NULL;
379: #else
380: PetscInt *colmap=NULL; /* local col number of off-diag col */
381: #endif
384: if (ctype == IS_COLORING_LOCAL) {
385: if (!map) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_INCOMP,"When using ghosted differencing matrix must have local to global mapping provided with MatSetLocalToGlobalMapping");
386: ISLocalToGlobalMappingGetIndices(map,<og);
387: }
389: MatGetBlockSize(mat,&bs);
390: PetscObjectTypeCompare((PetscObject)mat,MATMPIBAIJ,&isBAIJ);
391: PetscObjectTypeCompare((PetscObject)mat,MATMPISELL,&isSELL);
392: if (isBAIJ) {
393: Mat_MPIBAIJ *baij=(Mat_MPIBAIJ*)mat->data;
394: Mat_SeqBAIJ *spA,*spB;
395: A = baij->A; spA = (Mat_SeqBAIJ*)A->data; A_val = spA->a;
396: B = baij->B; spB = (Mat_SeqBAIJ*)B->data; B_val = spB->a;
397: nz = spA->nz + spB->nz; /* total nonzero entries of mat */
398: if (!baij->colmap) {
399: MatCreateColmap_MPIBAIJ_Private(mat);
400: }
401: colmap = baij->colmap;
402: MatGetColumnIJ_SeqBAIJ_Color(A,0,PETSC_FALSE,PETSC_FALSE,&ncols,&A_ci,&A_cj,&spidxA,NULL);
403: MatGetColumnIJ_SeqBAIJ_Color(B,0,PETSC_FALSE,PETSC_FALSE,&ncols,&B_ci,&B_cj,&spidxB,NULL);
405: if (ctype == IS_COLORING_GLOBAL && c->htype[0] == 'd') { /* create vscale for storing dx */
406: PetscInt *garray;
407: PetscMalloc1(B->cmap->n,&garray);
408: for (i=0; i<baij->B->cmap->n/bs; i++) {
409: for (j=0; j<bs; j++) {
410: garray[i*bs+j] = bs*baij->garray[i]+j;
411: }
412: }
413: VecCreateGhost(PetscObjectComm((PetscObject)mat),mat->cmap->n,PETSC_DETERMINE,B->cmap->n,garray,&c->vscale);
414: PetscFree(garray);
415: }
416: } else if (isSELL) {
417: Mat_MPISELL *sell=(Mat_MPISELL*)mat->data;
418: Mat_SeqSELL *spA,*spB;
419: A = sell->A; spA = (Mat_SeqSELL*)A->data; A_val = spA->val;
420: B = sell->B; spB = (Mat_SeqSELL*)B->data; B_val = spB->val;
421: nz = spA->nz + spB->nz; /* total nonzero entries of mat */
422: if (!sell->colmap) {
423: /* Allow access to data structures of local part of matrix
424: - creates aij->colmap which maps global column number to local number in part B */
425: MatCreateColmap_MPISELL_Private(mat);
426: }
427: colmap = sell->colmap;
428: MatGetColumnIJ_SeqSELL_Color(A,0,PETSC_FALSE,PETSC_FALSE,&ncols,&A_ci,&A_cj,&spidxA,NULL);
429: MatGetColumnIJ_SeqSELL_Color(B,0,PETSC_FALSE,PETSC_FALSE,&ncols,&B_ci,&B_cj,&spidxB,NULL);
431: bs = 1; /* only bs=1 is supported for non MPIBAIJ matrix */
433: if (ctype == IS_COLORING_GLOBAL && c->htype[0] == 'd') { /* create vscale for storing dx */
434: VecCreateGhost(PetscObjectComm((PetscObject)mat),mat->cmap->n,PETSC_DETERMINE,B->cmap->n,sell->garray,&c->vscale);
435: }
436: } else {
437: Mat_MPIAIJ *aij=(Mat_MPIAIJ*)mat->data;
438: Mat_SeqAIJ *spA,*spB;
439: A = aij->A; spA = (Mat_SeqAIJ*)A->data; A_val = spA->a;
440: B = aij->B; spB = (Mat_SeqAIJ*)B->data; B_val = spB->a;
441: nz = spA->nz + spB->nz; /* total nonzero entries of mat */
442: if (!aij->colmap) {
443: /* Allow access to data structures of local part of matrix
444: - creates aij->colmap which maps global column number to local number in part B */
445: MatCreateColmap_MPIAIJ_Private(mat);
446: }
447: colmap = aij->colmap;
448: MatGetColumnIJ_SeqAIJ_Color(A,0,PETSC_FALSE,PETSC_FALSE,&ncols,&A_ci,&A_cj,&spidxA,NULL);
449: MatGetColumnIJ_SeqAIJ_Color(B,0,PETSC_FALSE,PETSC_FALSE,&ncols,&B_ci,&B_cj,&spidxB,NULL);
451: bs = 1; /* only bs=1 is supported for non MPIBAIJ matrix */
453: if (ctype == IS_COLORING_GLOBAL && c->htype[0] == 'd') { /* create vscale for storing dx */
454: VecCreateGhost(PetscObjectComm((PetscObject)mat),mat->cmap->n,PETSC_DETERMINE,B->cmap->n,aij->garray,&c->vscale);
455: }
456: }
458: m = mat->rmap->n/bs;
459: cstart = mat->cmap->rstart/bs;
460: cend = mat->cmap->rend/bs;
462: PetscMalloc1(nis,&c->ncolumns);
463: PetscMalloc1(nis,&c->columns);
464: PetscCalloc1(nis,&c->nrows);
465: PetscLogObjectMemory((PetscObject)c,3*nis*sizeof(PetscInt));
467: if (c->htype[0] == 'd') {
468: PetscMalloc1(nz,&Jentry);
469: PetscLogObjectMemory((PetscObject)c,nz*sizeof(MatEntry));
470: c->matentry = Jentry;
471: } else if (c->htype[0] == 'w') {
472: PetscMalloc1(nz,&Jentry2);
473: PetscLogObjectMemory((PetscObject)c,nz*sizeof(MatEntry2));
474: c->matentry2 = Jentry2;
475: } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"htype is not supported");
477: PetscMalloc2(m+1,&rowhit,m+1,&valaddrhit);
478: nz = 0;
479: ISColoringGetIS(iscoloring,PETSC_IGNORE,&isa);
480: for (i=0; i<nis; i++) { /* for each local color */
481: ISGetLocalSize(isa[i],&n);
482: ISGetIndices(isa[i],&is);
484: c->ncolumns[i] = n; /* local number of columns of this color on this process */
485: if (n) {
486: PetscMalloc1(n,&c->columns[i]);
487: PetscLogObjectMemory((PetscObject)c,n*sizeof(PetscInt));
488: PetscMemcpy(c->columns[i],is,n*sizeof(PetscInt));
489: } else {
490: c->columns[i] = 0;
491: }
493: if (ctype == IS_COLORING_GLOBAL) {
494: /* Determine nctot, the total (parallel) number of columns of this color */
495: MPI_Comm_size(PetscObjectComm((PetscObject)mat),&size);
496: PetscMalloc2(size,&ncolsonproc,size,&disp);
498: /* ncolsonproc[j]: local ncolumns on proc[j] of this color */
499: PetscMPIIntCast(n,&nn);
500: MPI_Allgather(&nn,1,MPI_INT,ncolsonproc,1,MPI_INT,PetscObjectComm((PetscObject)mat));
501: nctot = 0; for (j=0; j<size; j++) nctot += ncolsonproc[j];
502: if (!nctot) {
503: PetscInfo(mat,"Coloring of matrix has some unneeded colors with no corresponding rows\n");
504: }
506: disp[0] = 0;
507: for (j=1; j<size; j++) {
508: disp[j] = disp[j-1] + ncolsonproc[j-1];
509: }
511: /* Get cols, the complete list of columns for this color on each process */
512: PetscMalloc1(nctot+1,&cols);
513: MPI_Allgatherv((void*)is,n,MPIU_INT,cols,ncolsonproc,disp,MPIU_INT,PetscObjectComm((PetscObject)mat));
514: PetscFree2(ncolsonproc,disp);
515: } else if (ctype == IS_COLORING_LOCAL) {
516: /* Determine local number of columns of this color on this process, including ghost points */
517: nctot = n;
518: PetscMalloc1(nctot+1,&cols);
519: PetscMemcpy(cols,is,n*sizeof(PetscInt));
520: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not provided for this MatFDColoring type");
522: /* Mark all rows affect by these columns */
523: PetscMemzero(rowhit,m*sizeof(PetscInt));
524: bs2 = bs*bs;
525: nrows_i = 0;
526: for (j=0; j<nctot; j++) { /* loop over columns*/
527: if (ctype == IS_COLORING_LOCAL) {
528: col = ltog[cols[j]];
529: } else {
530: col = cols[j];
531: }
532: if (col >= cstart && col < cend) { /* column is in A, diagonal block of mat */
533: row = A_cj + A_ci[col-cstart];
534: nrows = A_ci[col-cstart+1] - A_ci[col-cstart];
535: nrows_i += nrows;
536: /* loop over columns of A marking them in rowhit */
537: for (k=0; k<nrows; k++) {
538: /* set valaddrhit for part A */
539: spidx = bs2*spidxA[A_ci[col-cstart] + k];
540: valaddrhit[*row] = &A_val[spidx];
541: rowhit[*row++] = col - cstart + 1; /* local column index */
542: }
543: } else { /* column is in B, off-diagonal block of mat */
544: #if defined(PETSC_USE_CTABLE)
545: PetscTableFind(colmap,col+1,&colb);
546: colb--;
547: #else
548: colb = colmap[col] - 1; /* local column index */
549: #endif
550: if (colb == -1) {
551: nrows = 0;
552: } else {
553: colb = colb/bs;
554: row = B_cj + B_ci[colb];
555: nrows = B_ci[colb+1] - B_ci[colb];
556: }
557: nrows_i += nrows;
558: /* loop over columns of B marking them in rowhit */
559: for (k=0; k<nrows; k++) {
560: /* set valaddrhit for part B */
561: spidx = bs2*spidxB[B_ci[colb] + k];
562: valaddrhit[*row] = &B_val[spidx];
563: rowhit[*row++] = colb + 1 + cend - cstart; /* local column index */
564: }
565: }
566: }
567: c->nrows[i] = nrows_i;
569: if (c->htype[0] == 'd') {
570: for (j=0; j<m; j++) {
571: if (rowhit[j]) {
572: Jentry[nz].row = j; /* local row index */
573: Jentry[nz].col = rowhit[j] - 1; /* local column index */
574: Jentry[nz].valaddr = valaddrhit[j]; /* address of mat value for this entry */
575: nz++;
576: }
577: }
578: } else { /* c->htype == 'wp' */
579: for (j=0; j<m; j++) {
580: if (rowhit[j]) {
581: Jentry2[nz].row = j; /* local row index */
582: Jentry2[nz].valaddr = valaddrhit[j]; /* address of mat value for this entry */
583: nz++;
584: }
585: }
586: }
587: PetscFree(cols);
588: }
590: if (bcols > 1) { /* reorder Jentry for faster MatFDColoringApply() */
591: MatFDColoringSetUpBlocked_AIJ_Private(mat,c,nz);
592: }
594: if (isBAIJ) {
595: MatRestoreColumnIJ_SeqBAIJ_Color(A,0,PETSC_FALSE,PETSC_FALSE,&ncols,&A_ci,&A_cj,&spidxA,NULL);
596: MatRestoreColumnIJ_SeqBAIJ_Color(B,0,PETSC_FALSE,PETSC_FALSE,&ncols,&B_ci,&B_cj,&spidxB,NULL);
597: PetscMalloc1(bs*mat->rmap->n,&c->dy);
598: } else if (isSELL) {
599: MatRestoreColumnIJ_SeqSELL_Color(A,0,PETSC_FALSE,PETSC_FALSE,&ncols,&A_ci,&A_cj,&spidxA,NULL);
600: MatRestoreColumnIJ_SeqSELL_Color(B,0,PETSC_FALSE,PETSC_FALSE,&ncols,&B_ci,&B_cj,&spidxB,NULL);
601: }else {
602: MatRestoreColumnIJ_SeqAIJ_Color(A,0,PETSC_FALSE,PETSC_FALSE,&ncols,&A_ci,&A_cj,&spidxA,NULL);
603: MatRestoreColumnIJ_SeqAIJ_Color(B,0,PETSC_FALSE,PETSC_FALSE,&ncols,&B_ci,&B_cj,&spidxB,NULL);
604: }
606: ISColoringRestoreIS(iscoloring,&isa);
607: PetscFree2(rowhit,valaddrhit);
609: if (ctype == IS_COLORING_LOCAL) {
610: ISLocalToGlobalMappingRestoreIndices(map,<og);
611: }
612: PetscInfo3(c,"ncolors %D, brows %D and bcols %D are used.\n",c->ncolors,c->brows,c->bcols);
613: return(0);
614: }
616: PetscErrorCode MatFDColoringCreate_MPIXAIJ(Mat mat,ISColoring iscoloring,MatFDColoring c)
617: {
619: PetscInt bs,nis=iscoloring->n,m=mat->rmap->n;
620: PetscBool isBAIJ,isSELL;
623: /* set default brows and bcols for speedup inserting the dense matrix into sparse Jacobian;
624: bcols is chosen s.t. dy-array takes 50% of memory space as mat */
625: MatGetBlockSize(mat,&bs);
626: PetscObjectTypeCompare((PetscObject)mat,MATMPIBAIJ,&isBAIJ);
627: PetscObjectTypeCompare((PetscObject)mat,MATMPISELL,&isSELL);
628: if (isBAIJ || m == 0) {
629: c->brows = m;
630: c->bcols = 1;
631: } else if (isSELL) {
632: /* bcols is chosen s.t. dy-array takes 50% of local memory space as mat */
633: Mat_MPISELL *sell=(Mat_MPISELL*)mat->data;
634: Mat_SeqSELL *spA,*spB;
635: Mat A,B;
636: PetscInt nz,brows,bcols;
637: PetscReal mem;
639: bs = 1; /* only bs=1 is supported for MPISELL matrix */
641: A = sell->A; spA = (Mat_SeqSELL*)A->data;
642: B = sell->B; spB = (Mat_SeqSELL*)B->data;
643: nz = spA->nz + spB->nz; /* total local nonzero entries of mat */
644: mem = nz*(sizeof(PetscScalar) + sizeof(PetscInt)) + 3*m*sizeof(PetscInt);
645: bcols = (PetscInt)(0.5*mem /(m*sizeof(PetscScalar)));
646: brows = 1000/bcols;
647: if (bcols > nis) bcols = nis;
648: if (brows == 0 || brows > m) brows = m;
649: c->brows = brows;
650: c->bcols = bcols;
651: } else { /* mpiaij matrix */
652: /* bcols is chosen s.t. dy-array takes 50% of local memory space as mat */
653: Mat_MPIAIJ *aij=(Mat_MPIAIJ*)mat->data;
654: Mat_SeqAIJ *spA,*spB;
655: Mat A,B;
656: PetscInt nz,brows,bcols;
657: PetscReal mem;
659: bs = 1; /* only bs=1 is supported for MPIAIJ matrix */
661: A = aij->A; spA = (Mat_SeqAIJ*)A->data;
662: B = aij->B; spB = (Mat_SeqAIJ*)B->data;
663: nz = spA->nz + spB->nz; /* total local nonzero entries of mat */
664: mem = nz*(sizeof(PetscScalar) + sizeof(PetscInt)) + 3*m*sizeof(PetscInt);
665: bcols = (PetscInt)(0.5*mem /(m*sizeof(PetscScalar)));
666: brows = 1000/bcols;
667: if (bcols > nis) bcols = nis;
668: if (brows == 0 || brows > m) brows = m;
669: c->brows = brows;
670: c->bcols = bcols;
671: }
673: c->M = mat->rmap->N/bs; /* set the global rows and columns and local rows */
674: c->N = mat->cmap->N/bs;
675: c->m = mat->rmap->n/bs;
676: c->rstart = mat->rmap->rstart/bs;
677: c->ncolors = nis;
678: return(0);
679: }