Actual source code: fdmpiaij.c
petsc-3.5.4 2015-05-23
2: #include <../src/mat/impls/aij/mpi/mpiaij.h>
3: #include <../src/mat/impls/baij/mpi/mpibaij.h>
7: PetscErrorCode MatFDColoringApply_BAIJ(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
8: {
9: PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void*))coloring->f;
11: PetscInt k,cstart,cend,l,row,col,nz,spidx,i,j;
12: PetscScalar dx=0.0,*xx,*w3_array,*dy_i,*dy=coloring->dy;
13: PetscScalar *vscale_array;
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: VecGetArray(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: VecRestoreArray(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: }
160: PetscErrorCode MatFDColoringApply_AIJ(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
161: {
162: PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void*))coloring->f;
164: PetscInt k,cstart,cend,l,row,col,nz;
165: PetscScalar dx=0.0,*y,*xx,*w3_array;
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: PetscInt ctype=coloring->ctype,nxloc,nrows_k;
171: MatEntry *Jentry=coloring->matentry;
172: MatEntry2 *Jentry2=coloring->matentry2;
173: const PetscInt ncolors=coloring->ncolors,*ncolumns=coloring->ncolumns,*nrows=coloring->nrows;
176: /* (1) Set w1 = F(x1) */
177: if (!coloring->fset) {
178: PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);
179: (*f)(sctx,x1,w1,fctx);
180: PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);
181: } else {
182: coloring->fset = PETSC_FALSE;
183: }
185: /* (2) Compute vscale = 1./dx - the local scale factors, including ghost points */
186: if (coloring->htype[0] == 'w') {
187: /* vscale = 1./dx is a constant scalar */
188: VecNorm(x1,NORM_2,&unorm);
189: dx = 1.0/(PetscSqrtReal(1.0 + unorm)*epsilon);
190: } else {
191: VecGetLocalSize(x1,&nxloc);
192: VecGetArray(x1,&xx);
193: VecGetArray(vscale,&vscale_array);
194: for (col=0; col<nxloc; col++) {
195: dx = xx[col];
196: if (PetscAbsScalar(dx) < umin) {
197: if (PetscRealPart(dx) >= 0.0) dx = umin;
198: else if (PetscRealPart(dx) < 0.0 ) dx = -umin;
199: }
200: dx *= epsilon;
201: vscale_array[col] = 1.0/dx;
202: }
203: VecRestoreArray(x1,&xx);
204: VecRestoreArray(vscale,&vscale_array);
205: }
206: if (ctype == IS_COLORING_GLOBAL && coloring->htype[0] == 'd') {
207: VecGhostUpdateBegin(vscale,INSERT_VALUES,SCATTER_FORWARD);
208: VecGhostUpdateEnd(vscale,INSERT_VALUES,SCATTER_FORWARD);
209: }
211: /* (3) Loop over each color */
212: if (!coloring->w3) {
213: VecDuplicate(x1,&coloring->w3);
214: PetscLogObjectParent((PetscObject)coloring,(PetscObject)coloring->w3);
215: }
216: w3 = coloring->w3;
218: VecGetOwnershipRange(x1,&cstart,&cend); /* used by ghosted vscale */
219: if (vscale) {
220: VecGetArray(vscale,&vscale_array);
221: }
222: nz = 0;
224: if (coloring->bcols > 1) { /* use blocked insertion of Jentry */
225: PetscInt i,m=J->rmap->n,nbcols,bcols=coloring->bcols;
226: PetscScalar *dy=coloring->dy,*dy_k;
228: nbcols = 0;
229: for (k=0; k<ncolors; k+=bcols) {
230: coloring->currentcolor = k;
232: /*
233: (3-1) Loop over each column associated with color
234: adding the perturbation to the vector w3 = x1 + dx.
235: */
237: dy_k = dy;
238: if (k + bcols > ncolors) bcols = ncolors - k;
239: for (i=0; i<bcols; i++) {
241: VecCopy(x1,w3);
242: VecGetArray(w3,&w3_array);
243: if (ctype == IS_COLORING_GLOBAL) w3_array -= cstart; /* shift pointer so global index can be used */
244: if (coloring->htype[0] == 'w') {
245: for (l=0; l<ncolumns[k+i]; l++) {
246: col = coloring->columns[k+i][l]; /* local column (in global index!) of the matrix we are probing for */
247: w3_array[col] += 1.0/dx;
248: }
249: } else { /* htype == 'ds' */
250: vscale_array -= cstart; /* shift pointer so global index can be used */
251: for (l=0; l<ncolumns[k+i]; l++) {
252: col = coloring->columns[k+i][l]; /* local column (in global index!) of the matrix we are probing for */
253: w3_array[col] += 1.0/vscale_array[col];
254: }
255: vscale_array += cstart;
256: }
257: if (ctype == IS_COLORING_GLOBAL) w3_array += cstart;
258: VecRestoreArray(w3,&w3_array);
260: /*
261: (3-2) Evaluate function at w3 = x1 + dx (here dx is a vector of perturbations)
262: w2 = F(x1 + dx) - F(x1)
263: */
264: PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);
265: VecPlaceArray(w2,dy_k); /* place w2 to the array dy_i */
266: (*f)(sctx,w3,w2,fctx);
267: PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);
268: VecAXPY(w2,-1.0,w1);
269: VecResetArray(w2);
270: dy_k += m; /* points to dy+i*nxloc */
271: }
273: /*
274: (3-3) Loop over block rows of vector, putting results into Jacobian matrix
275: */
276: nrows_k = nrows[nbcols++];
277: VecGetArray(w2,&y);
279: if (coloring->htype[0] == 'w') {
280: for (l=0; l<nrows_k; l++) {
281: row = Jentry2[nz].row; /* local row index */
282: *(Jentry2[nz++].valaddr) = dy[row]*dx;
283: }
284: } else { /* htype == 'ds' */
285: for (l=0; l<nrows_k; l++) {
286: row = Jentry[nz].row; /* local row index */
287: *(Jentry[nz].valaddr) = dy[row]*vscale_array[Jentry[nz].col];
288: nz++;
289: }
290: }
291: VecRestoreArray(w2,&y);
292: }
293: } else { /* bcols == 1 */
294: for (k=0; k<ncolors; k++) {
295: coloring->currentcolor = k;
297: /*
298: (3-1) Loop over each column associated with color
299: adding the perturbation to the vector w3 = x1 + dx.
300: */
301: VecCopy(x1,w3);
302: VecGetArray(w3,&w3_array);
303: if (ctype == IS_COLORING_GLOBAL) w3_array -= cstart; /* shift pointer so global index can be used */
304: if (coloring->htype[0] == 'w') {
305: for (l=0; l<ncolumns[k]; l++) {
306: col = coloring->columns[k][l]; /* local column (in global index!) of the matrix we are probing for */
307: w3_array[col] += 1.0/dx;
308: }
309: } else { /* htype == 'ds' */
310: vscale_array -= cstart; /* shift pointer so global index can be used */
311: for (l=0; l<ncolumns[k]; l++) {
312: col = coloring->columns[k][l]; /* local column (in global index!) of the matrix we are probing for */
313: w3_array[col] += 1.0/vscale_array[col];
314: }
315: vscale_array += cstart;
316: }
317: if (ctype == IS_COLORING_GLOBAL) w3_array += cstart;
318: VecRestoreArray(w3,&w3_array);
320: /*
321: (3-2) Evaluate function at w3 = x1 + dx (here dx is a vector of perturbations)
322: w2 = F(x1 + dx) - F(x1)
323: */
324: PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);
325: (*f)(sctx,w3,w2,fctx);
326: PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);
327: VecAXPY(w2,-1.0,w1);
329: /*
330: (3-3) Loop over rows of vector, putting results into Jacobian matrix
331: */
332: nrows_k = nrows[k];
333: VecGetArray(w2,&y);
334: if (coloring->htype[0] == 'w') {
335: for (l=0; l<nrows_k; l++) {
336: row = Jentry2[nz].row; /* local row index */
337: *(Jentry2[nz++].valaddr) = y[row]*dx;
338: }
339: } else { /* htype == 'ds' */
340: for (l=0; l<nrows_k; l++) {
341: row = Jentry[nz].row; /* local row index */
342: *(Jentry[nz].valaddr) = y[row]*vscale_array[Jentry[nz].col];
343: nz++;
344: }
345: }
346: VecRestoreArray(w2,&y);
347: }
348: }
350: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
351: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
352: if (vscale) {
353: VecRestoreArray(vscale,&vscale_array);
354: }
355: coloring->currentcolor = -1;
356: return(0);
357: }
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;
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_GHOSTED) {
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: if (isBAIJ) {
392: Mat_MPIBAIJ *baij=(Mat_MPIBAIJ*)mat->data;
393: Mat_SeqBAIJ *spA,*spB;
394: A = baij->A; spA = (Mat_SeqBAIJ*)A->data; A_val = spA->a;
395: B = baij->B; spB = (Mat_SeqBAIJ*)B->data; B_val = spB->a;
396: nz = spA->nz + spB->nz; /* total nonzero entries of mat */
397: if (!baij->colmap) {
398: MatCreateColmap_MPIBAIJ_Private(mat);
399: colmap = baij->colmap;
400: }
401: MatGetColumnIJ_SeqBAIJ_Color(A,0,PETSC_FALSE,PETSC_FALSE,&ncols,&A_ci,&A_cj,&spidxA,NULL);
402: MatGetColumnIJ_SeqBAIJ_Color(B,0,PETSC_FALSE,PETSC_FALSE,&ncols,&B_ci,&B_cj,&spidxB,NULL);
404: if (ctype == IS_COLORING_GLOBAL && c->htype[0] == 'd') { /* create vscale for storing dx */
405: PetscInt *garray;
406: PetscMalloc1(B->cmap->n,&garray);
407: for (i=0; i<baij->B->cmap->n/bs; i++) {
408: for (j=0; j<bs; j++) {
409: garray[i*bs+j] = bs*baij->garray[i]+j;
410: }
411: }
412: VecCreateGhost(PetscObjectComm((PetscObject)mat),mat->cmap->n,PETSC_DETERMINE,B->cmap->n,garray,&c->vscale);
413: PetscFree(garray);
414: }
415: } else {
416: Mat_MPIAIJ *aij=(Mat_MPIAIJ*)mat->data;
417: Mat_SeqAIJ *spA,*spB;
418: A = aij->A; spA = (Mat_SeqAIJ*)A->data; A_val = spA->a;
419: B = aij->B; spB = (Mat_SeqAIJ*)B->data; B_val = spB->a;
420: nz = spA->nz + spB->nz; /* total nonzero entries of mat */
421: if (!aij->colmap) {
422: /* Allow access to data structures of local part of matrix
423: - creates aij->colmap which maps global column number to local number in part B */
424: MatCreateColmap_MPIAIJ_Private(mat);
425: colmap = aij->colmap;
426: }
427: MatGetColumnIJ_SeqAIJ_Color(A,0,PETSC_FALSE,PETSC_FALSE,&ncols,&A_ci,&A_cj,&spidxA,NULL);
428: MatGetColumnIJ_SeqAIJ_Color(B,0,PETSC_FALSE,PETSC_FALSE,&ncols,&B_ci,&B_cj,&spidxB,NULL);
430: bs = 1; /* only bs=1 is supported for non MPIBAIJ matrix */
432: if (ctype == IS_COLORING_GLOBAL && c->htype[0] == 'd') { /* create vscale for storing dx */
433: VecCreateGhost(PetscObjectComm((PetscObject)mat),mat->cmap->n,PETSC_DETERMINE,B->cmap->n,aij->garray,&c->vscale);
434: }
435: }
437: m = mat->rmap->n/bs;
438: cstart = mat->cmap->rstart/bs;
439: cend = mat->cmap->rend/bs;
441: PetscMalloc1(nis,&c->ncolumns);
442: PetscMalloc1(nis,&c->columns);
443: PetscMalloc1(nis,&c->nrows);
444: PetscLogObjectMemory((PetscObject)c,3*nis*sizeof(PetscInt));
446: if (c->htype[0] == 'd') {
447: PetscMalloc1(nz,&Jentry);
448: PetscLogObjectMemory((PetscObject)c,nz*sizeof(MatEntry));
449: c->matentry = Jentry;
450: } else if (c->htype[0] == 'w') {
451: PetscMalloc1(nz,&Jentry2);
452: PetscLogObjectMemory((PetscObject)c,nz*sizeof(MatEntry2));
453: c->matentry2 = Jentry2;
454: } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"htype is not supported");
456: PetscMalloc2(m+1,&rowhit,m+1,&valaddrhit);
457: nz = 0;
458: ISColoringGetIS(iscoloring,PETSC_IGNORE,&isa);
459: for (i=0; i<nis; i++) { /* for each local color */
460: ISGetLocalSize(isa[i],&n);
461: ISGetIndices(isa[i],&is);
463: c->ncolumns[i] = n; /* local number of columns of this color on this process */
464: if (n) {
465: PetscMalloc1(n,&c->columns[i]);
466: PetscLogObjectMemory((PetscObject)c,n*sizeof(PetscInt));
467: PetscMemcpy(c->columns[i],is,n*sizeof(PetscInt));
468: } else {
469: c->columns[i] = 0;
470: }
472: if (ctype == IS_COLORING_GLOBAL) {
473: /* Determine nctot, the total (parallel) number of columns of this color */
474: MPI_Comm_size(PetscObjectComm((PetscObject)mat),&size);
475: PetscMalloc2(size,&ncolsonproc,size,&disp);
477: /* ncolsonproc[j]: local ncolumns on proc[j] of this color */
478: PetscMPIIntCast(n,&nn);
479: MPI_Allgather(&nn,1,MPI_INT,ncolsonproc,1,MPI_INT,PetscObjectComm((PetscObject)mat));
480: nctot = 0; for (j=0; j<size; j++) nctot += ncolsonproc[j];
481: if (!nctot) {
482: PetscInfo(mat,"Coloring of matrix has some unneeded colors with no corresponding rows\n");
483: }
485: disp[0] = 0;
486: for (j=1; j<size; j++) {
487: disp[j] = disp[j-1] + ncolsonproc[j-1];
488: }
490: /* Get cols, the complete list of columns for this color on each process */
491: PetscMalloc1((nctot+1),&cols);
492: MPI_Allgatherv((void*)is,n,MPIU_INT,cols,ncolsonproc,disp,MPIU_INT,PetscObjectComm((PetscObject)mat));
493: PetscFree2(ncolsonproc,disp);
494: } else if (ctype == IS_COLORING_GHOSTED) {
495: /* Determine local number of columns of this color on this process, including ghost points */
496: nctot = n;
497: PetscMalloc1((nctot+1),&cols);
498: PetscMemcpy(cols,is,n*sizeof(PetscInt));
499: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not provided for this MatFDColoring type");
501: /* Mark all rows affect by these columns */
502: PetscMemzero(rowhit,m*sizeof(PetscInt));
503: bs2 = bs*bs;
504: nrows_i = 0;
505: for (j=0; j<nctot; j++) { /* loop over columns*/
506: if (ctype == IS_COLORING_GHOSTED) {
507: col = ltog[cols[j]];
508: } else {
509: col = cols[j];
510: }
511: if (col >= cstart && col < cend) { /* column is in A, diagonal block of mat */
512: row = A_cj + A_ci[col-cstart];
513: nrows = A_ci[col-cstart+1] - A_ci[col-cstart];
514: nrows_i += nrows;
515: /* loop over columns of A marking them in rowhit */
516: for (k=0; k<nrows; k++) {
517: /* set valaddrhit for part A */
518: spidx = bs2*spidxA[A_ci[col-cstart] + k];
519: valaddrhit[*row] = &A_val[spidx];
520: rowhit[*row++] = col - cstart + 1; /* local column index */
521: }
522: } else { /* column is in B, off-diagonal block of mat */
523: #if defined(PETSC_USE_CTABLE)
524: PetscTableFind(colmap,col+1,&colb);
525: colb--;
526: #else
527: colb = colmap[col] - 1; /* local column index */
528: #endif
529: if (colb == -1) {
530: nrows = 0;
531: } else {
532: colb = colb/bs;
533: row = B_cj + B_ci[colb];
534: nrows = B_ci[colb+1] - B_ci[colb];
535: }
536: nrows_i += nrows;
537: /* loop over columns of B marking them in rowhit */
538: for (k=0; k<nrows; k++) {
539: /* set valaddrhit for part B */
540: spidx = bs2*spidxB[B_ci[colb] + k];
541: valaddrhit[*row] = &B_val[spidx];
542: rowhit[*row++] = colb + 1 + cend - cstart; /* local column index */
543: }
544: }
545: }
546: c->nrows[i] = nrows_i;
548: if (c->htype[0] == 'd') {
549: for (j=0; j<m; j++) {
550: if (rowhit[j]) {
551: Jentry[nz].row = j; /* local row index */
552: Jentry[nz].col = rowhit[j] - 1; /* local column index */
553: Jentry[nz].valaddr = valaddrhit[j]; /* address of mat value for this entry */
554: nz++;
555: }
556: }
557: } else { /* c->htype == 'wp' */
558: for (j=0; j<m; j++) {
559: if (rowhit[j]) {
560: Jentry2[nz].row = j; /* local row index */
561: Jentry2[nz].valaddr = valaddrhit[j]; /* address of mat value for this entry */
562: nz++;
563: }
564: }
565: }
566: PetscFree(cols);
567: }
569: if (bcols > 1) { /* reorder Jentry for faster MatFDColoringApply() */
570: MatFDColoringSetUpBlocked_AIJ_Private(mat,c,nz);
571: }
573: if (isBAIJ) {
574: MatRestoreColumnIJ_SeqBAIJ_Color(A,0,PETSC_FALSE,PETSC_FALSE,&ncols,&A_ci,&A_cj,&spidxA,NULL);
575: MatRestoreColumnIJ_SeqBAIJ_Color(B,0,PETSC_FALSE,PETSC_FALSE,&ncols,&B_ci,&B_cj,&spidxB,NULL);
576: PetscMalloc1(bs*mat->rmap->n,&c->dy);
577: } else {
578: MatRestoreColumnIJ_SeqAIJ_Color(A,0,PETSC_FALSE,PETSC_FALSE,&ncols,&A_ci,&A_cj,&spidxA,NULL);
579: MatRestoreColumnIJ_SeqAIJ_Color(B,0,PETSC_FALSE,PETSC_FALSE,&ncols,&B_ci,&B_cj,&spidxB,NULL);
580: }
582: ISColoringRestoreIS(iscoloring,&isa);
583: PetscFree2(rowhit,valaddrhit);
585: if (ctype == IS_COLORING_GHOSTED) {
586: ISLocalToGlobalMappingRestoreIndices(map,<og);
587: }
588: PetscInfo3(c,"ncolors %D, brows %D and bcols %D are used.\n",c->ncolors,c->brows,c->bcols);
589: return(0);
590: }
594: PetscErrorCode MatFDColoringCreate_MPIXAIJ(Mat mat,ISColoring iscoloring,MatFDColoring c)
595: {
597: PetscInt bs,nis=iscoloring->n,m=mat->rmap->n;
598: PetscBool isBAIJ;
601: /* set default brows and bcols for speedup inserting the dense matrix into sparse Jacobian;
602: bcols is chosen s.t. dy-array takes 50% of memory space as mat */
603: MatGetBlockSize(mat,&bs);
604: PetscObjectTypeCompare((PetscObject)mat,MATMPIBAIJ,&isBAIJ);
605: if (isBAIJ) {
606: c->brows = m;
607: c->bcols = 1;
608: } else { /* mpiaij matrix */
609: /* bcols is chosen s.t. dy-array takes 50% of local memory space as mat */
610: Mat_MPIAIJ *aij=(Mat_MPIAIJ*)mat->data;
611: Mat_SeqAIJ *spA,*spB;
612: Mat A,B;
613: PetscInt nz,brows,bcols;
614: PetscReal mem;
616: bs = 1; /* only bs=1 is supported for MPIAIJ matrix */
618: A = aij->A; spA = (Mat_SeqAIJ*)A->data;
619: B = aij->B; spB = (Mat_SeqAIJ*)B->data;
620: nz = spA->nz + spB->nz; /* total local nonzero entries of mat */
621: mem = nz*(sizeof(PetscScalar) + sizeof(PetscInt)) + 3*m*sizeof(PetscInt);
622: bcols = (PetscInt)(0.5*mem /(m*sizeof(PetscScalar)));
623: brows = 1000/bcols;
624: if (bcols > nis) bcols = nis;
625: if (brows == 0 || brows > m) brows = m;
626: c->brows = brows;
627: c->bcols = bcols;
628: }
630: c->M = mat->rmap->N/bs; /* set the global rows and columns and local rows */
631: c->N = mat->cmap->N/bs;
632: c->m = mat->rmap->n/bs;
633: c->rstart = mat->rmap->rstart/bs;
634: c->ncolors = nis;
635: return(0);
636: }