xref: /petsc/src/mat/impls/aij/seq/matmatmult.c (revision 3116ef58e11bc77b08297b681ddd79bc46c3f215) 
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__ "MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ"
193 PetscErrorCode MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C)
194 {
195   PetscErrorCode ierr;
196   PetscFreeSpaceList free_space=PETSC_NULL,current_space=PETSC_NULL;
197   Mat_SeqAIJ         *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
198   PetscInt           *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,*ci,*cj,*acol,*bcol;
199   PetscInt           am=A->rmap->N,bm=B->rmap->N;
200   PetscInt           i,j,anzi,bnzj,cnzi,nlnk,*lnk,nspacedouble=0,ka,kb,index[1];
201   MatScalar          *ca;
202   PetscBT            lnkbt;
203   PetscReal          afill;
204 
205   PetscFunctionBegin;
206   /* Allocate row pointer array ci  */
207   ierr = PetscMalloc(((am+1)+1)*sizeof(PetscInt),&ci);CHKERRQ(ierr);
208   ci[0] = 0;
209 
210   /* Create and initialize a linked list for C columns */
211   nlnk = bm+1;
212   ierr = PetscLLCreate(bm,bm,nlnk,lnk,lnkbt);CHKERRQ(ierr);
213 
214   /* Initial FreeSpace with size fill*(nnz(A)+nnz(B)) */
215   ierr = PetscFreeSpaceGet((PetscInt)(fill*(ai[am]+bi[bm])),&free_space);CHKERRQ(ierr);
216   current_space = free_space;
217 
218   /* Determine symbolic info for each row of the product A*B^T: */
219   for (i=0; i<am; i++) {
220     anzi = ai[i+1] - ai[i];
221     cnzi = 0;
222     acol = aj + ai[i];
223     for (j=0; j<bm; j++){
224       bnzj = bi[j+1] - bi[j];
225       bcol= bj + bi[j];
226       /* sparse inner-product A[i,:] * B[j,:]^T */
227       ka = 0; kb = 0;
228       while (ka < anzi && kb < bnzj){
229         while (acol[ka] < bcol[kb] && ka < anzi){
230           ka++;
231         }
232         while (acol[ka] > bcol[kb] && kb < bnzj){
233           kb++;
234         }
235         if (acol[ka] == bcol[kb]){ /* add nonzero c(i,j) to lnk */
236           index[0] = j;
237           ierr = PetscLLAdd(1,index,bm,nlnk,lnk,lnkbt);CHKERRQ(ierr);
238           cnzi++;
239           break;
240         }
241       }
242     }
243 
244     /* If free space is not available, make more free space */
245     /* Double the amount of total space in the list */
246     if (current_space->local_remaining<cnzi) {
247       ierr = PetscFreeSpaceGet(cnzi+current_space->total_array_size,&current_space);CHKERRQ(ierr);
248       nspacedouble++;
249     }
250 
251     /* Copy data into free space, then initialize lnk */
252     ierr = PetscLLClean(bm,bm,cnzi,lnk,current_space->array,lnkbt);CHKERRQ(ierr);
253     current_space->array           += cnzi;
254     current_space->local_used      += cnzi;
255     current_space->local_remaining -= cnzi;
256 
257     ci[i+1] = ci[i] + cnzi;
258   }
259 
260 
261   /* Column indices are in the list of free space.
262      Allocate array cj, copy column indices to cj, and destroy list of free space */
263   ierr = PetscMalloc((ci[am]+1)*sizeof(PetscInt),&cj);CHKERRQ(ierr);
264   ierr = PetscFreeSpaceContiguous(&free_space,cj);CHKERRQ(ierr);
265   ierr = PetscLLDestroy(lnk,lnkbt);CHKERRQ(ierr);
266 
267   /* Allocate space for ca */
268   ierr = PetscMalloc((ci[am]+1)*sizeof(MatScalar),&ca);CHKERRQ(ierr);
269   ierr = PetscMemzero(ca,(ci[am]+1)*sizeof(MatScalar));CHKERRQ(ierr);
270 
271   /* put together the new symbolic matrix */
272   ierr = MatCreateSeqAIJWithArrays(((PetscObject)A)->comm,am,bm,ci,cj,ca,C);CHKERRQ(ierr);
273 
274   /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
275   /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
276   c = (Mat_SeqAIJ *)((*C)->data);
277   c->free_a   = PETSC_TRUE;
278   c->free_ij  = PETSC_TRUE;
279   c->nonew    = 0;
280 
281   /* set MatInfo */
282   afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;
283   if (afill < 1.0) afill = 1.0;
284   c->maxnz                     = ci[am];
285   c->nz                        = ci[am];
286   (*C)->info.mallocs           = nspacedouble;
287   (*C)->info.fill_ratio_given  = fill;
288   (*C)->info.fill_ratio_needed = afill;
289 
290 #if defined(PETSC_USE_INFO)
291   if (ci[am]) {
292     ierr = PetscInfo3((*C),"Reallocs %D; Fill ratio: given %G needed %G.\n",nspacedouble,fill,afill);CHKERRQ(ierr);
293     ierr = PetscInfo1((*C),"Use MatMatMultTranspose(A,B,MatReuse,%G,&C) for best performance.;\n",afill);CHKERRQ(ierr);
294   } else {
295     ierr = PetscInfo((*C),"Empty matrix product\n");CHKERRQ(ierr);
296   }
297 #endif
298   PetscFunctionReturn(0);
299 }
300 
301 #undef __FUNCT__
302 #define __FUNCT__ "MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ"
303 PetscErrorCode MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C)
304 {
305 #if defined(TMP)
306   PetscErrorCode ierr;
307   Mat_SeqAIJ     *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c=(Mat_SeqAIJ*)C->data;
308   PetscInt       am=A->rmap->n,anzi,*ai=a->i,*aj=a->j,*bi=b->i,*bj,bnzi,nextb;
309   PetscInt       cm=C->rmap->n,*ci=c->i,*cj=c->j,i,j,k,cnzi,*ccol;
310   PetscLogDouble flops=0.0;
311   MatScalar      *aa=a->a,*ba,*ca=c->a,*caj;
312 #endif
313 
314   PetscFunctionBegin;
315 #if defined(TMP)
316   ierr = MatView(A, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
317   ierr = MatView(B, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
318 
319   for (i=0; i<cm; i++) {
320     printf("C row %d\n",i);
321     cnzi = ci[i+1] - ci[i];
322     ccol = cj + ci[i];
323     for (j=0; j<cnzi; j++){
324       printf(" %d, ",ccol[j]);
325     }
326     printf("   cnzi %d\n",cnzi);
327   }
328 #endif
329   SETERRQ(PETSC_COMM_SELF,0,"MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ Not done yet");
330   PetscFunctionReturn(0);
331 }
332 
333 #undef __FUNCT__
334 #define __FUNCT__ "MatMatTransposeMult_SeqAIJ_SeqAIJ"
335 PetscErrorCode MatMatTransposeMult_SeqAIJ_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) {
336   PetscErrorCode ierr;
337 
338   PetscFunctionBegin;
339   if (scall == MAT_INITIAL_MATRIX){
340     ierr = MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ(A,B,fill,C);CHKERRQ(ierr);
341   }
342   ierr = MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(A,B,*C);CHKERRQ(ierr);
343   PetscFunctionReturn(0);
344 }
345 
346 #undef __FUNCT__
347 #define __FUNCT__ "MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ"
348 PetscErrorCode MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C)
349 {
350   PetscErrorCode ierr;
351   Mat            At;
352   PetscInt       *ati,*atj;
353 
354   PetscFunctionBegin;
355   /* create symbolic At */
356   ierr = MatGetSymbolicTranspose_SeqAIJ(A,&ati,&atj);CHKERRQ(ierr);
357   ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,A->cmap->n,A->rmap->n,ati,atj,PETSC_NULL,&At);CHKERRQ(ierr);
358 
359   /* get symbolic C=At*B */
360   ierr = MatMatMultSymbolic_SeqAIJ_SeqAIJ(At,B,fill,C);CHKERRQ(ierr);
361 
362   /* clean up */
363   ierr = MatDestroy(&At);CHKERRQ(ierr);
364   ierr = MatRestoreSymbolicTranspose_SeqAIJ(A,&ati,&atj);CHKERRQ(ierr);
365   PetscFunctionReturn(0);
366 }
367 
368 #undef __FUNCT__
369 #define __FUNCT__ "MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ"
370 PetscErrorCode MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C)
371 {
372   PetscErrorCode ierr;
373   Mat_SeqAIJ     *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c=(Mat_SeqAIJ*)C->data;
374   PetscInt       am=A->rmap->n,anzi,*ai=a->i,*aj=a->j,*bi=b->i,*bj,bnzi,nextb;
375   PetscInt       cm=C->rmap->n,*ci=c->i,*cj=c->j,crow,*cjj,i,j,k;
376   PetscLogDouble flops=0.0;
377   MatScalar      *aa=a->a,*ba,*ca=c->a,*caj;
378 
379   PetscFunctionBegin;
380   /* clear old values in C */
381   ierr = PetscMemzero(ca,ci[cm]*sizeof(MatScalar));CHKERRQ(ierr);
382 
383   /* compute A^T*B using outer product (A^T)[:,i]*B[i,:] */
384   for (i=0;i<am;i++) {
385     bj   = b->j + bi[i];
386     ba   = b->a + bi[i];
387     bnzi = bi[i+1] - bi[i];
388     anzi = ai[i+1] - ai[i];
389     for (j=0; j<anzi; j++) {
390       nextb = 0;
391       crow  = *aj++;
392       cjj   = cj + ci[crow];
393       caj   = ca + ci[crow];
394       /* perform sparse axpy operation.  Note cjj includes bj. */
395       for (k=0; nextb<bnzi; k++) {
396         if (cjj[k] == *(bj+nextb)) { /* ccol == bcol */
397           caj[k] += (*aa)*(*(ba+nextb));
398           nextb++;
399         }
400       }
401       flops += 2*bnzi;
402       aa++;
403     }
404   }
405 
406   /* Assemble the final matrix and clean up */
407   ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
408   ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
409   ierr = PetscLogFlops(flops);CHKERRQ(ierr);
410   PetscFunctionReturn(0);
411 }
412 
413 EXTERN_C_BEGIN
414 #undef __FUNCT__
415 #define __FUNCT__ "MatMatMult_SeqAIJ_SeqDense"
416 PetscErrorCode MatMatMult_SeqAIJ_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
417 {
418   PetscErrorCode ierr;
419 
420   PetscFunctionBegin;
421   if (scall == MAT_INITIAL_MATRIX){
422     ierr = MatMatMultSymbolic_SeqAIJ_SeqDense(A,B,fill,C);CHKERRQ(ierr);
423   }
424   ierr = MatMatMultNumeric_SeqAIJ_SeqDense(A,B,*C);CHKERRQ(ierr);
425   PetscFunctionReturn(0);
426 }
427 EXTERN_C_END
428 
429 #undef __FUNCT__
430 #define __FUNCT__ "MatMatMultSymbolic_SeqAIJ_SeqDense"
431 PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
432 {
433   PetscErrorCode ierr;
434 
435   PetscFunctionBegin;
436   ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,0.0,C);CHKERRQ(ierr);
437   PetscFunctionReturn(0);
438 }
439 
440 #undef __FUNCT__
441 #define __FUNCT__ "MatMatMultNumeric_SeqAIJ_SeqDense"
442 PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqDense(Mat A,Mat B,Mat C)
443 {
444   Mat_SeqAIJ     *a=(Mat_SeqAIJ*)A->data;
445   PetscErrorCode ierr;
446   PetscScalar    *b,*c,r1,r2,r3,r4,*b1,*b2,*b3,*b4;
447   MatScalar      *aa;
448   PetscInt       cm=C->rmap->n, cn=B->cmap->n, bm=B->rmap->n, col, i,j,n,*aj, am = A->rmap->n;
449   PetscInt       am2 = 2*am, am3 = 3*am,  bm4 = 4*bm,colam;
450 
451   PetscFunctionBegin;
452   if (!cm || !cn) PetscFunctionReturn(0);
453   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);
454   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);
455   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);
456   ierr = MatGetArray(B,&b);CHKERRQ(ierr);
457   ierr = MatGetArray(C,&c);CHKERRQ(ierr);
458   b1 = b; b2 = b1 + bm; b3 = b2 + bm; b4 = b3 + bm;
459   for (col=0; col<cn-4; col += 4){  /* over columns of C */
460     colam = col*am;
461     for (i=0; i<am; i++) {        /* over rows of C in those columns */
462       r1 = r2 = r3 = r4 = 0.0;
463       n   = a->i[i+1] - a->i[i];
464       aj  = a->j + a->i[i];
465       aa  = a->a + a->i[i];
466       for (j=0; j<n; j++) {
467         r1 += (*aa)*b1[*aj];
468         r2 += (*aa)*b2[*aj];
469         r3 += (*aa)*b3[*aj];
470         r4 += (*aa++)*b4[*aj++];
471       }
472       c[colam + i]       = r1;
473       c[colam + am + i]  = r2;
474       c[colam + am2 + i] = r3;
475       c[colam + am3 + i] = r4;
476     }
477     b1 += bm4;
478     b2 += bm4;
479     b3 += bm4;
480     b4 += bm4;
481   }
482   for (;col<cn; col++){     /* over extra columns of C */
483     for (i=0; i<am; i++) {  /* over rows of C in those columns */
484       r1 = 0.0;
485       n   = a->i[i+1] - a->i[i];
486       aj  = a->j + a->i[i];
487       aa  = a->a + a->i[i];
488 
489       for (j=0; j<n; j++) {
490         r1 += (*aa++)*b1[*aj++];
491       }
492       c[col*am + i]     = r1;
493     }
494     b1 += bm;
495   }
496   ierr = PetscLogFlops(cn*(2.0*a->nz));CHKERRQ(ierr);
497   ierr = MatRestoreArray(B,&b);CHKERRQ(ierr);
498   ierr = MatRestoreArray(C,&c);CHKERRQ(ierr);
499   ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
500   ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
501   PetscFunctionReturn(0);
502 }
503 
504 /*
505    Note very similar to MatMult_SeqAIJ(), should generate both codes from same base
506 */
507 #undef __FUNCT__
508 #define __FUNCT__ "MatMatMultNumericAdd_SeqAIJ_SeqDense"
509 PetscErrorCode MatMatMultNumericAdd_SeqAIJ_SeqDense(Mat A,Mat B,Mat C)
510 {
511   Mat_SeqAIJ     *a=(Mat_SeqAIJ*)A->data;
512   PetscErrorCode ierr;
513   PetscScalar    *b,*c,r1,r2,r3,r4,*b1,*b2,*b3,*b4;
514   MatScalar      *aa;
515   PetscInt       cm=C->rmap->n, cn=B->cmap->n, bm=B->rmap->n, col, i,j,n,*aj, am = A->rmap->n,*ii,arm;
516   PetscInt       am2 = 2*am, am3 = 3*am,  bm4 = 4*bm,colam,*ridx;
517 
518   PetscFunctionBegin;
519   if (!cm || !cn) PetscFunctionReturn(0);
520   ierr = MatGetArray(B,&b);CHKERRQ(ierr);
521   ierr = MatGetArray(C,&c);CHKERRQ(ierr);
522   b1 = b; b2 = b1 + bm; b3 = b2 + bm; b4 = b3 + bm;
523 
524   if (a->compressedrow.use){ /* use compressed row format */
525     for (col=0; col<cn-4; col += 4){  /* over columns of C */
526       colam = col*am;
527       arm   = a->compressedrow.nrows;
528       ii    = a->compressedrow.i;
529       ridx  = a->compressedrow.rindex;
530       for (i=0; i<arm; i++) {        /* over rows of C in those columns */
531 	r1 = r2 = r3 = r4 = 0.0;
532 	n   = ii[i+1] - ii[i];
533 	aj  = a->j + ii[i];
534 	aa  = a->a + ii[i];
535 	for (j=0; j<n; j++) {
536 	  r1 += (*aa)*b1[*aj];
537 	  r2 += (*aa)*b2[*aj];
538 	  r3 += (*aa)*b3[*aj];
539 	  r4 += (*aa++)*b4[*aj++];
540 	}
541 	c[colam       + ridx[i]] += r1;
542 	c[colam + am  + ridx[i]] += r2;
543 	c[colam + am2 + ridx[i]] += r3;
544 	c[colam + am3 + ridx[i]] += r4;
545       }
546       b1 += bm4;
547       b2 += bm4;
548       b3 += bm4;
549       b4 += bm4;
550     }
551     for (;col<cn; col++){     /* over extra columns of C */
552       colam = col*am;
553       arm   = a->compressedrow.nrows;
554       ii    = a->compressedrow.i;
555       ridx  = a->compressedrow.rindex;
556       for (i=0; i<arm; i++) {  /* over rows of C in those columns */
557 	r1 = 0.0;
558 	n   = ii[i+1] - ii[i];
559 	aj  = a->j + ii[i];
560 	aa  = a->a + ii[i];
561 
562 	for (j=0; j<n; j++) {
563 	  r1 += (*aa++)*b1[*aj++];
564 	}
565 	c[col*am + ridx[i]] += r1;
566       }
567       b1 += bm;
568     }
569   } else {
570     for (col=0; col<cn-4; col += 4){  /* over columns of C */
571       colam = col*am;
572       for (i=0; i<am; i++) {        /* over rows of C in those columns */
573 	r1 = r2 = r3 = r4 = 0.0;
574 	n   = a->i[i+1] - a->i[i];
575 	aj  = a->j + a->i[i];
576 	aa  = a->a + a->i[i];
577 	for (j=0; j<n; j++) {
578 	  r1 += (*aa)*b1[*aj];
579 	  r2 += (*aa)*b2[*aj];
580 	  r3 += (*aa)*b3[*aj];
581 	  r4 += (*aa++)*b4[*aj++];
582 	}
583 	c[colam + i]       += r1;
584 	c[colam + am + i]  += r2;
585 	c[colam + am2 + i] += r3;
586 	c[colam + am3 + i] += r4;
587       }
588       b1 += bm4;
589       b2 += bm4;
590       b3 += bm4;
591       b4 += bm4;
592     }
593     for (;col<cn; col++){     /* over extra columns of C */
594       for (i=0; i<am; i++) {  /* over rows of C in those columns */
595 	r1 = 0.0;
596 	n   = a->i[i+1] - a->i[i];
597 	aj  = a->j + a->i[i];
598 	aa  = a->a + a->i[i];
599 
600 	for (j=0; j<n; j++) {
601 	  r1 += (*aa++)*b1[*aj++];
602 	}
603 	c[col*am + i]     += r1;
604       }
605       b1 += bm;
606     }
607   }
608   ierr = PetscLogFlops(cn*2.0*a->nz);CHKERRQ(ierr);
609   ierr = MatRestoreArray(B,&b);CHKERRQ(ierr);
610   ierr = MatRestoreArray(C,&c);CHKERRQ(ierr);
611   PetscFunctionReturn(0);
612 }
613