xref: /petsc/src/mat/impls/aij/seq/matmatmult.c (revision d6cbdb99fb73f50d903551d9fc83deb7cc6e3d53) 
1 
2 /*
3   Defines matrix-matrix product routines for pairs of SeqAIJ matrices
4           C = A * B
5 */
6 
7 #include <../src/mat/impls/aij/seq/aij.h> /*I "petscmat.h" I*/
8 #include <../src/mat/utils/freespace.h>
9 #include <petscbt.h>
10 #include <../src/mat/impls/dense/seq/dense.h> /*I "petscmat.h" I*/
11 
12 EXTERN_C_BEGIN
13 #undef __FUNCT__
14 #define __FUNCT__ "MatMatMult_SeqAIJ_SeqAIJ"
15 PetscErrorCode MatMatMult_SeqAIJ_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
16 {
17   PetscErrorCode ierr;
18 
19   PetscFunctionBegin;
20   if (scall == MAT_INITIAL_MATRIX){
21     ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(A,B,fill,C);CHKERRQ(ierr);
22   }
23   ierr = MatMatMultNumeric_SeqAIJ_SeqAIJ(A,B,*C);CHKERRQ(ierr);
24   PetscFunctionReturn(0);
25 }
26 EXTERN_C_END
27 
28 #undef __FUNCT__
29 #define __FUNCT__ "MatMatMultSymbolic_SeqAIJ_SeqAIJ"
30 PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C)
31 {
32   PetscErrorCode     ierr;
33   PetscFreeSpaceList free_space=PETSC_NULL,current_space=PETSC_NULL;
34   Mat_SeqAIJ         *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
35   PetscInt           *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,*bjj,*ci,*cj;
36   PetscInt           am=A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;
37   PetscInt           i,j,anzi,brow,bnzj,cnzi,nlnk,*lnk,nspacedouble=0;
38   MatScalar          *ca;
39   PetscBT            lnkbt;
40   PetscReal          afill;
41 
42   PetscFunctionBegin;
43   /* Allocate ci array, arrays for fill computation and */
44   /* free space for accumulating nonzero column info */
45   ierr = PetscMalloc(((am+1)+1)*sizeof(PetscInt),&ci);CHKERRQ(ierr);
46   ci[0] = 0;
47 
48   /* create and initialize a linked list */
49   nlnk = bn+1;
50   ierr = PetscLLCreate(bn,bn,nlnk,lnk,lnkbt);CHKERRQ(ierr);
51 
52   /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */
53   ierr = PetscFreeSpaceGet((PetscInt)(fill*(ai[am]+bi[bm])),&free_space);CHKERRQ(ierr);
54   current_space = free_space;
55 
56   /* Determine symbolic info for each row of the product: */
57   for (i=0;i<am;i++) {
58     anzi = ai[i+1] - ai[i];
59     cnzi = 0;
60     j    = anzi;
61     aj   = a->j + ai[i];
62     while (j){/* assume cols are almost in increasing order, starting from its end saves computation */
63       j--;
64       brow = *(aj + j);
65       bnzj = bi[brow+1] - bi[brow];
66       bjj  = bj + bi[brow];
67       /* add non-zero cols of B into the sorted linked list lnk */
68       ierr = PetscLLAdd(bnzj,bjj,bn,nlnk,lnk,lnkbt);CHKERRQ(ierr);
69       cnzi += nlnk;
70     }
71 
72     /* If free space is not available, make more free space */
73     /* Double the amount of total space in the list */
74     if (current_space->local_remaining<cnzi) {
75       ierr = PetscFreeSpaceGet(cnzi+current_space->total_array_size,&current_space);CHKERRQ(ierr);
76       nspacedouble++;
77     }
78 
79     /* Copy data into free space, then initialize lnk */
80     ierr = PetscLLClean(bn,bn,cnzi,lnk,current_space->array,lnkbt);CHKERRQ(ierr);
81     current_space->array           += cnzi;
82     current_space->local_used      += cnzi;
83     current_space->local_remaining -= cnzi;
84 
85     ci[i+1] = ci[i] + cnzi;
86   }
87 
88   /* Column indices are in the list of free space */
89   /* Allocate space for cj, initialize cj, and */
90   /* destroy list of free space and other temporary array(s) */
91   ierr = PetscMalloc((ci[am]+1)*sizeof(PetscInt),&cj);CHKERRQ(ierr);
92   ierr = PetscFreeSpaceContiguous(&free_space,cj);CHKERRQ(ierr);
93   ierr = PetscLLDestroy(lnk,lnkbt);CHKERRQ(ierr);
94 
95   /* Allocate space for ca */
96   ierr = PetscMalloc((ci[am]+1)*sizeof(MatScalar),&ca);CHKERRQ(ierr);
97   ierr = PetscMemzero(ca,(ci[am]+1)*sizeof(MatScalar));CHKERRQ(ierr);
98 
99   /* put together the new symbolic matrix */
100   ierr = MatCreateSeqAIJWithArrays(((PetscObject)A)->comm,am,bn,ci,cj,ca,C);CHKERRQ(ierr);
101 
102   /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
103   /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
104   c = (Mat_SeqAIJ *)((*C)->data);
105   c->free_a   = PETSC_TRUE;
106   c->free_ij  = PETSC_TRUE;
107   c->nonew    = 0;
108 
109   /* set MatInfo */
110   afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;
111   if (afill < 1.0) afill = 1.0;
112   c->maxnz                     = ci[am];
113   c->nz                        = ci[am];
114   (*C)->info.mallocs           = nspacedouble;
115   (*C)->info.fill_ratio_given  = fill;
116   (*C)->info.fill_ratio_needed = afill;
117 
118 #if defined(PETSC_USE_INFO)
119   if (ci[am]) {
120     ierr = PetscInfo3((*C),"Reallocs %D; Fill ratio: given %G needed %G.\n",nspacedouble,fill,afill);CHKERRQ(ierr);
121     ierr = PetscInfo1((*C),"Use MatMatMult(A,B,MatReuse,%G,&C) for best performance.;\n",afill);CHKERRQ(ierr);
122   } else {
123     ierr = PetscInfo((*C),"Empty matrix product\n");CHKERRQ(ierr);
124   }
125 #endif
126   PetscFunctionReturn(0);
127 }
128 
129 
130 #undef __FUNCT__
131 #define __FUNCT__ "MatMatMultNumeric_SeqAIJ_SeqAIJ"
132 PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C)
133 {
134   PetscErrorCode ierr;
135   PetscLogDouble flops=0.0;
136   Mat_SeqAIJ     *a = (Mat_SeqAIJ *)A->data;
137   Mat_SeqAIJ     *b = (Mat_SeqAIJ *)B->data;
138   Mat_SeqAIJ     *c = (Mat_SeqAIJ *)C->data;
139   PetscInt       *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,*bjj,*ci=c->i,*cj=c->j;
140   PetscInt       am=A->rmap->N,cm=C->rmap->N;
141   PetscInt       i,j,k,anzi,bnzi,cnzi,brow,nextb;
142   MatScalar      *aa=a->a,*ba=b->a,*baj,*ca=c->a;
143 
144   PetscFunctionBegin;
145   /* clean old values in C */
146   ierr = PetscMemzero(ca,ci[cm]*sizeof(MatScalar));CHKERRQ(ierr);
147   /* Traverse A row-wise. */
148   /* Build the ith row in C by summing over nonzero columns in A, */
149   /* the rows of B corresponding to nonzeros of A. */
150   for (i=0;i<am;i++) {
151     anzi = ai[i+1] - ai[i];
152     for (j=0;j<anzi;j++) {
153       brow = *aj++;
154       bnzi = bi[brow+1] - bi[brow];
155       bjj  = bj + bi[brow];
156       baj  = ba + bi[brow];
157       nextb = 0;
158       for (k=0; nextb<bnzi; k++) {
159         if (cj[k] == bjj[nextb]){ /* ccol == bcol */
160           ca[k] += (*aa)*baj[nextb++];
161         }
162       }
163       flops += 2*bnzi;
164       aa++;
165     }
166     cnzi = ci[i+1] - ci[i];
167     ca += cnzi;
168     cj += cnzi;
169   }
170   ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
171   ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
172 
173   ierr = PetscLogFlops(flops);CHKERRQ(ierr);
174   PetscFunctionReturn(0);
175 }
176 
177 #undef __FUNCT__
178 #define __FUNCT__ "MatMatMultTranspose_SeqAIJ_SeqAIJ"
179 PetscErrorCode MatMatMultTranspose_SeqAIJ_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
180 {
181   PetscErrorCode ierr;
182 
183   PetscFunctionBegin;
184   if (scall == MAT_INITIAL_MATRIX){
185     ierr = MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ(A,B,fill,C);CHKERRQ(ierr);
186   }
187   ierr = MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ(A,B,*C);CHKERRQ(ierr);
188   PetscFunctionReturn(0);
189 }
190 
191 #undef __FUNCT__
192 #define __FUNCT__ "PetscContainerDestroy_Mat_MatMatMultTrans"
193 PetscErrorCode PetscContainerDestroy_Mat_MatMatMultTrans(void *ptr)
194 {
195   PetscErrorCode      ierr;
196   Mat_MatMatMultTrans *multtrans=(Mat_MatMatMultTrans*)ptr;
197 
198   PetscFunctionBegin;
199   ierr = MatMultTransposeColoringDestroy(&multtrans->matcoloring);CHKERRQ(ierr);
200   ierr = MatDestroy(&multtrans->Bt_den);CHKERRQ(ierr);
201   ierr = MatDestroy(&multtrans->ABt_den);CHKERRQ(ierr);
202   ierr = PetscFree(multtrans);CHKERRQ(ierr);
203   PetscFunctionReturn(0);
204 }
205 
206 #undef __FUNCT__
207 #define __FUNCT__ "MatDestroy_SeqAIJ_MatMatMultTrans"
208 PetscErrorCode MatDestroy_SeqAIJ_MatMatMultTrans(Mat A)
209 {
210   PetscErrorCode      ierr;
211   PetscContainer      container;
212   Mat_MatMatMultTrans *multtrans=PETSC_NULL;
213 
214   PetscFunctionBegin;
215   ierr = PetscObjectQuery((PetscObject)A,"Mat_MatMatMultTrans",(PetscObject *)&container);CHKERRQ(ierr);
216   if (!container) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Container does not exit");
217   ierr = PetscContainerGetPointer(container,(void **)&multtrans);CHKERRQ(ierr);
218   A->ops->destroy   = multtrans->destroy;
219   if (A->ops->destroy) {
220     ierr = (*A->ops->destroy)(A);CHKERRQ(ierr);
221   }
222   ierr = PetscObjectCompose((PetscObject)A,"Mat_MatMatMultTrans",0);CHKERRQ(ierr);
223   PetscFunctionReturn(0);
224 }
225 
226 #undef __FUNCT__
227 #define __FUNCT__ "MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ"
228 PetscErrorCode MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C)
229 {
230   PetscErrorCode      ierr;
231   Mat                 Bt;
232   PetscInt            *bti,*btj;
233   Mat_MatMatMultTrans *multtrans;
234   PetscContainer      container;
235 
236   PetscFunctionBegin;
237    /* create symbolic Bt */
238   ierr = MatGetSymbolicTranspose_SeqAIJ(B,&bti,&btj);CHKERRQ(ierr);
239   ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,B->cmap->n,B->rmap->n,bti,btj,PETSC_NULL,&Bt);CHKERRQ(ierr);
240 
241   /* get symbolic C=A*Bt */
242   ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(A,Bt,fill,C);CHKERRQ(ierr);
243 
244   /* clean up */
245   ierr = MatDestroy(&Bt);CHKERRQ(ierr);
246   ierr = MatRestoreSymbolicTranspose_SeqAIJ(B,&bti,&btj);CHKERRQ(ierr);
247 
248   /* create a supporting struct for reuse intermidiate dense matrices with matcoloring */
249   ierr = PetscNew(Mat_MatMatMultTrans,&multtrans);CHKERRQ(ierr);
250 
251   /* attach the supporting struct to C */
252   ierr = PetscContainerCreate(PETSC_COMM_SELF,&container);CHKERRQ(ierr);
253   ierr = PetscContainerSetPointer(container,multtrans);CHKERRQ(ierr);
254   ierr = PetscContainerSetUserDestroy(container,PetscContainerDestroy_Mat_MatMatMultTrans);CHKERRQ(ierr);
255   ierr = PetscObjectCompose((PetscObject)(*C),"Mat_MatMatMultTrans",(PetscObject)container);CHKERRQ(ierr);
256   ierr = PetscContainerDestroy(&container);CHKERRQ(ierr);
257 
258   multtrans->usecoloring = PETSC_FALSE;
259   multtrans->destroy = (*C)->ops->destroy;
260   (*C)->ops->destroy = MatDestroy_SeqAIJ_MatMatMultTrans;
261 
262   ierr = PetscOptionsGetBool(PETSC_NULL,"-matmulttrans_color",&multtrans->usecoloring,PETSC_NULL);CHKERRQ(ierr);
263   if (multtrans->usecoloring){
264     /* Create MatMultTransposeColoring from symbolic C=A*B^T */
265     MatMultTransposeColoring  matcoloring;
266     ISColoring                iscoloring;
267     Mat                       Bt_dense,C_dense;
268     PetscLogDouble            t0,tf,etime0=0.0,etime1=0.0;
269 
270     ierr = PetscGetTime(&t0);CHKERRQ(ierr);
271     ierr = MatGetColoring(*C,MATCOLORINGSL,&iscoloring);CHKERRQ(ierr);
272     ierr = MatMultTransposeColoringCreate(*C,iscoloring,&matcoloring);CHKERRQ(ierr);
273     multtrans->matcoloring = matcoloring;
274     ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr);
275     ierr = PetscGetTime(&tf);CHKERRQ(ierr);
276     etime0 += tf - t0;
277 
278     ierr = PetscGetTime(&t0);CHKERRQ(ierr);
279     /* Create Bt_dense and C_dense = A*Bt_dense */
280     ierr = MatCreate(PETSC_COMM_SELF,&Bt_dense);CHKERRQ(ierr);
281     ierr = MatSetSizes(Bt_dense,A->cmap->n,matcoloring->ncolors,A->cmap->n,matcoloring->ncolors);CHKERRQ(ierr);
282     ierr = MatSetType(Bt_dense,MATSEQDENSE);CHKERRQ(ierr);
283     ierr = MatSeqDenseSetPreallocation(Bt_dense,PETSC_NULL);CHKERRQ(ierr);
284     Bt_dense->assembled = PETSC_TRUE;
285     multtrans->Bt_den = Bt_dense;
286 
287     ierr = MatCreate(PETSC_COMM_SELF,&C_dense);CHKERRQ(ierr);
288     ierr = MatSetSizes(C_dense,A->rmap->n,matcoloring->ncolors,A->rmap->n,matcoloring->ncolors);CHKERRQ(ierr);
289     ierr = MatSetType(C_dense,MATSEQDENSE);CHKERRQ(ierr);
290     ierr = MatSeqDenseSetPreallocation(C_dense,PETSC_NULL);CHKERRQ(ierr);
291     Bt_dense->assembled = PETSC_TRUE;
292     multtrans->ABt_den = C_dense;
293     ierr = PetscGetTime(&tf);CHKERRQ(ierr);
294     etime1 += tf - t0;
295     printf("MatMultTransColorCreate %g, MatDenseCreate %g\n",etime0,etime1);
296   }
297 
298 #if defined(INEFFICIENT_ALGORITHM)
299   /* The algorithm below computes am*bm sparse inner-product - inefficient! It will be deleted later. */
300   PetscFreeSpaceList free_space=PETSC_NULL,current_space=PETSC_NULL;
301   Mat_SeqAIJ         *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
302   PetscInt           *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,*ci,*cj,*acol,*bcol;
303   PetscInt           am=A->rmap->N,bm=B->rmap->N;
304   PetscInt           i,j,anzi,bnzj,cnzi,nlnk,*lnk,nspacedouble=0,ka,kb,index[1];
305   MatScalar          *ca;
306   PetscBT            lnkbt;
307   PetscReal          afill;
308 
309   /* Allocate row pointer array ci  */
310   ierr = PetscMalloc(((am+1)+1)*sizeof(PetscInt),&ci);CHKERRQ(ierr);
311   ci[0] = 0;
312 
313   /* Create and initialize a linked list for C columns */
314   nlnk = bm+1;
315   ierr = PetscLLCreate(bm,bm,nlnk,lnk,lnkbt);CHKERRQ(ierr);
316 
317   /* Initial FreeSpace with size fill*(nnz(A)+nnz(B)) */
318   ierr = PetscFreeSpaceGet((PetscInt)(fill*(ai[am]+bi[bm])),&free_space);CHKERRQ(ierr);
319   current_space = free_space;
320 
321   /* Determine symbolic info for each row of the product A*B^T: */
322   for (i=0; i<am; i++) {
323     anzi = ai[i+1] - ai[i];
324     cnzi = 0;
325     acol = aj + ai[i];
326     for (j=0; j<bm; j++){
327       bnzj = bi[j+1] - bi[j];
328       bcol= bj + bi[j];
329       /* sparse inner-product c(i,j)=A[i,:]*B[j,:]^T */
330       ka = 0; kb = 0;
331       while (ka < anzi && kb < bnzj){
332         while (acol[ka] < bcol[kb] && ka < anzi) ka++;
333         if (ka == anzi) break;
334         while (acol[ka] > bcol[kb] && kb < bnzj) kb++;
335         if (kb == bnzj) break;
336         if (acol[ka] == bcol[kb]){ /* add nonzero c(i,j) to lnk */
337           index[0] = j;
338           ierr = PetscLLAdd(1,index,bm,nlnk,lnk,lnkbt);CHKERRQ(ierr);
339           cnzi++;
340           break;
341         }
342       }
343     }
344 
345     /* If free space is not available, make more free space */
346     /* Double the amount of total space in the list */
347     if (current_space->local_remaining<cnzi) {
348       ierr = PetscFreeSpaceGet(cnzi+current_space->total_array_size,&current_space);CHKERRQ(ierr);
349       nspacedouble++;
350     }
351 
352     /* Copy data into free space, then initialize lnk */
353     ierr = PetscLLClean(bm,bm,cnzi,lnk,current_space->array,lnkbt);CHKERRQ(ierr);
354     current_space->array           += cnzi;
355     current_space->local_used      += cnzi;
356     current_space->local_remaining -= cnzi;
357 
358     ci[i+1] = ci[i] + cnzi;
359   }
360 
361 
362   /* Column indices are in the list of free space.
363      Allocate array cj, copy column indices to cj, and destroy list of free space */
364   ierr = PetscMalloc((ci[am]+1)*sizeof(PetscInt),&cj);CHKERRQ(ierr);
365   ierr = PetscFreeSpaceContiguous(&free_space,cj);CHKERRQ(ierr);
366   ierr = PetscLLDestroy(lnk,lnkbt);CHKERRQ(ierr);
367 
368   /* Allocate space for ca */
369   ierr = PetscMalloc((ci[am]+1)*sizeof(MatScalar),&ca);CHKERRQ(ierr);
370   ierr = PetscMemzero(ca,(ci[am]+1)*sizeof(MatScalar));CHKERRQ(ierr);
371 
372   /* put together the new symbolic matrix */
373   ierr = MatCreateSeqAIJWithArrays(((PetscObject)A)->comm,am,bm,ci,cj,ca,C);CHKERRQ(ierr);
374 
375   /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
376   /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
377   c = (Mat_SeqAIJ *)((*C)->data);
378   c->free_a   = PETSC_TRUE;
379   c->free_ij  = PETSC_TRUE;
380   c->nonew    = 0;
381 
382   /* set MatInfo */
383   afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;
384   if (afill < 1.0) afill = 1.0;
385   c->maxnz                     = ci[am];
386   c->nz                        = ci[am];
387   (*C)->info.mallocs           = nspacedouble;
388   (*C)->info.fill_ratio_given  = fill;
389   (*C)->info.fill_ratio_needed = afill;
390 
391 #if defined(PETSC_USE_INFO)
392   if (ci[am]) {
393     ierr = PetscInfo3((*C),"Reallocs %D; Fill ratio: given %G needed %G.\n",nspacedouble,fill,afill);CHKERRQ(ierr);
394     ierr = PetscInfo1((*C),"Use MatMatMultTranspose(A,B,MatReuse,%G,&C) for best performance.;\n",afill);CHKERRQ(ierr);
395   } else {
396     ierr = PetscInfo((*C),"Empty matrix product\n");CHKERRQ(ierr);
397   }
398 #endif
399 #endif
400   PetscFunctionReturn(0);
401 }
402 
403 /* #define USE_ARRAY - for sparse dot product. Slower than !USE_ARRAY */
404 #undef __FUNCT__
405 #define __FUNCT__ "MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ"
406 PetscErrorCode MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C)
407 {
408   PetscErrorCode ierr;
409   Mat_SeqAIJ     *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c=(Mat_SeqAIJ*)C->data;
410   PetscInt       *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,anzi,bnzj,nexta,nextb,*acol,*bcol,brow;
411   PetscInt       cm=C->rmap->n,*ci=c->i,*cj=c->j,i,j,cnzi,*ccol;
412   PetscLogDouble flops=0.0;
413   MatScalar      *aa=a->a,*aval,*ba=b->a,*bval,*ca=c->a,*cval;
414   Mat_MatMatMultTrans *multtrans;
415   PetscContainer      container;
416 #if defined(USE_ARRAY)
417   MatScalar      *spdot;
418 #endif
419 
420   PetscFunctionBegin;
421   ierr = PetscObjectQuery((PetscObject)C,"Mat_MatMatMultTrans",(PetscObject *)&container);CHKERRQ(ierr);
422   if (!container) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Container does not exit");
423   ierr  = PetscContainerGetPointer(container,(void **)&multtrans);CHKERRQ(ierr);
424   if (multtrans->usecoloring){
425     MatMultTransposeColoring  matcoloring = multtrans->matcoloring;
426     Mat                       Bt_dense;
427     PetscInt                  m,n;
428     PetscLogDouble t0,tf,etime0=0.0,etime1=0.0,etime2=0.0;
429     Mat C_dense = multtrans->ABt_den;
430 
431     Bt_dense = multtrans->Bt_den;
432     ierr = MatGetLocalSize(Bt_dense,&m,&n);CHKERRQ(ierr);
433     printf("Bt_dense: %d,%d\n",m,n);
434 
435     /* Get Bt_dense by Apply MatMultTransposeColoring to B */
436     ierr = PetscGetTime(&t0);CHKERRQ(ierr);
437     ierr = MatMultTransposeColoringApply(B,Bt_dense,matcoloring);CHKERRQ(ierr);
438     ierr = PetscGetTime(&tf);CHKERRQ(ierr);
439     etime0 += tf - t0;
440 
441     /* C_dense = A*Bt_dense */
442     ierr = PetscGetTime(&t0);CHKERRQ(ierr);
443     ierr = MatMatMultNumeric_SeqAIJ_SeqDense(A,Bt_dense,C_dense);CHKERRQ(ierr);
444     ierr = PetscGetTime(&tf);CHKERRQ(ierr);
445     etime2 += tf - t0;
446 
447     /* Recover C from C_dense */
448     ierr = PetscGetTime(&t0);CHKERRQ(ierr);
449     ierr = MatMultTransColoringApplyDenToSp(matcoloring,C_dense,C);CHKERRQ(ierr);
450     ierr = PetscGetTime(&tf);CHKERRQ(ierr);
451     etime1 += tf - t0;
452     printf("etime ColoringApply: %g %g; MatMatMultNumeric_sp_dense: %g\n",etime0,etime1,etime2);
453     PetscFunctionReturn(0);
454   }
455 
456 #if defined(USE_ARRAY)
457   /* allocate an array for implementing sparse inner-product */
458   ierr = PetscMalloc((A->cmap->n+1)*sizeof(MatScalar),&spdot);CHKERRQ(ierr);
459   ierr = PetscMemzero(spdot,(A->cmap->n+1)*sizeof(MatScalar));CHKERRQ(ierr);
460 #endif
461 
462   /* clear old values in C */
463   ierr = PetscMemzero(ca,ci[cm]*sizeof(MatScalar));CHKERRQ(ierr);
464 
465   for (i=0; i<cm; i++) {
466     anzi = ai[i+1] - ai[i];
467     acol = aj + ai[i];
468     aval = aa + ai[i];
469     cnzi = ci[i+1] - ci[i];
470     ccol = cj + ci[i];
471     cval = ca + ci[i];
472     for (j=0; j<cnzi; j++){
473       brow = ccol[j];
474       bnzj = bi[brow+1] - bi[brow];
475       bcol = bj + bi[brow];
476       bval = ba + bi[brow];
477 
478       /* perform sparse inner-product c(i,j)=A[i,:]*B[j,:]^T */
479 #if defined(USE_ARRAY)
480       /* put ba to spdot array */
481       for (nextb=0; nextb<bnzj; nextb++) spdot[bcol[nextb]] = bval[nextb];
482       /* c(i,j)=A[i,:]*B[j,:]^T */
483       for (nexta=0; nexta<anzi; nexta++){
484         cval[j] += spdot[acol[nexta]]*aval[nexta];
485       }
486       /* zero spdot array */
487       for (nextb=0; nextb<bnzj; nextb++) spdot[bcol[nextb]] = 0.0;
488 #else
489       nexta = 0; nextb = 0;
490       while (nexta<anzi && nextb<bnzj){
491         while (acol[nexta] < bcol[nextb] && nexta < anzi) nexta++;
492         if (nexta == anzi) break;
493         while (acol[nexta] > bcol[nextb] && nextb < bnzj) nextb++;
494         if (nextb == bnzj) break;
495         if (acol[nexta] == bcol[nextb]){
496           cval[j] += aval[nexta]*bval[nextb];
497           nexta++; nextb++;
498           flops += 2;
499         }
500       }
501 #endif
502     }
503   }
504   ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
505   ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
506   ierr = PetscLogFlops(flops);CHKERRQ(ierr);
507 #if defined(USE_ARRAY)
508   ierr = PetscFree(spdot);
509 #endif
510   PetscFunctionReturn(0);
511 }
512 
513 #undef __FUNCT__
514 #define __FUNCT__ "MatMatTransposeMult_SeqAIJ_SeqAIJ"
515 PetscErrorCode MatMatTransposeMult_SeqAIJ_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) {
516   PetscErrorCode ierr;
517 
518   PetscFunctionBegin;
519   if (scall == MAT_INITIAL_MATRIX){
520     ierr = MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ(A,B,fill,C);CHKERRQ(ierr);
521   }
522   ierr = MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(A,B,*C);CHKERRQ(ierr);
523   PetscFunctionReturn(0);
524 }
525 
526 #undef __FUNCT__
527 #define __FUNCT__ "MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ"
528 PetscErrorCode MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C)
529 {
530   PetscErrorCode ierr;
531   Mat            At;
532   PetscInt       *ati,*atj;
533 
534   PetscFunctionBegin;
535   /* create symbolic At */
536   ierr = MatGetSymbolicTranspose_SeqAIJ(A,&ati,&atj);CHKERRQ(ierr);
537   ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,A->cmap->n,A->rmap->n,ati,atj,PETSC_NULL,&At);CHKERRQ(ierr);
538 
539   /* get symbolic C=At*B */
540   ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(At,B,fill,C);CHKERRQ(ierr);
541 
542   /* clean up */
543   ierr = MatDestroy(&At);CHKERRQ(ierr);
544   ierr = MatRestoreSymbolicTranspose_SeqAIJ(A,&ati,&atj);CHKERRQ(ierr);
545   PetscFunctionReturn(0);
546 }
547 
548 #undef __FUNCT__
549 #define __FUNCT__ "MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ"
550 PetscErrorCode MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C)
551 {
552   PetscErrorCode ierr;
553   Mat_SeqAIJ     *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c=(Mat_SeqAIJ*)C->data;
554   PetscInt       am=A->rmap->n,anzi,*ai=a->i,*aj=a->j,*bi=b->i,*bj,bnzi,nextb;
555   PetscInt       cm=C->rmap->n,*ci=c->i,*cj=c->j,crow,*cjj,i,j,k;
556   PetscLogDouble flops=0.0;
557   MatScalar      *aa=a->a,*ba,*ca=c->a,*caj;
558 
559   PetscFunctionBegin;
560   /* clear old values in C */
561   ierr = PetscMemzero(ca,ci[cm]*sizeof(MatScalar));CHKERRQ(ierr);
562 
563   /* compute A^T*B using outer product (A^T)[:,i]*B[i,:] */
564   for (i=0;i<am;i++) {
565     bj   = b->j + bi[i];
566     ba   = b->a + bi[i];
567     bnzi = bi[i+1] - bi[i];
568     anzi = ai[i+1] - ai[i];
569     for (j=0; j<anzi; j++) {
570       nextb = 0;
571       crow  = *aj++;
572       cjj   = cj + ci[crow];
573       caj   = ca + ci[crow];
574       /* perform sparse axpy operation.  Note cjj includes bj. */
575       for (k=0; nextb<bnzi; k++) {
576         if (cjj[k] == *(bj+nextb)) { /* ccol == bcol */
577           caj[k] += (*aa)*(*(ba+nextb));
578           nextb++;
579         }
580       }
581       flops += 2*bnzi;
582       aa++;
583     }
584   }
585 
586   /* Assemble the final matrix and clean up */
587   ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
588   ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
589   ierr = PetscLogFlops(flops);CHKERRQ(ierr);
590   PetscFunctionReturn(0);
591 }
592 
593 EXTERN_C_BEGIN
594 #undef __FUNCT__
595 #define __FUNCT__ "MatMatMult_SeqAIJ_SeqDense"
596 PetscErrorCode MatMatMult_SeqAIJ_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
597 {
598   PetscErrorCode ierr;
599 
600   PetscFunctionBegin;
601   if (scall == MAT_INITIAL_MATRIX){
602     ierr = MatMatMultSymbolic_SeqAIJ_SeqDense(A,B,fill,C);CHKERRQ(ierr);
603   }
604   ierr = MatMatMultNumeric_SeqAIJ_SeqDense(A,B,*C);CHKERRQ(ierr);
605   PetscFunctionReturn(0);
606 }
607 EXTERN_C_END
608 
609 #undef __FUNCT__
610 #define __FUNCT__ "MatMatMultSymbolic_SeqAIJ_SeqDense"
611 PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
612 {
613   PetscErrorCode ierr;
614 
615   PetscFunctionBegin;
616   ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,0.0,C);CHKERRQ(ierr);
617   PetscFunctionReturn(0);
618 }
619 
620 #undef __FUNCT__
621 #define __FUNCT__ "MatMatMultNumeric_SeqAIJ_SeqDense"
622 PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqDense(Mat A,Mat B,Mat C)
623 {
624   Mat_SeqAIJ     *a=(Mat_SeqAIJ*)A->data;
625   PetscErrorCode ierr;
626   PetscScalar    *b,*c,r1,r2,r3,r4,*b1,*b2,*b3,*b4;
627   MatScalar      *aa;
628   PetscInt       cm=C->rmap->n, cn=B->cmap->n, bm=B->rmap->n, col, i,j,n,*aj, am = A->rmap->n;
629   PetscInt       am2 = 2*am, am3 = 3*am,  bm4 = 4*bm,colam;
630 
631   PetscFunctionBegin;
632   if (!cm || !cn) PetscFunctionReturn(0);
633   if (bm != A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number columns in A %D not equal rows in B %D\n",A->cmap->n,bm);
634   if (A->rmap->n != C->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number rows in C %D not equal rows in A %D\n",C->rmap->n,A->rmap->n);
635   if (B->cmap->n != C->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number columns in B %D not equal columns in C %D\n",B->cmap->n,C->cmap->n);
636   ierr = MatGetArray(B,&b);CHKERRQ(ierr);
637   ierr = MatGetArray(C,&c);CHKERRQ(ierr);
638   b1 = b; b2 = b1 + bm; b3 = b2 + bm; b4 = b3 + bm;
639   for (col=0; col<cn-4; col += 4){  /* over columns of C */
640     colam = col*am;
641     for (i=0; i<am; i++) {        /* over rows of C in those columns */
642       r1 = r2 = r3 = r4 = 0.0;
643       n   = a->i[i+1] - a->i[i];
644       aj  = a->j + a->i[i];
645       aa  = a->a + a->i[i];
646       for (j=0; j<n; j++) {
647         r1 += (*aa)*b1[*aj];
648         r2 += (*aa)*b2[*aj];
649         r3 += (*aa)*b3[*aj];
650         r4 += (*aa++)*b4[*aj++];
651       }
652       c[colam + i]       = r1;
653       c[colam + am + i]  = r2;
654       c[colam + am2 + i] = r3;
655       c[colam + am3 + i] = r4;
656     }
657     b1 += bm4;
658     b2 += bm4;
659     b3 += bm4;
660     b4 += bm4;
661   }
662   for (;col<cn; col++){     /* over extra columns of C */
663     for (i=0; i<am; i++) {  /* over rows of C in those columns */
664       r1 = 0.0;
665       n   = a->i[i+1] - a->i[i];
666       aj  = a->j + a->i[i];
667       aa  = a->a + a->i[i];
668 
669       for (j=0; j<n; j++) {
670         r1 += (*aa++)*b1[*aj++];
671       }
672       c[col*am + i]     = r1;
673     }
674     b1 += bm;
675   }
676   ierr = PetscLogFlops(cn*(2.0*a->nz));CHKERRQ(ierr);
677   ierr = MatRestoreArray(B,&b);CHKERRQ(ierr);
678   ierr = MatRestoreArray(C,&c);CHKERRQ(ierr);
679   ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
680   ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
681   PetscFunctionReturn(0);
682 }
683 
684 /*
685    Note very similar to MatMult_SeqAIJ(), should generate both codes from same base
686 */
687 #undef __FUNCT__
688 #define __FUNCT__ "MatMatMultNumericAdd_SeqAIJ_SeqDense"
689 PetscErrorCode MatMatMultNumericAdd_SeqAIJ_SeqDense(Mat A,Mat B,Mat C)
690 {
691   Mat_SeqAIJ     *a=(Mat_SeqAIJ*)A->data;
692   PetscErrorCode ierr;
693   PetscScalar    *b,*c,r1,r2,r3,r4,*b1,*b2,*b3,*b4;
694   MatScalar      *aa;
695   PetscInt       cm=C->rmap->n, cn=B->cmap->n, bm=B->rmap->n, col, i,j,n,*aj, am = A->rmap->n,*ii,arm;
696   PetscInt       am2 = 2*am, am3 = 3*am,  bm4 = 4*bm,colam,*ridx;
697 
698   PetscFunctionBegin;
699   if (!cm || !cn) PetscFunctionReturn(0);
700   ierr = MatGetArray(B,&b);CHKERRQ(ierr);
701   ierr = MatGetArray(C,&c);CHKERRQ(ierr);
702   b1 = b; b2 = b1 + bm; b3 = b2 + bm; b4 = b3 + bm;
703 
704   if (a->compressedrow.use){ /* use compressed row format */
705     for (col=0; col<cn-4; col += 4){  /* over columns of C */
706       colam = col*am;
707       arm   = a->compressedrow.nrows;
708       ii    = a->compressedrow.i;
709       ridx  = a->compressedrow.rindex;
710       for (i=0; i<arm; i++) {        /* over rows of C in those columns */
711 	r1 = r2 = r3 = r4 = 0.0;
712 	n   = ii[i+1] - ii[i];
713 	aj  = a->j + ii[i];
714 	aa  = a->a + ii[i];
715 	for (j=0; j<n; j++) {
716 	  r1 += (*aa)*b1[*aj];
717 	  r2 += (*aa)*b2[*aj];
718 	  r3 += (*aa)*b3[*aj];
719 	  r4 += (*aa++)*b4[*aj++];
720 	}
721 	c[colam       + ridx[i]] += r1;
722 	c[colam + am  + ridx[i]] += r2;
723 	c[colam + am2 + ridx[i]] += r3;
724 	c[colam + am3 + ridx[i]] += r4;
725       }
726       b1 += bm4;
727       b2 += bm4;
728       b3 += bm4;
729       b4 += bm4;
730     }
731     for (;col<cn; col++){     /* over extra columns of C */
732       colam = col*am;
733       arm   = a->compressedrow.nrows;
734       ii    = a->compressedrow.i;
735       ridx  = a->compressedrow.rindex;
736       for (i=0; i<arm; i++) {  /* over rows of C in those columns */
737 	r1 = 0.0;
738 	n   = ii[i+1] - ii[i];
739 	aj  = a->j + ii[i];
740 	aa  = a->a + ii[i];
741 
742 	for (j=0; j<n; j++) {
743 	  r1 += (*aa++)*b1[*aj++];
744 	}
745 	c[col*am + ridx[i]] += r1;
746       }
747       b1 += bm;
748     }
749   } else {
750     for (col=0; col<cn-4; col += 4){  /* over columns of C */
751       colam = col*am;
752       for (i=0; i<am; i++) {        /* over rows of C in those columns */
753 	r1 = r2 = r3 = r4 = 0.0;
754 	n   = a->i[i+1] - a->i[i];
755 	aj  = a->j + a->i[i];
756 	aa  = a->a + a->i[i];
757 	for (j=0; j<n; j++) {
758 	  r1 += (*aa)*b1[*aj];
759 	  r2 += (*aa)*b2[*aj];
760 	  r3 += (*aa)*b3[*aj];
761 	  r4 += (*aa++)*b4[*aj++];
762 	}
763 	c[colam + i]       += r1;
764 	c[colam + am + i]  += r2;
765 	c[colam + am2 + i] += r3;
766 	c[colam + am3 + i] += r4;
767       }
768       b1 += bm4;
769       b2 += bm4;
770       b3 += bm4;
771       b4 += bm4;
772     }
773     for (;col<cn; col++){     /* over extra columns of C */
774       for (i=0; i<am; i++) {  /* over rows of C in those columns */
775 	r1 = 0.0;
776 	n   = a->i[i+1] - a->i[i];
777 	aj  = a->j + a->i[i];
778 	aa  = a->a + a->i[i];
779 
780 	for (j=0; j<n; j++) {
781 	  r1 += (*aa++)*b1[*aj++];
782 	}
783 	c[col*am + i]     += r1;
784       }
785       b1 += bm;
786     }
787   }
788   ierr = PetscLogFlops(cn*2.0*a->nz);CHKERRQ(ierr);
789   ierr = MatRestoreArray(B,&b);CHKERRQ(ierr);
790   ierr = MatRestoreArray(C,&c);CHKERRQ(ierr);
791   PetscFunctionReturn(0);
792 }
793 
794 #undef __FUNCT__
795 #define __FUNCT__ "MatMultTransposeColoringApply_SeqAIJ"
796 PetscErrorCode  MatMultTransposeColoringApply_SeqAIJ(Mat B,Mat Btdense,MatMultTransposeColoring coloring)
797 {
798   PetscErrorCode ierr;
799   Mat_SeqAIJ     *a = (Mat_SeqAIJ*)B->data;
800   Mat_SeqDense   *atdense = (Mat_SeqDense*)Btdense->data;
801   PetscInt       m=Btdense->rmap->n,n=Btdense->cmap->n,j,k,l,col,anz,*atcol,brow,bcol;
802   MatScalar      *atval,*bval;
803 
804   PetscFunctionBegin;
805   ierr = PetscMemzero(atdense->v,(m*n)*sizeof(MatScalar));CHKERRQ(ierr);
806   for (k=0; k<coloring->ncolors; k++) {
807     for (l=0; l<coloring->ncolumns[k]; l++) { /* insert a row of B to a column of Btdense */
808       col = coloring->columns[k][l];   /* =row of B */
809       anz = a->i[col+1] - a->i[col];
810       for (j=0; j<anz; j++){
811         atcol = a->j + a->i[col];
812         atval = a->a + a->i[col];
813         bval  = atdense->v;
814         brow  = atcol[j];
815         bcol  = k;
816         bval[bcol*m+brow] = atval[j];
817       }
818     }
819   }
820   PetscFunctionReturn(0);
821 }
822 
823 #undef __FUNCT__
824 #define __FUNCT__ "MatMultTransColoringApplyDenToSp_SeqAIJ"
825 PetscErrorCode MatMultTransColoringApplyDenToSp_SeqAIJ(MatMultTransposeColoring matcoloring,Mat Cden,Mat Csp)
826 {
827   PetscErrorCode ierr;
828   Mat_SeqAIJ     *csp = (Mat_SeqAIJ*)Csp->data;
829   PetscInt       k,l,row,col,m;
830   PetscScalar    *ca_den,*cp_den,*ca=csp->a;
831   PetscInt       *ci=csp->i,*cj=csp->j,*rp,low,high,*cilen = csp->ilen,t,i;
832 #if defined(PETSC_USE_DEBUG)
833   PetscBool      found;
834 #endif
835 
836   PetscFunctionBegin;
837   ierr = MatGetLocalSize(Csp,&m,PETSC_NULL);CHKERRQ(ierr);
838   ierr = MatGetArray(Cden,&ca_den);CHKERRQ(ierr);
839   cp_den = ca_den;
840   for (k=0; k<matcoloring->ncolors; k++) {
841     for (l=0; l<matcoloring->nrows[k]; l++){
842       row   = matcoloring->rows[k][l];             /* local row index */
843       col   = matcoloring->columnsforrow[k][l];    /* global column index */
844 
845       /* below is optimized from  MatSetValues(Csp,1,&row,1,&col,cval+row,INSERT_VALUES); */
846       rp   = cj + ci[row];
847       low  = 0; high = cilen[row];
848       while (high-low > 5) {
849         t = (low+high)/2;
850         if (rp[t] > col) high = t;
851         else             low  = t;
852       }
853 #if defined(PETSC_USE_DEBUG)
854       found = PETSC_FALSE;
855 #endif
856       for (i=low; i<high; i++) {
857         if (rp[i] == col) {
858           *(ca + ci[row] + i) = cp_den[row];
859 #if defined(PETSC_USE_DEBUG)
860           found = PETSC_TRUE;
861 #endif
862           break;
863         }
864       }
865 #if defined(PETSC_USE_DEBUG)
866       if (!found) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"%d %d is not inserted",row,col);
867 #endif
868     }
869     cp_den += m;
870   }
871   ierr = MatRestoreArray(Cden,&ca_den);CHKERRQ(ierr);
872   PetscFunctionReturn(0);
873 }
874 
875 #undef __FUNCT__
876 #define __FUNCT__ "MatMultTransposeColoringCreate_SeqAIJ"
877 PetscErrorCode MatMultTransposeColoringCreate_SeqAIJ(Mat mat,ISColoring iscoloring,MatMultTransposeColoring c)
878 {
879   PetscErrorCode ierr;
880   PetscInt       i,n,nrows,N,j,k,m,*rows,*ci,*cj,ncols,col,cm;
881   const PetscInt *is;
882   PetscInt       nis = iscoloring->n,*rowhit,*columnsforrow,bs = 1;
883   IS             *isa;
884   PetscBool      done;
885   PetscBool      flg1,flg2;
886 
887   PetscFunctionBegin;
888   ierr = ISColoringGetIS(iscoloring,PETSC_IGNORE,&isa);CHKERRQ(ierr);
889   /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */
890   ierr = PetscTypeCompare((PetscObject)mat,MATSEQBAIJ,&flg1);CHKERRQ(ierr);
891   ierr = PetscTypeCompare((PetscObject)mat,MATMPIBAIJ,&flg2);CHKERRQ(ierr);
892   if (flg1 || flg2) {
893     ierr = MatGetBlockSize(mat,&bs);CHKERRQ(ierr);
894   }
895 
896   N          = mat->cmap->N/bs;
897   c->M       = mat->rmap->N/bs;  /* set total rows, columns and local rows */
898   c->N       = mat->cmap->N/bs;
899   c->m       = mat->rmap->N/bs;
900   c->rstart  = 0;
901 
902   c->ncolors = nis;
903   cm         = c->m;
904   ierr       = PetscMalloc(nis*sizeof(PetscInt),&c->ncolumns);CHKERRQ(ierr);
905   ierr       = PetscMalloc(nis*sizeof(PetscInt*),&c->columns);CHKERRQ(ierr);
906   ierr       = PetscMalloc(cm*sizeof(PetscInt),&c->nrows);CHKERRQ(ierr);
907   ierr       = PetscMalloc(cm*sizeof(PetscInt*),&c->rows);CHKERRQ(ierr);
908   ierr       = PetscMalloc(cm*sizeof(PetscInt*),&c->columnsforrow);CHKERRQ(ierr);
909 
910   ierr = MatGetColumnIJ(mat,0,PETSC_FALSE,PETSC_FALSE,&ncols,&ci,&cj,&done);CHKERRQ(ierr);
911   if (!done) SETERRQ1(((PetscObject)mat)->comm,PETSC_ERR_SUP,"MatGetColumnIJ() not supported for matrix type %s",((PetscObject)mat)->type_name);
912 
913   ierr = PetscMalloc2(cm+1,PetscInt,&rowhit,cm+1,PetscInt,&columnsforrow);CHKERRQ(ierr);
914   for (i=0; i<nis; i++) {
915     ierr = ISGetLocalSize(isa[i],&n);CHKERRQ(ierr);
916     ierr = ISGetIndices(isa[i],&is);CHKERRQ(ierr);
917     c->ncolumns[i] = n;
918     if (n) {
919       ierr = PetscMalloc(n*sizeof(PetscInt),&c->columns[i]);CHKERRQ(ierr);
920       ierr = PetscMemcpy(c->columns[i],is,n*sizeof(PetscInt));CHKERRQ(ierr);
921     } else {
922       c->columns[i]  = 0;
923     }
924 
925     /* fast, crude version requires O(N*N) work */
926     ierr = PetscMemzero(rowhit,cm*sizeof(PetscInt));CHKERRQ(ierr);
927     /* loop over columns*/
928     for (j=0; j<n; j++) {
929       col  = is[j];
930       rows = cj + ci[col];
931       m    = ci[col+1] - ci[col];
932       /* loop over columns marking them in rowhit */
933       for (k=0; k<m; k++) {
934         rowhit[*rows++] = col + 1;
935       }
936     }
937     /* count the number of hits */
938     nrows = 0;
939     for (j=0; j<cm; j++) {
940       if (rowhit[j]) nrows++;
941     }
942     c->nrows[i] = nrows;
943     ierr        = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->rows[i]);CHKERRQ(ierr);
944     ierr        = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->columnsforrow[i]);CHKERRQ(ierr);
945     nrows       = 0;
946     for (j=0; j<cm; j++) {
947       if (rowhit[j]) {
948         c->rows[i][nrows]          = j;
949         c->columnsforrow[i][nrows] = rowhit[j] - 1;
950         nrows++;
951       }
952     }
953     ierr = ISRestoreIndices(isa[i],&is);CHKERRQ(ierr);
954   }
955   ierr = MatRestoreColumnIJ(mat,0,PETSC_FALSE,PETSC_FALSE,&ncols,&ci,&cj,&done);CHKERRQ(ierr);
956 
957   ierr = PetscFree2(rowhit,columnsforrow);CHKERRQ(ierr);
958   ierr = ISColoringRestoreIS(iscoloring,&isa);CHKERRQ(ierr);
959   PetscFunctionReturn(0);
960 }
961