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