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