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