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