xref: /petsc/src/mat/impls/aij/mpi/fdmpiaij.c (revision 9e917edb6df560653d8ab07c604d4ed93450162d) 
1 
2 #include <../src/mat/impls/aij/mpi/mpiaij.h>
3 
4 extern PetscErrorCode MatCreateColmap_MPIAIJ_Private(Mat);
5 
6 #undef __FUNCT__
7 #define __FUNCT__ "MatFDColoringApply_MPIAIJ"
8 PetscErrorCode  MatFDColoringApply_MPIAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx)
9 {
10   PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void*))coloring->f;
11   PetscErrorCode ierr;
12   PetscInt       k,cstart,cend,l,row,col,colb,nz;
13   PetscScalar    dx=0.0,*y,*xx,*w3_array;
14   PetscScalar    *vscale_array;
15   PetscReal      epsilon=coloring->error_rel,umin=coloring->umin,unorm;
16   Vec            w1=coloring->w1,w2=coloring->w2,w3,vscale=coloring->vscale;
17   void           *fctx=coloring->fctx;
18   PetscBool      flg=PETSC_FALSE;
19   PetscInt       ctype=coloring->ctype,nxloc;
20   Mat_MPIAIJ     *aij=(Mat_MPIAIJ*)J->data;
21   MatEntry       *Jentry=coloring->matentry;
22   const PetscInt ncolors=coloring->ncolors,*ncolumns=coloring->ncolumns,*nrows=coloring->nrows;
23 
24   PetscFunctionBegin;
25   ierr = MatSetUnfactored(J);CHKERRQ(ierr);
26   ierr = PetscOptionsGetBool(NULL,"-mat_fd_coloring_dont_rezero",&flg,NULL);CHKERRQ(ierr);
27   if (flg) {
28     ierr = PetscInfo(coloring,"Not calling MatZeroEntries()\n");CHKERRQ(ierr);
29   } else {
30     PetscBool assembled;
31     ierr = MatAssembled(J,&assembled);CHKERRQ(ierr);
32     if (assembled) {
33       ierr = MatZeroEntries(J);CHKERRQ(ierr);
34     }
35   }
36 
37   /* create vscale for storing dx */
38   if (!vscale) {
39     if (ctype == IS_COLORING_GLOBAL && coloring->htype[0] == 'd') {
40       ierr = VecCreateGhost(PetscObjectComm((PetscObject)J),J->cmap->n,PETSC_DETERMINE,aij->B->cmap->n,aij->garray,&vscale);CHKERRQ(ierr);
41     } else if (ctype == IS_COLORING_GHOSTED) {
42       ierr = VecDuplicate(x1,&vscale);CHKERRQ(ierr);
43     }
44     coloring->vscale = vscale;
45   }
46 
47   /* (1) Set w1 = F(x1) */
48   if (!coloring->fset) {
49     ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr);
50     ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr);
51     ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr);
52   } else {
53     coloring->fset = PETSC_FALSE;
54   }
55 
56   /* (2) Compute vscale = 1./dx - the local scale factors, including ghost points */
57   if (coloring->htype[0] == 'w') {
58     /* vscale = dx is a constant scalar */
59     ierr = VecNorm(x1,NORM_2,&unorm);CHKERRQ(ierr);
60     dx = 1.0/((1.0 + unorm)*epsilon);
61   } else {
62     ierr = VecGetLocalSize(x1,&nxloc);CHKERRQ(ierr);
63     ierr = VecGetArray(x1,&xx);CHKERRQ(ierr);
64     ierr = VecGetArray(vscale,&vscale_array);CHKERRQ(ierr);
65     for (col=0; col<nxloc; col++) {
66       dx = xx[col];
67       if (PetscAbsScalar(dx) < umin) {
68         if (PetscRealPart(dx) >= 0.0)      dx = umin;
69         else if (PetscRealPart(dx) < 0.0 ) dx = -umin;
70       }
71       dx               *= epsilon;
72       vscale_array[col] = 1.0/dx;
73     }
74     ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr);
75     ierr = VecRestoreArray(vscale,&vscale_array);CHKERRQ(ierr);
76   }
77   if (ctype == IS_COLORING_GLOBAL && coloring->htype[0] != 'w') {
78     ierr = VecGhostUpdateBegin(vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
79     ierr = VecGhostUpdateEnd(vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
80   }
81 
82   /* (3) Loop over each color */
83   if (!coloring->w3) {
84     ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr);
85     ierr = PetscLogObjectParent((PetscObject)coloring,(PetscObject)coloring->w3);CHKERRQ(ierr);
86   }
87   w3 = coloring->w3;
88 
89   ierr = VecGetOwnershipRange(x1,&cstart,&cend);CHKERRQ(ierr); /* used by ghosted vscale */
90   if (vscale) {
91     ierr = VecGetArray(vscale,&vscale_array);CHKERRQ(ierr);
92   }
93   nz   = 0;
94   for (k=0; k<ncolors; k++) {
95     coloring->currentcolor = k;
96 
97     /*
98       (3-1) Loop over each column associated with color
99       adding the perturbation to the vector w3 = x1 + dx.
100     */
101     ierr = VecCopy(x1,w3);CHKERRQ(ierr);
102     ierr = VecGetArray(w3,&w3_array);CHKERRQ(ierr);
103     if (ctype == IS_COLORING_GLOBAL) w3_array = w3_array - cstart; /* shift pointer so global index can be used */
104     for (l=0; l<ncolumns[k]; l++) {
105       col = coloring->columns[k][l]; /* local column (in global index!) of the matrix we are probing for */
106       if (coloring->htype[0] == 'w') {
107         w3_array[col] += 1.0/dx;
108         continue;
109       }
110       /* convert global col to local colb */
111       if (col >= cstart && col < cend) {
112         colb = col - cstart;
113       } else {
114 #if defined(PETSC_USE_CTABLE)
115         ierr = PetscTableFind(aij->colmap,col+1,&colb);CHKERRQ(ierr);
116         colb--;
117 #else
118         colb = aij->colmap[col] - 1; /* local column index */
119 #endif
120         colb += cend - cstart;
121       }
122       w3_array[col] += 1.0/vscale_array[colb];
123     }
124     if (ctype == IS_COLORING_GLOBAL) w3_array = w3_array + cstart;
125     ierr = VecRestoreArray(w3,&w3_array);CHKERRQ(ierr);
126 
127     /*
128       (3-2) Evaluate function at w3 = x1 + dx (here dx is a vector of perturbations)
129                            w2 = F(x1 + dx) - F(x1)
130     */
131     ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr);
132     ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr);
133     ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr);
134     ierr = VecAXPY(w2,-1.0,w1);CHKERRQ(ierr);
135 
136     /*
137      (3-3) Loop over rows of vector, putting results into Jacobian matrix
138     */
139     ierr = VecGetArray(w2,&y);CHKERRQ(ierr);
140     for (l=0; l<nrows[k]; l++) {
141       row                   = Jentry[nz].row;   /* local row index */
142       if (coloring->htype[0] == 'w') {
143         *(Jentry[nz].valaddr) = y[row]*dx;
144       } else {
145         *(Jentry[nz].valaddr) = y[row]*vscale_array[Jentry[nz].col];
146       }
147       nz++;
148     }
149     ierr = VecRestoreArray(w2,&y);CHKERRQ(ierr);
150   }
151   ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
152   ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
153   if (vscale) {
154     ierr = VecRestoreArray(vscale,&vscale_array);CHKERRQ(ierr);
155   }
156 
157   coloring->currentcolor = -1;
158   PetscFunctionReturn(0);
159 }
160 
161 #undef __FUNCT__
162 #define __FUNCT__ "MatFDColoringCreate_MPIAIJ_new"
163 PetscErrorCode MatFDColoringCreate_MPIAIJ_new(Mat mat,ISColoring iscoloring,MatFDColoring c)
164 {
165   Mat_MPIAIJ             *aij=(Mat_MPIAIJ*)mat->data;
166   PetscErrorCode         ierr;
167   PetscMPIInt            size,*ncolsonproc,*disp,nn;
168   PetscInt               i,n,nrows,j,k,m,ncols,col;
169   const PetscInt         *is,*A_ci,*A_cj,*B_ci,*B_cj,*row = NULL,*ltog=NULL;
170   PetscInt               nis=iscoloring->n,nctot,*cols;
171   PetscInt               *rowhit,cstart,cend,colb;
172   IS                     *isa;
173   ISLocalToGlobalMapping map=mat->cmap->mapping;
174   PetscInt               ctype=c->ctype;
175   Mat                    A=aij->A,B=aij->B;
176   Mat_SeqAIJ             *spA=(Mat_SeqAIJ*)A->data,*spB=(Mat_SeqAIJ*)B->data;
177   PetscScalar            *A_val=spA->a,*B_val=spB->a;
178   PetscInt               spidx;
179   PetscInt               *spidxA,*spidxB,nz;
180   PetscScalar            **valaddrhit;
181   MatEntry               *Jentry;
182 
183   PetscFunctionBegin;
184   if (ctype == IS_COLORING_GHOSTED) {
185     if(!map) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_INCOMP,"When using ghosted differencing matrix must have local to global mapping provided with MatSetLocalToGlobalMapping");
186     ierr = ISLocalToGlobalMappingGetIndices(map,&ltog);CHKERRQ(ierr);
187   }
188 
189   m         = mat->rmap->n;
190   cstart    = mat->cmap->rstart;
191   cend      = mat->cmap->rend;
192   c->M      = mat->rmap->N;         /* set the global rows and columns and local rows */
193   c->N      = mat->cmap->N;
194   c->m      = m;
195   c->rstart = mat->rmap->rstart;
196 
197   c->ncolors = nis;
198   ierr       = PetscMalloc(nis*sizeof(PetscInt),&c->ncolumns);CHKERRQ(ierr);
199   ierr       = PetscMalloc(nis*sizeof(PetscInt*),&c->columns);CHKERRQ(ierr);
200   ierr       = PetscMalloc(nis*sizeof(PetscInt),&c->nrows);CHKERRQ(ierr);
201   ierr       = PetscLogObjectMemory((PetscObject)c,3*nis*sizeof(PetscInt));CHKERRQ(ierr);
202 
203   nz         = spA->nz + spB->nz; /* total nonzero entries of mat */
204   ierr       = PetscMalloc(nz*sizeof(MatEntry),&Jentry);CHKERRQ(ierr);
205   ierr       = PetscLogObjectMemory((PetscObject)c,nz*sizeof(MatEntry));CHKERRQ(ierr);
206   c->matentry = Jentry;
207 
208   /* Allow access to data structures of local part of matrix
209    - creates aij->colmap which maps global column number to local number in part B */
210   if (!aij->colmap) {
211     ierr = MatCreateColmap_MPIAIJ_Private(mat);CHKERRQ(ierr);
212   }
213   ierr = MatGetColumnIJ_SeqAIJ_Color(aij->A,0,PETSC_FALSE,PETSC_FALSE,&ncols,&A_ci,&A_cj,&spidxA,NULL);CHKERRQ(ierr);
214   ierr = MatGetColumnIJ_SeqAIJ_Color(aij->B,0,PETSC_FALSE,PETSC_FALSE,&ncols,&B_ci,&B_cj,&spidxB,NULL);CHKERRQ(ierr);
215 
216   ierr = PetscMalloc2(m+1,PetscInt,&rowhit,m+1,PetscScalar*,&valaddrhit);CHKERRQ(ierr);
217   nz = 0;
218   ierr = ISColoringGetIS(iscoloring,PETSC_IGNORE,&isa);CHKERRQ(ierr);
219   for (i=0; i<nis; i++) { /* for each local color */
220     ierr = ISGetLocalSize(isa[i],&n);CHKERRQ(ierr);
221     ierr = ISGetIndices(isa[i],&is);CHKERRQ(ierr);
222 
223     c->ncolumns[i] = n; /* local number of columns of this color on this process */
224     if (n) {
225       ierr = PetscMalloc(n*sizeof(PetscInt),&c->columns[i]);CHKERRQ(ierr);
226       ierr = PetscLogObjectMemory((PetscObject)c,n*sizeof(PetscInt));CHKERRQ(ierr);
227       ierr = PetscMemcpy(c->columns[i],is,n*sizeof(PetscInt));CHKERRQ(ierr);
228     } else {
229       c->columns[i] = 0;
230     }
231 
232     if (ctype == IS_COLORING_GLOBAL) {
233       /* Determine nctot, the total (parallel) number of columns of this color */
234       ierr = MPI_Comm_size(PetscObjectComm((PetscObject)mat),&size);CHKERRQ(ierr);
235       ierr = PetscMalloc2(size,PetscMPIInt,&ncolsonproc,size,PetscMPIInt,&disp);CHKERRQ(ierr);
236 
237       /* ncolsonproc[j]: local ncolumns on proc[j] of this color */
238       ierr  = PetscMPIIntCast(n,&nn);CHKERRQ(ierr);
239       ierr  = MPI_Allgather(&nn,1,MPI_INT,ncolsonproc,1,MPI_INT,PetscObjectComm((PetscObject)mat));CHKERRQ(ierr);
240       nctot = 0; for (j=0; j<size; j++) nctot += ncolsonproc[j];
241       if (!nctot) {
242         ierr = PetscInfo(mat,"Coloring of matrix has some unneeded colors with no corresponding rows\n");CHKERRQ(ierr);
243       }
244 
245       disp[0] = 0;
246       for (j=1; j<size; j++) {
247         disp[j] = disp[j-1] + ncolsonproc[j-1];
248       }
249 
250       /* Get cols, the complete list of columns for this color on each process */
251       ierr = PetscMalloc((nctot+1)*sizeof(PetscInt),&cols);CHKERRQ(ierr);
252       ierr = MPI_Allgatherv((void*)is,n,MPIU_INT,cols,ncolsonproc,disp,MPIU_INT,PetscObjectComm((PetscObject)mat));CHKERRQ(ierr);
253       ierr = PetscFree2(ncolsonproc,disp);CHKERRQ(ierr);
254     } else if (ctype == IS_COLORING_GHOSTED) {
255       /* Determine local number of columns of this color on this process, including ghost points */
256       nctot = n;
257       ierr  = PetscMalloc((nctot+1)*sizeof(PetscInt),&cols);CHKERRQ(ierr);
258       ierr  = PetscMemcpy(cols,is,n*sizeof(PetscInt));CHKERRQ(ierr);
259     } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not provided for this MatFDColoring type");
260 
261     /* Mark all rows affect by these columns */
262     ierr = PetscMemzero(rowhit,m*sizeof(PetscInt));CHKERRQ(ierr);
263     for (j=0; j<nctot; j++) { /* loop over columns*/
264       if (ctype == IS_COLORING_GHOSTED) {
265         col = ltog[cols[j]];
266       } else {
267         col = cols[j];
268       }
269       if (col >= cstart && col < cend) { /* column is in diagonal block of matrix A */
270         row = A_cj + A_ci[col-cstart];
271         nrows = A_ci[col-cstart+1] - A_ci[col-cstart];
272         /* loop over columns of A marking them in rowhit */
273         for (k=0; k<nrows; k++) {
274           /* set valaddrhit for part A */
275           spidx            = spidxA[A_ci[col-cstart] + k];
276           valaddrhit[*row] = &A_val[spidx];
277           rowhit[*row++]   = col - cstart + 1; /* local column index */
278         }
279       } else { /* column is in off-diagonal block of matrix B */
280 #if defined(PETSC_USE_CTABLE)
281         ierr = PetscTableFind(aij->colmap,col+1,&colb);CHKERRQ(ierr);
282         colb--;
283 #else
284         colb = aij->colmap[col] - 1; /* local column index */
285 #endif
286         if (colb == -1) {
287           nrows = 0;
288         } else {
289           row   = B_cj + B_ci[colb];
290           nrows = B_ci[colb+1] - B_ci[colb];
291         }
292         /* loop over columns of B marking them in rowhit */
293         for (k=0; k<nrows; k++) {
294           /* set valaddrhit for part B */
295           spidx            = spidxB[B_ci[colb] + k];
296           valaddrhit[*row] = &B_val[spidx];
297           rowhit[*row++]   = colb + 1 + cend - cstart; /* local column index */
298         }
299       }
300     }
301 
302     /* count the number of hits */
303     nrows = 0;
304     for (j=0; j<m; j++) {
305       if (rowhit[j]) nrows++;
306     }
307     c->nrows[i] = nrows;
308 
309     for (j=0; j<m; j++) {
310       if (rowhit[j]) {
311         Jentry[nz].row     = j;              /* local row index */
312         Jentry[nz].col     = rowhit[j] - 1;  /* local column index */
313         Jentry[nz].valaddr = valaddrhit[j];  /* address of mat value for this entry */
314         nz++;
315       }
316     }
317     ierr = PetscFree(cols);CHKERRQ(ierr);
318   }
319   ierr = ISColoringRestoreIS(iscoloring,&isa);CHKERRQ(ierr);
320   if (nz != spA->nz + spB->nz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"nz %d != mat->nz %d\n",nz,spA->nz+spB->nz);
321 
322   ierr = PetscFree2(rowhit,valaddrhit);CHKERRQ(ierr);
323   ierr = MatRestoreColumnIJ_SeqAIJ_Color(aij->A,0,PETSC_FALSE,PETSC_FALSE,&ncols,&A_ci,&A_cj,&spidxA,NULL);CHKERRQ(ierr);
324   ierr = MatRestoreColumnIJ_SeqAIJ_Color(aij->B,0,PETSC_FALSE,PETSC_FALSE,&ncols,&B_ci,&B_cj,&spidxB,NULL);CHKERRQ(ierr);
325   if (ctype == IS_COLORING_GHOSTED) {
326     ierr = ISLocalToGlobalMappingRestoreIndices(map,&ltog);CHKERRQ(ierr);
327   }
328 
329   mat->ops->fdcoloringapply = MatFDColoringApply_MPIAIJ;
330   PetscFunctionReturn(0);
331 }
332 
333 /*------------------------------------------------------*/
334 #undef __FUNCT__
335 #define __FUNCT__ "MatFDColoringCreate_MPIAIJ"
336 PetscErrorCode MatFDColoringCreate_MPIAIJ(Mat mat,ISColoring iscoloring,MatFDColoring c)
337 {
338   Mat_MPIAIJ             *aij = (Mat_MPIAIJ*)mat->data;
339   PetscErrorCode         ierr;
340   PetscMPIInt            size,*ncolsonproc,*disp,nn;
341   PetscInt               i,n,nrows,j,k,m,ncols,col;
342   const PetscInt         *is,*A_ci,*A_cj,*B_ci,*B_cj,*rows = 0,*ltog;
343   PetscInt               nis = iscoloring->n,nctot,*cols;
344   PetscInt               *rowhit,M,cstart,cend,colb;
345   PetscInt               *columnsforrow,l;
346   IS                     *isa;
347   PetscBool              done,flg;
348   ISLocalToGlobalMapping map   = mat->cmap->mapping;
349   PetscInt               ctype=c->ctype;
350   PetscBool              new_impl=PETSC_FALSE;
351 
352   PetscFunctionBegin;
353   ierr = PetscOptionsName("-new","using new impls","",&new_impl);CHKERRQ(ierr);
354   if (new_impl){
355     ierr =  MatFDColoringCreate_MPIAIJ_new(mat,iscoloring,c);CHKERRQ(ierr);
356     PetscFunctionReturn(0);
357   }
358   if (ctype == IS_COLORING_GHOSTED && !map) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_INCOMP,"When using ghosted differencing matrix must have local to global mapping provided with MatSetLocalToGlobalMapping");
359 
360   if (map) {ierr = ISLocalToGlobalMappingGetIndices(map,&ltog);CHKERRQ(ierr);}
361   else     ltog = NULL;
362   ierr = ISColoringGetIS(iscoloring,PETSC_IGNORE,&isa);CHKERRQ(ierr);
363 
364   M         = mat->rmap->n;
365   cstart    = mat->cmap->rstart;
366   cend      = mat->cmap->rend;
367   c->M      = mat->rmap->N;         /* set the global rows and columns and local rows */
368   c->N      = mat->cmap->N;
369   c->m      = mat->rmap->n;
370   c->rstart = mat->rmap->rstart;
371 
372   c->ncolors = nis;
373   ierr       = PetscMalloc(nis*sizeof(PetscInt),&c->ncolumns);CHKERRQ(ierr);
374   ierr       = PetscMalloc(nis*sizeof(PetscInt*),&c->columns);CHKERRQ(ierr);
375   ierr       = PetscMalloc(nis*sizeof(PetscInt),&c->nrows);CHKERRQ(ierr);
376   ierr       = PetscMalloc(nis*sizeof(PetscInt*),&c->rows);CHKERRQ(ierr);
377   ierr       = PetscMalloc(nis*sizeof(PetscInt*),&c->columnsforrow);CHKERRQ(ierr);
378   ierr       = PetscLogObjectMemory((PetscObject)c,5*nis*sizeof(PetscInt));CHKERRQ(ierr);
379 
380   /* Allow access to data structures of local part of matrix */
381   if (!aij->colmap) {
382     ierr = MatCreateColmap_MPIAIJ_Private(mat);CHKERRQ(ierr);
383   }
384   ierr = MatGetColumnIJ(aij->A,0,PETSC_FALSE,PETSC_FALSE,&ncols,&A_ci,&A_cj,&done);CHKERRQ(ierr);
385   ierr = MatGetColumnIJ(aij->B,0,PETSC_FALSE,PETSC_FALSE,&ncols,&B_ci,&B_cj,&done);CHKERRQ(ierr);
386 
387   ierr = PetscMalloc((M+1)*sizeof(PetscInt),&rowhit);CHKERRQ(ierr);
388   ierr = PetscMalloc((M+1)*sizeof(PetscInt),&columnsforrow);CHKERRQ(ierr);
389 
390   for (i=0; i<nis; i++) {
391     ierr = ISGetLocalSize(isa[i],&n);CHKERRQ(ierr);
392     ierr = ISGetIndices(isa[i],&is);CHKERRQ(ierr);
393 
394     c->ncolumns[i] = n; /* local number of columns of this color on this process */
395     if (n) {
396       ierr = PetscMalloc(n*sizeof(PetscInt),&c->columns[i]);CHKERRQ(ierr);
397       ierr = PetscLogObjectMemory((PetscObject)c,n*sizeof(PetscInt));CHKERRQ(ierr);
398       ierr = PetscMemcpy(c->columns[i],is,n*sizeof(PetscInt));CHKERRQ(ierr);
399     } else {
400       c->columns[i] = 0;
401     }
402 
403     if (ctype == IS_COLORING_GLOBAL) {
404       /* Determine the total (parallel) number of columns of this color */
405       ierr = MPI_Comm_size(PetscObjectComm((PetscObject)mat),&size);CHKERRQ(ierr);
406       ierr = PetscMalloc2(size,PetscMPIInt,&ncolsonproc,size,PetscMPIInt,&disp);CHKERRQ(ierr);
407 
408       ierr  = PetscMPIIntCast(n,&nn);CHKERRQ(ierr);
409       ierr  = MPI_Allgather(&nn,1,MPI_INT,ncolsonproc,1,MPI_INT,PetscObjectComm((PetscObject)mat));CHKERRQ(ierr);
410       nctot = 0; for (j=0; j<size; j++) nctot += ncolsonproc[j];
411       if (!nctot) {
412         ierr = PetscInfo(mat,"Coloring of matrix has some unneeded colors with no corresponding rows\n");CHKERRQ(ierr);
413       }
414 
415       disp[0] = 0;
416       for (j=1; j<size; j++) {
417         disp[j] = disp[j-1] + ncolsonproc[j-1];
418       }
419 
420       /* Get complete list of columns for color on each processor */
421       ierr = PetscMalloc((nctot+1)*sizeof(PetscInt),&cols);CHKERRQ(ierr);
422       ierr = MPI_Allgatherv((void*)is,n,MPIU_INT,cols,ncolsonproc,disp,MPIU_INT,PetscObjectComm((PetscObject)mat));CHKERRQ(ierr);
423       ierr = PetscFree2(ncolsonproc,disp);CHKERRQ(ierr);
424     } else if (ctype == IS_COLORING_GHOSTED) {
425       /* Determine local number of columns of this color on this process, including ghost points */
426       nctot = n;
427       ierr  = PetscMalloc((nctot+1)*sizeof(PetscInt),&cols);CHKERRQ(ierr);
428       ierr  = PetscMemcpy(cols,is,n*sizeof(PetscInt));CHKERRQ(ierr);
429     } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not provided for this MatFDColoring type");
430 
431     /*
432        Mark all rows affect by these columns
433     */
434     /* Temporary option to allow for debugging/testing */
435     flg  = PETSC_FALSE;
436     ierr = PetscOptionsGetBool(NULL,"-matfdcoloring_slow",&flg,NULL);CHKERRQ(ierr);
437     if (!flg) { /*-----------------------------------------------------------------------------*/
438       /* crude, fast version */
439       ierr = PetscMemzero(rowhit,M*sizeof(PetscInt));CHKERRQ(ierr);
440       /* loop over columns*/
441       for (j=0; j<nctot; j++) {
442         if (ctype == IS_COLORING_GHOSTED) {
443           col = ltog[cols[j]];
444         } else {
445           col = cols[j];
446         }
447         if (col >= cstart && col < cend) {
448           /* column is in diagonal block of matrix */
449           rows = A_cj + A_ci[col-cstart];
450           m    = A_ci[col-cstart+1] - A_ci[col-cstart];
451         } else {
452 #if defined(PETSC_USE_CTABLE)
453           ierr = PetscTableFind(aij->colmap,col+1,&colb);CHKERRQ(ierr);
454           colb--;
455 #else
456           colb = aij->colmap[col] - 1;
457 #endif
458           if (colb == -1) {
459             m = 0;
460           } else {
461             rows = B_cj + B_ci[colb];
462             m    = B_ci[colb+1] - B_ci[colb];
463           }
464         }
465         /* loop over columns marking them in rowhit */
466         for (k=0; k<m; k++) {
467           rowhit[*rows++] = col + 1;
468         }
469       }
470 
471       /* count the number of hits */
472       nrows = 0;
473       for (j=0; j<M; j++) {
474         if (rowhit[j]) nrows++;
475       }
476       c->nrows[i] = nrows;
477       ierr        = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->rows[i]);CHKERRQ(ierr);
478       ierr        = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->columnsforrow[i]);CHKERRQ(ierr);
479       ierr        = PetscLogObjectMemory((PetscObject)c,2*(nrows+1)*sizeof(PetscInt));CHKERRQ(ierr);
480       nrows       = 0;
481       for (j=0; j<M; j++) {
482         if (rowhit[j]) {
483           c->rows[i][nrows]          = j;              /* local row index */
484           c->columnsforrow[i][nrows] = rowhit[j] - 1;  /* global column index */
485           nrows++;
486         }
487       }
488     } else { /*-------------------------------------------------------------------------------*/
489       /* slow version, using rowhit as a linked list */
490       PetscInt currentcol,fm,mfm;
491       rowhit[M] = M;
492       nrows     = 0;
493       /* loop over columns*/
494       for (j=0; j<nctot; j++) {
495         if (ctype == IS_COLORING_GHOSTED) {
496           col = ltog[cols[j]];
497         } else {
498           col = cols[j];
499         }
500         if (col >= cstart && col < cend) {
501           /* column is in diagonal block of matrix */
502           rows = A_cj + A_ci[col-cstart];
503           m    = A_ci[col-cstart+1] - A_ci[col-cstart];
504         } else {
505 #if defined(PETSC_USE_CTABLE)
506           ierr = PetscTableFind(aij->colmap,col+1,&colb);CHKERRQ(ierr);
507           colb--;
508 #else
509           colb = aij->colmap[col] - 1;
510 #endif
511           if (colb == -1) {
512             m = 0;
513           } else {
514             rows = B_cj + B_ci[colb];
515             m    = B_ci[colb+1] - B_ci[colb];
516           }
517         }
518 
519         /* loop over columns marking them in rowhit */
520         fm = M;    /* fm points to first entry in linked list */
521         for (k=0; k<m; k++) {
522           currentcol = *rows++;
523           /* is it already in the list? */
524           do {
525             mfm = fm;
526             fm  = rowhit[fm];
527           } while (fm < currentcol);
528           /* not in list so add it */
529           if (fm != currentcol) {
530             nrows++;
531             columnsforrow[currentcol] = col;
532             /* next three lines insert new entry into linked list */
533             rowhit[mfm]        = currentcol;
534             rowhit[currentcol] = fm;
535             fm                 = currentcol;
536             /* fm points to present position in list since we know the columns are sorted */
537           } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Invalid coloring of matrix detected");
538         }
539       }
540       c->nrows[i] = nrows;
541 
542       ierr = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->rows[i]);CHKERRQ(ierr);
543       ierr = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->columnsforrow[i]);CHKERRQ(ierr);
544       ierr = PetscLogObjectMemory((PetscObject)c,(nrows+1)*sizeof(PetscInt));CHKERRQ(ierr);
545       /* now store the linked list of rows into c->rows[i] */
546       nrows = 0;
547       fm    = rowhit[M];
548       do {
549         c->rows[i][nrows]            = fm;
550         c->columnsforrow[i][nrows++] = columnsforrow[fm];
551         fm                           = rowhit[fm];
552       } while (fm < M);
553     } /* ---------------------------------------------------------------------------------------*/
554     ierr = PetscFree(cols);CHKERRQ(ierr);
555   }
556 
557   /* Optimize by adding the vscale, and scaleforrow[][] fields */
558   /*
559        vscale will contain the "diagonal" on processor scalings followed by the off processor
560   */
561   if (ctype == IS_COLORING_GLOBAL) {
562     ierr = VecCreateGhost(PetscObjectComm((PetscObject)mat),aij->A->rmap->n,PETSC_DETERMINE,aij->B->cmap->n,aij->garray,&c->vscale);CHKERRQ(ierr);
563     ierr = PetscMalloc(c->ncolors*sizeof(PetscInt*),&c->vscaleforrow);CHKERRQ(ierr);
564     for (k=0; k<c->ncolors; k++) {
565       ierr = PetscMalloc((c->nrows[k]+1)*sizeof(PetscInt),&c->vscaleforrow[k]);CHKERRQ(ierr);
566       for (l=0; l<c->nrows[k]; l++) {
567         col = c->columnsforrow[k][l];
568         if (col >= cstart && col < cend) {
569           /* column is in diagonal block of matrix */
570           colb = col - cstart;
571         } else {
572           /* column  is in "off-processor" part */
573 #if defined(PETSC_USE_CTABLE)
574           ierr = PetscTableFind(aij->colmap,col+1,&colb);CHKERRQ(ierr);
575           colb--;
576 #else
577           colb = aij->colmap[col] - 1;
578 #endif
579           colb += cend - cstart;
580         }
581         c->vscaleforrow[k][l] = colb;
582       }
583     }
584   } else if (ctype == IS_COLORING_GHOSTED) {
585     /* Get gtol mapping */
586     PetscInt N = mat->cmap->N,nlocal,*gtol;
587     ierr = PetscMalloc((N+1)*sizeof(PetscInt),&gtol);CHKERRQ(ierr);
588     for (i=0; i<N; i++) gtol[i] = -1;
589     ierr = ISLocalToGlobalMappingGetSize(map,&nlocal);CHKERRQ(ierr);
590     for (i=0; i<nlocal; i++) gtol[ltog[i]] = i;
591 
592     c->vscale = 0; /* will be created in MatFDColoringApply() */
593     ierr      = PetscMalloc(c->ncolors*sizeof(PetscInt*),&c->vscaleforrow);CHKERRQ(ierr);
594     for (k=0; k<c->ncolors; k++) {
595       ierr = PetscMalloc((c->nrows[k]+1)*sizeof(PetscInt),&c->vscaleforrow[k]);CHKERRQ(ierr);
596       for (l=0; l<c->nrows[k]; l++) {
597         col = c->columnsforrow[k][l];      /* global column index */
598         c->vscaleforrow[k][l] = gtol[col]; /* local column index */
599       }
600     }
601     ierr = PetscFree(gtol);CHKERRQ(ierr);
602   }
603   ierr = ISColoringRestoreIS(iscoloring,&isa);CHKERRQ(ierr);
604 
605   ierr = PetscFree(rowhit);CHKERRQ(ierr);
606   ierr = PetscFree(columnsforrow);CHKERRQ(ierr);
607   ierr = MatRestoreColumnIJ(aij->A,0,PETSC_FALSE,PETSC_FALSE,&ncols,&A_ci,&A_cj,&done);CHKERRQ(ierr);
608   ierr = MatRestoreColumnIJ(aij->B,0,PETSC_FALSE,PETSC_FALSE,&ncols,&B_ci,&B_cj,&done);CHKERRQ(ierr);
609   if (map) {ierr = ISLocalToGlobalMappingRestoreIndices(map,&ltog);CHKERRQ(ierr);}
610   PetscFunctionReturn(0);
611 }
612 
613 
614 
615 
616 
617 
618