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