Actual source code: matmatmult.c
1: /*
2: Defines matrix-matrix product routines for pairs of SeqAIJ matrices
3: C = A * B
4: */
6: #include <../src/mat/impls/aij/seq/aij.h>
7: #include <../src/mat/utils/freespace.h>
8: #include <petscbt.h>
9: #include <petsc/private/isimpl.h>
10: #include <../src/mat/impls/dense/seq/dense.h>
12: PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ(Mat A, Mat B, Mat C)
13: {
14: PetscFunctionBegin;
15: if (C->ops->matmultnumeric) PetscCall((*C->ops->matmultnumeric)(A, B, C));
16: else PetscCall(MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted(A, B, C));
17: PetscFunctionReturn(PETSC_SUCCESS);
18: }
20: /* Modified from MatCreateSeqAIJWithArrays() */
21: PETSC_INTERN PetscErrorCode MatSetSeqAIJWithArrays_private(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt i[], PetscInt j[], PetscScalar a[], MatType mtype, Mat mat)
22: {
23: PetscInt ii;
24: Mat_SeqAIJ *aij;
25: PetscBool isseqaij, osingle, ofree_a, ofree_ij;
27: PetscFunctionBegin;
28: PetscCheck(m <= 0 || !i[0], PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "i (row indices) must start with 0");
29: PetscCall(MatSetSizes(mat, m, n, m, n));
31: if (!mtype) {
32: PetscCall(PetscObjectBaseTypeCompare((PetscObject)mat, MATSEQAIJ, &isseqaij));
33: if (!isseqaij) PetscCall(MatSetType(mat, MATSEQAIJ));
34: } else {
35: PetscCall(MatSetType(mat, mtype));
36: }
38: aij = (Mat_SeqAIJ *)(mat)->data;
39: osingle = aij->singlemalloc;
40: ofree_a = aij->free_a;
41: ofree_ij = aij->free_ij;
42: /* changes the free flags */
43: PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(mat, MAT_SKIP_ALLOCATION, NULL));
45: PetscCall(PetscFree(aij->ilen));
46: PetscCall(PetscFree(aij->imax));
47: PetscCall(PetscMalloc1(m, &aij->imax));
48: PetscCall(PetscMalloc1(m, &aij->ilen));
49: for (ii = 0, aij->nonzerorowcnt = 0, aij->rmax = 0; ii < m; ii++) {
50: const PetscInt rnz = i[ii + 1] - i[ii];
51: aij->nonzerorowcnt += !!rnz;
52: aij->rmax = PetscMax(aij->rmax, rnz);
53: aij->ilen[ii] = aij->imax[ii] = i[ii + 1] - i[ii];
54: }
55: aij->maxnz = i[m];
56: aij->nz = i[m];
58: if (osingle) {
59: PetscCall(PetscFree3(aij->a, aij->j, aij->i));
60: } else {
61: if (ofree_a) PetscCall(PetscFree(aij->a));
62: if (ofree_ij) PetscCall(PetscFree(aij->j));
63: if (ofree_ij) PetscCall(PetscFree(aij->i));
64: }
65: aij->i = i;
66: aij->j = j;
67: aij->a = a;
68: aij->nonew = -1; /* this indicates that inserting a new value in the matrix that generates a new nonzero is an error */
69: /* default to not retain ownership */
70: aij->singlemalloc = PETSC_FALSE;
71: aij->free_a = PETSC_FALSE;
72: aij->free_ij = PETSC_FALSE;
73: PetscCall(MatCheckCompressedRow(mat, aij->nonzerorowcnt, &aij->compressedrow, aij->i, m, 0.6));
74: // Always build the diag info when i, j are set
75: PetscCall(MatMarkDiagonal_SeqAIJ(mat));
76: PetscFunctionReturn(PETSC_SUCCESS);
77: }
79: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ(Mat A, Mat B, PetscReal fill, Mat C)
80: {
81: Mat_Product *product = C->product;
82: MatProductAlgorithm alg;
83: PetscBool flg;
85: PetscFunctionBegin;
86: if (product) {
87: alg = product->alg;
88: } else {
89: alg = "sorted";
90: }
91: /* sorted */
92: PetscCall(PetscStrcmp(alg, "sorted", &flg));
93: if (flg) {
94: PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ_Sorted(A, B, fill, C));
95: PetscFunctionReturn(PETSC_SUCCESS);
96: }
98: /* scalable */
99: PetscCall(PetscStrcmp(alg, "scalable", &flg));
100: if (flg) {
101: PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable(A, B, fill, C));
102: PetscFunctionReturn(PETSC_SUCCESS);
103: }
105: /* scalable_fast */
106: PetscCall(PetscStrcmp(alg, "scalable_fast", &flg));
107: if (flg) {
108: PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable_fast(A, B, fill, C));
109: PetscFunctionReturn(PETSC_SUCCESS);
110: }
112: /* heap */
113: PetscCall(PetscStrcmp(alg, "heap", &flg));
114: if (flg) {
115: PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ_Heap(A, B, fill, C));
116: PetscFunctionReturn(PETSC_SUCCESS);
117: }
119: /* btheap */
120: PetscCall(PetscStrcmp(alg, "btheap", &flg));
121: if (flg) {
122: PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ_BTHeap(A, B, fill, C));
123: PetscFunctionReturn(PETSC_SUCCESS);
124: }
126: /* llcondensed */
127: PetscCall(PetscStrcmp(alg, "llcondensed", &flg));
128: if (flg) {
129: PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ_LLCondensed(A, B, fill, C));
130: PetscFunctionReturn(PETSC_SUCCESS);
131: }
133: /* rowmerge */
134: PetscCall(PetscStrcmp(alg, "rowmerge", &flg));
135: if (flg) {
136: PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ_RowMerge(A, B, fill, C));
137: PetscFunctionReturn(PETSC_SUCCESS);
138: }
140: #if defined(PETSC_HAVE_HYPRE)
141: PetscCall(PetscStrcmp(alg, "hypre", &flg));
142: if (flg) {
143: PetscCall(MatMatMultSymbolic_AIJ_AIJ_wHYPRE(A, B, fill, C));
144: PetscFunctionReturn(PETSC_SUCCESS);
145: }
146: #endif
148: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Mat Product Algorithm is not supported");
149: }
151: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_LLCondensed(Mat A, Mat B, PetscReal fill, Mat C)
152: {
153: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data, *c;
154: PetscInt *ai = a->i, *bi = b->i, *ci, *cj;
155: PetscInt am = A->rmap->N, bn = B->cmap->N, bm = B->rmap->N;
156: PetscReal afill;
157: PetscInt i, j, anzi, brow, bnzj, cnzi, *bj, *aj, *lnk, ndouble = 0, Crmax;
158: PetscHMapI ta;
159: PetscBT lnkbt;
160: PetscFreeSpaceList free_space = NULL, current_space = NULL;
162: PetscFunctionBegin;
163: /* Get ci and cj */
164: /* Allocate ci array, arrays for fill computation and */
165: /* free space for accumulating nonzero column info */
166: PetscCall(PetscMalloc1(am + 2, &ci));
167: ci[0] = 0;
169: /* create and initialize a linked list */
170: PetscCall(PetscHMapICreateWithSize(bn, &ta));
171: MatRowMergeMax_SeqAIJ(b, bm, ta);
172: PetscCall(PetscHMapIGetSize(ta, &Crmax));
173: PetscCall(PetscHMapIDestroy(&ta));
175: PetscCall(PetscLLCondensedCreate(Crmax, bn, &lnk, &lnkbt));
177: /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */
178: PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill, PetscIntSumTruncate(ai[am], bi[bm])), &free_space));
180: current_space = free_space;
182: /* Determine ci and cj */
183: for (i = 0; i < am; i++) {
184: anzi = ai[i + 1] - ai[i];
185: aj = a->j + ai[i];
186: for (j = 0; j < anzi; j++) {
187: brow = aj[j];
188: bnzj = bi[brow + 1] - bi[brow];
189: bj = b->j + bi[brow];
190: /* add non-zero cols of B into the sorted linked list lnk */
191: PetscCall(PetscLLCondensedAddSorted(bnzj, bj, lnk, lnkbt));
192: }
193: /* add possible missing diagonal entry */
194: if (C->force_diagonals) PetscCall(PetscLLCondensedAddSorted(1, &i, lnk, lnkbt));
195: cnzi = lnk[0];
197: /* If free space is not available, make more free space */
198: /* Double the amount of total space in the list */
199: if (current_space->local_remaining < cnzi) {
200: PetscCall(PetscFreeSpaceGet(PetscIntSumTruncate(cnzi, current_space->total_array_size), ¤t_space));
201: ndouble++;
202: }
204: /* Copy data into free space, then initialize lnk */
205: PetscCall(PetscLLCondensedClean(bn, cnzi, current_space->array, lnk, lnkbt));
207: current_space->array += cnzi;
208: current_space->local_used += cnzi;
209: current_space->local_remaining -= cnzi;
211: ci[i + 1] = ci[i] + cnzi;
212: }
214: /* Column indices are in the list of free space */
215: /* Allocate space for cj, initialize cj, and */
216: /* destroy list of free space and other temporary array(s) */
217: PetscCall(PetscMalloc1(ci[am] + 1, &cj));
218: PetscCall(PetscFreeSpaceContiguous(&free_space, cj));
219: PetscCall(PetscLLCondensedDestroy(lnk, lnkbt));
221: /* put together the new symbolic matrix */
222: PetscCall(MatSetSeqAIJWithArrays_private(PetscObjectComm((PetscObject)A), am, bn, ci, cj, NULL, ((PetscObject)A)->type_name, C));
223: PetscCall(MatSetBlockSizesFromMats(C, A, B));
225: /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
226: /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
227: c = (Mat_SeqAIJ *)C->data;
228: c->free_a = PETSC_FALSE;
229: c->free_ij = PETSC_TRUE;
230: c->nonew = 0;
232: /* fast, needs non-scalable O(bn) array 'abdense' */
233: C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted;
235: /* set MatInfo */
236: afill = (PetscReal)ci[am] / (ai[am] + bi[bm]) + 1.e-5;
237: if (afill < 1.0) afill = 1.0;
238: C->info.mallocs = ndouble;
239: C->info.fill_ratio_given = fill;
240: C->info.fill_ratio_needed = afill;
242: #if defined(PETSC_USE_INFO)
243: if (ci[am]) {
244: PetscCall(PetscInfo(C, "Reallocs %" PetscInt_FMT "; Fill ratio: given %g needed %g.\n", ndouble, (double)fill, (double)afill));
245: PetscCall(PetscInfo(C, "Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n", (double)afill));
246: } else {
247: PetscCall(PetscInfo(C, "Empty matrix product\n"));
248: }
249: #endif
250: PetscFunctionReturn(PETSC_SUCCESS);
251: }
253: PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted(Mat A, Mat B, Mat C)
254: {
255: PetscLogDouble flops = 0.0;
256: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
257: Mat_SeqAIJ *b = (Mat_SeqAIJ *)B->data;
258: Mat_SeqAIJ *c = (Mat_SeqAIJ *)C->data;
259: PetscInt *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j, *bjj, *ci = c->i, *cj = c->j;
260: PetscInt am = A->rmap->n, cm = C->rmap->n;
261: PetscInt i, j, k, anzi, bnzi, cnzi, brow;
262: PetscScalar *ca, valtmp;
263: PetscScalar *ab_dense;
264: PetscContainer cab_dense;
265: const PetscScalar *aa, *ba, *baj;
267: PetscFunctionBegin;
268: PetscCall(MatSeqAIJGetArrayRead(A, &aa));
269: PetscCall(MatSeqAIJGetArrayRead(B, &ba));
270: if (!c->a) { /* first call of MatMatMultNumeric_SeqAIJ_SeqAIJ, allocate ca and matmult_abdense */
271: PetscCall(PetscMalloc1(ci[cm] + 1, &ca));
272: c->a = ca;
273: c->free_a = PETSC_TRUE;
274: } else ca = c->a;
276: /* TODO this should be done in the symbolic phase */
277: /* However, this function is so heavily used (sometimes in an hidden way through multnumeric function pointers
278: that is hard to eradicate) */
279: PetscCall(PetscObjectQuery((PetscObject)C, "__PETSc__ab_dense", (PetscObject *)&cab_dense));
280: if (!cab_dense) {
281: PetscCall(PetscMalloc1(B->cmap->N, &ab_dense));
282: PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &cab_dense));
283: PetscCall(PetscContainerSetPointer(cab_dense, ab_dense));
284: PetscCall(PetscContainerSetUserDestroy(cab_dense, PetscContainerUserDestroyDefault));
285: PetscCall(PetscObjectCompose((PetscObject)C, "__PETSc__ab_dense", (PetscObject)cab_dense));
286: PetscCall(PetscObjectDereference((PetscObject)cab_dense));
287: }
288: PetscCall(PetscContainerGetPointer(cab_dense, (void **)&ab_dense));
289: PetscCall(PetscArrayzero(ab_dense, B->cmap->N));
291: /* clean old values in C */
292: PetscCall(PetscArrayzero(ca, ci[cm]));
293: /* Traverse A row-wise. */
294: /* Build the ith row in C by summing over nonzero columns in A, */
295: /* the rows of B corresponding to nonzeros of A. */
296: for (i = 0; i < am; i++) {
297: anzi = ai[i + 1] - ai[i];
298: for (j = 0; j < anzi; j++) {
299: brow = aj[j];
300: bnzi = bi[brow + 1] - bi[brow];
301: bjj = PetscSafePointerPlusOffset(bj, bi[brow]);
302: baj = PetscSafePointerPlusOffset(ba, bi[brow]);
303: /* perform dense axpy */
304: valtmp = aa[j];
305: for (k = 0; k < bnzi; k++) ab_dense[bjj[k]] += valtmp * baj[k];
306: flops += 2 * bnzi;
307: }
308: aj = PetscSafePointerPlusOffset(aj, anzi);
309: aa = PetscSafePointerPlusOffset(aa, anzi);
311: cnzi = ci[i + 1] - ci[i];
312: for (k = 0; k < cnzi; k++) {
313: ca[k] += ab_dense[cj[k]];
314: ab_dense[cj[k]] = 0.0; /* zero ab_dense */
315: }
316: flops += cnzi;
317: cj = PetscSafePointerPlusOffset(cj, cnzi);
318: ca += cnzi;
319: }
320: #if defined(PETSC_HAVE_DEVICE)
321: if (C->offloadmask != PETSC_OFFLOAD_UNALLOCATED) C->offloadmask = PETSC_OFFLOAD_CPU;
322: #endif
323: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
324: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
325: PetscCall(PetscLogFlops(flops));
326: PetscCall(MatSeqAIJRestoreArrayRead(A, &aa));
327: PetscCall(MatSeqAIJRestoreArrayRead(B, &ba));
328: PetscFunctionReturn(PETSC_SUCCESS);
329: }
331: PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable(Mat A, Mat B, Mat C)
332: {
333: PetscLogDouble flops = 0.0;
334: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
335: Mat_SeqAIJ *b = (Mat_SeqAIJ *)B->data;
336: Mat_SeqAIJ *c = (Mat_SeqAIJ *)C->data;
337: PetscInt *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j, *bjj, *ci = c->i, *cj = c->j;
338: PetscInt am = A->rmap->N, cm = C->rmap->N;
339: PetscInt i, j, k, anzi, bnzi, cnzi, brow;
340: PetscScalar *ca = c->a, valtmp;
341: const PetscScalar *aa, *ba, *baj;
342: PetscInt nextb;
344: PetscFunctionBegin;
345: PetscCall(MatSeqAIJGetArrayRead(A, &aa));
346: PetscCall(MatSeqAIJGetArrayRead(B, &ba));
347: if (!ca) { /* first call of MatMatMultNumeric_SeqAIJ_SeqAIJ, allocate ca and matmult_abdense */
348: PetscCall(PetscMalloc1(ci[cm] + 1, &ca));
349: c->a = ca;
350: c->free_a = PETSC_TRUE;
351: }
353: /* clean old values in C */
354: PetscCall(PetscArrayzero(ca, ci[cm]));
355: /* Traverse A row-wise. */
356: /* Build the ith row in C by summing over nonzero columns in A, */
357: /* the rows of B corresponding to nonzeros of A. */
358: for (i = 0; i < am; i++) {
359: anzi = ai[i + 1] - ai[i];
360: cnzi = ci[i + 1] - ci[i];
361: for (j = 0; j < anzi; j++) {
362: brow = aj[j];
363: bnzi = bi[brow + 1] - bi[brow];
364: bjj = bj + bi[brow];
365: baj = ba + bi[brow];
366: /* perform sparse axpy */
367: valtmp = aa[j];
368: nextb = 0;
369: for (k = 0; nextb < bnzi; k++) {
370: if (cj[k] == bjj[nextb]) { /* ccol == bcol */
371: ca[k] += valtmp * baj[nextb++];
372: }
373: }
374: flops += 2 * bnzi;
375: }
376: aj += anzi;
377: aa += anzi;
378: cj += cnzi;
379: ca += cnzi;
380: }
381: #if defined(PETSC_HAVE_DEVICE)
382: if (C->offloadmask != PETSC_OFFLOAD_UNALLOCATED) C->offloadmask = PETSC_OFFLOAD_CPU;
383: #endif
384: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
385: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
386: PetscCall(PetscLogFlops(flops));
387: PetscCall(MatSeqAIJRestoreArrayRead(A, &aa));
388: PetscCall(MatSeqAIJRestoreArrayRead(B, &ba));
389: PetscFunctionReturn(PETSC_SUCCESS);
390: }
392: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable_fast(Mat A, Mat B, PetscReal fill, Mat C)
393: {
394: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data, *c;
395: PetscInt *ai = a->i, *bi = b->i, *ci, *cj;
396: PetscInt am = A->rmap->N, bn = B->cmap->N, bm = B->rmap->N;
397: MatScalar *ca;
398: PetscReal afill;
399: PetscInt i, j, anzi, brow, bnzj, cnzi, *bj, *aj, *lnk, ndouble = 0, Crmax;
400: PetscHMapI ta;
401: PetscFreeSpaceList free_space = NULL, current_space = NULL;
403: PetscFunctionBegin;
404: /* Get ci and cj - same as MatMatMultSymbolic_SeqAIJ_SeqAIJ except using PetscLLxxx_fast() */
405: /* Allocate arrays for fill computation and free space for accumulating nonzero column */
406: PetscCall(PetscMalloc1(am + 2, &ci));
407: ci[0] = 0;
409: /* create and initialize a linked list */
410: PetscCall(PetscHMapICreateWithSize(bn, &ta));
411: MatRowMergeMax_SeqAIJ(b, bm, ta);
412: PetscCall(PetscHMapIGetSize(ta, &Crmax));
413: PetscCall(PetscHMapIDestroy(&ta));
415: PetscCall(PetscLLCondensedCreate_fast(Crmax, &lnk));
417: /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */
418: PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill, PetscIntSumTruncate(ai[am], bi[bm])), &free_space));
419: current_space = free_space;
421: /* Determine ci and cj */
422: for (i = 0; i < am; i++) {
423: anzi = ai[i + 1] - ai[i];
424: aj = a->j + ai[i];
425: for (j = 0; j < anzi; j++) {
426: brow = aj[j];
427: bnzj = bi[brow + 1] - bi[brow];
428: bj = b->j + bi[brow];
429: /* add non-zero cols of B into the sorted linked list lnk */
430: PetscCall(PetscLLCondensedAddSorted_fast(bnzj, bj, lnk));
431: }
432: /* add possible missing diagonal entry */
433: if (C->force_diagonals) PetscCall(PetscLLCondensedAddSorted_fast(1, &i, lnk));
434: cnzi = lnk[1];
436: /* If free space is not available, make more free space */
437: /* Double the amount of total space in the list */
438: if (current_space->local_remaining < cnzi) {
439: PetscCall(PetscFreeSpaceGet(PetscIntSumTruncate(cnzi, current_space->total_array_size), ¤t_space));
440: ndouble++;
441: }
443: /* Copy data into free space, then initialize lnk */
444: PetscCall(PetscLLCondensedClean_fast(cnzi, current_space->array, lnk));
446: current_space->array += cnzi;
447: current_space->local_used += cnzi;
448: current_space->local_remaining -= cnzi;
450: ci[i + 1] = ci[i] + cnzi;
451: }
453: /* Column indices are in the list of free space */
454: /* Allocate space for cj, initialize cj, and */
455: /* destroy list of free space and other temporary array(s) */
456: PetscCall(PetscMalloc1(ci[am] + 1, &cj));
457: PetscCall(PetscFreeSpaceContiguous(&free_space, cj));
458: PetscCall(PetscLLCondensedDestroy_fast(lnk));
460: /* Allocate space for ca */
461: PetscCall(PetscCalloc1(ci[am] + 1, &ca));
463: /* put together the new symbolic matrix */
464: PetscCall(MatSetSeqAIJWithArrays_private(PetscObjectComm((PetscObject)A), am, bn, ci, cj, ca, ((PetscObject)A)->type_name, C));
465: PetscCall(MatSetBlockSizesFromMats(C, A, B));
467: /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
468: /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
469: c = (Mat_SeqAIJ *)C->data;
470: c->free_a = PETSC_TRUE;
471: c->free_ij = PETSC_TRUE;
472: c->nonew = 0;
474: /* slower, less memory */
475: C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable;
477: /* set MatInfo */
478: afill = (PetscReal)ci[am] / (ai[am] + bi[bm]) + 1.e-5;
479: if (afill < 1.0) afill = 1.0;
480: C->info.mallocs = ndouble;
481: C->info.fill_ratio_given = fill;
482: C->info.fill_ratio_needed = afill;
484: #if defined(PETSC_USE_INFO)
485: if (ci[am]) {
486: PetscCall(PetscInfo(C, "Reallocs %" PetscInt_FMT "; Fill ratio: given %g needed %g.\n", ndouble, (double)fill, (double)afill));
487: PetscCall(PetscInfo(C, "Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n", (double)afill));
488: } else {
489: PetscCall(PetscInfo(C, "Empty matrix product\n"));
490: }
491: #endif
492: PetscFunctionReturn(PETSC_SUCCESS);
493: }
495: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_Scalable(Mat A, Mat B, PetscReal fill, Mat C)
496: {
497: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data, *c;
498: PetscInt *ai = a->i, *bi = b->i, *ci, *cj;
499: PetscInt am = A->rmap->N, bn = B->cmap->N, bm = B->rmap->N;
500: MatScalar *ca;
501: PetscReal afill;
502: PetscInt i, j, anzi, brow, bnzj, cnzi, *bj, *aj, *lnk, ndouble = 0, Crmax;
503: PetscHMapI ta;
504: PetscFreeSpaceList free_space = NULL, current_space = NULL;
506: PetscFunctionBegin;
507: /* Get ci and cj - same as MatMatMultSymbolic_SeqAIJ_SeqAIJ except using PetscLLxxx_Scalalbe() */
508: /* Allocate arrays for fill computation and free space for accumulating nonzero column */
509: PetscCall(PetscMalloc1(am + 2, &ci));
510: ci[0] = 0;
512: /* create and initialize a linked list */
513: PetscCall(PetscHMapICreateWithSize(bn, &ta));
514: MatRowMergeMax_SeqAIJ(b, bm, ta);
515: PetscCall(PetscHMapIGetSize(ta, &Crmax));
516: PetscCall(PetscHMapIDestroy(&ta));
517: PetscCall(PetscLLCondensedCreate_Scalable(Crmax, &lnk));
519: /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */
520: PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill, PetscIntSumTruncate(ai[am], bi[bm])), &free_space));
521: current_space = free_space;
523: /* Determine ci and cj */
524: for (i = 0; i < am; i++) {
525: anzi = ai[i + 1] - ai[i];
526: aj = a->j + ai[i];
527: for (j = 0; j < anzi; j++) {
528: brow = aj[j];
529: bnzj = bi[brow + 1] - bi[brow];
530: bj = b->j + bi[brow];
531: /* add non-zero cols of B into the sorted linked list lnk */
532: PetscCall(PetscLLCondensedAddSorted_Scalable(bnzj, bj, lnk));
533: }
534: /* add possible missing diagonal entry */
535: if (C->force_diagonals) PetscCall(PetscLLCondensedAddSorted_Scalable(1, &i, lnk));
537: cnzi = lnk[0];
539: /* If free space is not available, make more free space */
540: /* Double the amount of total space in the list */
541: if (current_space->local_remaining < cnzi) {
542: PetscCall(PetscFreeSpaceGet(PetscIntSumTruncate(cnzi, current_space->total_array_size), ¤t_space));
543: ndouble++;
544: }
546: /* Copy data into free space, then initialize lnk */
547: PetscCall(PetscLLCondensedClean_Scalable(cnzi, current_space->array, lnk));
549: current_space->array += cnzi;
550: current_space->local_used += cnzi;
551: current_space->local_remaining -= cnzi;
553: ci[i + 1] = ci[i] + cnzi;
554: }
556: /* Column indices are in the list of free space */
557: /* Allocate space for cj, initialize cj, and */
558: /* destroy list of free space and other temporary array(s) */
559: PetscCall(PetscMalloc1(ci[am] + 1, &cj));
560: PetscCall(PetscFreeSpaceContiguous(&free_space, cj));
561: PetscCall(PetscLLCondensedDestroy_Scalable(lnk));
563: /* Allocate space for ca */
564: PetscCall(PetscCalloc1(ci[am] + 1, &ca));
566: /* put together the new symbolic matrix */
567: PetscCall(MatSetSeqAIJWithArrays_private(PetscObjectComm((PetscObject)A), am, bn, ci, cj, ca, ((PetscObject)A)->type_name, C));
568: PetscCall(MatSetBlockSizesFromMats(C, A, B));
570: /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
571: /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
572: c = (Mat_SeqAIJ *)C->data;
573: c->free_a = PETSC_TRUE;
574: c->free_ij = PETSC_TRUE;
575: c->nonew = 0;
577: /* slower, less memory */
578: C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Scalable;
580: /* set MatInfo */
581: afill = (PetscReal)ci[am] / (ai[am] + bi[bm]) + 1.e-5;
582: if (afill < 1.0) afill = 1.0;
583: C->info.mallocs = ndouble;
584: C->info.fill_ratio_given = fill;
585: C->info.fill_ratio_needed = afill;
587: #if defined(PETSC_USE_INFO)
588: if (ci[am]) {
589: PetscCall(PetscInfo(C, "Reallocs %" PetscInt_FMT "; Fill ratio: given %g needed %g.\n", ndouble, (double)fill, (double)afill));
590: PetscCall(PetscInfo(C, "Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n", (double)afill));
591: } else {
592: PetscCall(PetscInfo(C, "Empty matrix product\n"));
593: }
594: #endif
595: PetscFunctionReturn(PETSC_SUCCESS);
596: }
598: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_Heap(Mat A, Mat B, PetscReal fill, Mat C)
599: {
600: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data, *c;
601: const PetscInt *ai = a->i, *bi = b->i, *aj = a->j, *bj = b->j;
602: PetscInt *ci, *cj, *bb;
603: PetscInt am = A->rmap->N, bn = B->cmap->N, bm = B->rmap->N;
604: PetscReal afill;
605: PetscInt i, j, col, ndouble = 0;
606: PetscFreeSpaceList free_space = NULL, current_space = NULL;
607: PetscHeap h;
609: PetscFunctionBegin;
610: /* Get ci and cj - by merging sorted rows using a heap */
611: /* Allocate arrays for fill computation and free space for accumulating nonzero column */
612: PetscCall(PetscMalloc1(am + 2, &ci));
613: ci[0] = 0;
615: /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */
616: PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill, PetscIntSumTruncate(ai[am], bi[bm])), &free_space));
617: current_space = free_space;
619: PetscCall(PetscHeapCreate(a->rmax, &h));
620: PetscCall(PetscMalloc1(a->rmax, &bb));
622: /* Determine ci and cj */
623: for (i = 0; i < am; i++) {
624: 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 */
625: const PetscInt *acol = aj + ai[i]; /* column indices of nonzero entries in this row */
626: ci[i + 1] = ci[i];
627: /* Populate the min heap */
628: for (j = 0; j < anzi; j++) {
629: bb[j] = bi[acol[j]]; /* bb points at the start of the row */
630: if (bb[j] < bi[acol[j] + 1]) { /* Add if row is nonempty */
631: PetscCall(PetscHeapAdd(h, j, bj[bb[j]++]));
632: }
633: }
634: /* Pick off the min element, adding it to free space */
635: PetscCall(PetscHeapPop(h, &j, &col));
636: while (j >= 0) {
637: if (current_space->local_remaining < 1) { /* double the size, but don't exceed 16 MiB */
638: PetscCall(PetscFreeSpaceGet(PetscMin(PetscIntMultTruncate(2, current_space->total_array_size), 16 << 20), ¤t_space));
639: ndouble++;
640: }
641: *(current_space->array++) = col;
642: current_space->local_used++;
643: current_space->local_remaining--;
644: ci[i + 1]++;
646: /* stash if anything else remains in this row of B */
647: if (bb[j] < bi[acol[j] + 1]) PetscCall(PetscHeapStash(h, j, bj[bb[j]++]));
648: while (1) { /* pop and stash any other rows of B that also had an entry in this column */
649: PetscInt j2, col2;
650: PetscCall(PetscHeapPeek(h, &j2, &col2));
651: if (col2 != col) break;
652: PetscCall(PetscHeapPop(h, &j2, &col2));
653: if (bb[j2] < bi[acol[j2] + 1]) PetscCall(PetscHeapStash(h, j2, bj[bb[j2]++]));
654: }
655: /* Put any stashed elements back into the min heap */
656: PetscCall(PetscHeapUnstash(h));
657: PetscCall(PetscHeapPop(h, &j, &col));
658: }
659: }
660: PetscCall(PetscFree(bb));
661: PetscCall(PetscHeapDestroy(&h));
663: /* Column indices are in the list of free space */
664: /* Allocate space for cj, initialize cj, and */
665: /* destroy list of free space and other temporary array(s) */
666: PetscCall(PetscMalloc1(ci[am], &cj));
667: PetscCall(PetscFreeSpaceContiguous(&free_space, cj));
669: /* put together the new symbolic matrix */
670: PetscCall(MatSetSeqAIJWithArrays_private(PetscObjectComm((PetscObject)A), am, bn, ci, cj, NULL, ((PetscObject)A)->type_name, C));
671: PetscCall(MatSetBlockSizesFromMats(C, A, B));
673: /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
674: /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
675: c = (Mat_SeqAIJ *)C->data;
676: c->free_a = PETSC_TRUE;
677: c->free_ij = PETSC_TRUE;
678: c->nonew = 0;
680: C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted;
682: /* set MatInfo */
683: afill = (PetscReal)ci[am] / (ai[am] + bi[bm]) + 1.e-5;
684: if (afill < 1.0) afill = 1.0;
685: C->info.mallocs = ndouble;
686: C->info.fill_ratio_given = fill;
687: C->info.fill_ratio_needed = afill;
689: #if defined(PETSC_USE_INFO)
690: if (ci[am]) {
691: PetscCall(PetscInfo(C, "Reallocs %" PetscInt_FMT "; Fill ratio: given %g needed %g.\n", ndouble, (double)fill, (double)afill));
692: PetscCall(PetscInfo(C, "Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n", (double)afill));
693: } else {
694: PetscCall(PetscInfo(C, "Empty matrix product\n"));
695: }
696: #endif
697: PetscFunctionReturn(PETSC_SUCCESS);
698: }
700: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_BTHeap(Mat A, Mat B, PetscReal fill, Mat C)
701: {
702: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data, *c;
703: const PetscInt *ai = a->i, *bi = b->i, *aj = a->j, *bj = b->j;
704: PetscInt *ci, *cj, *bb;
705: PetscInt am = A->rmap->N, bn = B->cmap->N, bm = B->rmap->N;
706: PetscReal afill;
707: PetscInt i, j, col, ndouble = 0;
708: PetscFreeSpaceList free_space = NULL, current_space = NULL;
709: PetscHeap h;
710: PetscBT bt;
712: PetscFunctionBegin;
713: /* Get ci and cj - using a heap for the sorted rows, but use BT so that each index is only added once */
714: /* Allocate arrays for fill computation and free space for accumulating nonzero column */
715: PetscCall(PetscMalloc1(am + 2, &ci));
716: ci[0] = 0;
718: /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */
719: PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(fill, PetscIntSumTruncate(ai[am], bi[bm])), &free_space));
721: current_space = free_space;
723: PetscCall(PetscHeapCreate(a->rmax, &h));
724: PetscCall(PetscMalloc1(a->rmax, &bb));
725: PetscCall(PetscBTCreate(bn, &bt));
727: /* Determine ci and cj */
728: for (i = 0; i < am; i++) {
729: 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 */
730: const PetscInt *acol = aj + ai[i]; /* column indices of nonzero entries in this row */
731: const PetscInt *fptr = current_space->array; /* Save beginning of the row so we can clear the BT later */
732: ci[i + 1] = ci[i];
733: /* Populate the min heap */
734: for (j = 0; j < anzi; j++) {
735: PetscInt brow = acol[j];
736: for (bb[j] = bi[brow]; bb[j] < bi[brow + 1]; bb[j]++) {
737: PetscInt bcol = bj[bb[j]];
738: if (!PetscBTLookupSet(bt, bcol)) { /* new entry */
739: PetscCall(PetscHeapAdd(h, j, bcol));
740: bb[j]++;
741: break;
742: }
743: }
744: }
745: /* Pick off the min element, adding it to free space */
746: PetscCall(PetscHeapPop(h, &j, &col));
747: while (j >= 0) {
748: if (current_space->local_remaining < 1) { /* double the size, but don't exceed 16 MiB */
749: fptr = NULL; /* need PetscBTMemzero */
750: PetscCall(PetscFreeSpaceGet(PetscMin(PetscIntMultTruncate(2, current_space->total_array_size), 16 << 20), ¤t_space));
751: ndouble++;
752: }
753: *(current_space->array++) = col;
754: current_space->local_used++;
755: current_space->local_remaining--;
756: ci[i + 1]++;
758: /* stash if anything else remains in this row of B */
759: for (; bb[j] < bi[acol[j] + 1]; bb[j]++) {
760: PetscInt bcol = bj[bb[j]];
761: if (!PetscBTLookupSet(bt, bcol)) { /* new entry */
762: PetscCall(PetscHeapAdd(h, j, bcol));
763: bb[j]++;
764: break;
765: }
766: }
767: PetscCall(PetscHeapPop(h, &j, &col));
768: }
769: if (fptr) { /* Clear the bits for this row */
770: for (; fptr < current_space->array; fptr++) PetscCall(PetscBTClear(bt, *fptr));
771: } else { /* We reallocated so we don't remember (easily) how to clear only the bits we changed */
772: PetscCall(PetscBTMemzero(bn, bt));
773: }
774: }
775: PetscCall(PetscFree(bb));
776: PetscCall(PetscHeapDestroy(&h));
777: PetscCall(PetscBTDestroy(&bt));
779: /* Column indices are in the list of free space */
780: /* Allocate space for cj, initialize cj, and */
781: /* destroy list of free space and other temporary array(s) */
782: PetscCall(PetscMalloc1(ci[am], &cj));
783: PetscCall(PetscFreeSpaceContiguous(&free_space, cj));
785: /* put together the new symbolic matrix */
786: PetscCall(MatSetSeqAIJWithArrays_private(PetscObjectComm((PetscObject)A), am, bn, ci, cj, NULL, ((PetscObject)A)->type_name, C));
787: PetscCall(MatSetBlockSizesFromMats(C, A, B));
789: /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
790: /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
791: c = (Mat_SeqAIJ *)C->data;
792: c->free_a = PETSC_TRUE;
793: c->free_ij = PETSC_TRUE;
794: c->nonew = 0;
796: C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted;
798: /* set MatInfo */
799: afill = (PetscReal)ci[am] / (ai[am] + bi[bm]) + 1.e-5;
800: if (afill < 1.0) afill = 1.0;
801: C->info.mallocs = ndouble;
802: C->info.fill_ratio_given = fill;
803: C->info.fill_ratio_needed = afill;
805: #if defined(PETSC_USE_INFO)
806: if (ci[am]) {
807: PetscCall(PetscInfo(C, "Reallocs %" PetscInt_FMT "; Fill ratio: given %g needed %g.\n", ndouble, (double)fill, (double)afill));
808: PetscCall(PetscInfo(C, "Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n", (double)afill));
809: } else {
810: PetscCall(PetscInfo(C, "Empty matrix product\n"));
811: }
812: #endif
813: PetscFunctionReturn(PETSC_SUCCESS);
814: }
816: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_RowMerge(Mat A, Mat B, PetscReal fill, Mat C)
817: {
818: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data, *c;
819: const PetscInt *ai = a->i, *bi = b->i, *aj = a->j, *bj = b->j, *inputi, *inputj, *inputcol, *inputcol_L1;
820: PetscInt *ci, *cj, *outputj, worki_L1[9], worki_L2[9];
821: PetscInt c_maxmem, a_maxrownnz = 0, a_rownnz;
822: const PetscInt workcol[8] = {0, 1, 2, 3, 4, 5, 6, 7};
823: const PetscInt am = A->rmap->N, bn = B->cmap->N, bm = B->rmap->N;
824: const PetscInt *brow_ptr[8], *brow_end[8];
825: PetscInt window[8];
826: PetscInt window_min, old_window_min, ci_nnz, outputi_nnz = 0, L1_nrows, L2_nrows;
827: PetscInt i, k, ndouble = 0, L1_rowsleft, rowsleft;
828: PetscReal afill;
829: PetscInt *workj_L1, *workj_L2, *workj_L3;
830: PetscInt L1_nnz, L2_nnz;
832: /* Step 1: Get upper bound on memory required for allocation.
833: Because of the way virtual memory works,
834: only the memory pages that are actually needed will be physically allocated. */
835: PetscFunctionBegin;
836: PetscCall(PetscMalloc1(am + 1, &ci));
837: for (i = 0; i < am; i++) {
838: 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 */
839: const PetscInt *acol = aj + ai[i]; /* column indices of nonzero entries in this row */
840: a_rownnz = 0;
841: for (k = 0; k < anzi; ++k) {
842: a_rownnz += bi[acol[k] + 1] - bi[acol[k]];
843: if (a_rownnz > bn) {
844: a_rownnz = bn;
845: break;
846: }
847: }
848: a_maxrownnz = PetscMax(a_maxrownnz, a_rownnz);
849: }
850: /* temporary work areas for merging rows */
851: PetscCall(PetscMalloc1(a_maxrownnz * 8, &workj_L1));
852: PetscCall(PetscMalloc1(a_maxrownnz * 8, &workj_L2));
853: PetscCall(PetscMalloc1(a_maxrownnz, &workj_L3));
855: /* This should be enough for almost all matrices. If not, memory is reallocated later. */
856: c_maxmem = 8 * (ai[am] + bi[bm]);
857: /* Step 2: Populate pattern for C */
858: PetscCall(PetscMalloc1(c_maxmem, &cj));
860: ci_nnz = 0;
861: ci[0] = 0;
862: worki_L1[0] = 0;
863: worki_L2[0] = 0;
864: for (i = 0; i < am; i++) {
865: 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 */
866: const PetscInt *acol = aj + ai[i]; /* column indices of nonzero entries in this row */
867: rowsleft = anzi;
868: inputcol_L1 = acol;
869: L2_nnz = 0;
870: L2_nrows = 1; /* Number of rows to be merged on Level 3. output of L3 already exists -> initial value 1 */
871: worki_L2[1] = 0;
872: outputi_nnz = 0;
874: /* If the number of indices in C so far + the max number of columns in the next row > c_maxmem -> allocate more memory */
875: while (ci_nnz + a_maxrownnz > c_maxmem) {
876: c_maxmem *= 2;
877: ndouble++;
878: PetscCall(PetscRealloc(sizeof(PetscInt) * c_maxmem, &cj));
879: }
881: while (rowsleft) {
882: L1_rowsleft = PetscMin(64, rowsleft); /* In the inner loop max 64 rows of B can be merged */
883: L1_nrows = 0;
884: L1_nnz = 0;
885: inputcol = inputcol_L1;
886: inputi = bi;
887: inputj = bj;
889: /* The following macro is used to specialize for small rows in A.
890: This helps with compiler unrolling, improving performance substantially.
891: Input: inputj inputi inputcol bn
892: Output: outputj outputi_nnz */
893: #define MatMatMultSymbolic_RowMergeMacro(ANNZ) \
894: do { \
895: window_min = bn; \
896: outputi_nnz = 0; \
897: for (k = 0; k < ANNZ; ++k) { \
898: brow_ptr[k] = inputj + inputi[inputcol[k]]; \
899: brow_end[k] = inputj + inputi[inputcol[k] + 1]; \
900: window[k] = (brow_ptr[k] != brow_end[k]) ? *brow_ptr[k] : bn; \
901: window_min = PetscMin(window[k], window_min); \
902: } \
903: while (window_min < bn) { \
904: outputj[outputi_nnz++] = window_min; \
905: /* advance front and compute new minimum */ \
906: old_window_min = window_min; \
907: window_min = bn; \
908: for (k = 0; k < ANNZ; ++k) { \
909: if (window[k] == old_window_min) { \
910: brow_ptr[k]++; \
911: window[k] = (brow_ptr[k] != brow_end[k]) ? *brow_ptr[k] : bn; \
912: } \
913: window_min = PetscMin(window[k], window_min); \
914: } \
915: } \
916: } while (0)
918: /************** L E V E L 1 ***************/
919: /* Merge up to 8 rows of B to L1 work array*/
920: while (L1_rowsleft) {
921: outputi_nnz = 0;
922: if (anzi > 8) outputj = workj_L1 + L1_nnz; /* Level 1 rowmerge*/
923: else outputj = cj + ci_nnz; /* Merge directly to C */
925: switch (L1_rowsleft) {
926: case 1:
927: brow_ptr[0] = inputj + inputi[inputcol[0]];
928: brow_end[0] = inputj + inputi[inputcol[0] + 1];
929: for (; brow_ptr[0] != brow_end[0]; ++brow_ptr[0]) outputj[outputi_nnz++] = *brow_ptr[0]; /* copy row in b over */
930: inputcol += L1_rowsleft;
931: rowsleft -= L1_rowsleft;
932: L1_rowsleft = 0;
933: break;
934: case 2:
935: MatMatMultSymbolic_RowMergeMacro(2);
936: inputcol += L1_rowsleft;
937: rowsleft -= L1_rowsleft;
938: L1_rowsleft = 0;
939: break;
940: case 3:
941: MatMatMultSymbolic_RowMergeMacro(3);
942: inputcol += L1_rowsleft;
943: rowsleft -= L1_rowsleft;
944: L1_rowsleft = 0;
945: break;
946: case 4:
947: MatMatMultSymbolic_RowMergeMacro(4);
948: inputcol += L1_rowsleft;
949: rowsleft -= L1_rowsleft;
950: L1_rowsleft = 0;
951: break;
952: case 5:
953: MatMatMultSymbolic_RowMergeMacro(5);
954: inputcol += L1_rowsleft;
955: rowsleft -= L1_rowsleft;
956: L1_rowsleft = 0;
957: break;
958: case 6:
959: MatMatMultSymbolic_RowMergeMacro(6);
960: inputcol += L1_rowsleft;
961: rowsleft -= L1_rowsleft;
962: L1_rowsleft = 0;
963: break;
964: case 7:
965: MatMatMultSymbolic_RowMergeMacro(7);
966: inputcol += L1_rowsleft;
967: rowsleft -= L1_rowsleft;
968: L1_rowsleft = 0;
969: break;
970: default:
971: MatMatMultSymbolic_RowMergeMacro(8);
972: inputcol += 8;
973: rowsleft -= 8;
974: L1_rowsleft -= 8;
975: break;
976: }
977: inputcol_L1 = inputcol;
978: L1_nnz += outputi_nnz;
979: worki_L1[++L1_nrows] = L1_nnz;
980: }
982: /********************** L E V E L 2 ************************/
983: /* Merge from L1 work array to either C or to L2 work array */
984: if (anzi > 8) {
985: inputi = worki_L1;
986: inputj = workj_L1;
987: inputcol = workcol;
988: outputi_nnz = 0;
990: if (anzi <= 64) outputj = cj + ci_nnz; /* Merge from L1 work array to C */
991: else outputj = workj_L2 + L2_nnz; /* Merge from L1 work array to L2 work array */
993: switch (L1_nrows) {
994: case 1:
995: brow_ptr[0] = inputj + inputi[inputcol[0]];
996: brow_end[0] = inputj + inputi[inputcol[0] + 1];
997: for (; brow_ptr[0] != brow_end[0]; ++brow_ptr[0]) outputj[outputi_nnz++] = *brow_ptr[0]; /* copy row in b over */
998: break;
999: case 2:
1000: MatMatMultSymbolic_RowMergeMacro(2);
1001: break;
1002: case 3:
1003: MatMatMultSymbolic_RowMergeMacro(3);
1004: break;
1005: case 4:
1006: MatMatMultSymbolic_RowMergeMacro(4);
1007: break;
1008: case 5:
1009: MatMatMultSymbolic_RowMergeMacro(5);
1010: break;
1011: case 6:
1012: MatMatMultSymbolic_RowMergeMacro(6);
1013: break;
1014: case 7:
1015: MatMatMultSymbolic_RowMergeMacro(7);
1016: break;
1017: case 8:
1018: MatMatMultSymbolic_RowMergeMacro(8);
1019: break;
1020: default:
1021: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "MatMatMult logic error: Not merging 1-8 rows from L1 work array!");
1022: }
1023: L2_nnz += outputi_nnz;
1024: worki_L2[++L2_nrows] = L2_nnz;
1026: /************************ L E V E L 3 **********************/
1027: /* Merge from L2 work array to either C or to L2 work array */
1028: if (anzi > 64 && (L2_nrows == 8 || rowsleft == 0)) {
1029: inputi = worki_L2;
1030: inputj = workj_L2;
1031: inputcol = workcol;
1032: outputi_nnz = 0;
1033: if (rowsleft) outputj = workj_L3;
1034: else outputj = cj + ci_nnz;
1035: switch (L2_nrows) {
1036: case 1:
1037: brow_ptr[0] = inputj + inputi[inputcol[0]];
1038: brow_end[0] = inputj + inputi[inputcol[0] + 1];
1039: for (; brow_ptr[0] != brow_end[0]; ++brow_ptr[0]) outputj[outputi_nnz++] = *brow_ptr[0]; /* copy row in b over */
1040: break;
1041: case 2:
1042: MatMatMultSymbolic_RowMergeMacro(2);
1043: break;
1044: case 3:
1045: MatMatMultSymbolic_RowMergeMacro(3);
1046: break;
1047: case 4:
1048: MatMatMultSymbolic_RowMergeMacro(4);
1049: break;
1050: case 5:
1051: MatMatMultSymbolic_RowMergeMacro(5);
1052: break;
1053: case 6:
1054: MatMatMultSymbolic_RowMergeMacro(6);
1055: break;
1056: case 7:
1057: MatMatMultSymbolic_RowMergeMacro(7);
1058: break;
1059: case 8:
1060: MatMatMultSymbolic_RowMergeMacro(8);
1061: break;
1062: default:
1063: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "MatMatMult logic error: Not merging 1-8 rows from L2 work array!");
1064: }
1065: L2_nrows = 1;
1066: L2_nnz = outputi_nnz;
1067: worki_L2[1] = outputi_nnz;
1068: /* Copy to workj_L2 */
1069: if (rowsleft) {
1070: for (k = 0; k < outputi_nnz; ++k) workj_L2[k] = outputj[k];
1071: }
1072: }
1073: }
1074: } /* while (rowsleft) */
1075: #undef MatMatMultSymbolic_RowMergeMacro
1077: /* terminate current row */
1078: ci_nnz += outputi_nnz;
1079: ci[i + 1] = ci_nnz;
1080: }
1082: /* Step 3: Create the new symbolic matrix */
1083: PetscCall(MatSetSeqAIJWithArrays_private(PetscObjectComm((PetscObject)A), am, bn, ci, cj, NULL, ((PetscObject)A)->type_name, C));
1084: PetscCall(MatSetBlockSizesFromMats(C, A, B));
1086: /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
1087: /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
1088: c = (Mat_SeqAIJ *)C->data;
1089: c->free_a = PETSC_TRUE;
1090: c->free_ij = PETSC_TRUE;
1091: c->nonew = 0;
1093: C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted;
1095: /* set MatInfo */
1096: afill = (PetscReal)ci[am] / (ai[am] + bi[bm]) + 1.e-5;
1097: if (afill < 1.0) afill = 1.0;
1098: C->info.mallocs = ndouble;
1099: C->info.fill_ratio_given = fill;
1100: C->info.fill_ratio_needed = afill;
1102: #if defined(PETSC_USE_INFO)
1103: if (ci[am]) {
1104: PetscCall(PetscInfo(C, "Reallocs %" PetscInt_FMT "; Fill ratio: given %g needed %g.\n", ndouble, (double)fill, (double)afill));
1105: PetscCall(PetscInfo(C, "Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n", (double)afill));
1106: } else {
1107: PetscCall(PetscInfo(C, "Empty matrix product\n"));
1108: }
1109: #endif
1111: /* Step 4: Free temporary work areas */
1112: PetscCall(PetscFree(workj_L1));
1113: PetscCall(PetscFree(workj_L2));
1114: PetscCall(PetscFree(workj_L3));
1115: PetscFunctionReturn(PETSC_SUCCESS);
1116: }
1118: /* concatenate unique entries and then sort */
1119: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqAIJ_Sorted(Mat A, Mat B, PetscReal fill, Mat C)
1120: {
1121: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data, *c;
1122: const PetscInt *ai = a->i, *bi = b->i, *aj = a->j, *bj = b->j;
1123: PetscInt *ci, *cj, bcol;
1124: PetscInt am = A->rmap->N, bn = B->cmap->N, bm = B->rmap->N;
1125: PetscReal afill;
1126: PetscInt i, j, ndouble = 0;
1127: PetscSegBuffer seg, segrow;
1128: char *seen;
1130: PetscFunctionBegin;
1131: PetscCall(PetscMalloc1(am + 1, &ci));
1132: ci[0] = 0;
1134: /* Initial FreeSpace size is fill*(nnz(A)+nnz(B)) */
1135: PetscCall(PetscSegBufferCreate(sizeof(PetscInt), (PetscInt)(fill * (ai[am] + bi[bm])), &seg));
1136: PetscCall(PetscSegBufferCreate(sizeof(PetscInt), 100, &segrow));
1137: PetscCall(PetscCalloc1(bn, &seen));
1139: /* Determine ci and cj */
1140: for (i = 0; i < am; i++) {
1141: 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 */
1142: const PetscInt *acol = PetscSafePointerPlusOffset(aj, ai[i]); /* column indices of nonzero entries in this row */
1143: PetscInt packlen = 0, *PETSC_RESTRICT crow;
1145: /* Pack segrow */
1146: for (j = 0; j < anzi; j++) {
1147: PetscInt brow = acol[j], bjstart = bi[brow], bjend = bi[brow + 1], k;
1148: for (k = bjstart; k < bjend; k++) {
1149: bcol = bj[k];
1150: if (!seen[bcol]) { /* new entry */
1151: PetscInt *PETSC_RESTRICT slot;
1152: PetscCall(PetscSegBufferGetInts(segrow, 1, &slot));
1153: *slot = bcol;
1154: seen[bcol] = 1;
1155: packlen++;
1156: }
1157: }
1158: }
1160: /* Check i-th diagonal entry */
1161: if (C->force_diagonals && !seen[i]) {
1162: PetscInt *PETSC_RESTRICT slot;
1163: PetscCall(PetscSegBufferGetInts(segrow, 1, &slot));
1164: *slot = i;
1165: seen[i] = 1;
1166: packlen++;
1167: }
1169: PetscCall(PetscSegBufferGetInts(seg, packlen, &crow));
1170: PetscCall(PetscSegBufferExtractTo(segrow, crow));
1171: PetscCall(PetscSortInt(packlen, crow));
1172: ci[i + 1] = ci[i] + packlen;
1173: for (j = 0; j < packlen; j++) seen[crow[j]] = 0;
1174: }
1175: PetscCall(PetscSegBufferDestroy(&segrow));
1176: PetscCall(PetscFree(seen));
1178: /* Column indices are in the segmented buffer */
1179: PetscCall(PetscSegBufferExtractAlloc(seg, &cj));
1180: PetscCall(PetscSegBufferDestroy(&seg));
1182: /* put together the new symbolic matrix */
1183: PetscCall(MatSetSeqAIJWithArrays_private(PetscObjectComm((PetscObject)A), am, bn, ci, cj, NULL, ((PetscObject)A)->type_name, C));
1184: PetscCall(MatSetBlockSizesFromMats(C, A, B));
1186: /* MatCreateSeqAIJWithArrays flags matrix so PETSc doesn't free the user's arrays. */
1187: /* These are PETSc arrays, so change flags so arrays can be deleted by PETSc */
1188: c = (Mat_SeqAIJ *)C->data;
1189: c->free_a = PETSC_TRUE;
1190: c->free_ij = PETSC_TRUE;
1191: c->nonew = 0;
1193: C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted;
1195: /* set MatInfo */
1196: afill = (PetscReal)ci[am] / PetscMax(ai[am] + bi[bm], 1) + 1.e-5;
1197: if (afill < 1.0) afill = 1.0;
1198: C->info.mallocs = ndouble;
1199: C->info.fill_ratio_given = fill;
1200: C->info.fill_ratio_needed = afill;
1202: #if defined(PETSC_USE_INFO)
1203: if (ci[am]) {
1204: PetscCall(PetscInfo(C, "Reallocs %" PetscInt_FMT "; Fill ratio: given %g needed %g.\n", ndouble, (double)fill, (double)afill));
1205: PetscCall(PetscInfo(C, "Use MatMatMult(A,B,MatReuse,%g,&C) for best performance.;\n", (double)afill));
1206: } else {
1207: PetscCall(PetscInfo(C, "Empty matrix product\n"));
1208: }
1209: #endif
1210: PetscFunctionReturn(PETSC_SUCCESS);
1211: }
1213: static PetscErrorCode MatDestroy_SeqAIJ_MatMatMultTrans(void *data)
1214: {
1215: Mat_MatMatTransMult *abt = (Mat_MatMatTransMult *)data;
1217: PetscFunctionBegin;
1218: PetscCall(MatTransposeColoringDestroy(&abt->matcoloring));
1219: PetscCall(MatDestroy(&abt->Bt_den));
1220: PetscCall(MatDestroy(&abt->ABt_den));
1221: PetscCall(PetscFree(abt));
1222: PetscFunctionReturn(PETSC_SUCCESS);
1223: }
1225: PetscErrorCode MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ(Mat A, Mat B, PetscReal fill, Mat C)
1226: {
1227: Mat Bt;
1228: Mat_MatMatTransMult *abt;
1229: Mat_Product *product = C->product;
1230: char *alg;
1232: PetscFunctionBegin;
1233: PetscCheck(product, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing product struct");
1234: PetscCheck(!product->data, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Extra product struct not empty");
1236: /* create symbolic Bt */
1237: PetscCall(MatTransposeSymbolic(B, &Bt));
1238: PetscCall(MatSetBlockSizes(Bt, PetscAbs(A->cmap->bs), PetscAbs(B->cmap->bs)));
1239: PetscCall(MatSetType(Bt, ((PetscObject)A)->type_name));
1241: /* get symbolic C=A*Bt */
1242: PetscCall(PetscStrallocpy(product->alg, &alg));
1243: PetscCall(MatProductSetAlgorithm(C, "sorted")); /* set algorithm for C = A*Bt */
1244: PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ(A, Bt, fill, C));
1245: PetscCall(MatProductSetAlgorithm(C, alg)); /* resume original algorithm for ABt product */
1246: PetscCall(PetscFree(alg));
1248: /* create a supporting struct for reuse intermediate dense matrices with matcoloring */
1249: PetscCall(PetscNew(&abt));
1251: product->data = abt;
1252: product->destroy = MatDestroy_SeqAIJ_MatMatMultTrans;
1254: C->ops->mattransposemultnumeric = MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ;
1256: abt->usecoloring = PETSC_FALSE;
1257: PetscCall(PetscStrcmp(product->alg, "color", &abt->usecoloring));
1258: if (abt->usecoloring) {
1259: /* Create MatTransposeColoring from symbolic C=A*B^T */
1260: MatTransposeColoring matcoloring;
1261: MatColoring coloring;
1262: ISColoring iscoloring;
1263: Mat Bt_dense, C_dense;
1265: /* inode causes memory problem */
1266: PetscCall(MatSetOption(C, MAT_USE_INODES, PETSC_FALSE));
1268: PetscCall(MatColoringCreate(C, &coloring));
1269: PetscCall(MatColoringSetDistance(coloring, 2));
1270: PetscCall(MatColoringSetType(coloring, MATCOLORINGSL));
1271: PetscCall(MatColoringSetFromOptions(coloring));
1272: PetscCall(MatColoringApply(coloring, &iscoloring));
1273: PetscCall(MatColoringDestroy(&coloring));
1274: PetscCall(MatTransposeColoringCreate(C, iscoloring, &matcoloring));
1276: abt->matcoloring = matcoloring;
1278: PetscCall(ISColoringDestroy(&iscoloring));
1280: /* Create Bt_dense and C_dense = A*Bt_dense */
1281: PetscCall(MatCreate(PETSC_COMM_SELF, &Bt_dense));
1282: PetscCall(MatSetSizes(Bt_dense, A->cmap->n, matcoloring->ncolors, A->cmap->n, matcoloring->ncolors));
1283: PetscCall(MatSetType(Bt_dense, MATSEQDENSE));
1284: PetscCall(MatSeqDenseSetPreallocation(Bt_dense, NULL));
1286: Bt_dense->assembled = PETSC_TRUE;
1287: abt->Bt_den = Bt_dense;
1289: PetscCall(MatCreate(PETSC_COMM_SELF, &C_dense));
1290: PetscCall(MatSetSizes(C_dense, A->rmap->n, matcoloring->ncolors, A->rmap->n, matcoloring->ncolors));
1291: PetscCall(MatSetType(C_dense, MATSEQDENSE));
1292: PetscCall(MatSeqDenseSetPreallocation(C_dense, NULL));
1294: Bt_dense->assembled = PETSC_TRUE;
1295: abt->ABt_den = C_dense;
1297: #if defined(PETSC_USE_INFO)
1298: {
1299: Mat_SeqAIJ *c = (Mat_SeqAIJ *)C->data;
1300: PetscCall(PetscInfo(C, "Use coloring of C=A*B^T; B^T: %" PetscInt_FMT " %" PetscInt_FMT ", Bt_dense: %" PetscInt_FMT ",%" PetscInt_FMT "; Cnz %" PetscInt_FMT " / (cm*ncolors %" PetscInt_FMT ") = %g\n", B->cmap->n, B->rmap->n, Bt_dense->rmap->n,
1301: Bt_dense->cmap->n, c->nz, A->rmap->n * matcoloring->ncolors, (double)(((PetscReal)c->nz) / ((PetscReal)(A->rmap->n * matcoloring->ncolors)))));
1302: }
1303: #endif
1304: }
1305: /* clean up */
1306: PetscCall(MatDestroy(&Bt));
1307: PetscFunctionReturn(PETSC_SUCCESS);
1308: }
1310: PetscErrorCode MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ(Mat A, Mat B, Mat C)
1311: {
1312: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data, *c = (Mat_SeqAIJ *)C->data;
1313: PetscInt *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j, anzi, bnzj, nexta, nextb, *acol, *bcol, brow;
1314: PetscInt cm = C->rmap->n, *ci = c->i, *cj = c->j, i, j, cnzi, *ccol;
1315: PetscLogDouble flops = 0.0;
1316: MatScalar *aa = a->a, *aval, *ba = b->a, *bval, *ca, *cval;
1317: Mat_MatMatTransMult *abt;
1318: Mat_Product *product = C->product;
1320: PetscFunctionBegin;
1321: PetscCheck(product, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing product struct");
1322: abt = (Mat_MatMatTransMult *)product->data;
1323: PetscCheck(abt, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing product struct");
1324: /* clear old values in C */
1325: if (!c->a) {
1326: PetscCall(PetscCalloc1(ci[cm] + 1, &ca));
1327: c->a = ca;
1328: c->free_a = PETSC_TRUE;
1329: } else {
1330: ca = c->a;
1331: PetscCall(PetscArrayzero(ca, ci[cm] + 1));
1332: }
1334: if (abt->usecoloring) {
1335: MatTransposeColoring matcoloring = abt->matcoloring;
1336: Mat Bt_dense, C_dense = abt->ABt_den;
1338: /* Get Bt_dense by Apply MatTransposeColoring to B */
1339: Bt_dense = abt->Bt_den;
1340: PetscCall(MatTransColoringApplySpToDen(matcoloring, B, Bt_dense));
1342: /* C_dense = A*Bt_dense */
1343: PetscCall(MatMatMultNumeric_SeqAIJ_SeqDense(A, Bt_dense, C_dense));
1345: /* Recover C from C_dense */
1346: PetscCall(MatTransColoringApplyDenToSp(matcoloring, C_dense, C));
1347: PetscFunctionReturn(PETSC_SUCCESS);
1348: }
1350: for (i = 0; i < cm; i++) {
1351: anzi = ai[i + 1] - ai[i];
1352: acol = PetscSafePointerPlusOffset(aj, ai[i]);
1353: aval = PetscSafePointerPlusOffset(aa, ai[i]);
1354: cnzi = ci[i + 1] - ci[i];
1355: ccol = PetscSafePointerPlusOffset(cj, ci[i]);
1356: cval = ca + ci[i];
1357: for (j = 0; j < cnzi; j++) {
1358: brow = ccol[j];
1359: bnzj = bi[brow + 1] - bi[brow];
1360: bcol = bj + bi[brow];
1361: bval = ba + bi[brow];
1363: /* perform sparse inner-product c(i,j)=A[i,:]*B[j,:]^T */
1364: nexta = 0;
1365: nextb = 0;
1366: while (nexta < anzi && nextb < bnzj) {
1367: while (nexta < anzi && acol[nexta] < bcol[nextb]) nexta++;
1368: if (nexta == anzi) break;
1369: while (nextb < bnzj && acol[nexta] > bcol[nextb]) nextb++;
1370: if (nextb == bnzj) break;
1371: if (acol[nexta] == bcol[nextb]) {
1372: cval[j] += aval[nexta] * bval[nextb];
1373: nexta++;
1374: nextb++;
1375: flops += 2;
1376: }
1377: }
1378: }
1379: }
1380: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
1381: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
1382: PetscCall(PetscLogFlops(flops));
1383: PetscFunctionReturn(PETSC_SUCCESS);
1384: }
1386: PetscErrorCode MatDestroy_SeqAIJ_MatTransMatMult(void *data)
1387: {
1388: Mat_MatTransMatMult *atb = (Mat_MatTransMatMult *)data;
1390: PetscFunctionBegin;
1391: PetscCall(MatDestroy(&atb->At));
1392: if (atb->destroy) PetscCall((*atb->destroy)(atb->data));
1393: PetscCall(PetscFree(atb));
1394: PetscFunctionReturn(PETSC_SUCCESS);
1395: }
1397: PetscErrorCode MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ(Mat A, Mat B, PetscReal fill, Mat C)
1398: {
1399: Mat At = NULL;
1400: Mat_Product *product = C->product;
1401: PetscBool flg, def, square;
1403: PetscFunctionBegin;
1404: MatCheckProduct(C, 4);
1405: square = (PetscBool)(A == B && A->symmetric == PETSC_BOOL3_TRUE);
1406: /* outerproduct */
1407: PetscCall(PetscStrcmp(product->alg, "outerproduct", &flg));
1408: if (flg) {
1409: /* create symbolic At */
1410: if (!square) {
1411: PetscCall(MatTransposeSymbolic(A, &At));
1412: PetscCall(MatSetBlockSizes(At, PetscAbs(A->cmap->bs), PetscAbs(B->cmap->bs)));
1413: PetscCall(MatSetType(At, ((PetscObject)A)->type_name));
1414: }
1415: /* get symbolic C=At*B */
1416: PetscCall(MatProductSetAlgorithm(C, "sorted"));
1417: PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ(square ? A : At, B, fill, C));
1419: /* clean up */
1420: if (!square) PetscCall(MatDestroy(&At));
1422: C->ops->mattransposemultnumeric = MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ; /* outerproduct */
1423: PetscCall(MatProductSetAlgorithm(C, "outerproduct"));
1424: PetscFunctionReturn(PETSC_SUCCESS);
1425: }
1427: /* matmatmult */
1428: PetscCall(PetscStrcmp(product->alg, "default", &def));
1429: PetscCall(PetscStrcmp(product->alg, "at*b", &flg));
1430: if (flg || def) {
1431: Mat_MatTransMatMult *atb;
1433: PetscCheck(!product->data, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Extra product struct not empty");
1434: PetscCall(PetscNew(&atb));
1435: if (!square) PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &At));
1436: PetscCall(MatProductSetAlgorithm(C, "sorted"));
1437: PetscCall(MatMatMultSymbolic_SeqAIJ_SeqAIJ(square ? A : At, B, fill, C));
1438: PetscCall(MatProductSetAlgorithm(C, "at*b"));
1439: product->data = atb;
1440: product->destroy = MatDestroy_SeqAIJ_MatTransMatMult;
1441: atb->At = At;
1443: C->ops->mattransposemultnumeric = NULL; /* see MatProductNumeric_AtB_SeqAIJ_SeqAIJ */
1444: PetscFunctionReturn(PETSC_SUCCESS);
1445: }
1447: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Mat Product Algorithm is not supported");
1448: }
1450: PetscErrorCode MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ(Mat A, Mat B, Mat C)
1451: {
1452: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data, *c = (Mat_SeqAIJ *)C->data;
1453: PetscInt am = A->rmap->n, anzi, *ai = a->i, *aj = a->j, *bi = b->i, *bj, bnzi, nextb;
1454: PetscInt cm = C->rmap->n, *ci = c->i, *cj = c->j, crow, *cjj, i, j, k;
1455: PetscLogDouble flops = 0.0;
1456: MatScalar *aa = a->a, *ba, *ca, *caj;
1458: PetscFunctionBegin;
1459: if (!c->a) {
1460: PetscCall(PetscCalloc1(ci[cm] + 1, &ca));
1462: c->a = ca;
1463: c->free_a = PETSC_TRUE;
1464: } else {
1465: ca = c->a;
1466: PetscCall(PetscArrayzero(ca, ci[cm]));
1467: }
1469: /* compute A^T*B using outer product (A^T)[:,i]*B[i,:] */
1470: for (i = 0; i < am; i++) {
1471: bj = b->j + bi[i];
1472: ba = b->a + bi[i];
1473: bnzi = bi[i + 1] - bi[i];
1474: anzi = ai[i + 1] - ai[i];
1475: for (j = 0; j < anzi; j++) {
1476: nextb = 0;
1477: crow = *aj++;
1478: cjj = cj + ci[crow];
1479: caj = ca + ci[crow];
1480: /* perform sparse axpy operation. Note cjj includes bj. */
1481: for (k = 0; nextb < bnzi; k++) {
1482: if (cjj[k] == *(bj + nextb)) { /* ccol == bcol */
1483: caj[k] += (*aa) * (*(ba + nextb));
1484: nextb++;
1485: }
1486: }
1487: flops += 2 * bnzi;
1488: aa++;
1489: }
1490: }
1492: /* Assemble the final matrix and clean up */
1493: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
1494: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
1495: PetscCall(PetscLogFlops(flops));
1496: PetscFunctionReturn(PETSC_SUCCESS);
1497: }
1499: PetscErrorCode MatMatMultSymbolic_SeqAIJ_SeqDense(Mat A, Mat B, PetscReal fill, Mat C)
1500: {
1501: PetscFunctionBegin;
1502: PetscCall(MatMatMultSymbolic_SeqDense_SeqDense(A, B, 0.0, C));
1503: C->ops->matmultnumeric = MatMatMultNumeric_SeqAIJ_SeqDense;
1504: PetscFunctionReturn(PETSC_SUCCESS);
1505: }
1507: PETSC_INTERN PetscErrorCode MatMatMultNumericAdd_SeqAIJ_SeqDense(Mat A, Mat B, Mat C, const PetscBool add)
1508: {
1509: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
1510: PetscScalar *c, r1, r2, r3, r4, *c1, *c2, *c3, *c4;
1511: const PetscScalar *aa, *b, *b1, *b2, *b3, *b4, *av;
1512: const PetscInt *aj;
1513: PetscInt cm = C->rmap->n, cn = B->cmap->n, bm, am = A->rmap->n;
1514: PetscInt clda;
1515: PetscInt am4, bm4, col, i, j, n;
1517: PetscFunctionBegin;
1518: if (!cm || !cn) PetscFunctionReturn(PETSC_SUCCESS);
1519: PetscCall(MatSeqAIJGetArrayRead(A, &av));
1520: if (add) {
1521: PetscCall(MatDenseGetArray(C, &c));
1522: } else {
1523: PetscCall(MatDenseGetArrayWrite(C, &c));
1524: }
1525: PetscCall(MatDenseGetArrayRead(B, &b));
1526: PetscCall(MatDenseGetLDA(B, &bm));
1527: PetscCall(MatDenseGetLDA(C, &clda));
1528: am4 = 4 * clda;
1529: bm4 = 4 * bm;
1530: if (b) {
1531: b1 = b;
1532: b2 = b1 + bm;
1533: b3 = b2 + bm;
1534: b4 = b3 + bm;
1535: } else b1 = b2 = b3 = b4 = NULL;
1536: c1 = c;
1537: c2 = c1 + clda;
1538: c3 = c2 + clda;
1539: c4 = c3 + clda;
1540: for (col = 0; col < (cn / 4) * 4; col += 4) { /* over columns of C */
1541: for (i = 0; i < am; i++) { /* over rows of A in those columns */
1542: r1 = r2 = r3 = r4 = 0.0;
1543: n = a->i[i + 1] - a->i[i];
1544: aj = PetscSafePointerPlusOffset(a->j, a->i[i]);
1545: aa = PetscSafePointerPlusOffset(av, a->i[i]);
1546: for (j = 0; j < n; j++) {
1547: const PetscScalar aatmp = aa[j];
1548: const PetscInt ajtmp = aj[j];
1549: r1 += aatmp * b1[ajtmp];
1550: r2 += aatmp * b2[ajtmp];
1551: r3 += aatmp * b3[ajtmp];
1552: r4 += aatmp * b4[ajtmp];
1553: }
1554: if (add) {
1555: c1[i] += r1;
1556: c2[i] += r2;
1557: c3[i] += r3;
1558: c4[i] += r4;
1559: } else {
1560: c1[i] = r1;
1561: c2[i] = r2;
1562: c3[i] = r3;
1563: c4[i] = r4;
1564: }
1565: }
1566: if (b) {
1567: b1 += bm4;
1568: b2 += bm4;
1569: b3 += bm4;
1570: b4 += bm4;
1571: }
1572: c1 += am4;
1573: c2 += am4;
1574: c3 += am4;
1575: c4 += am4;
1576: }
1577: /* process remaining columns */
1578: if (col != cn) {
1579: PetscInt rc = cn - col;
1581: if (rc == 1) {
1582: for (i = 0; i < am; i++) {
1583: r1 = 0.0;
1584: n = a->i[i + 1] - a->i[i];
1585: aj = PetscSafePointerPlusOffset(a->j, a->i[i]);
1586: aa = PetscSafePointerPlusOffset(av, a->i[i]);
1587: for (j = 0; j < n; j++) r1 += aa[j] * b1[aj[j]];
1588: if (add) c1[i] += r1;
1589: else c1[i] = r1;
1590: }
1591: } else if (rc == 2) {
1592: for (i = 0; i < am; i++) {
1593: r1 = r2 = 0.0;
1594: n = a->i[i + 1] - a->i[i];
1595: aj = PetscSafePointerPlusOffset(a->j, a->i[i]);
1596: aa = PetscSafePointerPlusOffset(av, a->i[i]);
1597: for (j = 0; j < n; j++) {
1598: const PetscScalar aatmp = aa[j];
1599: const PetscInt ajtmp = aj[j];
1600: r1 += aatmp * b1[ajtmp];
1601: r2 += aatmp * b2[ajtmp];
1602: }
1603: if (add) {
1604: c1[i] += r1;
1605: c2[i] += r2;
1606: } else {
1607: c1[i] = r1;
1608: c2[i] = r2;
1609: }
1610: }
1611: } else {
1612: for (i = 0; i < am; i++) {
1613: r1 = r2 = r3 = 0.0;
1614: n = a->i[i + 1] - a->i[i];
1615: aj = PetscSafePointerPlusOffset(a->j, a->i[i]);
1616: aa = PetscSafePointerPlusOffset(av, a->i[i]);
1617: for (j = 0; j < n; j++) {
1618: const PetscScalar aatmp = aa[j];
1619: const PetscInt ajtmp = aj[j];
1620: r1 += aatmp * b1[ajtmp];
1621: r2 += aatmp * b2[ajtmp];
1622: r3 += aatmp * b3[ajtmp];
1623: }
1624: if (add) {
1625: c1[i] += r1;
1626: c2[i] += r2;
1627: c3[i] += r3;
1628: } else {
1629: c1[i] = r1;
1630: c2[i] = r2;
1631: c3[i] = r3;
1632: }
1633: }
1634: }
1635: }
1636: PetscCall(PetscLogFlops(cn * (2.0 * a->nz)));
1637: if (add) {
1638: PetscCall(MatDenseRestoreArray(C, &c));
1639: } else {
1640: PetscCall(MatDenseRestoreArrayWrite(C, &c));
1641: }
1642: PetscCall(MatDenseRestoreArrayRead(B, &b));
1643: PetscCall(MatSeqAIJRestoreArrayRead(A, &av));
1644: PetscFunctionReturn(PETSC_SUCCESS);
1645: }
1647: PetscErrorCode MatMatMultNumeric_SeqAIJ_SeqDense(Mat A, Mat B, Mat C)
1648: {
1649: PetscFunctionBegin;
1650: PetscCheck(B->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number columns in A %" PetscInt_FMT " not equal rows in B %" PetscInt_FMT, A->cmap->n, B->rmap->n);
1651: PetscCheck(A->rmap->n == C->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number rows in C %" PetscInt_FMT " not equal rows in A %" PetscInt_FMT, C->rmap->n, A->rmap->n);
1652: PetscCheck(B->cmap->n == C->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number columns in B %" PetscInt_FMT " not equal columns in C %" PetscInt_FMT, B->cmap->n, C->cmap->n);
1654: PetscCall(MatMatMultNumericAdd_SeqAIJ_SeqDense(A, B, C, PETSC_FALSE));
1655: PetscFunctionReturn(PETSC_SUCCESS);
1656: }
1658: static PetscErrorCode MatProductSetFromOptions_SeqAIJ_SeqDense_AB(Mat C)
1659: {
1660: PetscFunctionBegin;
1661: C->ops->matmultsymbolic = MatMatMultSymbolic_SeqAIJ_SeqDense;
1662: C->ops->productsymbolic = MatProductSymbolic_AB;
1663: PetscFunctionReturn(PETSC_SUCCESS);
1664: }
1666: PETSC_INTERN PetscErrorCode MatTMatTMultSymbolic_SeqAIJ_SeqDense(Mat, Mat, PetscReal, Mat);
1668: static PetscErrorCode MatProductSetFromOptions_SeqAIJ_SeqDense_AtB(Mat C)
1669: {
1670: PetscFunctionBegin;
1671: C->ops->transposematmultsymbolic = MatTMatTMultSymbolic_SeqAIJ_SeqDense;
1672: C->ops->productsymbolic = MatProductSymbolic_AtB;
1673: PetscFunctionReturn(PETSC_SUCCESS);
1674: }
1676: static PetscErrorCode MatProductSetFromOptions_SeqAIJ_SeqDense_ABt(Mat C)
1677: {
1678: PetscFunctionBegin;
1679: C->ops->mattransposemultsymbolic = MatTMatTMultSymbolic_SeqAIJ_SeqDense;
1680: C->ops->productsymbolic = MatProductSymbolic_ABt;
1681: PetscFunctionReturn(PETSC_SUCCESS);
1682: }
1684: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqAIJ_SeqDense(Mat C)
1685: {
1686: Mat_Product *product = C->product;
1688: PetscFunctionBegin;
1689: switch (product->type) {
1690: case MATPRODUCT_AB:
1691: PetscCall(MatProductSetFromOptions_SeqAIJ_SeqDense_AB(C));
1692: break;
1693: case MATPRODUCT_AtB:
1694: PetscCall(MatProductSetFromOptions_SeqAIJ_SeqDense_AtB(C));
1695: break;
1696: case MATPRODUCT_ABt:
1697: PetscCall(MatProductSetFromOptions_SeqAIJ_SeqDense_ABt(C));
1698: break;
1699: default:
1700: break;
1701: }
1702: PetscFunctionReturn(PETSC_SUCCESS);
1703: }
1705: static PetscErrorCode MatProductSetFromOptions_SeqXBAIJ_SeqDense_AB(Mat C)
1706: {
1707: Mat_Product *product = C->product;
1708: Mat A = product->A;
1709: PetscBool baij;
1711: PetscFunctionBegin;
1712: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQBAIJ, &baij));
1713: if (!baij) { /* A is seqsbaij */
1714: PetscBool sbaij;
1715: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQSBAIJ, &sbaij));
1716: PetscCheck(sbaij, PetscObjectComm((PetscObject)C), PETSC_ERR_ARG_WRONGSTATE, "Mat must be either seqbaij or seqsbaij format");
1718: C->ops->matmultsymbolic = MatMatMultSymbolic_SeqSBAIJ_SeqDense;
1719: } else { /* A is seqbaij */
1720: C->ops->matmultsymbolic = MatMatMultSymbolic_SeqBAIJ_SeqDense;
1721: }
1723: C->ops->productsymbolic = MatProductSymbolic_AB;
1724: PetscFunctionReturn(PETSC_SUCCESS);
1725: }
1727: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqXBAIJ_SeqDense(Mat C)
1728: {
1729: Mat_Product *product = C->product;
1731: PetscFunctionBegin;
1732: MatCheckProduct(C, 1);
1733: PetscCheck(product->A, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing A");
1734: if (product->type == MATPRODUCT_AB || (product->type == MATPRODUCT_AtB && product->A->symmetric == PETSC_BOOL3_TRUE)) PetscCall(MatProductSetFromOptions_SeqXBAIJ_SeqDense_AB(C));
1735: PetscFunctionReturn(PETSC_SUCCESS);
1736: }
1738: static PetscErrorCode MatProductSetFromOptions_SeqDense_SeqAIJ_AB(Mat C)
1739: {
1740: PetscFunctionBegin;
1741: C->ops->matmultsymbolic = MatMatMultSymbolic_SeqDense_SeqAIJ;
1742: C->ops->productsymbolic = MatProductSymbolic_AB;
1743: PetscFunctionReturn(PETSC_SUCCESS);
1744: }
1746: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_SeqDense_SeqAIJ(Mat C)
1747: {
1748: Mat_Product *product = C->product;
1750: PetscFunctionBegin;
1751: if (product->type == MATPRODUCT_AB) PetscCall(MatProductSetFromOptions_SeqDense_SeqAIJ_AB(C));
1752: PetscFunctionReturn(PETSC_SUCCESS);
1753: }
1755: PetscErrorCode MatTransColoringApplySpToDen_SeqAIJ(MatTransposeColoring coloring, Mat B, Mat Btdense)
1756: {
1757: Mat_SeqAIJ *b = (Mat_SeqAIJ *)B->data;
1758: Mat_SeqDense *btdense = (Mat_SeqDense *)Btdense->data;
1759: PetscInt *bi = b->i, *bj = b->j;
1760: PetscInt m = Btdense->rmap->n, n = Btdense->cmap->n, j, k, l, col, anz, *btcol, brow, ncolumns;
1761: MatScalar *btval, *btval_den, *ba = b->a;
1762: PetscInt *columns = coloring->columns, *colorforcol = coloring->colorforcol, ncolors = coloring->ncolors;
1764: PetscFunctionBegin;
1765: btval_den = btdense->v;
1766: PetscCall(PetscArrayzero(btval_den, m * n));
1767: for (k = 0; k < ncolors; k++) {
1768: ncolumns = coloring->ncolumns[k];
1769: for (l = 0; l < ncolumns; l++) { /* insert a row of B to a column of Btdense */
1770: col = *(columns + colorforcol[k] + l);
1771: btcol = bj + bi[col];
1772: btval = ba + bi[col];
1773: anz = bi[col + 1] - bi[col];
1774: for (j = 0; j < anz; j++) {
1775: brow = btcol[j];
1776: btval_den[brow] = btval[j];
1777: }
1778: }
1779: btval_den += m;
1780: }
1781: PetscFunctionReturn(PETSC_SUCCESS);
1782: }
1784: PetscErrorCode MatTransColoringApplyDenToSp_SeqAIJ(MatTransposeColoring matcoloring, Mat Cden, Mat Csp)
1785: {
1786: Mat_SeqAIJ *csp = (Mat_SeqAIJ *)Csp->data;
1787: const PetscScalar *ca_den, *ca_den_ptr;
1788: PetscScalar *ca = csp->a;
1789: PetscInt k, l, m = Cden->rmap->n, ncolors = matcoloring->ncolors;
1790: PetscInt brows = matcoloring->brows, *den2sp = matcoloring->den2sp;
1791: PetscInt nrows, *row, *idx;
1792: PetscInt *rows = matcoloring->rows, *colorforrow = matcoloring->colorforrow;
1794: PetscFunctionBegin;
1795: PetscCall(MatDenseGetArrayRead(Cden, &ca_den));
1797: if (brows > 0) {
1798: PetscInt *lstart, row_end, row_start;
1799: lstart = matcoloring->lstart;
1800: PetscCall(PetscArrayzero(lstart, ncolors));
1802: row_end = brows;
1803: if (row_end > m) row_end = m;
1804: for (row_start = 0; row_start < m; row_start += brows) { /* loop over row blocks of Csp */
1805: ca_den_ptr = ca_den;
1806: for (k = 0; k < ncolors; k++) { /* loop over colors (columns of Cden) */
1807: nrows = matcoloring->nrows[k];
1808: row = rows + colorforrow[k];
1809: idx = den2sp + colorforrow[k];
1810: for (l = lstart[k]; l < nrows; l++) {
1811: if (row[l] >= row_end) {
1812: lstart[k] = l;
1813: break;
1814: } else {
1815: ca[idx[l]] = ca_den_ptr[row[l]];
1816: }
1817: }
1818: ca_den_ptr += m;
1819: }
1820: row_end += brows;
1821: if (row_end > m) row_end = m;
1822: }
1823: } else { /* non-blocked impl: loop over columns of Csp - slow if Csp is large */
1824: ca_den_ptr = ca_den;
1825: for (k = 0; k < ncolors; k++) {
1826: nrows = matcoloring->nrows[k];
1827: row = rows + colorforrow[k];
1828: idx = den2sp + colorforrow[k];
1829: for (l = 0; l < nrows; l++) ca[idx[l]] = ca_den_ptr[row[l]];
1830: ca_den_ptr += m;
1831: }
1832: }
1834: PetscCall(MatDenseRestoreArrayRead(Cden, &ca_den));
1835: #if defined(PETSC_USE_INFO)
1836: if (matcoloring->brows > 0) {
1837: PetscCall(PetscInfo(Csp, "Loop over %" PetscInt_FMT " row blocks for den2sp\n", brows));
1838: } else {
1839: PetscCall(PetscInfo(Csp, "Loop over colors/columns of Cden, inefficient for large sparse matrix product \n"));
1840: }
1841: #endif
1842: PetscFunctionReturn(PETSC_SUCCESS);
1843: }
1845: PetscErrorCode MatTransposeColoringCreate_SeqAIJ(Mat mat, ISColoring iscoloring, MatTransposeColoring c)
1846: {
1847: PetscInt i, n, nrows, Nbs, j, k, m, ncols, col, cm;
1848: const PetscInt *is, *ci, *cj, *row_idx;
1849: PetscInt nis = iscoloring->n, *rowhit, bs = 1;
1850: IS *isa;
1851: Mat_SeqAIJ *csp = (Mat_SeqAIJ *)mat->data;
1852: PetscInt *colorforrow, *rows, *rows_i, *idxhit, *spidx, *den2sp, *den2sp_i;
1853: PetscInt *colorforcol, *columns, *columns_i, brows;
1854: PetscBool flg;
1856: PetscFunctionBegin;
1857: PetscCall(ISColoringGetIS(iscoloring, PETSC_USE_POINTER, PETSC_IGNORE, &isa));
1859: /* bs >1 is not being tested yet! */
1860: Nbs = mat->cmap->N / bs;
1861: c->M = mat->rmap->N / bs; /* set total rows, columns and local rows */
1862: c->N = Nbs;
1863: c->m = c->M;
1864: c->rstart = 0;
1865: c->brows = 100;
1867: c->ncolors = nis;
1868: PetscCall(PetscMalloc3(nis, &c->ncolumns, nis, &c->nrows, nis + 1, &colorforrow));
1869: PetscCall(PetscMalloc1(csp->nz + 1, &rows));
1870: PetscCall(PetscMalloc1(csp->nz + 1, &den2sp));
1872: brows = c->brows;
1873: PetscCall(PetscOptionsGetInt(NULL, NULL, "-matden2sp_brows", &brows, &flg));
1874: if (flg) c->brows = brows;
1875: if (brows > 0) PetscCall(PetscMalloc1(nis + 1, &c->lstart));
1877: colorforrow[0] = 0;
1878: rows_i = rows;
1879: den2sp_i = den2sp;
1881: PetscCall(PetscMalloc1(nis + 1, &colorforcol));
1882: PetscCall(PetscMalloc1(Nbs + 1, &columns));
1884: colorforcol[0] = 0;
1885: columns_i = columns;
1887: /* get column-wise storage of mat */
1888: PetscCall(MatGetColumnIJ_SeqAIJ_Color(mat, 0, PETSC_FALSE, PETSC_FALSE, &ncols, &ci, &cj, &spidx, NULL));
1890: cm = c->m;
1891: PetscCall(PetscMalloc1(cm + 1, &rowhit));
1892: PetscCall(PetscMalloc1(cm + 1, &idxhit));
1893: for (i = 0; i < nis; i++) { /* loop over color */
1894: PetscCall(ISGetLocalSize(isa[i], &n));
1895: PetscCall(ISGetIndices(isa[i], &is));
1897: c->ncolumns[i] = n;
1898: if (n) PetscCall(PetscArraycpy(columns_i, is, n));
1899: colorforcol[i + 1] = colorforcol[i] + n;
1900: columns_i += n;
1902: /* fast, crude version requires O(N*N) work */
1903: PetscCall(PetscArrayzero(rowhit, cm));
1905: for (j = 0; j < n; j++) { /* loop over columns*/
1906: col = is[j];
1907: row_idx = cj + ci[col];
1908: m = ci[col + 1] - ci[col];
1909: for (k = 0; k < m; k++) { /* loop over columns marking them in rowhit */
1910: idxhit[*row_idx] = spidx[ci[col] + k];
1911: rowhit[*row_idx++] = col + 1;
1912: }
1913: }
1914: /* count the number of hits */
1915: nrows = 0;
1916: for (j = 0; j < cm; j++) {
1917: if (rowhit[j]) nrows++;
1918: }
1919: c->nrows[i] = nrows;
1920: colorforrow[i + 1] = colorforrow[i] + nrows;
1922: nrows = 0;
1923: for (j = 0; j < cm; j++) { /* loop over rows */
1924: if (rowhit[j]) {
1925: rows_i[nrows] = j;
1926: den2sp_i[nrows] = idxhit[j];
1927: nrows++;
1928: }
1929: }
1930: den2sp_i += nrows;
1932: PetscCall(ISRestoreIndices(isa[i], &is));
1933: rows_i += nrows;
1934: }
1935: PetscCall(MatRestoreColumnIJ_SeqAIJ_Color(mat, 0, PETSC_FALSE, PETSC_FALSE, &ncols, &ci, &cj, &spidx, NULL));
1936: PetscCall(PetscFree(rowhit));
1937: PetscCall(ISColoringRestoreIS(iscoloring, PETSC_USE_POINTER, &isa));
1938: PetscCheck(csp->nz == colorforrow[nis], PETSC_COMM_SELF, PETSC_ERR_PLIB, "csp->nz %" PetscInt_FMT " != colorforrow[nis] %" PetscInt_FMT, csp->nz, colorforrow[nis]);
1940: c->colorforrow = colorforrow;
1941: c->rows = rows;
1942: c->den2sp = den2sp;
1943: c->colorforcol = colorforcol;
1944: c->columns = columns;
1946: PetscCall(PetscFree(idxhit));
1947: PetscFunctionReturn(PETSC_SUCCESS);
1948: }
1950: static PetscErrorCode MatProductNumeric_AtB_SeqAIJ_SeqAIJ(Mat C)
1951: {
1952: Mat_Product *product = C->product;
1953: Mat A = product->A, B = product->B;
1955: PetscFunctionBegin;
1956: if (C->ops->mattransposemultnumeric) {
1957: /* Alg: "outerproduct" */
1958: PetscCall((*C->ops->mattransposemultnumeric)(A, B, C));
1959: } else {
1960: /* Alg: "matmatmult" -- C = At*B */
1961: Mat_MatTransMatMult *atb = (Mat_MatTransMatMult *)product->data;
1963: PetscCheck(atb, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing product struct");
1964: if (atb->At) {
1965: /* At is computed in MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ();
1966: user may have called MatProductReplaceMats() to get this A=product->A */
1967: PetscCall(MatTransposeSetPrecursor(A, atb->At));
1968: PetscCall(MatTranspose(A, MAT_REUSE_MATRIX, &atb->At));
1969: }
1970: PetscCall(MatMatMultNumeric_SeqAIJ_SeqAIJ(atb->At ? atb->At : A, B, C));
1971: }
1972: PetscFunctionReturn(PETSC_SUCCESS);
1973: }
1975: static PetscErrorCode MatProductSymbolic_AtB_SeqAIJ_SeqAIJ(Mat C)
1976: {
1977: Mat_Product *product = C->product;
1978: Mat A = product->A, B = product->B;
1979: PetscReal fill = product->fill;
1981: PetscFunctionBegin;
1982: PetscCall(MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ(A, B, fill, C));
1984: C->ops->productnumeric = MatProductNumeric_AtB_SeqAIJ_SeqAIJ;
1985: PetscFunctionReturn(PETSC_SUCCESS);
1986: }
1988: static PetscErrorCode MatProductSetFromOptions_SeqAIJ_AB(Mat C)
1989: {
1990: Mat_Product *product = C->product;
1991: PetscInt alg = 0; /* default algorithm */
1992: PetscBool flg = PETSC_FALSE;
1993: #if !defined(PETSC_HAVE_HYPRE)
1994: const char *algTypes[7] = {"sorted", "scalable", "scalable_fast", "heap", "btheap", "llcondensed", "rowmerge"};
1995: PetscInt nalg = 7;
1996: #else
1997: const char *algTypes[8] = {"sorted", "scalable", "scalable_fast", "heap", "btheap", "llcondensed", "rowmerge", "hypre"};
1998: PetscInt nalg = 8;
1999: #endif
2001: PetscFunctionBegin;
2002: /* Set default algorithm */
2003: PetscCall(PetscStrcmp(C->product->alg, "default", &flg));
2004: if (flg) PetscCall(MatProductSetAlgorithm(C, (MatProductAlgorithm)algTypes[alg]));
2006: /* Get runtime option */
2007: if (product->api_user) {
2008: PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatMatMult", "Mat");
2009: PetscCall(PetscOptionsEList("-matmatmult_via", "Algorithmic approach", "MatMatMult", algTypes, nalg, algTypes[0], &alg, &flg));
2010: PetscOptionsEnd();
2011: } else {
2012: PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatProduct_AB", "Mat");
2013: PetscCall(PetscOptionsEList("-mat_product_algorithm", "Algorithmic approach", "MatProduct_AB", algTypes, nalg, algTypes[0], &alg, &flg));
2014: PetscOptionsEnd();
2015: }
2016: if (flg) PetscCall(MatProductSetAlgorithm(C, (MatProductAlgorithm)algTypes[alg]));
2018: C->ops->productsymbolic = MatProductSymbolic_AB;
2019: C->ops->matmultsymbolic = MatMatMultSymbolic_SeqAIJ_SeqAIJ;
2020: PetscFunctionReturn(PETSC_SUCCESS);
2021: }
2023: static PetscErrorCode MatProductSetFromOptions_SeqAIJ_AtB(Mat C)
2024: {
2025: Mat_Product *product = C->product;
2026: PetscInt alg = 0; /* default algorithm */
2027: PetscBool flg = PETSC_FALSE;
2028: const char *algTypes[3] = {"default", "at*b", "outerproduct"};
2029: PetscInt nalg = 3;
2031: PetscFunctionBegin;
2032: /* Get runtime option */
2033: if (product->api_user) {
2034: PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatTransposeMatMult", "Mat");
2035: PetscCall(PetscOptionsEList("-mattransposematmult_via", "Algorithmic approach", "MatTransposeMatMult", algTypes, nalg, algTypes[alg], &alg, &flg));
2036: PetscOptionsEnd();
2037: } else {
2038: PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatProduct_AtB", "Mat");
2039: PetscCall(PetscOptionsEList("-mat_product_algorithm", "Algorithmic approach", "MatProduct_AtB", algTypes, nalg, algTypes[alg], &alg, &flg));
2040: PetscOptionsEnd();
2041: }
2042: if (flg) PetscCall(MatProductSetAlgorithm(C, (MatProductAlgorithm)algTypes[alg]));
2044: C->ops->productsymbolic = MatProductSymbolic_AtB_SeqAIJ_SeqAIJ;
2045: PetscFunctionReturn(PETSC_SUCCESS);
2046: }
2048: static PetscErrorCode MatProductSetFromOptions_SeqAIJ_ABt(Mat C)
2049: {
2050: Mat_Product *product = C->product;
2051: PetscInt alg = 0; /* default algorithm */
2052: PetscBool flg = PETSC_FALSE;
2053: const char *algTypes[2] = {"default", "color"};
2054: PetscInt nalg = 2;
2056: PetscFunctionBegin;
2057: /* Set default algorithm */
2058: PetscCall(PetscStrcmp(C->product->alg, "default", &flg));
2059: if (!flg) {
2060: alg = 1;
2061: PetscCall(MatProductSetAlgorithm(C, (MatProductAlgorithm)algTypes[alg]));
2062: }
2064: /* Get runtime option */
2065: if (product->api_user) {
2066: PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatMatTransposeMult", "Mat");
2067: PetscCall(PetscOptionsEList("-matmattransmult_via", "Algorithmic approach", "MatMatTransposeMult", algTypes, nalg, algTypes[alg], &alg, &flg));
2068: PetscOptionsEnd();
2069: } else {
2070: PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatProduct_ABt", "Mat");
2071: PetscCall(PetscOptionsEList("-mat_product_algorithm", "Algorithmic approach", "MatProduct_ABt", algTypes, nalg, algTypes[alg], &alg, &flg));
2072: PetscOptionsEnd();
2073: }
2074: if (flg) PetscCall(MatProductSetAlgorithm(C, (MatProductAlgorithm)algTypes[alg]));
2076: C->ops->mattransposemultsymbolic = MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ;
2077: C->ops->productsymbolic = MatProductSymbolic_ABt;
2078: PetscFunctionReturn(PETSC_SUCCESS);
2079: }
2081: static PetscErrorCode MatProductSetFromOptions_SeqAIJ_PtAP(Mat C)
2082: {
2083: Mat_Product *product = C->product;
2084: PetscBool flg = PETSC_FALSE;
2085: PetscInt alg = 0; /* default algorithm -- alg=1 should be default!!! */
2086: #if !defined(PETSC_HAVE_HYPRE)
2087: const char *algTypes[2] = {"scalable", "rap"};
2088: PetscInt nalg = 2;
2089: #else
2090: const char *algTypes[3] = {"scalable", "rap", "hypre"};
2091: PetscInt nalg = 3;
2092: #endif
2094: PetscFunctionBegin;
2095: /* Set default algorithm */
2096: PetscCall(PetscStrcmp(product->alg, "default", &flg));
2097: if (flg) PetscCall(MatProductSetAlgorithm(C, (MatProductAlgorithm)algTypes[alg]));
2099: /* Get runtime option */
2100: if (product->api_user) {
2101: PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatPtAP", "Mat");
2102: PetscCall(PetscOptionsEList("-matptap_via", "Algorithmic approach", "MatPtAP", algTypes, nalg, algTypes[0], &alg, &flg));
2103: PetscOptionsEnd();
2104: } else {
2105: PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatProduct_PtAP", "Mat");
2106: PetscCall(PetscOptionsEList("-mat_product_algorithm", "Algorithmic approach", "MatProduct_PtAP", algTypes, nalg, algTypes[0], &alg, &flg));
2107: PetscOptionsEnd();
2108: }
2109: if (flg) PetscCall(MatProductSetAlgorithm(C, (MatProductAlgorithm)algTypes[alg]));
2111: C->ops->productsymbolic = MatProductSymbolic_PtAP_SeqAIJ_SeqAIJ;
2112: PetscFunctionReturn(PETSC_SUCCESS);
2113: }
2115: static PetscErrorCode MatProductSetFromOptions_SeqAIJ_RARt(Mat C)
2116: {
2117: Mat_Product *product = C->product;
2118: PetscBool flg = PETSC_FALSE;
2119: PetscInt alg = 0; /* default algorithm */
2120: const char *algTypes[3] = {"r*a*rt", "r*art", "coloring_rart"};
2121: PetscInt nalg = 3;
2123: PetscFunctionBegin;
2124: /* Set default algorithm */
2125: PetscCall(PetscStrcmp(product->alg, "default", &flg));
2126: if (flg) PetscCall(MatProductSetAlgorithm(C, (MatProductAlgorithm)algTypes[alg]));
2128: /* Get runtime option */
2129: if (product->api_user) {
2130: PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatRARt", "Mat");
2131: PetscCall(PetscOptionsEList("-matrart_via", "Algorithmic approach", "MatRARt", algTypes, nalg, algTypes[0], &alg, &flg));
2132: PetscOptionsEnd();
2133: } else {
2134: PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatProduct_RARt", "Mat");
2135: PetscCall(PetscOptionsEList("-mat_product_algorithm", "Algorithmic approach", "MatProduct_RARt", algTypes, nalg, algTypes[0], &alg, &flg));
2136: PetscOptionsEnd();
2137: }
2138: if (flg) PetscCall(MatProductSetAlgorithm(C, (MatProductAlgorithm)algTypes[alg]));
2140: C->ops->productsymbolic = MatProductSymbolic_RARt_SeqAIJ_SeqAIJ;
2141: PetscFunctionReturn(PETSC_SUCCESS);
2142: }
2144: /* ABC = A*B*C = A*(B*C); ABC's algorithm must be chosen from AB's algorithm */
2145: static PetscErrorCode MatProductSetFromOptions_SeqAIJ_ABC(Mat C)
2146: {
2147: Mat_Product *product = C->product;
2148: PetscInt alg = 0; /* default algorithm */
2149: PetscBool flg = PETSC_FALSE;
2150: const char *algTypes[7] = {"sorted", "scalable", "scalable_fast", "heap", "btheap", "llcondensed", "rowmerge"};
2151: PetscInt nalg = 7;
2153: PetscFunctionBegin;
2154: /* Set default algorithm */
2155: PetscCall(PetscStrcmp(product->alg, "default", &flg));
2156: if (flg) PetscCall(MatProductSetAlgorithm(C, (MatProductAlgorithm)algTypes[alg]));
2158: /* Get runtime option */
2159: if (product->api_user) {
2160: PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatMatMatMult", "Mat");
2161: PetscCall(PetscOptionsEList("-matmatmatmult_via", "Algorithmic approach", "MatMatMatMult", algTypes, nalg, algTypes[alg], &alg, &flg));
2162: PetscOptionsEnd();
2163: } else {
2164: PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatProduct_ABC", "Mat");
2165: PetscCall(PetscOptionsEList("-mat_product_algorithm", "Algorithmic approach", "MatProduct_ABC", algTypes, nalg, algTypes[alg], &alg, &flg));
2166: PetscOptionsEnd();
2167: }
2168: if (flg) PetscCall(MatProductSetAlgorithm(C, (MatProductAlgorithm)algTypes[alg]));
2170: C->ops->matmatmultsymbolic = MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ;
2171: C->ops->productsymbolic = MatProductSymbolic_ABC;
2172: PetscFunctionReturn(PETSC_SUCCESS);
2173: }
2175: PetscErrorCode MatProductSetFromOptions_SeqAIJ(Mat C)
2176: {
2177: Mat_Product *product = C->product;
2179: PetscFunctionBegin;
2180: switch (product->type) {
2181: case MATPRODUCT_AB:
2182: PetscCall(MatProductSetFromOptions_SeqAIJ_AB(C));
2183: break;
2184: case MATPRODUCT_AtB:
2185: PetscCall(MatProductSetFromOptions_SeqAIJ_AtB(C));
2186: break;
2187: case MATPRODUCT_ABt:
2188: PetscCall(MatProductSetFromOptions_SeqAIJ_ABt(C));
2189: break;
2190: case MATPRODUCT_PtAP:
2191: PetscCall(MatProductSetFromOptions_SeqAIJ_PtAP(C));
2192: break;
2193: case MATPRODUCT_RARt:
2194: PetscCall(MatProductSetFromOptions_SeqAIJ_RARt(C));
2195: break;
2196: case MATPRODUCT_ABC:
2197: PetscCall(MatProductSetFromOptions_SeqAIJ_ABC(C));
2198: break;
2199: default:
2200: break;
2201: }
2202: PetscFunctionReturn(PETSC_SUCCESS);
2203: }