Actual source code: matmatmult.c
petsc-3.11.4 2019-09-28
2: /*
3: Defines matrix-matrix product routines for pairs of SeqAIJ matrices
4: C = A * B
5: */
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 <petsc/private/isimpl.h>
11: #include <../src/mat/impls/dense/seq/dense.h>
14: PETSC_INTERN PetscErrorCode MatMatMult_SeqAIJ_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
15: {
19: if (scall == MAT_INITIAL_MATRIX) {
20: PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);
21: MatMatMultSymbolic_SeqAIJ_SeqAIJ(A,B,fill,C);
22: PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);
23: }
25: PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);
26: MatMatMultNumeric_SeqAIJ_SeqAIJ(A,B,*C);
27: PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);
28: return(0);
29: }
31: PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C)
32: {
36: if (C->ops->matmultnumeric) {
37: (*C->ops->matmultnumeric)(A,B,C);
38: } else {
39: MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted(A,B,C);
40: }
41: return(0);
42: }
44: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C)
45: {
47: #if !defined(PETSC_HAVE_HYPRE)
48: const char *algTypes[8] = {"sorted","scalable","scalable_fast","heap","btheap","llcondensed","combined","rowmerge"};
49: PetscInt nalg = 8;
50: #else
51: const char *algTypes[9] = {"sorted","scalable","scalable_fast","heap","btheap","llcondensed","combined","rowmerge","hypre"};
52: PetscInt nalg = 9;
53: #endif
54: PetscInt alg = 0; /* set default algorithm */
57: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MatMatMult","Mat");
58: PetscOptionsEList("-matmatmult_via","Algorithmic approach","MatMatMult",algTypes,nalg,algTypes[0],&alg,NULL);
59: PetscOptionsEnd();
60: switch (alg) {
61: case 1:
62: MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable(A,B,fill,C);
63: break;
64: case 2:
65: MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable_fast(A,B,fill,C);
66: break;
67: case 3:
68: MatMatMultSymbolic_SeqAIJ_SeqAIJ_Heap(A,B,fill,C);
69: break;
70: case 4:
71: MatMatMultSymbolic_SeqAIJ_SeqAIJ_BTHeap(A,B,fill,C);
72: break;
73: case 5:
74: MatMatMultSymbolic_SeqAIJ_SeqAIJ_LLCondensed(A,B,fill,C);
75: break;
76: case 6:
77: MatMatMult_SeqAIJ_SeqAIJ_Combined(A,B,fill,C);
78: break;
79: case 7:
80: MatMatMultSymbolic_SeqAIJ_SeqAIJ_RowMerge(A,B,fill,C);
81: break;
82: #if defined(PETSC_HAVE_HYPRE)
83: case 8:
84: MatMatMultSymbolic_AIJ_AIJ_wHYPRE(A,B,fill,C);
85: break;
86: #endif
87: default:
88: MatMatMultSymbolic_SeqAIJ_SeqAIJ_Sorted(A,B,fill,C);
89: break;
90: }
91: return(0);
92: }
94: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_LLCondensed(Mat A,Mat B,PetscReal fill,Mat *C)
95: {
96: PetscErrorCode ierr;
97: Mat_SeqAIJ *a =(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
98: PetscInt *ai=a->i,*bi=b->i,*ci,*cj;
99: PetscInt am =A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;
100: PetscReal afill;
101: PetscInt i,j,anzi,brow,bnzj,cnzi,*bj,*aj,*lnk,ndouble=0,Crmax;
102: PetscTable ta;
103: PetscBT lnkbt;
104: PetscFreeSpaceList free_space=NULL,current_space=NULL;
107: /* Get ci and cj */
108: /*---------------*/
109: /* Allocate ci array, arrays for fill computation and */
110: /* free space for accumulating nonzero column info */
111: PetscMalloc1(am+2,&ci);
112: ci[0] = 0;
114: /* create and initialize a linked list */
115: PetscTableCreate(bn,bn,&ta);
116: MatRowMergeMax_SeqAIJ(b,bm,ta);
117: PetscTableGetCount(ta,&Crmax);
118: PetscTableDestroy(&ta);
120: PetscLLCondensedCreate(Crmax,bn,&lnk,&lnkbt);
122: /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */
123: PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ai[am],bi[bm])),&free_space);
125: current_space = free_space;
127: /* Determine ci and cj */
128: for (i=0; i<am; i++) {
129: anzi = ai[i+1] - ai[i];
130: aj = a->j + ai[i];
131: for (j=0; j<anzi; j++) {
132: brow = aj[j];
133: bnzj = bi[brow+1] - bi[brow];
134: bj = b->j + bi[brow];
135: /* add non-zero cols of B into the sorted linked list lnk */
136: PetscLLCondensedAddSorted(bnzj,bj,lnk,lnkbt);
137: }
138: cnzi = lnk[0];
140: /* If free space is not available, make more free space */
141: /* Double the amount of total space in the list */
142: if (current_space->local_remaining<cnzi) {
143: PetscFreeSpaceGet(PetscIntSumTruncate(cnzi,current_space->total_array_size),¤t_space);
144: ndouble++;
145: }
147: /* Copy data into free space, then initialize lnk */
148: PetscLLCondensedClean(bn,cnzi,current_space->array,lnk,lnkbt);
150: current_space->array += cnzi;
151: current_space->local_used += cnzi;
152: current_space->local_remaining -= cnzi;
154: ci[i+1] = ci[i] + cnzi;
155: }
157: /* Column indices are in the list of free space */
158: /* Allocate space for cj, initialize cj, and */
159: /* destroy list of free space and other temporary array(s) */
160: PetscMalloc1(ci[am]+1,&cj);
161: PetscFreeSpaceContiguous(&free_space,cj);
162: PetscLLCondensedDestroy(lnk,lnkbt);
164: /* put together the new symbolic matrix */
165: MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),am,bn,ci,cj,NULL,C);
166: MatSetBlockSizesFromMats(*C,A,B);
167: MatSetType(*C,((PetscObject)A)->type_name);
169: /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
170: /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
171: c = (Mat_SeqAIJ*)((*C)->data);
172: c->free_a = PETSC_FALSE;
173: c->free_ij = PETSC_TRUE;
174: c->nonew = 0;
175: (*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted; /* fast, needs non-scalable O(bn) array 'abdense' */
177: /* set MatInfo */
178: afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;
179: if (afill < 1.0) afill = 1.0;
180: c->maxnz = ci[am];
181: c->nz = ci[am];
182: (*C)->info.mallocs = ndouble;
183: (*C)->info.fill_ratio_given = fill;
184: (*C)->info.fill_ratio_needed = afill;
186: #if defined(PETSC_USE_INFO)
187: if (ci[am]) {
188: PetscInfo3((*C),"Reallocs %D; Fill ratio: given %g needed %g.\n",ndouble,(double)fill,(double)afill);
189: PetscInfo1((*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n",(double)afill);
190: } else {
191: PetscInfo((*C),"Empty matrix product\n");
192: }
193: #endif
194: return(0);
195: }
197: PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted(Mat A,Mat B,Mat C)
198: {
200: PetscLogDouble flops=0.0;
201: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
202: Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data;
203: Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data;
204: PetscInt *ai =a->i,*aj=a->j,*bi=b->i,*bj=b->j,*bjj,*ci=c->i,*cj=c->j;
205: PetscInt am =A->rmap->n,cm=C->rmap->n;
206: PetscInt i,j,k,anzi,bnzi,cnzi,brow;
207: PetscScalar *aa=a->a,*ba=b->a,*baj,*ca,valtmp;
208: PetscScalar *ab_dense;
211: if (!c->a) { /* first call of MatMatMultNumeric_SeqAIJ_SeqAIJ, allocate ca and matmult_abdense */
212: PetscMalloc1(ci[cm]+1,&ca);
213: c->a = ca;
214: c->free_a = PETSC_TRUE;
215: } else {
216: ca = c->a;
217: }
218: if (!c->matmult_abdense) {
219: PetscCalloc1(B->cmap->N,&ab_dense);
220: c->matmult_abdense = ab_dense;
221: } else {
222: ab_dense = c->matmult_abdense;
223: }
225: /* clean old values in C */
226: PetscMemzero(ca,ci[cm]*sizeof(MatScalar));
227: /* Traverse A row-wise. */
228: /* Build the ith row in C by summing over nonzero columns in A, */
229: /* the rows of B corresponding to nonzeros of A. */
230: for (i=0; i<am; i++) {
231: anzi = ai[i+1] - ai[i];
232: for (j=0; j<anzi; j++) {
233: brow = aj[j];
234: bnzi = bi[brow+1] - bi[brow];
235: bjj = bj + bi[brow];
236: baj = ba + bi[brow];
237: /* perform dense axpy */
238: valtmp = aa[j];
239: for (k=0; k<bnzi; k++) {
240: ab_dense[bjj[k]] += valtmp*baj[k];
241: }
242: flops += 2*bnzi;
243: }
244: aj += anzi; aa += anzi;
246: cnzi = ci[i+1] - ci[i];
247: for (k=0; k<cnzi; k++) {
248: ca[k] += ab_dense[cj[k]];
249: ab_dense[cj[k]] = 0.0; /* zero ab_dense */
250: }
251: flops += cnzi;
252: cj += cnzi; ca += cnzi;
253: }
254: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
255: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
256: PetscLogFlops(flops);
257: return(0);
258: }
260: PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable(Mat A,Mat B,Mat C)
261: {
263: PetscLogDouble flops=0.0;
264: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
265: Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data;
266: Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data;
267: PetscInt *ai = a->i,*aj=a->j,*bi=b->i,*bj=b->j,*bjj,*ci=c->i,*cj=c->j;
268: PetscInt am = A->rmap->N,cm=C->rmap->N;
269: PetscInt i,j,k,anzi,bnzi,cnzi,brow;
270: PetscScalar *aa=a->a,*ba=b->a,*baj,*ca=c->a,valtmp;
271: PetscInt nextb;
274: if (!ca) { /* first call of MatMatMultNumeric_SeqAIJ_SeqAIJ, allocate ca and matmult_abdense */
275: PetscMalloc1(ci[cm]+1,&ca);
276: c->a = ca;
277: c->free_a = PETSC_TRUE;
278: }
280: /* clean old values in C */
281: PetscMemzero(ca,ci[cm]*sizeof(MatScalar));
282: /* Traverse A row-wise. */
283: /* Build the ith row in C by summing over nonzero columns in A, */
284: /* the rows of B corresponding to nonzeros of A. */
285: for (i=0; i<am; i++) {
286: anzi = ai[i+1] - ai[i];
287: cnzi = ci[i+1] - ci[i];
288: for (j=0; j<anzi; j++) {
289: brow = aj[j];
290: bnzi = bi[brow+1] - bi[brow];
291: bjj = bj + bi[brow];
292: baj = ba + bi[brow];
293: /* perform sparse axpy */
294: valtmp = aa[j];
295: nextb = 0;
296: for (k=0; nextb<bnzi; k++) {
297: if (cj[k] == bjj[nextb]) { /* ccol == bcol */
298: ca[k] += valtmp*baj[nextb++];
299: }
300: }
301: flops += 2*bnzi;
302: }
303: aj += anzi; aa += anzi;
304: cj += cnzi; ca += cnzi;
305: }
307: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
308: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
309: PetscLogFlops(flops);
310: return(0);
311: }
313: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable_fast(Mat A,Mat B,PetscReal fill,Mat *C)
314: {
315: PetscErrorCode ierr;
316: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
317: PetscInt *ai = a->i,*bi=b->i,*ci,*cj;
318: PetscInt am = A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;
319: MatScalar *ca;
320: PetscReal afill;
321: PetscInt i,j,anzi,brow,bnzj,cnzi,*bj,*aj,*lnk,ndouble=0,Crmax;
322: PetscTable ta;
323: PetscFreeSpaceList free_space=NULL,current_space=NULL;
326: /* Get ci and cj - same as MatMatMultSymbolic_SeqAIJ_SeqAIJ except using PetscLLxxx_fast() */
327: /*-----------------------------------------------------------------------------------------*/
328: /* Allocate arrays for fill computation and free space for accumulating nonzero column */
329: PetscMalloc1(am+2,&ci);
330: ci[0] = 0;
332: /* create and initialize a linked list */
333: PetscTableCreate(bn,bn,&ta);
334: MatRowMergeMax_SeqAIJ(b,bm,ta);
335: PetscTableGetCount(ta,&Crmax);
336: PetscTableDestroy(&ta);
338: PetscLLCondensedCreate_fast(Crmax,&lnk);
340: /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */
341: PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ai[am],bi[bm])),&free_space);
342: current_space = free_space;
344: /* Determine ci and cj */
345: for (i=0; i<am; i++) {
346: anzi = ai[i+1] - ai[i];
347: aj = a->j + ai[i];
348: for (j=0; j<anzi; j++) {
349: brow = aj[j];
350: bnzj = bi[brow+1] - bi[brow];
351: bj = b->j + bi[brow];
352: /* add non-zero cols of B into the sorted linked list lnk */
353: PetscLLCondensedAddSorted_fast(bnzj,bj,lnk);
354: }
355: cnzi = lnk[1];
357: /* If free space is not available, make more free space */
358: /* Double the amount of total space in the list */
359: if (current_space->local_remaining<cnzi) {
360: PetscFreeSpaceGet(PetscIntSumTruncate(cnzi,current_space->total_array_size),¤t_space);
361: ndouble++;
362: }
364: /* Copy data into free space, then initialize lnk */
365: PetscLLCondensedClean_fast(cnzi,current_space->array,lnk);
367: current_space->array += cnzi;
368: current_space->local_used += cnzi;
369: current_space->local_remaining -= cnzi;
371: ci[i+1] = ci[i] + cnzi;
372: }
374: /* Column indices are in the list of free space */
375: /* Allocate space for cj, initialize cj, and */
376: /* destroy list of free space and other temporary array(s) */
377: PetscMalloc1(ci[am]+1,&cj);
378: PetscFreeSpaceContiguous(&free_space,cj);
379: PetscLLCondensedDestroy_fast(lnk);
381: /* Allocate space for ca */
382: PetscMalloc1(ci[am]+1,&ca);
383: PetscMemzero(ca,(ci[am]+1)*sizeof(MatScalar));
385: /* put together the new symbolic matrix */
386: MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),am,bn,ci,cj,ca,C);
387: MatSetBlockSizesFromMats(*C,A,B);
388: MatSetType(*C,((PetscObject)A)->type_name);
390: /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
391: /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
392: c = (Mat_SeqAIJ*)((*C)->data);
393: c->free_a = PETSC_TRUE;
394: c->free_ij = PETSC_TRUE;
395: c->nonew = 0;
397: (*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable; /* slower, less memory */
399: /* set MatInfo */
400: afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;
401: if (afill < 1.0) afill = 1.0;
402: c->maxnz = ci[am];
403: c->nz = ci[am];
404: (*C)->info.mallocs = ndouble;
405: (*C)->info.fill_ratio_given = fill;
406: (*C)->info.fill_ratio_needed = afill;
408: #if defined(PETSC_USE_INFO)
409: if (ci[am]) {
410: PetscInfo3((*C),"Reallocs %D; Fill ratio: given %g needed %g.\n",ndouble,(double)fill,(double)afill);
411: PetscInfo1((*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n",(double)afill);
412: } else {
413: PetscInfo((*C),"Empty matrix product\n");
414: }
415: #endif
416: return(0);
417: }
420: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable(Mat A,Mat B,PetscReal fill,Mat *C)
421: {
422: PetscErrorCode ierr;
423: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
424: PetscInt *ai = a->i,*bi=b->i,*ci,*cj;
425: PetscInt am = A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;
426: MatScalar *ca;
427: PetscReal afill;
428: PetscInt i,j,anzi,brow,bnzj,cnzi,*bj,*aj,*lnk,ndouble=0,Crmax;
429: PetscTable ta;
430: PetscFreeSpaceList free_space=NULL,current_space=NULL;
433: /* Get ci and cj - same as MatMatMultSymbolic_SeqAIJ_SeqAIJ except using PetscLLxxx_Scalalbe() */
434: /*---------------------------------------------------------------------------------------------*/
435: /* Allocate arrays for fill computation and free space for accumulating nonzero column */
436: PetscMalloc1(am+2,&ci);
437: ci[0] = 0;
439: /* create and initialize a linked list */
440: PetscTableCreate(bn,bn,&ta);
441: MatRowMergeMax_SeqAIJ(b,bm,ta);
442: PetscTableGetCount(ta,&Crmax);
443: PetscTableDestroy(&ta);
444: PetscLLCondensedCreate_Scalable(Crmax,&lnk);
446: /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */
447: PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ai[am],bi[bm])),&free_space);
448: current_space = free_space;
450: /* Determine ci and cj */
451: for (i=0; i<am; i++) {
452: anzi = ai[i+1] - ai[i];
453: aj = a->j + ai[i];
454: for (j=0; j<anzi; j++) {
455: brow = aj[j];
456: bnzj = bi[brow+1] - bi[brow];
457: bj = b->j + bi[brow];
458: /* add non-zero cols of B into the sorted linked list lnk */
459: PetscLLCondensedAddSorted_Scalable(bnzj,bj,lnk);
460: }
461: cnzi = lnk[0];
463: /* If free space is not available, make more free space */
464: /* Double the amount of total space in the list */
465: if (current_space->local_remaining<cnzi) {
466: PetscFreeSpaceGet(PetscIntSumTruncate(cnzi,current_space->total_array_size),¤t_space);
467: ndouble++;
468: }
470: /* Copy data into free space, then initialize lnk */
471: PetscLLCondensedClean_Scalable(cnzi,current_space->array,lnk);
473: current_space->array += cnzi;
474: current_space->local_used += cnzi;
475: current_space->local_remaining -= cnzi;
477: ci[i+1] = ci[i] + cnzi;
478: }
480: /* Column indices are in the list of free space */
481: /* Allocate space for cj, initialize cj, and */
482: /* destroy list of free space and other temporary array(s) */
483: PetscMalloc1(ci[am]+1,&cj);
484: PetscFreeSpaceContiguous(&free_space,cj);
485: PetscLLCondensedDestroy_Scalable(lnk);
487: /* Allocate space for ca */
488: /*-----------------------*/
489: PetscMalloc1(ci[am]+1,&ca);
490: PetscMemzero(ca,(ci[am]+1)*sizeof(MatScalar));
492: /* put together the new symbolic matrix */
493: MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),am,bn,ci,cj,ca,C);
494: MatSetBlockSizesFromMats(*C,A,B);
495: MatSetType(*C,((PetscObject)A)->type_name);
497: /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
498: /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
499: c = (Mat_SeqAIJ*)((*C)->data);
500: c->free_a = PETSC_TRUE;
501: c->free_ij = PETSC_TRUE;
502: c->nonew = 0;
504: (*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable; /* slower, less memory */
506: /* set MatInfo */
507: afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;
508: if (afill < 1.0) afill = 1.0;
509: c->maxnz = ci[am];
510: c->nz = ci[am];
511: (*C)->info.mallocs = ndouble;
512: (*C)->info.fill_ratio_given = fill;
513: (*C)->info.fill_ratio_needed = afill;
515: #if defined(PETSC_USE_INFO)
516: if (ci[am]) {
517: PetscInfo3((*C),"Reallocs %D; Fill ratio: given %g needed %g.\n",ndouble,(double)fill,(double)afill);
518: PetscInfo1((*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n",(double)afill);
519: } else {
520: PetscInfo((*C),"Empty matrix product\n");
521: }
522: #endif
523: return(0);
524: }
526: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_Heap(Mat A,Mat B,PetscReal fill,Mat *C)
527: {
528: PetscErrorCode ierr;
529: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
530: const PetscInt *ai=a->i,*bi=b->i,*aj=a->j,*bj=b->j;
531: PetscInt *ci,*cj,*bb;
532: PetscInt am=A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;
533: PetscReal afill;
534: PetscInt i,j,col,ndouble = 0;
535: PetscFreeSpaceList free_space=NULL,current_space=NULL;
536: PetscHeap h;
539: /* Get ci and cj - by merging sorted rows using a heap */
540: /*---------------------------------------------------------------------------------------------*/
541: /* Allocate arrays for fill computation and free space for accumulating nonzero column */
542: PetscMalloc1(am+2,&ci);
543: ci[0] = 0;
545: /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */
546: PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ai[am],bi[bm])),&free_space);
547: current_space = free_space;
549: PetscHeapCreate(a->rmax,&h);
550: PetscMalloc1(a->rmax,&bb);
552: /* Determine ci and cj */
553: for (i=0; i<am; i++) {
554: const PetscInt anzi = ai[i+1] - ai[i]; /* number of nonzeros in this row of A, this is the number of rows of B that we merge */
555: const PetscInt *acol = aj + ai[i]; /* column indices of nonzero entries in this row */
556: ci[i+1] = ci[i];
557: /* Populate the min heap */
558: for (j=0; j<anzi; j++) {
559: bb[j] = bi[acol[j]]; /* bb points at the start of the row */
560: if (bb[j] < bi[acol[j]+1]) { /* Add if row is nonempty */
561: PetscHeapAdd(h,j,bj[bb[j]++]);
562: }
563: }
564: /* Pick off the min element, adding it to free space */
565: PetscHeapPop(h,&j,&col);
566: while (j >= 0) {
567: if (current_space->local_remaining < 1) { /* double the size, but don't exceed 16 MiB */
568: PetscFreeSpaceGet(PetscMin(PetscIntMultTruncate(2,current_space->total_array_size),16 << 20),¤t_space);
569: ndouble++;
570: }
571: *(current_space->array++) = col;
572: current_space->local_used++;
573: current_space->local_remaining--;
574: ci[i+1]++;
576: /* stash if anything else remains in this row of B */
577: if (bb[j] < bi[acol[j]+1]) {PetscHeapStash(h,j,bj[bb[j]++]);}
578: while (1) { /* pop and stash any other rows of B that also had an entry in this column */
579: PetscInt j2,col2;
580: PetscHeapPeek(h,&j2,&col2);
581: if (col2 != col) break;
582: PetscHeapPop(h,&j2,&col2);
583: if (bb[j2] < bi[acol[j2]+1]) {PetscHeapStash(h,j2,bj[bb[j2]++]);}
584: }
585: /* Put any stashed elements back into the min heap */
586: PetscHeapUnstash(h);
587: PetscHeapPop(h,&j,&col);
588: }
589: }
590: PetscFree(bb);
591: PetscHeapDestroy(&h);
593: /* Column indices are in the list of free space */
594: /* Allocate space for cj, initialize cj, and */
595: /* destroy list of free space and other temporary array(s) */
596: PetscMalloc1(ci[am],&cj);
597: PetscFreeSpaceContiguous(&free_space,cj);
599: /* put together the new symbolic matrix */
600: MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),am,bn,ci,cj,NULL,C);
601: MatSetBlockSizesFromMats(*C,A,B);
602: MatSetType(*C,((PetscObject)A)->type_name);
604: /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
605: /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
606: c = (Mat_SeqAIJ*)((*C)->data);
607: c->free_a = PETSC_TRUE;
608: c->free_ij = PETSC_TRUE;
609: c->nonew = 0;
611: (*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted;
613: /* set MatInfo */
614: afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;
615: if (afill < 1.0) afill = 1.0;
616: c->maxnz = ci[am];
617: c->nz = ci[am];
618: (*C)->info.mallocs = ndouble;
619: (*C)->info.fill_ratio_given = fill;
620: (*C)->info.fill_ratio_needed = afill;
622: #if defined(PETSC_USE_INFO)
623: if (ci[am]) {
624: PetscInfo3((*C),"Reallocs %D; Fill ratio: given %g needed %g.\n",ndouble,(double)fill,(double)afill);
625: PetscInfo1((*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n",(double)afill);
626: } else {
627: PetscInfo((*C),"Empty matrix product\n");
628: }
629: #endif
630: return(0);
631: }
633: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_BTHeap(Mat A,Mat B,PetscReal fill,Mat *C)
634: {
635: PetscErrorCode ierr;
636: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
637: const PetscInt *ai = a->i,*bi=b->i,*aj=a->j,*bj=b->j;
638: PetscInt *ci,*cj,*bb;
639: PetscInt am=A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;
640: PetscReal afill;
641: PetscInt i,j,col,ndouble = 0;
642: PetscFreeSpaceList free_space=NULL,current_space=NULL;
643: PetscHeap h;
644: PetscBT bt;
647: /* Get ci and cj - using a heap for the sorted rows, but use BT so that each index is only added once */
648: /*---------------------------------------------------------------------------------------------*/
649: /* Allocate arrays for fill computation and free space for accumulating nonzero column */
650: PetscMalloc1(am+2,&ci);
651: ci[0] = 0;
653: /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */
654: PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ai[am],bi[bm])),&free_space);
656: current_space = free_space;
658: PetscHeapCreate(a->rmax,&h);
659: PetscMalloc1(a->rmax,&bb);
660: PetscBTCreate(bn,&bt);
662: /* Determine ci and cj */
663: for (i=0; i<am; i++) {
664: const PetscInt anzi = ai[i+1] - ai[i]; /* number of nonzeros in this row of A, this is the number of rows of B that we merge */
665: const PetscInt *acol = aj + ai[i]; /* column indices of nonzero entries in this row */
666: const PetscInt *fptr = current_space->array; /* Save beginning of the row so we can clear the BT later */
667: ci[i+1] = ci[i];
668: /* Populate the min heap */
669: for (j=0; j<anzi; j++) {
670: PetscInt brow = acol[j];
671: for (bb[j] = bi[brow]; bb[j] < bi[brow+1]; bb[j]++) {
672: PetscInt bcol = bj[bb[j]];
673: if (!PetscBTLookupSet(bt,bcol)) { /* new entry */
674: PetscHeapAdd(h,j,bcol);
675: bb[j]++;
676: break;
677: }
678: }
679: }
680: /* Pick off the min element, adding it to free space */
681: PetscHeapPop(h,&j,&col);
682: while (j >= 0) {
683: if (current_space->local_remaining < 1) { /* double the size, but don't exceed 16 MiB */
684: fptr = NULL; /* need PetscBTMemzero */
685: PetscFreeSpaceGet(PetscMin(PetscIntMultTruncate(2,current_space->total_array_size),16 << 20),¤t_space);
686: ndouble++;
687: }
688: *(current_space->array++) = col;
689: current_space->local_used++;
690: current_space->local_remaining--;
691: ci[i+1]++;
693: /* stash if anything else remains in this row of B */
694: for (; bb[j] < bi[acol[j]+1]; bb[j]++) {
695: PetscInt bcol = bj[bb[j]];
696: if (!PetscBTLookupSet(bt,bcol)) { /* new entry */
697: PetscHeapAdd(h,j,bcol);
698: bb[j]++;
699: break;
700: }
701: }
702: PetscHeapPop(h,&j,&col);
703: }
704: if (fptr) { /* Clear the bits for this row */
705: for (; fptr<current_space->array; fptr++) {PetscBTClear(bt,*fptr);}
706: } else { /* We reallocated so we don't remember (easily) how to clear only the bits we changed */
707: PetscBTMemzero(bn,bt);
708: }
709: }
710: PetscFree(bb);
711: PetscHeapDestroy(&h);
712: PetscBTDestroy(&bt);
714: /* Column indices are in the list of free space */
715: /* Allocate space for cj, initialize cj, and */
716: /* destroy list of free space and other temporary array(s) */
717: PetscMalloc1(ci[am],&cj);
718: PetscFreeSpaceContiguous(&free_space,cj);
720: /* put together the new symbolic matrix */
721: MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),am,bn,ci,cj,NULL,C);
722: MatSetBlockSizesFromMats(*C,A,B);
723: MatSetType(*C,((PetscObject)A)->type_name);
725: /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
726: /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
727: c = (Mat_SeqAIJ*)((*C)->data);
728: c->free_a = PETSC_TRUE;
729: c->free_ij = PETSC_TRUE;
730: c->nonew = 0;
732: (*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted;
734: /* set MatInfo */
735: afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;
736: if (afill < 1.0) afill = 1.0;
737: c->maxnz = ci[am];
738: c->nz = ci[am];
739: (*C)->info.mallocs = ndouble;
740: (*C)->info.fill_ratio_given = fill;
741: (*C)->info.fill_ratio_needed = afill;
743: #if defined(PETSC_USE_INFO)
744: if (ci[am]) {
745: PetscInfo3((*C),"Reallocs %D; Fill ratio: given %g needed %g.\n",ndouble,(double)fill,(double)afill);
746: PetscInfo1((*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n",(double)afill);
747: } else {
748: PetscInfo((*C),"Empty matrix product\n");
749: }
750: #endif
751: return(0);
752: }
755: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_RowMerge(Mat A,Mat B,PetscReal fill,Mat *C)
756: {
757: PetscErrorCode ierr;
758: Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
759: const PetscInt *ai=a->i,*bi=b->i,*aj=a->j,*bj=b->j,*inputi,*inputj,*inputcol,*inputcol_L1;
760: PetscInt *ci,*cj,*outputj,worki_L1[9],worki_L2[9];
761: PetscInt c_maxmem,a_maxrownnz=0,a_rownnz;
762: const PetscInt workcol[8]={0,1,2,3,4,5,6,7};
763: const PetscInt am=A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;
764: const PetscInt *brow_ptr[8],*brow_end[8];
765: PetscInt window[8];
766: PetscInt window_min,old_window_min,ci_nnz,outputi_nnz=0,L1_nrows,L2_nrows;
767: PetscInt i,k,ndouble=0,L1_rowsleft,rowsleft;
768: PetscReal afill;
769: PetscInt *workj_L1,*workj_L2,*workj_L3;
770: PetscInt L1_nnz,L2_nnz;
772: /* Step 1: Get upper bound on memory required for allocation.
773: Because of the way virtual memory works,
774: only the memory pages that are actually needed will be physically allocated. */
776: PetscMalloc1(am+1,&ci);
777: for (i=0; i<am; i++) {
778: const PetscInt anzi = ai[i+1] - ai[i]; /* number of nonzeros in this row of A, this is the number of rows of B that we merge */
779: const PetscInt *acol = aj + ai[i]; /* column indices of nonzero entries in this row */
780: a_rownnz = 0;
781: for (k=0; k<anzi; ++k) {
782: a_rownnz += bi[acol[k]+1] - bi[acol[k]];
783: if (a_rownnz > bn) {
784: a_rownnz = bn;
785: break;
786: }
787: }
788: a_maxrownnz = PetscMax(a_maxrownnz, a_rownnz);
789: }
790: /* temporary work areas for merging rows */
791: PetscMalloc1(a_maxrownnz*8,&workj_L1);
792: PetscMalloc1(a_maxrownnz*8,&workj_L2);
793: PetscMalloc1(a_maxrownnz,&workj_L3);
795: /* This should be enough for almost all matrices. If not, memory is reallocated later. */
796: c_maxmem = 8*(ai[am]+bi[bm]);
797: /* Step 2: Populate pattern for C */
798: PetscMalloc1(c_maxmem,&cj);
800: ci_nnz = 0;
801: ci[0] = 0;
802: worki_L1[0] = 0;
803: worki_L2[0] = 0;
804: for (i=0; i<am; i++) {
805: const PetscInt anzi = ai[i+1] - ai[i]; /* number of nonzeros in this row of A, this is the number of rows of B that we merge */
806: const PetscInt *acol = aj + ai[i]; /* column indices of nonzero entries in this row */
807: rowsleft = anzi;
808: inputcol_L1 = acol;
809: L2_nnz = 0;
810: L2_nrows = 1; /* Number of rows to be merged on Level 3. output of L3 already exists -> initial value 1 */
811: worki_L2[1] = 0;
812: outputi_nnz = 0;
814: /* If the number of indices in C so far + the max number of columns in the next row > c_maxmem -> allocate more memory */
815: while (ci_nnz+a_maxrownnz > c_maxmem) {
816: c_maxmem *= 2;
817: ndouble++;
818: PetscRealloc(sizeof(PetscInt)*c_maxmem,&cj);
819: }
821: while (rowsleft) {
822: L1_rowsleft = PetscMin(64, rowsleft); /* In the inner loop max 64 rows of B can be merged */
823: L1_nrows = 0;
824: L1_nnz = 0;
825: inputcol = inputcol_L1;
826: inputi = bi;
827: inputj = bj;
829: /* The following macro is used to specialize for small rows in A.
830: This helps with compiler unrolling, improving performance substantially.
831: Input: inputj inputi inputcol bn
832: Output: outputj outputi_nnz */
833: #define MatMatMultSymbolic_RowMergeMacro(ANNZ) \
834: window_min = bn; \
835: outputi_nnz = 0; \
836: for (k=0; k<ANNZ; ++k) { \
837: brow_ptr[k] = inputj + inputi[inputcol[k]]; \
838: brow_end[k] = inputj + inputi[inputcol[k]+1]; \
839: window[k] = (brow_ptr[k] != brow_end[k]) ? *brow_ptr[k] : bn; \
840: window_min = PetscMin(window[k], window_min); \
841: } \
842: while (window_min < bn) { \
843: outputj[outputi_nnz++] = window_min; \
844: /* advance front and compute new minimum */ \
845: old_window_min = window_min; \
846: window_min = bn; \
847: for (k=0; k<ANNZ; ++k) { \
848: if (window[k] == old_window_min) { \
849: brow_ptr[k]++; \
850: window[k] = (brow_ptr[k] != brow_end[k]) ? *brow_ptr[k] : bn; \
851: } \
852: window_min = PetscMin(window[k], window_min); \
853: } \
854: }
856: /************** L E V E L 1 ***************/
857: /* Merge up to 8 rows of B to L1 work array*/
858: while (L1_rowsleft) {
859: outputi_nnz = 0;
860: if (anzi > 8) outputj = workj_L1 + L1_nnz; /* Level 1 rowmerge*/
861: else outputj = cj + ci_nnz; /* Merge directly to C */
863: switch (L1_rowsleft) {
864: case 1: brow_ptr[0] = inputj + inputi[inputcol[0]];
865: brow_end[0] = inputj + inputi[inputcol[0]+1];
866: for (; brow_ptr[0] != brow_end[0]; ++brow_ptr[0]) outputj[outputi_nnz++] = *brow_ptr[0]; /* copy row in b over */
867: inputcol += L1_rowsleft;
868: rowsleft -= L1_rowsleft;
869: L1_rowsleft = 0;
870: break;
871: case 2: MatMatMultSymbolic_RowMergeMacro(2);
872: inputcol += L1_rowsleft;
873: rowsleft -= L1_rowsleft;
874: L1_rowsleft = 0;
875: break;
876: case 3: MatMatMultSymbolic_RowMergeMacro(3);
877: inputcol += L1_rowsleft;
878: rowsleft -= L1_rowsleft;
879: L1_rowsleft = 0;
880: break;
881: case 4: MatMatMultSymbolic_RowMergeMacro(4);
882: inputcol += L1_rowsleft;
883: rowsleft -= L1_rowsleft;
884: L1_rowsleft = 0;
885: break;
886: case 5: MatMatMultSymbolic_RowMergeMacro(5);
887: inputcol += L1_rowsleft;
888: rowsleft -= L1_rowsleft;
889: L1_rowsleft = 0;
890: break;
891: case 6: MatMatMultSymbolic_RowMergeMacro(6);
892: inputcol += L1_rowsleft;
893: rowsleft -= L1_rowsleft;
894: L1_rowsleft = 0;
895: break;
896: case 7: MatMatMultSymbolic_RowMergeMacro(7);
897: inputcol += L1_rowsleft;
898: rowsleft -= L1_rowsleft;
899: L1_rowsleft = 0;
900: break;
901: default: MatMatMultSymbolic_RowMergeMacro(8);
902: inputcol += 8;
903: rowsleft -= 8;
904: L1_rowsleft -= 8;
905: break;
906: }
907: inputcol_L1 = inputcol;
908: L1_nnz += outputi_nnz;
909: worki_L1[++L1_nrows] = L1_nnz;
910: }
912: /********************** L E V E L 2 ************************/
913: /* Merge from L1 work array to either C or to L2 work array */
914: if (anzi > 8) {
915: inputi = worki_L1;
916: inputj = workj_L1;
917: inputcol = workcol;
918: outputi_nnz = 0;
920: if (anzi <= 64) outputj = cj + ci_nnz; /* Merge from L1 work array to C */
921: else outputj = workj_L2 + L2_nnz; /* Merge from L1 work array to L2 work array */
923: switch (L1_nrows) {
924: case 1: brow_ptr[0] = inputj + inputi[inputcol[0]];
925: brow_end[0] = inputj + inputi[inputcol[0]+1];
926: for (; brow_ptr[0] != brow_end[0]; ++brow_ptr[0]) outputj[outputi_nnz++] = *brow_ptr[0]; /* copy row in b over */
927: break;
928: case 2: MatMatMultSymbolic_RowMergeMacro(2); break;
929: case 3: MatMatMultSymbolic_RowMergeMacro(3); break;
930: case 4: MatMatMultSymbolic_RowMergeMacro(4); break;
931: case 5: MatMatMultSymbolic_RowMergeMacro(5); break;
932: case 6: MatMatMultSymbolic_RowMergeMacro(6); break;
933: case 7: MatMatMultSymbolic_RowMergeMacro(7); break;
934: case 8: MatMatMultSymbolic_RowMergeMacro(8); break;
935: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatMatMult logic error: Not merging 1-8 rows from L1 work array!");
936: }
937: L2_nnz += outputi_nnz;
938: worki_L2[++L2_nrows] = L2_nnz;
940: /************************ L E V E L 3 **********************/
941: /* Merge from L2 work array to either C or to L2 work array */
942: if (anzi > 64 && (L2_nrows == 8 || rowsleft == 0)) {
943: inputi = worki_L2;
944: inputj = workj_L2;
945: inputcol = workcol;
946: outputi_nnz = 0;
947: if (rowsleft) outputj = workj_L3;
948: else outputj = cj + ci_nnz;
949: switch (L2_nrows) {
950: case 1: brow_ptr[0] = inputj + inputi[inputcol[0]];
951: brow_end[0] = inputj + inputi[inputcol[0]+1];
952: for (; brow_ptr[0] != brow_end[0]; ++brow_ptr[0]) outputj[outputi_nnz++] = *brow_ptr[0]; /* copy row in b over */
953: break;
954: case 2: MatMatMultSymbolic_RowMergeMacro(2); break;
955: case 3: MatMatMultSymbolic_RowMergeMacro(3); break;
956: case 4: MatMatMultSymbolic_RowMergeMacro(4); break;
957: case 5: MatMatMultSymbolic_RowMergeMacro(5); break;
958: case 6: MatMatMultSymbolic_RowMergeMacro(6); break;
959: case 7: MatMatMultSymbolic_RowMergeMacro(7); break;
960: case 8: MatMatMultSymbolic_RowMergeMacro(8); break;
961: default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MatMatMult logic error: Not merging 1-8 rows from L2 work array!");
962: }
963: L2_nrows = 1;
964: L2_nnz = outputi_nnz;
965: worki_L2[1] = outputi_nnz;
966: /* Copy to workj_L2 */
967: if (rowsleft) {
968: for (k=0; k<outputi_nnz; ++k) workj_L2[k] = outputj[k];
969: }
970: }
971: }
972: } /* while (rowsleft) */
973: #undef MatMatMultSymbolic_RowMergeMacro
975: /* terminate current row */
976: ci_nnz += outputi_nnz;
977: ci[i+1] = ci_nnz;
978: }
980: /* Step 3: Create the new symbolic matrix */
981: MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),am,bn,ci,cj,NULL,C);
982: MatSetBlockSizesFromMats(*C,A,B);
983: MatSetType(*C,((PetscObject)A)->type_name);
985: /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
986: /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
987: c = (Mat_SeqAIJ*)((*C)->data);
988: c->free_a = PETSC_TRUE;
989: c->free_ij = PETSC_TRUE;
990: c->nonew = 0;
992: (*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted;
994: /* set MatInfo */
995: afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;
996: if (afill < 1.0) afill = 1.0;
997: c->maxnz = ci[am];
998: c->nz = ci[am];
999: (*C)->info.mallocs = ndouble;
1000: (*C)->info.fill_ratio_given = fill;
1001: (*C)->info.fill_ratio_needed = afill;
1003: #if defined(PETSC_USE_INFO)
1004: if (ci[am]) {
1005: PetscInfo3((*C),"Reallocs %D; Fill ratio: given %g needed %g.\n",ndouble,(double)fill,(double)afill);
1006: PetscInfo1((*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n",(double)afill);
1007: } else {
1008: PetscInfo((*C),"Empty matrix product\n");
1009: }
1010: #endif
1012: /* Step 4: Free temporary work areas */
1013: PetscFree(workj_L1);
1014: PetscFree(workj_L2);
1015: PetscFree(workj_L3);
1016: return(0);
1017: }
1019: /* concatenate unique entries and then sort */
1020: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_Sorted(Mat A,Mat B,PetscReal fill,Mat *C)
1021: {
1022: PetscErrorCode ierr;
1023: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
1024: const PetscInt *ai = a->i,*bi=b->i,*aj=a->j,*bj=b->j;
1025: PetscInt *ci,*cj;
1026: PetscInt am=A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;
1027: PetscReal afill;
1028: PetscInt i,j,ndouble = 0;
1029: PetscSegBuffer seg,segrow;
1030: char *seen;
1033: PetscMalloc1(am+1,&ci);
1034: ci[0] = 0;
1036: /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */
1037: PetscSegBufferCreate(sizeof(PetscInt),(PetscInt)(fill*(ai[am]+bi[bm])),&seg);
1038: PetscSegBufferCreate(sizeof(PetscInt),100,&segrow);
1039: PetscMalloc1(bn,&seen);
1040: PetscMemzero(seen,bn*sizeof(char));
1042: /* Determine ci and cj */
1043: for (i=0; i<am; i++) {
1044: const PetscInt anzi = ai[i+1] - ai[i]; /* number of nonzeros in this row of A, this is the number of rows of B that we merge */
1045: const PetscInt *acol = aj + ai[i]; /* column indices of nonzero entries in this row */
1046: PetscInt packlen = 0,*PETSC_RESTRICT crow;
1047: /* Pack segrow */
1048: for (j=0; j<anzi; j++) {
1049: PetscInt brow = acol[j],bjstart = bi[brow],bjend = bi[brow+1],k;
1050: for (k=bjstart; k<bjend; k++) {
1051: PetscInt bcol = bj[k];
1052: if (!seen[bcol]) { /* new entry */
1053: PetscInt *PETSC_RESTRICT slot;
1054: PetscSegBufferGetInts(segrow,1,&slot);
1055: *slot = bcol;
1056: seen[bcol] = 1;
1057: packlen++;
1058: }
1059: }
1060: }
1061: PetscSegBufferGetInts(seg,packlen,&crow);
1062: PetscSegBufferExtractTo(segrow,crow);
1063: PetscSortInt(packlen,crow);
1064: ci[i+1] = ci[i] + packlen;
1065: for (j=0; j<packlen; j++) seen[crow[j]] = 0;
1066: }
1067: PetscSegBufferDestroy(&segrow);
1068: PetscFree(seen);
1070: /* Column indices are in the segmented buffer */
1071: PetscSegBufferExtractAlloc(seg,&cj);
1072: PetscSegBufferDestroy(&seg);
1074: /* put together the new symbolic matrix */
1075: MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),am,bn,ci,cj,NULL,C);
1076: MatSetBlockSizesFromMats(*C,A,B);
1077: MatSetType(*C,((PetscObject)A)->type_name);
1079: /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
1080: /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
1081: c = (Mat_SeqAIJ*)((*C)->data);
1082: c->free_a = PETSC_TRUE;
1083: c->free_ij = PETSC_TRUE;
1084: c->nonew = 0;
1086: (*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted;
1088: /* set MatInfo */
1089: afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;
1090: if (afill < 1.0) afill = 1.0;
1091: c->maxnz = ci[am];
1092: c->nz = ci[am];
1093: (*C)->info.mallocs = ndouble;
1094: (*C)->info.fill_ratio_given = fill;
1095: (*C)->info.fill_ratio_needed = afill;
1097: #if defined(PETSC_USE_INFO)
1098: if (ci[am]) {
1099: PetscInfo3((*C),"Reallocs %D; Fill ratio: given %g needed %g.\n",ndouble,(double)fill,(double)afill);
1100: PetscInfo1((*C),"Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n",(double)afill);
1101: } else {
1102: PetscInfo((*C),"Empty matrix product\n");
1103: }
1104: #endif
1105: return(0);
1106: }
1108: /* This routine is not used. Should be removed! */
1109: PetscErrorCode MatMatTransposeMult_SeqAIJ_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
1110: {
1114: if (scall == MAT_INITIAL_MATRIX) {
1115: PetscLogEventBegin(MAT_MatTransposeMultSymbolic,A,B,0,0);
1116: MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ(A,B,fill,C);
1117: PetscLogEventEnd(MAT_MatTransposeMultSymbolic,A,B,0,0);
1118: }
1119: PetscLogEventBegin(MAT_MatTransposeMultNumeric,A,B,0,0);
1120: MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(A,B,*C);
1121: PetscLogEventEnd(MAT_MatTransposeMultNumeric,A,B,0,0);
1122: return(0);
1123: }
1125: PetscErrorCode MatDestroy_SeqAIJ_MatMatMultTrans(Mat A)
1126: {
1127: PetscErrorCode ierr;
1128: Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data;
1129: Mat_MatMatTransMult *abt=a->abt;
1132: (abt->destroy)(A);
1133: MatTransposeColoringDestroy(&abt->matcoloring);
1134: MatDestroy(&abt->Bt_den);
1135: MatDestroy(&abt->ABt_den);
1136: PetscFree(abt);
1137: return(0);
1138: }
1140: PetscErrorCode MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C)
1141: {
1142: PetscErrorCode ierr;
1143: Mat Bt;
1144: PetscInt *bti,*btj;
1145: Mat_MatMatTransMult *abt;
1146: Mat_SeqAIJ *c;
1149: /* create symbolic Bt */
1150: MatGetSymbolicTranspose_SeqAIJ(B,&bti,&btj);
1151: MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,B->cmap->n,B->rmap->n,bti,btj,NULL,&Bt);
1152: MatSetBlockSizes(Bt,PetscAbs(A->cmap->bs),PetscAbs(B->cmap->bs));
1153: MatSetType(Bt,((PetscObject)A)->type_name);
1155: /* get symbolic C=A*Bt */
1156: MatMatMultSymbolic_SeqAIJ_SeqAIJ(A,Bt,fill,C);
1158: /* create a supporting struct for reuse intermidiate dense matrices with matcoloring */
1159: PetscNew(&abt);
1160: c = (Mat_SeqAIJ*)(*C)->data;
1161: c->abt = abt;
1163: abt->usecoloring = PETSC_FALSE;
1164: abt->destroy = (*C)->ops->destroy;
1165: (*C)->ops->destroy = MatDestroy_SeqAIJ_MatMatMultTrans;
1167: PetscOptionsGetBool(((PetscObject)A)->options,NULL,"-matmattransmult_color",&abt->usecoloring,NULL);
1168: if (abt->usecoloring) {
1169: /* Create MatTransposeColoring from symbolic C=A*B^T */
1170: MatTransposeColoring matcoloring;
1171: MatColoring coloring;
1172: ISColoring iscoloring;
1173: Mat Bt_dense,C_dense;
1174: Mat_SeqAIJ *c=(Mat_SeqAIJ*)(*C)->data;
1175: /* inode causes memory problem, don't know why */
1176: if (c->inode.use) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"MAT_USE_INODES is not supported. Use '-mat_no_inode'");
1178: MatColoringCreate(*C,&coloring);
1179: MatColoringSetDistance(coloring,2);
1180: MatColoringSetType(coloring,MATCOLORINGSL);
1181: MatColoringSetFromOptions(coloring);
1182: MatColoringApply(coloring,&iscoloring);
1183: MatColoringDestroy(&coloring);
1184: MatTransposeColoringCreate(*C,iscoloring,&matcoloring);
1186: abt->matcoloring = matcoloring;
1188: ISColoringDestroy(&iscoloring);
1190: /* Create Bt_dense and C_dense = A*Bt_dense */
1191: MatCreate(PETSC_COMM_SELF,&Bt_dense);
1192: MatSetSizes(Bt_dense,A->cmap->n,matcoloring->ncolors,A->cmap->n,matcoloring->ncolors);
1193: MatSetType(Bt_dense,MATSEQDENSE);
1194: MatSeqDenseSetPreallocation(Bt_dense,NULL);
1196: Bt_dense->assembled = PETSC_TRUE;
1197: abt->Bt_den = Bt_dense;
1199: MatCreate(PETSC_COMM_SELF,&C_dense);
1200: MatSetSizes(C_dense,A->rmap->n,matcoloring->ncolors,A->rmap->n,matcoloring->ncolors);
1201: MatSetType(C_dense,MATSEQDENSE);
1202: MatSeqDenseSetPreallocation(C_dense,NULL);
1204: Bt_dense->assembled = PETSC_TRUE;
1205: abt->ABt_den = C_dense;
1207: #if defined(PETSC_USE_INFO)
1208: {
1209: Mat_SeqAIJ *c = (Mat_SeqAIJ*)(*C)->data;
1210: PetscInfo7(*C,"Use coloring of C=A*B^T; B^T: %D %D, Bt_dense: %D,%D; Cnz %D / (cm*ncolors %D) = %g\n",B->cmap->n,B->rmap->n,Bt_dense->rmap->n,Bt_dense->cmap->n,c->nz,A->rmap->n*matcoloring->ncolors,(PetscReal)(c->nz)/(A->rmap->n*matcoloring->ncolors));
1211: }
1212: #endif
1213: }
1214: /* clean up */
1215: MatDestroy(&Bt);
1216: MatRestoreSymbolicTranspose_SeqAIJ(B,&bti,&btj);
1217: return(0);
1218: }
1220: PetscErrorCode MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C)
1221: {
1222: PetscErrorCode ierr;
1223: Mat_SeqAIJ *a =(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c=(Mat_SeqAIJ*)C->data;
1224: PetscInt *ai =a->i,*aj=a->j,*bi=b->i,*bj=b->j,anzi,bnzj,nexta,nextb,*acol,*bcol,brow;
1225: PetscInt cm =C->rmap->n,*ci=c->i,*cj=c->j,i,j,cnzi,*ccol;
1226: PetscLogDouble flops=0.0;
1227: MatScalar *aa =a->a,*aval,*ba=b->a,*bval,*ca,*cval;
1228: Mat_MatMatTransMult *abt = c->abt;
1231: /* clear old values in C */
1232: if (!c->a) {
1233: PetscMalloc1(ci[cm]+1,&ca);
1234: c->a = ca;
1235: c->free_a = PETSC_TRUE;
1236: } else {
1237: ca = c->a;
1238: }
1239: PetscMemzero(ca,ci[cm]*sizeof(MatScalar));
1241: if (abt->usecoloring) {
1242: MatTransposeColoring matcoloring = abt->matcoloring;
1243: Mat Bt_dense,C_dense = abt->ABt_den;
1245: /* Get Bt_dense by Apply MatTransposeColoring to B */
1246: Bt_dense = abt->Bt_den;
1247: MatTransColoringApplySpToDen(matcoloring,B,Bt_dense);
1249: /* C_dense = A*Bt_dense */
1250: MatMatMultNumeric_SeqAIJ_SeqDense(A,Bt_dense,C_dense);
1252: /* Recover C from C_dense */
1253: MatTransColoringApplyDenToSp(matcoloring,C_dense,C);
1254: return(0);
1255: }
1257: for (i=0; i<cm; i++) {
1258: anzi = ai[i+1] - ai[i];
1259: acol = aj + ai[i];
1260: aval = aa + ai[i];
1261: cnzi = ci[i+1] - ci[i];
1262: ccol = cj + ci[i];
1263: cval = ca + ci[i];
1264: for (j=0; j<cnzi; j++) {
1265: brow = ccol[j];
1266: bnzj = bi[brow+1] - bi[brow];
1267: bcol = bj + bi[brow];
1268: bval = ba + bi[brow];
1270: /* perform sparse inner-product c(i,j)=A[i,:]*B[j,:]^T */
1271: nexta = 0; nextb = 0;
1272: while (nexta<anzi && nextb<bnzj) {
1273: while (nexta < anzi && acol[nexta] < bcol[nextb]) nexta++;
1274: if (nexta == anzi) break;
1275: while (nextb < bnzj && acol[nexta] > bcol[nextb]) nextb++;
1276: if (nextb == bnzj) break;
1277: if (acol[nexta] == bcol[nextb]) {
1278: cval[j] += aval[nexta]*bval[nextb];
1279: nexta++; nextb++;
1280: flops += 2;
1281: }
1282: }
1283: }
1284: }
1285: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
1286: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
1287: PetscLogFlops(flops);
1288: return(0);
1289: }
1291: PetscErrorCode MatDestroy_SeqAIJ_MatTransMatMult(Mat A)
1292: {
1293: PetscErrorCode ierr;
1294: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
1295: Mat_MatTransMatMult *atb = a->atb;
1298: if (atb) {
1299: MatDestroy(&atb->At);
1300: (*atb->destroy)(A);
1301: }
1302: PetscFree(atb);
1303: return(0);
1304: }
1306: PetscErrorCode MatTransposeMatMult_SeqAIJ_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
1307: {
1308: PetscErrorCode ierr;
1309: const char *algTypes[2] = {"matmatmult","outerproduct"};
1310: PetscInt alg=0; /* set default algorithm */
1311: Mat At;
1312: Mat_MatTransMatMult *atb;
1313: Mat_SeqAIJ *c;
1316: if (scall == MAT_INITIAL_MATRIX) {
1317: PetscOptionsBegin(PetscObjectComm((PetscObject)A),((PetscObject)A)->prefix,"MatTransposeMatMult","Mat");
1318: PetscOptionsEList("-mattransposematmult_via","Algorithmic approach","MatTransposeMatMult",algTypes,2,algTypes[0],&alg,NULL);
1319: PetscOptionsEnd();
1321: switch (alg) {
1322: case 1:
1323: MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ(A,B,fill,C);
1324: break;
1325: default:
1326: PetscNew(&atb);
1327: MatTranspose_SeqAIJ(A,MAT_INITIAL_MATRIX,&At);
1328: MatMatMult_SeqAIJ_SeqAIJ(At,B,MAT_INITIAL_MATRIX,fill,C);
1330: c = (Mat_SeqAIJ*)(*C)->data;
1331: c->atb = atb;
1332: atb->At = At;
1333: atb->destroy = (*C)->ops->destroy;
1334: (*C)->ops->destroy = MatDestroy_SeqAIJ_MatTransMatMult;
1336: break;
1337: }
1338: }
1339: if (alg) {
1340: (*(*C)->ops->mattransposemultnumeric)(A,B,*C);
1341: } else if (!alg && scall == MAT_REUSE_MATRIX) {
1342: c = (Mat_SeqAIJ*)(*C)->data;
1343: atb = c->atb;
1344: At = atb->At;
1345: MatTranspose_SeqAIJ(A,MAT_REUSE_MATRIX,&At);
1346: MatMatMult_SeqAIJ_SeqAIJ(At,B,MAT_REUSE_MATRIX,fill,C);
1347: }
1348: return(0);
1349: }
1351: PetscErrorCode MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C)
1352: {
1354: Mat At;
1355: PetscInt *ati,*atj;
1358: /* create symbolic At */
1359: MatGetSymbolicTranspose_SeqAIJ(A,&ati,&atj);
1360: MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,A->cmap->n,A->rmap->n,ati,atj,NULL,&At);
1361: MatSetBlockSizes(At,PetscAbs(A->cmap->bs),PetscAbs(B->cmap->bs));
1362: MatSetType(At,((PetscObject)A)->type_name);
1364: /* get symbolic C=At*B */
1365: MatMatMultSymbolic_SeqAIJ_SeqAIJ(At,B,fill,C);
1367: /* clean up */
1368: MatDestroy(&At);
1369: MatRestoreSymbolicTranspose_SeqAIJ(A,&ati,&atj);
1371: (*C)->ops->mattransposemultnumeric = MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ;
1372: return(0);
1373: }
1375: PetscErrorCode MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ(Mat A,Mat B,Mat C)
1376: {
1378: Mat_SeqAIJ *a =(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c=(Mat_SeqAIJ*)C->data;
1379: PetscInt am =A->rmap->n,anzi,*ai=a->i,*aj=a->j,*bi=b->i,*bj,bnzi,nextb;
1380: PetscInt cm =C->rmap->n,*ci=c->i,*cj=c->j,crow,*cjj,i,j,k;
1381: PetscLogDouble flops=0.0;
1382: MatScalar *aa =a->a,*ba,*ca,*caj;
1385: if (!c->a) {
1386: PetscMalloc1(ci[cm]+1,&ca);
1388: c->a = ca;
1389: c->free_a = PETSC_TRUE;
1390: } else {
1391: ca = c->a;
1392: }
1393: /* clear old values in C */
1394: PetscMemzero(ca,ci[cm]*sizeof(MatScalar));
1396: /* compute A^T*B using outer product (A^T)[:,i]*B[i,:] */
1397: for (i=0; i<am; i++) {
1398: bj = b->j + bi[i];
1399: ba = b->a + bi[i];
1400: bnzi = bi[i+1] - bi[i];
1401: anzi = ai[i+1] - ai[i];
1402: for (j=0; j<anzi; j++) {
1403: nextb = 0;
1404: crow = *aj++;
1405: cjj = cj + ci[crow];
1406: caj = ca + ci[crow];
1407: /* perform sparse axpy operation. Note cjj includes bj. */
1408: for (k=0; nextb<bnzi; k++) {
1409: if (cjj[k] == *(bj+nextb)) { /* ccol == bcol */
1410: caj[k] += (*aa)*(*(ba+nextb));
1411: nextb++;
1412: }
1413: }
1414: flops += 2*bnzi;
1415: aa++;
1416: }
1417: }
1419: /* Assemble the final matrix and clean up */
1420: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
1421: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
1422: PetscLogFlops(flops);
1423: return(0);
1424: }
1426: PETSC_INTERN PetscErrorCode MatMatMult_SeqAIJ_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
1427: {
1431: if (scall == MAT_INITIAL_MATRIX) {
1432: PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);
1433: MatMatMultSymbolic_SeqAIJ_SeqDense(A,B,fill,C);
1434: PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);
1435: }
1436: PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);
1437: MatMatMultNumeric_SeqAIJ_SeqDense(A,B,*C);
1438: PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);
1439: return(0);
1440: }
1442: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C)
1443: {
1447: MatMatMultSymbolic_SeqDense_SeqDense(A,B,0.0,C);
1449: (*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqDense;
1450: return(0);
1451: }
1453: PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqDense(Mat A,Mat B,Mat C)
1454: {
1455: Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data;
1456: Mat_SeqDense *bd = (Mat_SeqDense*)B->data;
1457: PetscErrorCode ierr;
1458: PetscScalar *c,*b,r1,r2,r3,r4,*c1,*c2,*c3,*c4,aatmp;
1459: const PetscScalar *aa,*b1,*b2,*b3,*b4;
1460: const PetscInt *aj;
1461: PetscInt cm=C->rmap->n,cn=B->cmap->n,bm=bd->lda,am=A->rmap->n;
1462: PetscInt am4=4*am,bm4=4*bm,col,i,j,n,ajtmp;
1465: if (!cm || !cn) return(0);
1466: if (B->rmap->n != 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,B->rmap->n);
1467: 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);
1468: 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);
1469: b = bd->v;
1470: MatDenseGetArray(C,&c);
1471: b1 = b; b2 = b1 + bm; b3 = b2 + bm; b4 = b3 + bm;
1472: c1 = c; c2 = c1 + am; c3 = c2 + am; c4 = c3 + am;
1473: for (col=0; col<cn-4; col += 4) { /* over columns of C */
1474: for (i=0; i<am; i++) { /* over rows of C in those columns */
1475: r1 = r2 = r3 = r4 = 0.0;
1476: n = a->i[i+1] - a->i[i];
1477: aj = a->j + a->i[i];
1478: aa = a->a + a->i[i];
1479: for (j=0; j<n; j++) {
1480: aatmp = aa[j]; ajtmp = aj[j];
1481: r1 += aatmp*b1[ajtmp];
1482: r2 += aatmp*b2[ajtmp];
1483: r3 += aatmp*b3[ajtmp];
1484: r4 += aatmp*b4[ajtmp];
1485: }
1486: c1[i] = r1;
1487: c2[i] = r2;
1488: c3[i] = r3;
1489: c4[i] = r4;
1490: }
1491: b1 += bm4; b2 += bm4; b3 += bm4; b4 += bm4;
1492: c1 += am4; c2 += am4; c3 += am4; c4 += am4;
1493: }
1494: for (; col<cn; col++) { /* over extra columns of C */
1495: for (i=0; i<am; i++) { /* over rows of C in those columns */
1496: r1 = 0.0;
1497: n = a->i[i+1] - a->i[i];
1498: aj = a->j + a->i[i];
1499: aa = a->a + a->i[i];
1500: for (j=0; j<n; j++) {
1501: r1 += aa[j]*b1[aj[j]];
1502: }
1503: c1[i] = r1;
1504: }
1505: b1 += bm;
1506: c1 += am;
1507: }
1508: PetscLogFlops(cn*(2.0*a->nz));
1509: MatDenseRestoreArray(C,&c);
1510: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
1511: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
1512: return(0);
1513: }
1515: /*
1516: Note very similar to MatMult_SeqAIJ(), should generate both codes from same base
1517: */
1518: PetscErrorCode MatMatMultNumericAdd_SeqAIJ_SeqDense(Mat A,Mat B,Mat C)
1519: {
1520: Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data;
1521: Mat_SeqDense *bd = (Mat_SeqDense*)B->data;
1523: PetscScalar *b,*c,r1,r2,r3,r4,*b1,*b2,*b3,*b4;
1524: MatScalar *aa;
1525: PetscInt cm = C->rmap->n, cn=B->cmap->n, bm=bd->lda, col, i,j,n,*aj, am = A->rmap->n,*ii,arm;
1526: PetscInt am2 = 2*am, am3 = 3*am, bm4 = 4*bm,colam,*ridx;
1529: if (!cm || !cn) return(0);
1530: b = bd->v;
1531: MatDenseGetArray(C,&c);
1532: b1 = b; b2 = b1 + bm; b3 = b2 + bm; b4 = b3 + bm;
1534: if (a->compressedrow.use) { /* use compressed row format */
1535: for (col=0; col<cn-4; col += 4) { /* over columns of C */
1536: colam = col*am;
1537: arm = a->compressedrow.nrows;
1538: ii = a->compressedrow.i;
1539: ridx = a->compressedrow.rindex;
1540: for (i=0; i<arm; i++) { /* over rows of C in those columns */
1541: r1 = r2 = r3 = r4 = 0.0;
1542: n = ii[i+1] - ii[i];
1543: aj = a->j + ii[i];
1544: aa = a->a + ii[i];
1545: for (j=0; j<n; j++) {
1546: r1 += (*aa)*b1[*aj];
1547: r2 += (*aa)*b2[*aj];
1548: r3 += (*aa)*b3[*aj];
1549: r4 += (*aa++)*b4[*aj++];
1550: }
1551: c[colam + ridx[i]] += r1;
1552: c[colam + am + ridx[i]] += r2;
1553: c[colam + am2 + ridx[i]] += r3;
1554: c[colam + am3 + ridx[i]] += r4;
1555: }
1556: b1 += bm4;
1557: b2 += bm4;
1558: b3 += bm4;
1559: b4 += bm4;
1560: }
1561: for (; col<cn; col++) { /* over extra columns of C */
1562: colam = col*am;
1563: arm = a->compressedrow.nrows;
1564: ii = a->compressedrow.i;
1565: ridx = a->compressedrow.rindex;
1566: for (i=0; i<arm; i++) { /* over rows of C in those columns */
1567: r1 = 0.0;
1568: n = ii[i+1] - ii[i];
1569: aj = a->j + ii[i];
1570: aa = a->a + ii[i];
1572: for (j=0; j<n; j++) {
1573: r1 += (*aa++)*b1[*aj++];
1574: }
1575: c[colam + ridx[i]] += r1;
1576: }
1577: b1 += bm;
1578: }
1579: } else {
1580: for (col=0; col<cn-4; col += 4) { /* over columns of C */
1581: colam = col*am;
1582: for (i=0; i<am; i++) { /* over rows of C in those columns */
1583: r1 = r2 = r3 = r4 = 0.0;
1584: n = a->i[i+1] - a->i[i];
1585: aj = a->j + a->i[i];
1586: aa = a->a + a->i[i];
1587: for (j=0; j<n; j++) {
1588: r1 += (*aa)*b1[*aj];
1589: r2 += (*aa)*b2[*aj];
1590: r3 += (*aa)*b3[*aj];
1591: r4 += (*aa++)*b4[*aj++];
1592: }
1593: c[colam + i] += r1;
1594: c[colam + am + i] += r2;
1595: c[colam + am2 + i] += r3;
1596: c[colam + am3 + i] += r4;
1597: }
1598: b1 += bm4;
1599: b2 += bm4;
1600: b3 += bm4;
1601: b4 += bm4;
1602: }
1603: for (; col<cn; col++) { /* over extra columns of C */
1604: colam = col*am;
1605: for (i=0; i<am; i++) { /* over rows of C in those columns */
1606: r1 = 0.0;
1607: n = a->i[i+1] - a->i[i];
1608: aj = a->j + a->i[i];
1609: aa = a->a + a->i[i];
1611: for (j=0; j<n; j++) {
1612: r1 += (*aa++)*b1[*aj++];
1613: }
1614: c[colam + i] += r1;
1615: }
1616: b1 += bm;
1617: }
1618: }
1619: PetscLogFlops(cn*2.0*a->nz);
1620: MatDenseRestoreArray(C,&c);
1621: return(0);
1622: }
1624: PetscErrorCode MatTransColoringApplySpToDen_SeqAIJ(MatTransposeColoring coloring,Mat B,Mat Btdense)
1625: {
1627: Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data;
1628: Mat_SeqDense *btdense = (Mat_SeqDense*)Btdense->data;
1629: PetscInt *bi = b->i,*bj=b->j;
1630: PetscInt m = Btdense->rmap->n,n=Btdense->cmap->n,j,k,l,col,anz,*btcol,brow,ncolumns;
1631: MatScalar *btval,*btval_den,*ba=b->a;
1632: PetscInt *columns=coloring->columns,*colorforcol=coloring->colorforcol,ncolors=coloring->ncolors;
1635: btval_den=btdense->v;
1636: PetscMemzero(btval_den,(m*n)*sizeof(MatScalar));
1637: for (k=0; k<ncolors; k++) {
1638: ncolumns = coloring->ncolumns[k];
1639: for (l=0; l<ncolumns; l++) { /* insert a row of B to a column of Btdense */
1640: col = *(columns + colorforcol[k] + l);
1641: btcol = bj + bi[col];
1642: btval = ba + bi[col];
1643: anz = bi[col+1] - bi[col];
1644: for (j=0; j<anz; j++) {
1645: brow = btcol[j];
1646: btval_den[brow] = btval[j];
1647: }
1648: }
1649: btval_den += m;
1650: }
1651: return(0);
1652: }
1654: PetscErrorCode MatTransColoringApplyDenToSp_SeqAIJ(MatTransposeColoring matcoloring,Mat Cden,Mat Csp)
1655: {
1657: Mat_SeqAIJ *csp = (Mat_SeqAIJ*)Csp->data;
1658: PetscScalar *ca_den,*ca_den_ptr,*ca=csp->a;
1659: PetscInt k,l,m=Cden->rmap->n,ncolors=matcoloring->ncolors;
1660: PetscInt brows=matcoloring->brows,*den2sp=matcoloring->den2sp;
1661: PetscInt nrows,*row,*idx;
1662: PetscInt *rows=matcoloring->rows,*colorforrow=matcoloring->colorforrow;
1665: MatDenseGetArray(Cden,&ca_den);
1667: if (brows > 0) {
1668: PetscInt *lstart,row_end,row_start;
1669: lstart = matcoloring->lstart;
1670: PetscMemzero(lstart,ncolors*sizeof(PetscInt));
1672: row_end = brows;
1673: if (row_end > m) row_end = m;
1674: for (row_start=0; row_start<m; row_start+=brows) { /* loop over row blocks of Csp */
1675: ca_den_ptr = ca_den;
1676: for (k=0; k<ncolors; k++) { /* loop over colors (columns of Cden) */
1677: nrows = matcoloring->nrows[k];
1678: row = rows + colorforrow[k];
1679: idx = den2sp + colorforrow[k];
1680: for (l=lstart[k]; l<nrows; l++) {
1681: if (row[l] >= row_end) {
1682: lstart[k] = l;
1683: break;
1684: } else {
1685: ca[idx[l]] = ca_den_ptr[row[l]];
1686: }
1687: }
1688: ca_den_ptr += m;
1689: }
1690: row_end += brows;
1691: if (row_end > m) row_end = m;
1692: }
1693: } else { /* non-blocked impl: loop over columns of Csp - slow if Csp is large */
1694: ca_den_ptr = ca_den;
1695: for (k=0; k<ncolors; k++) {
1696: nrows = matcoloring->nrows[k];
1697: row = rows + colorforrow[k];
1698: idx = den2sp + colorforrow[k];
1699: for (l=0; l<nrows; l++) {
1700: ca[idx[l]] = ca_den_ptr[row[l]];
1701: }
1702: ca_den_ptr += m;
1703: }
1704: }
1706: MatDenseRestoreArray(Cden,&ca_den);
1707: #if defined(PETSC_USE_INFO)
1708: if (matcoloring->brows > 0) {
1709: PetscInfo1(Csp,"Loop over %D row blocks for den2sp\n",brows);
1710: } else {
1711: PetscInfo(Csp,"Loop over colors/columns of Cden, inefficient for large sparse matrix product \n");
1712: }
1713: #endif
1714: return(0);
1715: }
1717: PetscErrorCode MatTransposeColoringCreate_SeqAIJ(Mat mat,ISColoring iscoloring,MatTransposeColoring c)
1718: {
1720: PetscInt i,n,nrows,Nbs,j,k,m,ncols,col,cm;
1721: const PetscInt *is,*ci,*cj,*row_idx;
1722: PetscInt nis = iscoloring->n,*rowhit,bs = 1;
1723: IS *isa;
1724: Mat_SeqAIJ *csp = (Mat_SeqAIJ*)mat->data;
1725: PetscInt *colorforrow,*rows,*rows_i,*idxhit,*spidx,*den2sp,*den2sp_i;
1726: PetscInt *colorforcol,*columns,*columns_i,brows;
1727: PetscBool flg;
1730: ISColoringGetIS(iscoloring,PETSC_IGNORE,&isa);
1732: /* bs >1 is not being tested yet! */
1733: Nbs = mat->cmap->N/bs;
1734: c->M = mat->rmap->N/bs; /* set total rows, columns and local rows */
1735: c->N = Nbs;
1736: c->m = c->M;
1737: c->rstart = 0;
1738: c->brows = 100;
1740: c->ncolors = nis;
1741: PetscMalloc3(nis,&c->ncolumns,nis,&c->nrows,nis+1,&colorforrow);
1742: PetscMalloc1(csp->nz+1,&rows);
1743: PetscMalloc1(csp->nz+1,&den2sp);
1745: brows = c->brows;
1746: PetscOptionsGetInt(NULL,NULL,"-matden2sp_brows",&brows,&flg);
1747: if (flg) c->brows = brows;
1748: if (brows > 0) {
1749: PetscMalloc1(nis+1,&c->lstart);
1750: }
1752: colorforrow[0] = 0;
1753: rows_i = rows;
1754: den2sp_i = den2sp;
1756: PetscMalloc1(nis+1,&colorforcol);
1757: PetscMalloc1(Nbs+1,&columns);
1759: colorforcol[0] = 0;
1760: columns_i = columns;
1762: /* get column-wise storage of mat */
1763: MatGetColumnIJ_SeqAIJ_Color(mat,0,PETSC_FALSE,PETSC_FALSE,&ncols,&ci,&cj,&spidx,NULL);
1765: cm = c->m;
1766: PetscMalloc1(cm+1,&rowhit);
1767: PetscMalloc1(cm+1,&idxhit);
1768: for (i=0; i<nis; i++) { /* loop over color */
1769: ISGetLocalSize(isa[i],&n);
1770: ISGetIndices(isa[i],&is);
1772: c->ncolumns[i] = n;
1773: if (n) {
1774: PetscMemcpy(columns_i,is,n*sizeof(PetscInt));
1775: }
1776: colorforcol[i+1] = colorforcol[i] + n;
1777: columns_i += n;
1779: /* fast, crude version requires O(N*N) work */
1780: PetscMemzero(rowhit,cm*sizeof(PetscInt));
1782: for (j=0; j<n; j++) { /* loop over columns*/
1783: col = is[j];
1784: row_idx = cj + ci[col];
1785: m = ci[col+1] - ci[col];
1786: for (k=0; k<m; k++) { /* loop over columns marking them in rowhit */
1787: idxhit[*row_idx] = spidx[ci[col] + k];
1788: rowhit[*row_idx++] = col + 1;
1789: }
1790: }
1791: /* count the number of hits */
1792: nrows = 0;
1793: for (j=0; j<cm; j++) {
1794: if (rowhit[j]) nrows++;
1795: }
1796: c->nrows[i] = nrows;
1797: colorforrow[i+1] = colorforrow[i] + nrows;
1799: nrows = 0;
1800: for (j=0; j<cm; j++) { /* loop over rows */
1801: if (rowhit[j]) {
1802: rows_i[nrows] = j;
1803: den2sp_i[nrows] = idxhit[j];
1804: nrows++;
1805: }
1806: }
1807: den2sp_i += nrows;
1809: ISRestoreIndices(isa[i],&is);
1810: rows_i += nrows;
1811: }
1812: MatRestoreColumnIJ_SeqAIJ_Color(mat,0,PETSC_FALSE,PETSC_FALSE,&ncols,&ci,&cj,&spidx,NULL);
1813: PetscFree(rowhit);
1814: ISColoringRestoreIS(iscoloring,&isa);
1815: if (csp->nz != colorforrow[nis]) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"csp->nz %d != colorforrow[nis] %d",csp->nz,colorforrow[nis]);
1817: c->colorforrow = colorforrow;
1818: c->rows = rows;
1819: c->den2sp = den2sp;
1820: c->colorforcol = colorforcol;
1821: c->columns = columns;
1823: PetscFree(idxhit);
1824: return(0);
1825: }
1827: /* The combine method has done the symbolic and numeric in the first phase, and so we just return here */
1828: PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ_Combined(Mat A,Mat B,Mat C)
1829: {
1831: /* Empty function */
1832: return(0);
1833: }
1835: /* This algorithm combines the symbolic and numeric phase of matrix-matrix multiplication. */
1836: PetscErrorCode MatMatMult_SeqAIJ_SeqAIJ_Combined(Mat A,Mat B,PetscReal fill,Mat *C)
1837: {
1838: PetscErrorCode ierr;
1839: PetscLogDouble flops=0.0;
1840: Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b=(Mat_SeqAIJ*)B->data,*c;
1841: const PetscInt *ai=a->i,*bi=b->i;
1842: PetscInt *ci,*cj,*cj_i;
1843: PetscScalar *ca,*ca_i;
1844: PetscInt b_maxmemrow,c_maxmem,a_col;
1845: PetscInt am=A->rmap->N,bn=B->cmap->N,bm=B->rmap->N;
1846: PetscInt i,k,ndouble=0;
1847: PetscReal afill;
1848: PetscScalar *c_row_val_dense;
1849: PetscBool *c_row_idx_flags;
1850: PetscInt *aj_i=a->j;
1851: PetscScalar *aa_i=a->a;
1855: /* Step 1: Determine upper bounds on memory for C and allocate memory */
1856: /* This should be enough for almost all matrices. If still more memory is needed, it is reallocated later. */
1857: c_maxmem = 8*(ai[am]+bi[bm]);
1858: b_maxmemrow = PetscMin(bi[bm],bn);
1859: PetscMalloc1(am+1,&ci);
1860: PetscMalloc1(bn,&c_row_val_dense);
1861: PetscMalloc1(bn,&c_row_idx_flags);
1862: PetscMalloc1(c_maxmem,&cj);
1863: PetscMalloc1(c_maxmem,&ca);
1864: ca_i = ca;
1865: cj_i = cj;
1866: ci[0] = 0;
1867: PetscMemzero(c_row_val_dense,bn*sizeof(PetscScalar));
1868: PetscMemzero(c_row_idx_flags,bn*sizeof(PetscBool));
1869: for (i=0; i<am; i++) {
1870: /* Step 2: Initialize the dense row vector for C */
1871: const PetscInt anzi = ai[i+1] - ai[i]; /* number of nonzeros in this row of A, this is the number of rows of B that we merge */
1872: PetscInt cnzi = 0;
1873: PetscInt *bj_i;
1874: PetscScalar *ba_i;
1875: /* If the number of indices in C so far + the max number of columns in the next row > c_maxmem -> allocate more memory
1876: Usually, there is enough memory in the first place, so this is not executed. */
1877: while (ci[i] + b_maxmemrow > c_maxmem) {
1878: c_maxmem *= 2;
1879: ndouble++;
1880: PetscRealloc(sizeof(PetscInt)*c_maxmem,&cj);
1881: PetscRealloc(sizeof(PetscScalar)*c_maxmem,&ca);
1882: }
1884: /* Step 3: Do the numerical calculations */
1885: for (a_col=0; a_col<anzi; a_col++) { /* iterate over all non zero values in a row of A */
1886: PetscInt a_col_index = aj_i[a_col];
1887: const PetscInt bnzi = bi[a_col_index+1] - bi[a_col_index];
1888: flops += 2*bnzi;
1889: bj_i = b->j + bi[a_col_index]; /* points to the current row in bj */
1890: ba_i = b->a + bi[a_col_index]; /* points to the current row in ba */
1891: for (k=0; k<bnzi; ++k) { /* iterate over all non zeros of this row in B */
1892: if (c_row_idx_flags[bj_i[k]] == PETSC_FALSE) {
1893: cj_i[cnzi++] = bj_i[k];
1894: c_row_idx_flags[bj_i[k]] = PETSC_TRUE;
1895: }
1896: c_row_val_dense[bj_i[k]] += aa_i[a_col] * ba_i[k];
1897: }
1898: }
1900: /* Sort array */
1901: PetscSortInt(cnzi,cj_i);
1902: /* Step 4 */
1903: for (k=0; k<cnzi; k++) {
1904: ca_i[k] = c_row_val_dense[cj_i[k]];
1905: c_row_val_dense[cj_i[k]] = 0.;
1906: c_row_idx_flags[cj_i[k]] = PETSC_FALSE;
1907: }
1908: /* terminate current row */
1909: aa_i += anzi;
1910: aj_i += anzi;
1911: ca_i += cnzi;
1912: cj_i += cnzi;
1913: ci[i+1] = ci[i] + cnzi;
1914: flops += cnzi;
1915: }
1917: /* Step 5 */
1918: /* Create the new matrix */
1919: MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),am,bn,ci,cj,NULL,C);
1920: MatSetBlockSizesFromMats(*C,A,B);
1921: MatSetType(*C,((PetscObject)A)->type_name);
1923: /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
1924: /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
1925: c = (Mat_SeqAIJ*)((*C)->data);
1926: c->a = ca;
1927: c->free_a = PETSC_TRUE;
1928: c->free_ij = PETSC_TRUE;
1929: c->nonew = 0;
1931: /* set MatInfo */
1932: afill = (PetscReal)ci[am]/(ai[am]+bi[bm]) + 1.e-5;
1933: if (afill < 1.0) afill = 1.0;
1934: c->maxnz = ci[am];
1935: c->nz = ci[am];
1936: (*C)->info.mallocs = ndouble;
1937: (*C)->info.fill_ratio_given = fill;
1938: (*C)->info.fill_ratio_needed = afill;
1939: (*C)->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Combined;
1941: PetscFree(c_row_val_dense);
1942: PetscFree(c_row_idx_flags);
1944: MatAssemblyBegin(*C,MAT_FINAL_ASSEMBLY);
1945: MatAssemblyEnd(*C,MAT_FINAL_ASSEMBLY);
1946: PetscLogFlops(flops);
1947: return(0);
1948: }