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