Actual source code: baij2.c
1: #include <../src/mat/impls/baij/seq/baij.h>
2: #include <../src/mat/impls/dense/seq/dense.h>
3: #include <petsc/private/kernels/blockinvert.h>
4: #include <petscbt.h>
5: #include <petscblaslapack.h>
7: #if PetscDefined(HAVE_IMMINTRIN_H) && defined(__AVX2__) && defined(__FMA__) && PetscDefined(USE_REAL_DOUBLE) && !PetscDefined(USE_COMPLEX) && !PetscDefined(USE_64BIT_INDICES)
8: #include <immintrin.h>
9: #elif PetscDefined(HAVE_XMMINTRIN_H)
10: #include <xmmintrin.h>
11: #endif
13: PetscErrorCode MatIncreaseOverlap_SeqBAIJ(Mat A, PetscInt is_max, IS is[], PetscInt ov)
14: {
15: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
16: PetscInt row, i, j, k, l, m, n, *nidx, isz, val, ival;
17: const PetscInt *idx;
18: PetscInt start, end, *ai, *aj, bs;
19: PetscBT table;
21: PetscFunctionBegin;
22: m = a->mbs;
23: ai = a->i;
24: aj = a->j;
25: bs = A->rmap->bs;
27: PetscCheck(ov >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative overlap specified");
29: PetscCall(PetscBTCreate(m, &table));
30: PetscCall(PetscMalloc1(m + 1, &nidx));
32: for (i = 0; i < is_max; i++) {
33: /* Initialise the two local arrays */
34: isz = 0;
35: PetscCall(PetscBTMemzero(m, table));
37: /* Extract the indices, assume there can be duplicate entries */
38: PetscCall(ISGetIndices(is[i], &idx));
39: PetscCall(ISGetLocalSize(is[i], &n));
41: /* Enter these into the temp arrays i.e mark table[row], enter row into new index */
42: for (j = 0; j < n; ++j) {
43: ival = idx[j] / bs; /* convert the indices into block indices */
44: PetscCheck(ival < m, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index greater than mat-dim");
45: if (!PetscBTLookupSet(table, ival)) nidx[isz++] = ival;
46: }
47: PetscCall(ISRestoreIndices(is[i], &idx));
48: PetscCall(ISDestroy(&is[i]));
50: k = 0;
51: for (j = 0; j < ov; j++) { /* for each overlap*/
52: n = isz;
53: for (; k < n; k++) { /* do only those rows in nidx[k], which are not done yet */
54: row = nidx[k];
55: start = ai[row];
56: end = ai[row + 1];
57: for (l = start; l < end; l++) {
58: val = aj[l];
59: if (!PetscBTLookupSet(table, val)) nidx[isz++] = val;
60: }
61: }
62: }
63: PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, isz, nidx, PETSC_COPY_VALUES, is + i));
64: }
65: PetscCall(PetscBTDestroy(&table));
66: PetscCall(PetscFree(nidx));
67: PetscFunctionReturn(PETSC_SUCCESS);
68: }
70: static PetscErrorCode MatCreateSubMatrix_SeqBAIJ_Private(Mat A, IS isrow, IS iscol, MatReuse scall, Mat *B)
71: {
72: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *c;
73: PetscInt *smap, i, k, kstart, kend, oldcols = a->nbs, *lens;
74: PetscInt row, mat_i, *mat_j, tcol, *mat_ilen;
75: const PetscInt *irow, *icol;
76: PetscInt nrows, ncols, *ssmap, bs = A->rmap->bs, bs2 = a->bs2;
77: PetscInt *aj = a->j, *ai = a->i;
78: MatScalar *mat_a;
79: Mat C;
80: PetscBool flag;
82: PetscFunctionBegin;
83: PetscCall(ISGetIndices(isrow, &irow));
84: PetscCall(ISGetIndices(iscol, &icol));
85: PetscCall(ISGetLocalSize(isrow, &nrows));
86: PetscCall(ISGetLocalSize(iscol, &ncols));
88: PetscCall(PetscCalloc1(1 + oldcols, &smap));
89: ssmap = smap;
90: PetscCall(PetscMalloc1(nrows, &lens));
91: for (i = 0; i < ncols; i++) smap[icol[i]] = i + 1;
92: /* determine lens of each row */
93: for (i = 0; i < nrows; i++) {
94: kstart = ai[irow[i]];
95: kend = kstart + a->ilen[irow[i]];
96: lens[i] = 0;
97: for (k = kstart; k < kend; k++) {
98: if (ssmap[aj[k]]) lens[i]++;
99: }
100: }
101: /* Create and fill new matrix */
102: if (scall == MAT_REUSE_MATRIX) {
103: c = (Mat_SeqBAIJ *)(*B)->data;
105: PetscCheck(c->mbs == nrows && c->nbs == ncols && (*B)->rmap->bs == bs, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Submatrix wrong size");
106: PetscCall(PetscArraycmp(c->ilen, lens, c->mbs, &flag));
107: PetscCheck(flag, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Cannot reuse matrix. wrong number of nonzeros");
108: PetscCall(PetscArrayzero(c->ilen, c->mbs));
109: C = *B;
110: } else {
111: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C));
112: PetscCall(MatSetSizes(C, nrows * bs, ncols * bs, PETSC_DETERMINE, PETSC_DETERMINE));
113: PetscCall(MatSetType(C, ((PetscObject)A)->type_name));
114: PetscCall(MatSetOption(C, MAT_STRUCTURE_ONLY, A->structure_only));
115: PetscCall(MatSeqBAIJSetPreallocation(C, bs, 0, lens));
116: }
117: c = (Mat_SeqBAIJ *)C->data;
118: for (i = 0; i < nrows; i++) {
119: row = irow[i];
120: kstart = ai[row];
121: kend = kstart + a->ilen[row];
122: mat_i = c->i[i];
123: mat_j = PetscSafePointerPlusOffset(c->j, mat_i);
124: mat_a = PetscSafePointerPlusOffset(c->a, mat_i * bs2);
125: mat_ilen = c->ilen + i;
126: for (k = kstart; k < kend; k++) {
127: if ((tcol = ssmap[a->j[k]])) {
128: *mat_j++ = tcol - 1;
129: if (!A->structure_only) {
130: PetscCall(PetscArraycpy(mat_a, a->a + k * bs2, bs2));
131: mat_a += bs2;
132: }
133: (*mat_ilen)++;
134: }
135: }
136: }
137: /* sort */
138: if (c->j) {
139: MatScalar *work = NULL;
141: if (!A->structure_only) PetscCall(PetscMalloc1(bs2, &work));
142: for (i = 0; i < nrows; i++) {
143: PetscInt ilen;
144: mat_i = c->i[i];
145: mat_j = c->j + mat_i;
146: mat_a = PetscSafePointerPlusOffset(c->a, mat_i * bs2);
147: ilen = c->ilen[i];
148: if (A->structure_only) PetscCall(PetscSortInt(ilen, mat_j));
149: else PetscCall(PetscSortIntWithDataArray(ilen, mat_j, mat_a, bs2 * sizeof(MatScalar), work));
150: }
151: PetscCall(PetscFree(work));
152: }
154: /* Free work space */
155: PetscCall(ISRestoreIndices(iscol, &icol));
156: PetscCall(PetscFree(smap));
157: PetscCall(PetscFree(lens));
158: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
159: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
161: PetscCall(ISRestoreIndices(isrow, &irow));
162: *B = C;
163: PetscFunctionReturn(PETSC_SUCCESS);
164: }
166: PetscErrorCode MatCreateSubMatrix_SeqBAIJ(Mat A, IS isrow, IS iscol, MatReuse scall, Mat *B)
167: {
168: IS is1, is2;
170: PetscFunctionBegin;
171: PetscCall(ISCompressIndicesGeneral(A->rmap->N, A->rmap->n, A->rmap->bs, 1, &isrow, &is1));
172: if (isrow == iscol) {
173: is2 = is1;
174: PetscCall(PetscObjectReference((PetscObject)is2));
175: } else PetscCall(ISCompressIndicesGeneral(A->cmap->N, A->cmap->n, A->cmap->bs, 1, &iscol, &is2));
176: PetscCall(MatCreateSubMatrix_SeqBAIJ_Private(A, is1, is2, scall, B));
177: PetscCall(ISDestroy(&is1));
178: PetscCall(ISDestroy(&is2));
179: PetscFunctionReturn(PETSC_SUCCESS);
180: }
182: PetscErrorCode MatDestroySubMatrix_SeqBAIJ(Mat C)
183: {
184: Mat_SeqBAIJ *c = (Mat_SeqBAIJ *)C->data;
185: Mat_SubSppt *submatj = c->submatis1;
187: PetscFunctionBegin;
188: PetscCall((*submatj->destroy)(C));
189: PetscCall(MatDestroySubMatrix_Private(submatj));
190: PetscFunctionReturn(PETSC_SUCCESS);
191: }
193: /* Note this has code duplication with MatDestroySubMatrices_SeqAIJ() */
194: PetscErrorCode MatDestroySubMatrices_SeqBAIJ(PetscInt n, Mat *mat[])
195: {
196: PetscInt i;
197: Mat C;
198: Mat_SeqBAIJ *c;
199: Mat_SubSppt *submatj;
201: PetscFunctionBegin;
202: for (i = 0; i < n; i++) {
203: C = (*mat)[i];
204: c = (Mat_SeqBAIJ *)C->data;
205: submatj = c->submatis1;
206: if (submatj) {
207: if (--((PetscObject)C)->refct <= 0) {
208: PetscCall(PetscFree(C->factorprefix));
209: PetscCall((*submatj->destroy)(C));
210: PetscCall(MatDestroySubMatrix_Private(submatj));
211: PetscCall(PetscFree(C->defaultvectype));
212: PetscCall(PetscFree(C->defaultrandtype));
213: PetscCall(PetscFree(C->solvertype));
214: PetscCall(PetscLayoutDestroy(&C->rmap));
215: PetscCall(PetscLayoutDestroy(&C->cmap));
216: PetscCall(PetscHeaderDestroy(&C));
217: }
218: } else {
219: PetscCall(MatDestroy(&C));
220: }
221: }
223: /* Destroy Dummy submatrices created for reuse */
224: PetscCall(MatDestroySubMatrices_Dummy(n, mat));
226: PetscCall(PetscFree(*mat));
227: PetscFunctionReturn(PETSC_SUCCESS);
228: }
230: PetscErrorCode MatCreateSubMatrices_SeqBAIJ(Mat A, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *B[])
231: {
232: PetscInt i;
234: PetscFunctionBegin;
235: if (scall == MAT_INITIAL_MATRIX) PetscCall(PetscCalloc1(n + 1, B));
237: for (i = 0; i < n; i++) PetscCall(MatCreateSubMatrix_SeqBAIJ(A, irow[i], icol[i], scall, &(*B)[i]));
238: PetscFunctionReturn(PETSC_SUCCESS);
239: }
241: /* Should check that shapes of vectors and matrices match */
242: PetscErrorCode MatMult_SeqBAIJ_1(Mat A, Vec xx, Vec zz)
243: {
244: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
245: PetscScalar *z, sum;
246: const PetscScalar *x;
247: const MatScalar *v;
248: PetscInt mbs, i, n;
249: const PetscInt *idx, *ii, *ridx = NULL;
250: PetscBool usecprow = a->compressedrow.use;
252: PetscFunctionBegin;
253: PetscCall(VecGetArrayRead(xx, &x));
254: PetscCall(VecGetArrayWrite(zz, &z));
256: if (usecprow) {
257: mbs = a->compressedrow.nrows;
258: ii = a->compressedrow.i;
259: ridx = a->compressedrow.rindex;
260: PetscCall(PetscArrayzero(z, a->mbs));
261: } else {
262: mbs = a->mbs;
263: ii = a->i;
264: }
266: for (i = 0; i < mbs; i++) {
267: n = ii[1] - ii[0];
268: v = a->a + ii[0];
269: idx = a->j + ii[0];
270: ii++;
271: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
272: PetscPrefetchBlock(v + 1 * n, 1 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
273: sum = 0.0;
274: PetscSparseDensePlusDot(sum, x, v, idx, n);
275: if (usecprow) {
276: z[ridx[i]] = sum;
277: } else {
278: z[i] = sum;
279: }
280: }
281: PetscCall(VecRestoreArrayRead(xx, &x));
282: PetscCall(VecRestoreArrayWrite(zz, &z));
283: PetscCall(PetscLogFlops(2.0 * a->nz - a->nonzerorowcnt));
284: PetscFunctionReturn(PETSC_SUCCESS);
285: }
287: PetscErrorCode MatMult_SeqBAIJ_2(Mat A, Vec xx, Vec zz)
288: {
289: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
290: PetscScalar *z = NULL, sum1, sum2, *zarray;
291: const PetscScalar *x, *xb;
292: PetscScalar x1, x2;
293: const MatScalar *v;
294: PetscInt mbs, i, *idx, *ii, j, n, *ridx = NULL;
295: PetscBool usecprow = a->compressedrow.use;
297: PetscFunctionBegin;
298: PetscCall(VecGetArrayRead(xx, &x));
299: PetscCall(VecGetArrayWrite(zz, &zarray));
301: idx = a->j;
302: v = a->a;
303: if (usecprow) {
304: mbs = a->compressedrow.nrows;
305: ii = a->compressedrow.i;
306: ridx = a->compressedrow.rindex;
307: PetscCall(PetscArrayzero(zarray, 2 * a->mbs));
308: } else {
309: mbs = a->mbs;
310: ii = a->i;
311: z = zarray;
312: }
314: for (i = 0; i < mbs; i++) {
315: n = ii[1] - ii[0];
316: ii++;
317: sum1 = 0.0;
318: sum2 = 0.0;
319: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
320: PetscPrefetchBlock(v + 4 * n, 4 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
321: for (j = 0; j < n; j++) {
322: xb = x + 2 * (*idx++);
323: x1 = xb[0];
324: x2 = xb[1];
325: sum1 += v[0] * x1 + v[2] * x2;
326: sum2 += v[1] * x1 + v[3] * x2;
327: v += 4;
328: }
329: if (usecprow) z = zarray + 2 * ridx[i];
330: z[0] = sum1;
331: z[1] = sum2;
332: if (!usecprow) z += 2;
333: }
334: PetscCall(VecRestoreArrayRead(xx, &x));
335: PetscCall(VecRestoreArrayWrite(zz, &zarray));
336: PetscCall(PetscLogFlops(8.0 * a->nz - 2.0 * a->nonzerorowcnt));
337: PetscFunctionReturn(PETSC_SUCCESS);
338: }
340: PetscErrorCode MatMult_SeqBAIJ_3(Mat A, Vec xx, Vec zz)
341: {
342: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
343: PetscScalar *z = NULL, sum1, sum2, sum3, x1, x2, x3, *zarray;
344: const PetscScalar *x, *xb;
345: const MatScalar *v;
346: PetscInt mbs, i, *idx, *ii, j, n, *ridx = NULL;
347: PetscBool usecprow = a->compressedrow.use;
349: #if PetscDefined(HAVE_PRAGMA_DISJOINT)
350: #pragma disjoint(*v, *z, *xb)
351: #endif
353: PetscFunctionBegin;
354: PetscCall(VecGetArrayRead(xx, &x));
355: PetscCall(VecGetArrayWrite(zz, &zarray));
357: idx = a->j;
358: v = a->a;
359: if (usecprow) {
360: mbs = a->compressedrow.nrows;
361: ii = a->compressedrow.i;
362: ridx = a->compressedrow.rindex;
363: PetscCall(PetscArrayzero(zarray, 3 * a->mbs));
364: } else {
365: mbs = a->mbs;
366: ii = a->i;
367: z = zarray;
368: }
370: for (i = 0; i < mbs; i++) {
371: n = ii[1] - ii[0];
372: ii++;
373: sum1 = 0.0;
374: sum2 = 0.0;
375: sum3 = 0.0;
376: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
377: PetscPrefetchBlock(v + 9 * n, 9 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
378: for (j = 0; j < n; j++) {
379: xb = x + 3 * (*idx++);
380: x1 = xb[0];
381: x2 = xb[1];
382: x3 = xb[2];
384: sum1 += v[0] * x1 + v[3] * x2 + v[6] * x3;
385: sum2 += v[1] * x1 + v[4] * x2 + v[7] * x3;
386: sum3 += v[2] * x1 + v[5] * x2 + v[8] * x3;
387: v += 9;
388: }
389: if (usecprow) z = zarray + 3 * ridx[i];
390: z[0] = sum1;
391: z[1] = sum2;
392: z[2] = sum3;
393: if (!usecprow) z += 3;
394: }
395: PetscCall(VecRestoreArrayRead(xx, &x));
396: PetscCall(VecRestoreArrayWrite(zz, &zarray));
397: PetscCall(PetscLogFlops(18.0 * a->nz - 3.0 * a->nonzerorowcnt));
398: PetscFunctionReturn(PETSC_SUCCESS);
399: }
401: PetscErrorCode MatMult_SeqBAIJ_4(Mat A, Vec xx, Vec zz)
402: {
403: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
404: PetscScalar *z = NULL, sum1, sum2, sum3, sum4, x1, x2, x3, x4, *zarray;
405: const PetscScalar *x, *xb;
406: const MatScalar *v;
407: PetscInt mbs, i, *idx, *ii, j, n, *ridx = NULL;
408: PetscBool usecprow = a->compressedrow.use;
410: PetscFunctionBegin;
411: PetscCall(VecGetArrayRead(xx, &x));
412: PetscCall(VecGetArrayWrite(zz, &zarray));
414: idx = a->j;
415: v = a->a;
416: if (usecprow) {
417: mbs = a->compressedrow.nrows;
418: ii = a->compressedrow.i;
419: ridx = a->compressedrow.rindex;
420: PetscCall(PetscArrayzero(zarray, 4 * a->mbs));
421: } else {
422: mbs = a->mbs;
423: ii = a->i;
424: z = zarray;
425: }
427: for (i = 0; i < mbs; i++) {
428: n = ii[1] - ii[0];
429: ii++;
430: sum1 = 0.0;
431: sum2 = 0.0;
432: sum3 = 0.0;
433: sum4 = 0.0;
435: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
436: PetscPrefetchBlock(v + 16 * n, 16 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
437: for (j = 0; j < n; j++) {
438: xb = x + 4 * (*idx++);
439: x1 = xb[0];
440: x2 = xb[1];
441: x3 = xb[2];
442: x4 = xb[3];
443: sum1 += v[0] * x1 + v[4] * x2 + v[8] * x3 + v[12] * x4;
444: sum2 += v[1] * x1 + v[5] * x2 + v[9] * x3 + v[13] * x4;
445: sum3 += v[2] * x1 + v[6] * x2 + v[10] * x3 + v[14] * x4;
446: sum4 += v[3] * x1 + v[7] * x2 + v[11] * x3 + v[15] * x4;
447: v += 16;
448: }
449: if (usecprow) z = zarray + 4 * ridx[i];
450: z[0] = sum1;
451: z[1] = sum2;
452: z[2] = sum3;
453: z[3] = sum4;
454: if (!usecprow) z += 4;
455: }
456: PetscCall(VecRestoreArrayRead(xx, &x));
457: PetscCall(VecRestoreArrayWrite(zz, &zarray));
458: PetscCall(PetscLogFlops(32.0 * a->nz - 4.0 * a->nonzerorowcnt));
459: PetscFunctionReturn(PETSC_SUCCESS);
460: }
462: PetscErrorCode MatMult_SeqBAIJ_5(Mat A, Vec xx, Vec zz)
463: {
464: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
465: PetscScalar *z = NULL, sum1, sum2, sum3, sum4, sum5, x1, x2, x3, x4, x5, *zarray;
466: const PetscScalar *xb, *x;
467: const MatScalar *v;
468: const PetscInt *idx, *ii, *ridx = NULL;
469: PetscInt mbs, i, j, n;
470: PetscBool usecprow = a->compressedrow.use;
472: PetscFunctionBegin;
473: PetscCall(VecGetArrayRead(xx, &x));
474: PetscCall(VecGetArrayWrite(zz, &zarray));
476: idx = a->j;
477: v = a->a;
478: if (usecprow) {
479: mbs = a->compressedrow.nrows;
480: ii = a->compressedrow.i;
481: ridx = a->compressedrow.rindex;
482: PetscCall(PetscArrayzero(zarray, 5 * a->mbs));
483: } else {
484: mbs = a->mbs;
485: ii = a->i;
486: z = zarray;
487: }
489: for (i = 0; i < mbs; i++) {
490: n = ii[1] - ii[0];
491: ii++;
492: sum1 = 0.0;
493: sum2 = 0.0;
494: sum3 = 0.0;
495: sum4 = 0.0;
496: sum5 = 0.0;
497: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
498: PetscPrefetchBlock(v + 25 * n, 25 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
499: for (j = 0; j < n; j++) {
500: xb = x + 5 * (*idx++);
501: x1 = xb[0];
502: x2 = xb[1];
503: x3 = xb[2];
504: x4 = xb[3];
505: x5 = xb[4];
506: sum1 += v[0] * x1 + v[5] * x2 + v[10] * x3 + v[15] * x4 + v[20] * x5;
507: sum2 += v[1] * x1 + v[6] * x2 + v[11] * x3 + v[16] * x4 + v[21] * x5;
508: sum3 += v[2] * x1 + v[7] * x2 + v[12] * x3 + v[17] * x4 + v[22] * x5;
509: sum4 += v[3] * x1 + v[8] * x2 + v[13] * x3 + v[18] * x4 + v[23] * x5;
510: sum5 += v[4] * x1 + v[9] * x2 + v[14] * x3 + v[19] * x4 + v[24] * x5;
511: v += 25;
512: }
513: if (usecprow) z = zarray + 5 * ridx[i];
514: z[0] = sum1;
515: z[1] = sum2;
516: z[2] = sum3;
517: z[3] = sum4;
518: z[4] = sum5;
519: if (!usecprow) z += 5;
520: }
521: PetscCall(VecRestoreArrayRead(xx, &x));
522: PetscCall(VecRestoreArrayWrite(zz, &zarray));
523: PetscCall(PetscLogFlops(50.0 * a->nz - 5.0 * a->nonzerorowcnt));
524: PetscFunctionReturn(PETSC_SUCCESS);
525: }
527: PetscErrorCode MatMult_SeqBAIJ_6(Mat A, Vec xx, Vec zz)
528: {
529: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
530: PetscScalar *z = NULL, sum1, sum2, sum3, sum4, sum5, sum6;
531: const PetscScalar *x, *xb;
532: PetscScalar x1, x2, x3, x4, x5, x6, *zarray;
533: const MatScalar *v;
534: PetscInt mbs, i, *idx, *ii, j, n, *ridx = NULL;
535: PetscBool usecprow = a->compressedrow.use;
537: PetscFunctionBegin;
538: PetscCall(VecGetArrayRead(xx, &x));
539: PetscCall(VecGetArrayWrite(zz, &zarray));
541: idx = a->j;
542: v = a->a;
543: if (usecprow) {
544: mbs = a->compressedrow.nrows;
545: ii = a->compressedrow.i;
546: ridx = a->compressedrow.rindex;
547: PetscCall(PetscArrayzero(zarray, 6 * a->mbs));
548: } else {
549: mbs = a->mbs;
550: ii = a->i;
551: z = zarray;
552: }
554: for (i = 0; i < mbs; i++) {
555: n = ii[1] - ii[0];
556: ii++;
557: sum1 = 0.0;
558: sum2 = 0.0;
559: sum3 = 0.0;
560: sum4 = 0.0;
561: sum5 = 0.0;
562: sum6 = 0.0;
564: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
565: PetscPrefetchBlock(v + 36 * n, 36 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
566: for (j = 0; j < n; j++) {
567: xb = x + 6 * (*idx++);
568: x1 = xb[0];
569: x2 = xb[1];
570: x3 = xb[2];
571: x4 = xb[3];
572: x5 = xb[4];
573: x6 = xb[5];
574: sum1 += v[0] * x1 + v[6] * x2 + v[12] * x3 + v[18] * x4 + v[24] * x5 + v[30] * x6;
575: sum2 += v[1] * x1 + v[7] * x2 + v[13] * x3 + v[19] * x4 + v[25] * x5 + v[31] * x6;
576: sum3 += v[2] * x1 + v[8] * x2 + v[14] * x3 + v[20] * x4 + v[26] * x5 + v[32] * x6;
577: sum4 += v[3] * x1 + v[9] * x2 + v[15] * x3 + v[21] * x4 + v[27] * x5 + v[33] * x6;
578: sum5 += v[4] * x1 + v[10] * x2 + v[16] * x3 + v[22] * x4 + v[28] * x5 + v[34] * x6;
579: sum6 += v[5] * x1 + v[11] * x2 + v[17] * x3 + v[23] * x4 + v[29] * x5 + v[35] * x6;
580: v += 36;
581: }
582: if (usecprow) z = zarray + 6 * ridx[i];
583: z[0] = sum1;
584: z[1] = sum2;
585: z[2] = sum3;
586: z[3] = sum4;
587: z[4] = sum5;
588: z[5] = sum6;
589: if (!usecprow) z += 6;
590: }
592: PetscCall(VecRestoreArrayRead(xx, &x));
593: PetscCall(VecRestoreArrayWrite(zz, &zarray));
594: PetscCall(PetscLogFlops(72.0 * a->nz - 6.0 * a->nonzerorowcnt));
595: PetscFunctionReturn(PETSC_SUCCESS);
596: }
598: PetscErrorCode MatMult_SeqBAIJ_7(Mat A, Vec xx, Vec zz)
599: {
600: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
601: PetscScalar *z = NULL, sum1, sum2, sum3, sum4, sum5, sum6, sum7;
602: const PetscScalar *x, *xb;
603: PetscScalar x1, x2, x3, x4, x5, x6, x7, *zarray;
604: const MatScalar *v;
605: PetscInt mbs, i, *idx, *ii, j, n, *ridx = NULL;
606: PetscBool usecprow = a->compressedrow.use;
608: PetscFunctionBegin;
609: PetscCall(VecGetArrayRead(xx, &x));
610: PetscCall(VecGetArrayWrite(zz, &zarray));
612: idx = a->j;
613: v = a->a;
614: if (usecprow) {
615: mbs = a->compressedrow.nrows;
616: ii = a->compressedrow.i;
617: ridx = a->compressedrow.rindex;
618: PetscCall(PetscArrayzero(zarray, 7 * a->mbs));
619: } else {
620: mbs = a->mbs;
621: ii = a->i;
622: z = zarray;
623: }
625: for (i = 0; i < mbs; i++) {
626: n = ii[1] - ii[0];
627: ii++;
628: sum1 = 0.0;
629: sum2 = 0.0;
630: sum3 = 0.0;
631: sum4 = 0.0;
632: sum5 = 0.0;
633: sum6 = 0.0;
634: sum7 = 0.0;
636: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
637: PetscPrefetchBlock(v + 49 * n, 49 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
638: for (j = 0; j < n; j++) {
639: xb = x + 7 * (*idx++);
640: x1 = xb[0];
641: x2 = xb[1];
642: x3 = xb[2];
643: x4 = xb[3];
644: x5 = xb[4];
645: x6 = xb[5];
646: x7 = xb[6];
647: sum1 += v[0] * x1 + v[7] * x2 + v[14] * x3 + v[21] * x4 + v[28] * x5 + v[35] * x6 + v[42] * x7;
648: sum2 += v[1] * x1 + v[8] * x2 + v[15] * x3 + v[22] * x4 + v[29] * x5 + v[36] * x6 + v[43] * x7;
649: sum3 += v[2] * x1 + v[9] * x2 + v[16] * x3 + v[23] * x4 + v[30] * x5 + v[37] * x6 + v[44] * x7;
650: sum4 += v[3] * x1 + v[10] * x2 + v[17] * x3 + v[24] * x4 + v[31] * x5 + v[38] * x6 + v[45] * x7;
651: sum5 += v[4] * x1 + v[11] * x2 + v[18] * x3 + v[25] * x4 + v[32] * x5 + v[39] * x6 + v[46] * x7;
652: sum6 += v[5] * x1 + v[12] * x2 + v[19] * x3 + v[26] * x4 + v[33] * x5 + v[40] * x6 + v[47] * x7;
653: sum7 += v[6] * x1 + v[13] * x2 + v[20] * x3 + v[27] * x4 + v[34] * x5 + v[41] * x6 + v[48] * x7;
654: v += 49;
655: }
656: if (usecprow) z = zarray + 7 * ridx[i];
657: z[0] = sum1;
658: z[1] = sum2;
659: z[2] = sum3;
660: z[3] = sum4;
661: z[4] = sum5;
662: z[5] = sum6;
663: z[6] = sum7;
664: if (!usecprow) z += 7;
665: }
667: PetscCall(VecRestoreArrayRead(xx, &x));
668: PetscCall(VecRestoreArrayWrite(zz, &zarray));
669: PetscCall(PetscLogFlops(98.0 * a->nz - 7.0 * a->nonzerorowcnt));
670: PetscFunctionReturn(PETSC_SUCCESS);
671: }
673: #if PetscDefined(HAVE_IMMINTRIN_H) && defined(__AVX2__) && defined(__FMA__) && PetscDefined(USE_REAL_DOUBLE) && !PetscDefined(USE_COMPLEX) && !PetscDefined(USE_64BIT_INDICES)
674: PetscErrorCode MatMult_SeqBAIJ_9_AVX2(Mat A, Vec xx, Vec zz)
675: {
676: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
677: PetscScalar *z = NULL, *work, *workt, *zarray;
678: const PetscScalar *x, *xb;
679: const MatScalar *v;
680: PetscInt mbs, i, bs = A->rmap->bs, j, n, bs2 = a->bs2;
681: const PetscInt *idx, *ii, *ridx = NULL;
682: PetscInt k;
683: PetscBool usecprow = a->compressedrow.use;
685: __m256d a0, a1, a2, a3, a4, a5;
686: __m256d w0, w1, w2, w3;
687: __m256d z0, z1, z2;
688: __m256i mask1 = _mm256_set_epi64x(0LL, 0LL, 0LL, 1LL << 63);
690: PetscFunctionBegin;
691: PetscCall(VecGetArrayRead(xx, &x));
692: PetscCall(VecGetArrayWrite(zz, &zarray));
694: idx = a->j;
695: v = a->a;
696: if (usecprow) {
697: mbs = a->compressedrow.nrows;
698: ii = a->compressedrow.i;
699: ridx = a->compressedrow.rindex;
700: PetscCall(PetscArrayzero(zarray, bs * a->mbs));
701: } else {
702: mbs = a->mbs;
703: ii = a->i;
704: z = zarray;
705: }
707: if (!a->mult_work) {
708: k = PetscMax(A->rmap->n, A->cmap->n);
709: PetscCall(PetscMalloc1(k + 1, &a->mult_work));
710: }
712: work = a->mult_work;
713: for (i = 0; i < mbs; i++) {
714: n = ii[1] - ii[0];
715: ii++;
716: workt = work;
717: for (j = 0; j < n; j++) {
718: xb = x + bs * (*idx++);
719: for (k = 0; k < bs; k++) workt[k] = xb[k];
720: workt += bs;
721: }
722: if (usecprow) z = zarray + bs * ridx[i];
724: z0 = _mm256_setzero_pd();
725: z1 = _mm256_setzero_pd();
726: z2 = _mm256_setzero_pd();
728: for (j = 0; j < n; j++) {
729: /* first column of a */
730: w0 = _mm256_set1_pd(work[j * 9]);
731: a0 = _mm256_loadu_pd(&v[j * 81]);
732: z0 = _mm256_fmadd_pd(a0, w0, z0);
733: a1 = _mm256_loadu_pd(&v[j * 81 + 4]);
734: z1 = _mm256_fmadd_pd(a1, w0, z1);
735: a2 = _mm256_loadu_pd(&v[j * 81 + 8]);
736: z2 = _mm256_fmadd_pd(a2, w0, z2);
738: /* second column of a */
739: w1 = _mm256_set1_pd(work[j * 9 + 1]);
740: a0 = _mm256_loadu_pd(&v[j * 81 + 9]);
741: z0 = _mm256_fmadd_pd(a0, w1, z0);
742: a1 = _mm256_loadu_pd(&v[j * 81 + 13]);
743: z1 = _mm256_fmadd_pd(a1, w1, z1);
744: a2 = _mm256_loadu_pd(&v[j * 81 + 17]);
745: z2 = _mm256_fmadd_pd(a2, w1, z2);
747: /* third column of a */
748: w2 = _mm256_set1_pd(work[j * 9 + 2]);
749: a3 = _mm256_loadu_pd(&v[j * 81 + 18]);
750: z0 = _mm256_fmadd_pd(a3, w2, z0);
751: a4 = _mm256_loadu_pd(&v[j * 81 + 22]);
752: z1 = _mm256_fmadd_pd(a4, w2, z1);
753: a5 = _mm256_loadu_pd(&v[j * 81 + 26]);
754: z2 = _mm256_fmadd_pd(a5, w2, z2);
756: /* fourth column of a */
757: w3 = _mm256_set1_pd(work[j * 9 + 3]);
758: a0 = _mm256_loadu_pd(&v[j * 81 + 27]);
759: z0 = _mm256_fmadd_pd(a0, w3, z0);
760: a1 = _mm256_loadu_pd(&v[j * 81 + 31]);
761: z1 = _mm256_fmadd_pd(a1, w3, z1);
762: a2 = _mm256_loadu_pd(&v[j * 81 + 35]);
763: z2 = _mm256_fmadd_pd(a2, w3, z2);
765: /* fifth column of a */
766: w0 = _mm256_set1_pd(work[j * 9 + 4]);
767: a3 = _mm256_loadu_pd(&v[j * 81 + 36]);
768: z0 = _mm256_fmadd_pd(a3, w0, z0);
769: a4 = _mm256_loadu_pd(&v[j * 81 + 40]);
770: z1 = _mm256_fmadd_pd(a4, w0, z1);
771: a5 = _mm256_loadu_pd(&v[j * 81 + 44]);
772: z2 = _mm256_fmadd_pd(a5, w0, z2);
774: /* sixth column of a */
775: w1 = _mm256_set1_pd(work[j * 9 + 5]);
776: a0 = _mm256_loadu_pd(&v[j * 81 + 45]);
777: z0 = _mm256_fmadd_pd(a0, w1, z0);
778: a1 = _mm256_loadu_pd(&v[j * 81 + 49]);
779: z1 = _mm256_fmadd_pd(a1, w1, z1);
780: a2 = _mm256_loadu_pd(&v[j * 81 + 53]);
781: z2 = _mm256_fmadd_pd(a2, w1, z2);
783: /* seventh column of a */
784: w2 = _mm256_set1_pd(work[j * 9 + 6]);
785: a0 = _mm256_loadu_pd(&v[j * 81 + 54]);
786: z0 = _mm256_fmadd_pd(a0, w2, z0);
787: a1 = _mm256_loadu_pd(&v[j * 81 + 58]);
788: z1 = _mm256_fmadd_pd(a1, w2, z1);
789: a2 = _mm256_loadu_pd(&v[j * 81 + 62]);
790: z2 = _mm256_fmadd_pd(a2, w2, z2);
792: /* eighth column of a */
793: w3 = _mm256_set1_pd(work[j * 9 + 7]);
794: a3 = _mm256_loadu_pd(&v[j * 81 + 63]);
795: z0 = _mm256_fmadd_pd(a3, w3, z0);
796: a4 = _mm256_loadu_pd(&v[j * 81 + 67]);
797: z1 = _mm256_fmadd_pd(a4, w3, z1);
798: a5 = _mm256_loadu_pd(&v[j * 81 + 71]);
799: z2 = _mm256_fmadd_pd(a5, w3, z2);
801: /* ninth column of a */
802: w0 = _mm256_set1_pd(work[j * 9 + 8]);
803: a0 = _mm256_loadu_pd(&v[j * 81 + 72]);
804: z0 = _mm256_fmadd_pd(a0, w0, z0);
805: a1 = _mm256_loadu_pd(&v[j * 81 + 76]);
806: z1 = _mm256_fmadd_pd(a1, w0, z1);
807: a2 = _mm256_maskload_pd(&v[j * 81 + 80], mask1);
808: z2 = _mm256_fmadd_pd(a2, w0, z2);
809: }
811: _mm256_storeu_pd(&z[0], z0);
812: _mm256_storeu_pd(&z[4], z1);
813: _mm256_maskstore_pd(&z[8], mask1, z2);
815: v += n * bs2;
816: if (!usecprow) z += bs;
817: }
818: PetscCall(VecRestoreArrayRead(xx, &x));
819: PetscCall(VecRestoreArrayWrite(zz, &zarray));
820: PetscCall(PetscLogFlops(2.0 * a->nz * bs2 - bs * a->nonzerorowcnt));
821: PetscFunctionReturn(PETSC_SUCCESS);
822: }
823: #endif
825: PetscErrorCode MatMult_SeqBAIJ_11(Mat A, Vec xx, Vec zz)
826: {
827: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
828: PetscScalar *z = NULL, sum1, sum2, sum3, sum4, sum5, sum6, sum7, sum8, sum9, sum10, sum11;
829: const PetscScalar *x, *xb;
830: PetscScalar *zarray, xv;
831: const MatScalar *v;
832: const PetscInt *ii, *ij = a->j, *idx;
833: PetscInt mbs, i, j, k, n, *ridx = NULL;
834: PetscBool usecprow = a->compressedrow.use;
836: PetscFunctionBegin;
837: PetscCall(VecGetArrayRead(xx, &x));
838: PetscCall(VecGetArrayWrite(zz, &zarray));
840: v = a->a;
841: if (usecprow) {
842: mbs = a->compressedrow.nrows;
843: ii = a->compressedrow.i;
844: ridx = a->compressedrow.rindex;
845: PetscCall(PetscArrayzero(zarray, 11 * a->mbs));
846: } else {
847: mbs = a->mbs;
848: ii = a->i;
849: z = zarray;
850: }
852: for (i = 0; i < mbs; i++) {
853: n = ii[i + 1] - ii[i];
854: idx = ij + ii[i];
855: sum1 = 0.0;
856: sum2 = 0.0;
857: sum3 = 0.0;
858: sum4 = 0.0;
859: sum5 = 0.0;
860: sum6 = 0.0;
861: sum7 = 0.0;
862: sum8 = 0.0;
863: sum9 = 0.0;
864: sum10 = 0.0;
865: sum11 = 0.0;
867: for (j = 0; j < n; j++) {
868: xb = x + 11 * idx[j];
870: for (k = 0; k < 11; k++) {
871: xv = xb[k];
872: sum1 += v[0] * xv;
873: sum2 += v[1] * xv;
874: sum3 += v[2] * xv;
875: sum4 += v[3] * xv;
876: sum5 += v[4] * xv;
877: sum6 += v[5] * xv;
878: sum7 += v[6] * xv;
879: sum8 += v[7] * xv;
880: sum9 += v[8] * xv;
881: sum10 += v[9] * xv;
882: sum11 += v[10] * xv;
883: v += 11;
884: }
885: }
886: if (usecprow) z = zarray + 11 * ridx[i];
887: z[0] = sum1;
888: z[1] = sum2;
889: z[2] = sum3;
890: z[3] = sum4;
891: z[4] = sum5;
892: z[5] = sum6;
893: z[6] = sum7;
894: z[7] = sum8;
895: z[8] = sum9;
896: z[9] = sum10;
897: z[10] = sum11;
899: if (!usecprow) z += 11;
900: }
902: PetscCall(VecRestoreArrayRead(xx, &x));
903: PetscCall(VecRestoreArrayWrite(zz, &zarray));
904: PetscCall(PetscLogFlops(242.0 * a->nz - 11.0 * a->nonzerorowcnt));
905: PetscFunctionReturn(PETSC_SUCCESS);
906: }
908: /* MatMult_SeqBAIJ_12 version 1: Columns in the block are accessed one at a time */
909: PetscErrorCode MatMult_SeqBAIJ_12_ver1(Mat A, Vec xx, Vec zz)
910: {
911: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
912: PetscScalar *z = NULL, sum1, sum2, sum3, sum4, sum5, sum6, sum7, sum8, sum9, sum10, sum11, sum12;
913: const PetscScalar *x, *xb;
914: PetscScalar *zarray, xv;
915: const MatScalar *v;
916: const PetscInt *ii, *ij = a->j, *idx;
917: PetscInt mbs, i, j, k, n, *ridx = NULL;
918: PetscBool usecprow = a->compressedrow.use;
920: PetscFunctionBegin;
921: PetscCall(VecGetArrayRead(xx, &x));
922: PetscCall(VecGetArrayWrite(zz, &zarray));
924: v = a->a;
925: if (usecprow) {
926: mbs = a->compressedrow.nrows;
927: ii = a->compressedrow.i;
928: ridx = a->compressedrow.rindex;
929: PetscCall(PetscArrayzero(zarray, 12 * a->mbs));
930: } else {
931: mbs = a->mbs;
932: ii = a->i;
933: z = zarray;
934: }
936: for (i = 0; i < mbs; i++) {
937: n = ii[i + 1] - ii[i];
938: idx = ij + ii[i];
939: sum1 = 0.0;
940: sum2 = 0.0;
941: sum3 = 0.0;
942: sum4 = 0.0;
943: sum5 = 0.0;
944: sum6 = 0.0;
945: sum7 = 0.0;
946: sum8 = 0.0;
947: sum9 = 0.0;
948: sum10 = 0.0;
949: sum11 = 0.0;
950: sum12 = 0.0;
952: for (j = 0; j < n; j++) {
953: xb = x + 12 * idx[j];
955: for (k = 0; k < 12; k++) {
956: xv = xb[k];
957: sum1 += v[0] * xv;
958: sum2 += v[1] * xv;
959: sum3 += v[2] * xv;
960: sum4 += v[3] * xv;
961: sum5 += v[4] * xv;
962: sum6 += v[5] * xv;
963: sum7 += v[6] * xv;
964: sum8 += v[7] * xv;
965: sum9 += v[8] * xv;
966: sum10 += v[9] * xv;
967: sum11 += v[10] * xv;
968: sum12 += v[11] * xv;
969: v += 12;
970: }
971: }
972: if (usecprow) z = zarray + 12 * ridx[i];
973: z[0] = sum1;
974: z[1] = sum2;
975: z[2] = sum3;
976: z[3] = sum4;
977: z[4] = sum5;
978: z[5] = sum6;
979: z[6] = sum7;
980: z[7] = sum8;
981: z[8] = sum9;
982: z[9] = sum10;
983: z[10] = sum11;
984: z[11] = sum12;
985: if (!usecprow) z += 12;
986: }
987: PetscCall(VecRestoreArrayRead(xx, &x));
988: PetscCall(VecRestoreArrayWrite(zz, &zarray));
989: PetscCall(PetscLogFlops(288.0 * a->nz - 12.0 * a->nonzerorowcnt));
990: PetscFunctionReturn(PETSC_SUCCESS);
991: }
993: PetscErrorCode MatMultAdd_SeqBAIJ_12_ver1(Mat A, Vec xx, Vec yy, Vec zz)
994: {
995: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
996: PetscScalar *z = NULL, *y = NULL, sum1, sum2, sum3, sum4, sum5, sum6, sum7, sum8, sum9, sum10, sum11, sum12;
997: const PetscScalar *x, *xb;
998: PetscScalar *zarray, *yarray, xv;
999: const MatScalar *v;
1000: const PetscInt *ii, *ij = a->j, *idx;
1001: PetscInt mbs = a->mbs, i, j, k, n, *ridx = NULL;
1002: PetscBool usecprow = a->compressedrow.use;
1004: PetscFunctionBegin;
1005: PetscCall(VecGetArrayRead(xx, &x));
1006: PetscCall(VecGetArrayPair(yy, zz, &yarray, &zarray));
1008: v = a->a;
1009: if (usecprow) {
1010: if (zz != yy) PetscCall(PetscArraycpy(zarray, yarray, 12 * mbs));
1011: mbs = a->compressedrow.nrows;
1012: ii = a->compressedrow.i;
1013: ridx = a->compressedrow.rindex;
1014: } else {
1015: ii = a->i;
1016: y = yarray;
1017: z = zarray;
1018: }
1020: for (i = 0; i < mbs; i++) {
1021: n = ii[i + 1] - ii[i];
1022: idx = ij + ii[i];
1024: if (usecprow) {
1025: y = yarray + 12 * ridx[i];
1026: z = zarray + 12 * ridx[i];
1027: }
1028: sum1 = y[0];
1029: sum2 = y[1];
1030: sum3 = y[2];
1031: sum4 = y[3];
1032: sum5 = y[4];
1033: sum6 = y[5];
1034: sum7 = y[6];
1035: sum8 = y[7];
1036: sum9 = y[8];
1037: sum10 = y[9];
1038: sum11 = y[10];
1039: sum12 = y[11];
1041: for (j = 0; j < n; j++) {
1042: xb = x + 12 * idx[j];
1044: for (k = 0; k < 12; k++) {
1045: xv = xb[k];
1046: sum1 += v[0] * xv;
1047: sum2 += v[1] * xv;
1048: sum3 += v[2] * xv;
1049: sum4 += v[3] * xv;
1050: sum5 += v[4] * xv;
1051: sum6 += v[5] * xv;
1052: sum7 += v[6] * xv;
1053: sum8 += v[7] * xv;
1054: sum9 += v[8] * xv;
1055: sum10 += v[9] * xv;
1056: sum11 += v[10] * xv;
1057: sum12 += v[11] * xv;
1058: v += 12;
1059: }
1060: }
1062: z[0] = sum1;
1063: z[1] = sum2;
1064: z[2] = sum3;
1065: z[3] = sum4;
1066: z[4] = sum5;
1067: z[5] = sum6;
1068: z[6] = sum7;
1069: z[7] = sum8;
1070: z[8] = sum9;
1071: z[9] = sum10;
1072: z[10] = sum11;
1073: z[11] = sum12;
1074: if (!usecprow) {
1075: y += 12;
1076: z += 12;
1077: }
1078: }
1079: PetscCall(VecRestoreArrayRead(xx, &x));
1080: PetscCall(VecRestoreArrayPair(yy, zz, &yarray, &zarray));
1081: PetscCall(PetscLogFlops(288.0 * a->nz - 12.0 * a->nonzerorowcnt));
1082: PetscFunctionReturn(PETSC_SUCCESS);
1083: }
1085: /* MatMult_SeqBAIJ_12_ver2 : Columns in the block are accessed in sets of 4,4,4 */
1086: PetscErrorCode MatMult_SeqBAIJ_12_ver2(Mat A, Vec xx, Vec zz)
1087: {
1088: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
1089: PetscScalar *z = NULL, sum1, sum2, sum3, sum4, sum5, sum6, sum7, sum8, sum9, sum10, sum11, sum12;
1090: const PetscScalar *x, *xb;
1091: PetscScalar x1, x2, x3, x4, *zarray;
1092: const MatScalar *v;
1093: const PetscInt *ii, *ij = a->j, *idx, *ridx = NULL;
1094: PetscInt mbs, i, j, n;
1095: PetscBool usecprow = a->compressedrow.use;
1097: PetscFunctionBegin;
1098: PetscCall(VecGetArrayRead(xx, &x));
1099: PetscCall(VecGetArrayWrite(zz, &zarray));
1101: v = a->a;
1102: if (usecprow) {
1103: mbs = a->compressedrow.nrows;
1104: ii = a->compressedrow.i;
1105: ridx = a->compressedrow.rindex;
1106: PetscCall(PetscArrayzero(zarray, 12 * a->mbs));
1107: } else {
1108: mbs = a->mbs;
1109: ii = a->i;
1110: z = zarray;
1111: }
1113: for (i = 0; i < mbs; i++) {
1114: n = ii[i + 1] - ii[i];
1115: idx = ij + ii[i];
1117: sum1 = sum2 = sum3 = sum4 = sum5 = sum6 = sum7 = sum8 = sum9 = sum10 = sum11 = sum12 = 0;
1118: for (j = 0; j < n; j++) {
1119: xb = x + 12 * idx[j];
1120: x1 = xb[0];
1121: x2 = xb[1];
1122: x3 = xb[2];
1123: x4 = xb[3];
1125: sum1 += v[0] * x1 + v[12] * x2 + v[24] * x3 + v[36] * x4;
1126: sum2 += v[1] * x1 + v[13] * x2 + v[25] * x3 + v[37] * x4;
1127: sum3 += v[2] * x1 + v[14] * x2 + v[26] * x3 + v[38] * x4;
1128: sum4 += v[3] * x1 + v[15] * x2 + v[27] * x3 + v[39] * x4;
1129: sum5 += v[4] * x1 + v[16] * x2 + v[28] * x3 + v[40] * x4;
1130: sum6 += v[5] * x1 + v[17] * x2 + v[29] * x3 + v[41] * x4;
1131: sum7 += v[6] * x1 + v[18] * x2 + v[30] * x3 + v[42] * x4;
1132: sum8 += v[7] * x1 + v[19] * x2 + v[31] * x3 + v[43] * x4;
1133: sum9 += v[8] * x1 + v[20] * x2 + v[32] * x3 + v[44] * x4;
1134: sum10 += v[9] * x1 + v[21] * x2 + v[33] * x3 + v[45] * x4;
1135: sum11 += v[10] * x1 + v[22] * x2 + v[34] * x3 + v[46] * x4;
1136: sum12 += v[11] * x1 + v[23] * x2 + v[35] * x3 + v[47] * x4;
1137: v += 48;
1139: x1 = xb[4];
1140: x2 = xb[5];
1141: x3 = xb[6];
1142: x4 = xb[7];
1144: sum1 += v[0] * x1 + v[12] * x2 + v[24] * x3 + v[36] * x4;
1145: sum2 += v[1] * x1 + v[13] * x2 + v[25] * x3 + v[37] * x4;
1146: sum3 += v[2] * x1 + v[14] * x2 + v[26] * x3 + v[38] * x4;
1147: sum4 += v[3] * x1 + v[15] * x2 + v[27] * x3 + v[39] * x4;
1148: sum5 += v[4] * x1 + v[16] * x2 + v[28] * x3 + v[40] * x4;
1149: sum6 += v[5] * x1 + v[17] * x2 + v[29] * x3 + v[41] * x4;
1150: sum7 += v[6] * x1 + v[18] * x2 + v[30] * x3 + v[42] * x4;
1151: sum8 += v[7] * x1 + v[19] * x2 + v[31] * x3 + v[43] * x4;
1152: sum9 += v[8] * x1 + v[20] * x2 + v[32] * x3 + v[44] * x4;
1153: sum10 += v[9] * x1 + v[21] * x2 + v[33] * x3 + v[45] * x4;
1154: sum11 += v[10] * x1 + v[22] * x2 + v[34] * x3 + v[46] * x4;
1155: sum12 += v[11] * x1 + v[23] * x2 + v[35] * x3 + v[47] * x4;
1156: v += 48;
1158: x1 = xb[8];
1159: x2 = xb[9];
1160: x3 = xb[10];
1161: x4 = xb[11];
1162: sum1 += v[0] * x1 + v[12] * x2 + v[24] * x3 + v[36] * x4;
1163: sum2 += v[1] * x1 + v[13] * x2 + v[25] * x3 + v[37] * x4;
1164: sum3 += v[2] * x1 + v[14] * x2 + v[26] * x3 + v[38] * x4;
1165: sum4 += v[3] * x1 + v[15] * x2 + v[27] * x3 + v[39] * x4;
1166: sum5 += v[4] * x1 + v[16] * x2 + v[28] * x3 + v[40] * x4;
1167: sum6 += v[5] * x1 + v[17] * x2 + v[29] * x3 + v[41] * x4;
1168: sum7 += v[6] * x1 + v[18] * x2 + v[30] * x3 + v[42] * x4;
1169: sum8 += v[7] * x1 + v[19] * x2 + v[31] * x3 + v[43] * x4;
1170: sum9 += v[8] * x1 + v[20] * x2 + v[32] * x3 + v[44] * x4;
1171: sum10 += v[9] * x1 + v[21] * x2 + v[33] * x3 + v[45] * x4;
1172: sum11 += v[10] * x1 + v[22] * x2 + v[34] * x3 + v[46] * x4;
1173: sum12 += v[11] * x1 + v[23] * x2 + v[35] * x3 + v[47] * x4;
1174: v += 48;
1175: }
1176: if (usecprow) z = zarray + 12 * ridx[i];
1177: z[0] = sum1;
1178: z[1] = sum2;
1179: z[2] = sum3;
1180: z[3] = sum4;
1181: z[4] = sum5;
1182: z[5] = sum6;
1183: z[6] = sum7;
1184: z[7] = sum8;
1185: z[8] = sum9;
1186: z[9] = sum10;
1187: z[10] = sum11;
1188: z[11] = sum12;
1189: if (!usecprow) z += 12;
1190: }
1191: PetscCall(VecRestoreArrayRead(xx, &x));
1192: PetscCall(VecRestoreArrayWrite(zz, &zarray));
1193: PetscCall(PetscLogFlops(288.0 * a->nz - 12.0 * a->nonzerorowcnt));
1194: PetscFunctionReturn(PETSC_SUCCESS);
1195: }
1197: /* MatMultAdd_SeqBAIJ_12_ver2 : Columns in the block are accessed in sets of 4,4,4 */
1198: PetscErrorCode MatMultAdd_SeqBAIJ_12_ver2(Mat A, Vec xx, Vec yy, Vec zz)
1199: {
1200: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
1201: PetscScalar *z = NULL, *y = NULL, sum1, sum2, sum3, sum4, sum5, sum6, sum7, sum8, sum9, sum10, sum11, sum12;
1202: const PetscScalar *x, *xb;
1203: PetscScalar x1, x2, x3, x4, *zarray, *yarray;
1204: const MatScalar *v;
1205: const PetscInt *ii, *ij = a->j, *idx, *ridx = NULL;
1206: PetscInt mbs = a->mbs, i, j, n;
1207: PetscBool usecprow = a->compressedrow.use;
1209: PetscFunctionBegin;
1210: PetscCall(VecGetArrayRead(xx, &x));
1211: PetscCall(VecGetArrayPair(yy, zz, &yarray, &zarray));
1213: v = a->a;
1214: if (usecprow) {
1215: if (zz != yy) PetscCall(PetscArraycpy(zarray, yarray, 12 * mbs));
1216: mbs = a->compressedrow.nrows;
1217: ii = a->compressedrow.i;
1218: ridx = a->compressedrow.rindex;
1219: } else {
1220: ii = a->i;
1221: y = yarray;
1222: z = zarray;
1223: }
1225: for (i = 0; i < mbs; i++) {
1226: n = ii[i + 1] - ii[i];
1227: idx = ij + ii[i];
1229: if (usecprow) {
1230: y = yarray + 12 * ridx[i];
1231: z = zarray + 12 * ridx[i];
1232: }
1233: sum1 = y[0];
1234: sum2 = y[1];
1235: sum3 = y[2];
1236: sum4 = y[3];
1237: sum5 = y[4];
1238: sum6 = y[5];
1239: sum7 = y[6];
1240: sum8 = y[7];
1241: sum9 = y[8];
1242: sum10 = y[9];
1243: sum11 = y[10];
1244: sum12 = y[11];
1246: for (j = 0; j < n; j++) {
1247: xb = x + 12 * idx[j];
1248: x1 = xb[0];
1249: x2 = xb[1];
1250: x3 = xb[2];
1251: x4 = xb[3];
1253: sum1 += v[0] * x1 + v[12] * x2 + v[24] * x3 + v[36] * x4;
1254: sum2 += v[1] * x1 + v[13] * x2 + v[25] * x3 + v[37] * x4;
1255: sum3 += v[2] * x1 + v[14] * x2 + v[26] * x3 + v[38] * x4;
1256: sum4 += v[3] * x1 + v[15] * x2 + v[27] * x3 + v[39] * x4;
1257: sum5 += v[4] * x1 + v[16] * x2 + v[28] * x3 + v[40] * x4;
1258: sum6 += v[5] * x1 + v[17] * x2 + v[29] * x3 + v[41] * x4;
1259: sum7 += v[6] * x1 + v[18] * x2 + v[30] * x3 + v[42] * x4;
1260: sum8 += v[7] * x1 + v[19] * x2 + v[31] * x3 + v[43] * x4;
1261: sum9 += v[8] * x1 + v[20] * x2 + v[32] * x3 + v[44] * x4;
1262: sum10 += v[9] * x1 + v[21] * x2 + v[33] * x3 + v[45] * x4;
1263: sum11 += v[10] * x1 + v[22] * x2 + v[34] * x3 + v[46] * x4;
1264: sum12 += v[11] * x1 + v[23] * x2 + v[35] * x3 + v[47] * x4;
1265: v += 48;
1267: x1 = xb[4];
1268: x2 = xb[5];
1269: x3 = xb[6];
1270: x4 = xb[7];
1272: sum1 += v[0] * x1 + v[12] * x2 + v[24] * x3 + v[36] * x4;
1273: sum2 += v[1] * x1 + v[13] * x2 + v[25] * x3 + v[37] * x4;
1274: sum3 += v[2] * x1 + v[14] * x2 + v[26] * x3 + v[38] * x4;
1275: sum4 += v[3] * x1 + v[15] * x2 + v[27] * x3 + v[39] * x4;
1276: sum5 += v[4] * x1 + v[16] * x2 + v[28] * x3 + v[40] * x4;
1277: sum6 += v[5] * x1 + v[17] * x2 + v[29] * x3 + v[41] * x4;
1278: sum7 += v[6] * x1 + v[18] * x2 + v[30] * x3 + v[42] * x4;
1279: sum8 += v[7] * x1 + v[19] * x2 + v[31] * x3 + v[43] * x4;
1280: sum9 += v[8] * x1 + v[20] * x2 + v[32] * x3 + v[44] * x4;
1281: sum10 += v[9] * x1 + v[21] * x2 + v[33] * x3 + v[45] * x4;
1282: sum11 += v[10] * x1 + v[22] * x2 + v[34] * x3 + v[46] * x4;
1283: sum12 += v[11] * x1 + v[23] * x2 + v[35] * x3 + v[47] * x4;
1284: v += 48;
1286: x1 = xb[8];
1287: x2 = xb[9];
1288: x3 = xb[10];
1289: x4 = xb[11];
1290: sum1 += v[0] * x1 + v[12] * x2 + v[24] * x3 + v[36] * x4;
1291: sum2 += v[1] * x1 + v[13] * x2 + v[25] * x3 + v[37] * x4;
1292: sum3 += v[2] * x1 + v[14] * x2 + v[26] * x3 + v[38] * x4;
1293: sum4 += v[3] * x1 + v[15] * x2 + v[27] * x3 + v[39] * x4;
1294: sum5 += v[4] * x1 + v[16] * x2 + v[28] * x3 + v[40] * x4;
1295: sum6 += v[5] * x1 + v[17] * x2 + v[29] * x3 + v[41] * x4;
1296: sum7 += v[6] * x1 + v[18] * x2 + v[30] * x3 + v[42] * x4;
1297: sum8 += v[7] * x1 + v[19] * x2 + v[31] * x3 + v[43] * x4;
1298: sum9 += v[8] * x1 + v[20] * x2 + v[32] * x3 + v[44] * x4;
1299: sum10 += v[9] * x1 + v[21] * x2 + v[33] * x3 + v[45] * x4;
1300: sum11 += v[10] * x1 + v[22] * x2 + v[34] * x3 + v[46] * x4;
1301: sum12 += v[11] * x1 + v[23] * x2 + v[35] * x3 + v[47] * x4;
1302: v += 48;
1303: }
1304: z[0] = sum1;
1305: z[1] = sum2;
1306: z[2] = sum3;
1307: z[3] = sum4;
1308: z[4] = sum5;
1309: z[5] = sum6;
1310: z[6] = sum7;
1311: z[7] = sum8;
1312: z[8] = sum9;
1313: z[9] = sum10;
1314: z[10] = sum11;
1315: z[11] = sum12;
1316: if (!usecprow) {
1317: y += 12;
1318: z += 12;
1319: }
1320: }
1321: PetscCall(VecRestoreArrayRead(xx, &x));
1322: PetscCall(VecRestoreArrayPair(yy, zz, &yarray, &zarray));
1323: PetscCall(PetscLogFlops(288.0 * a->nz - 12.0 * a->nonzerorowcnt));
1324: PetscFunctionReturn(PETSC_SUCCESS);
1325: }
1327: #if PetscDefined(HAVE_IMMINTRIN_H) && defined(__AVX2__) && defined(__FMA__) && PetscDefined(USE_REAL_DOUBLE) && !PetscDefined(USE_COMPLEX) && !PetscDefined(USE_64BIT_INDICES)
1328: PetscErrorCode MatMult_SeqBAIJ_12_AVX2(Mat A, Vec xx, Vec zz)
1329: {
1330: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
1331: PetscScalar *z = NULL, *zarray;
1332: const PetscScalar *x, *work;
1333: const MatScalar *v = a->a;
1334: PetscInt mbs, i, j, n;
1335: const PetscInt *idx = a->j, *ii, *ridx = NULL;
1336: PetscBool usecprow = a->compressedrow.use;
1337: const PetscInt bs = 12, bs2 = 144;
1339: __m256d a0, a1, a2, a3, a4, a5;
1340: __m256d w0, w1, w2, w3;
1341: __m256d z0, z1, z2;
1343: PetscFunctionBegin;
1344: PetscCall(VecGetArrayRead(xx, &x));
1345: PetscCall(VecGetArrayWrite(zz, &zarray));
1347: if (usecprow) {
1348: mbs = a->compressedrow.nrows;
1349: ii = a->compressedrow.i;
1350: ridx = a->compressedrow.rindex;
1351: PetscCall(PetscArrayzero(zarray, bs * a->mbs));
1352: } else {
1353: mbs = a->mbs;
1354: ii = a->i;
1355: z = zarray;
1356: }
1358: for (i = 0; i < mbs; i++) {
1359: z0 = _mm256_setzero_pd();
1360: z1 = _mm256_setzero_pd();
1361: z2 = _mm256_setzero_pd();
1363: n = ii[1] - ii[0];
1364: ii++;
1365: for (j = 0; j < n; j++) {
1366: work = x + bs * (*idx++);
1368: /* first column of a */
1369: w0 = _mm256_set1_pd(work[0]);
1370: a0 = _mm256_loadu_pd(v + 0);
1371: z0 = _mm256_fmadd_pd(a0, w0, z0);
1372: a1 = _mm256_loadu_pd(v + 4);
1373: z1 = _mm256_fmadd_pd(a1, w0, z1);
1374: a2 = _mm256_loadu_pd(v + 8);
1375: z2 = _mm256_fmadd_pd(a2, w0, z2);
1377: /* second column of a */
1378: w1 = _mm256_set1_pd(work[1]);
1379: a3 = _mm256_loadu_pd(v + 12);
1380: z0 = _mm256_fmadd_pd(a3, w1, z0);
1381: a4 = _mm256_loadu_pd(v + 16);
1382: z1 = _mm256_fmadd_pd(a4, w1, z1);
1383: a5 = _mm256_loadu_pd(v + 20);
1384: z2 = _mm256_fmadd_pd(a5, w1, z2);
1386: /* third column of a */
1387: w2 = _mm256_set1_pd(work[2]);
1388: a0 = _mm256_loadu_pd(v + 24);
1389: z0 = _mm256_fmadd_pd(a0, w2, z0);
1390: a1 = _mm256_loadu_pd(v + 28);
1391: z1 = _mm256_fmadd_pd(a1, w2, z1);
1392: a2 = _mm256_loadu_pd(v + 32);
1393: z2 = _mm256_fmadd_pd(a2, w2, z2);
1395: /* fourth column of a */
1396: w3 = _mm256_set1_pd(work[3]);
1397: a3 = _mm256_loadu_pd(v + 36);
1398: z0 = _mm256_fmadd_pd(a3, w3, z0);
1399: a4 = _mm256_loadu_pd(v + 40);
1400: z1 = _mm256_fmadd_pd(a4, w3, z1);
1401: a5 = _mm256_loadu_pd(v + 44);
1402: z2 = _mm256_fmadd_pd(a5, w3, z2);
1404: /* fifth column of a */
1405: w0 = _mm256_set1_pd(work[4]);
1406: a0 = _mm256_loadu_pd(v + 48);
1407: z0 = _mm256_fmadd_pd(a0, w0, z0);
1408: a1 = _mm256_loadu_pd(v + 52);
1409: z1 = _mm256_fmadd_pd(a1, w0, z1);
1410: a2 = _mm256_loadu_pd(v + 56);
1411: z2 = _mm256_fmadd_pd(a2, w0, z2);
1413: /* sixth column of a */
1414: w1 = _mm256_set1_pd(work[5]);
1415: a3 = _mm256_loadu_pd(v + 60);
1416: z0 = _mm256_fmadd_pd(a3, w1, z0);
1417: a4 = _mm256_loadu_pd(v + 64);
1418: z1 = _mm256_fmadd_pd(a4, w1, z1);
1419: a5 = _mm256_loadu_pd(v + 68);
1420: z2 = _mm256_fmadd_pd(a5, w1, z2);
1422: /* seventh column of a */
1423: w2 = _mm256_set1_pd(work[6]);
1424: a0 = _mm256_loadu_pd(v + 72);
1425: z0 = _mm256_fmadd_pd(a0, w2, z0);
1426: a1 = _mm256_loadu_pd(v + 76);
1427: z1 = _mm256_fmadd_pd(a1, w2, z1);
1428: a2 = _mm256_loadu_pd(v + 80);
1429: z2 = _mm256_fmadd_pd(a2, w2, z2);
1431: /* eighth column of a */
1432: w3 = _mm256_set1_pd(work[7]);
1433: a3 = _mm256_loadu_pd(v + 84);
1434: z0 = _mm256_fmadd_pd(a3, w3, z0);
1435: a4 = _mm256_loadu_pd(v + 88);
1436: z1 = _mm256_fmadd_pd(a4, w3, z1);
1437: a5 = _mm256_loadu_pd(v + 92);
1438: z2 = _mm256_fmadd_pd(a5, w3, z2);
1440: /* ninth column of a */
1441: w0 = _mm256_set1_pd(work[8]);
1442: a0 = _mm256_loadu_pd(v + 96);
1443: z0 = _mm256_fmadd_pd(a0, w0, z0);
1444: a1 = _mm256_loadu_pd(v + 100);
1445: z1 = _mm256_fmadd_pd(a1, w0, z1);
1446: a2 = _mm256_loadu_pd(v + 104);
1447: z2 = _mm256_fmadd_pd(a2, w0, z2);
1449: /* tenth column of a */
1450: w1 = _mm256_set1_pd(work[9]);
1451: a3 = _mm256_loadu_pd(v + 108);
1452: z0 = _mm256_fmadd_pd(a3, w1, z0);
1453: a4 = _mm256_loadu_pd(v + 112);
1454: z1 = _mm256_fmadd_pd(a4, w1, z1);
1455: a5 = _mm256_loadu_pd(v + 116);
1456: z2 = _mm256_fmadd_pd(a5, w1, z2);
1458: /* eleventh column of a */
1459: w2 = _mm256_set1_pd(work[10]);
1460: a0 = _mm256_loadu_pd(v + 120);
1461: z0 = _mm256_fmadd_pd(a0, w2, z0);
1462: a1 = _mm256_loadu_pd(v + 124);
1463: z1 = _mm256_fmadd_pd(a1, w2, z1);
1464: a2 = _mm256_loadu_pd(v + 128);
1465: z2 = _mm256_fmadd_pd(a2, w2, z2);
1467: /* twelveth column of a */
1468: w3 = _mm256_set1_pd(work[11]);
1469: a3 = _mm256_loadu_pd(v + 132);
1470: z0 = _mm256_fmadd_pd(a3, w3, z0);
1471: a4 = _mm256_loadu_pd(v + 136);
1472: z1 = _mm256_fmadd_pd(a4, w3, z1);
1473: a5 = _mm256_loadu_pd(v + 140);
1474: z2 = _mm256_fmadd_pd(a5, w3, z2);
1476: v += bs2;
1477: }
1478: if (usecprow) z = zarray + bs * ridx[i];
1479: _mm256_storeu_pd(&z[0], z0);
1480: _mm256_storeu_pd(&z[4], z1);
1481: _mm256_storeu_pd(&z[8], z2);
1482: if (!usecprow) z += bs;
1483: }
1484: PetscCall(VecRestoreArrayRead(xx, &x));
1485: PetscCall(VecRestoreArrayWrite(zz, &zarray));
1486: PetscCall(PetscLogFlops(2.0 * a->nz * bs2 - bs * a->nonzerorowcnt));
1487: PetscFunctionReturn(PETSC_SUCCESS);
1488: }
1489: #endif
1491: /* MatMult_SeqBAIJ_15 version 1: Columns in the block are accessed one at a time */
1492: /* Default MatMult for block size 15 */
1493: PetscErrorCode MatMult_SeqBAIJ_15_ver1(Mat A, Vec xx, Vec zz)
1494: {
1495: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
1496: PetscScalar *z = NULL, sum1, sum2, sum3, sum4, sum5, sum6, sum7, sum8, sum9, sum10, sum11, sum12, sum13, sum14, sum15;
1497: const PetscScalar *x, *xb;
1498: PetscScalar *zarray, xv;
1499: const MatScalar *v;
1500: const PetscInt *ii, *ij = a->j, *idx;
1501: PetscInt mbs, i, j, k, n, *ridx = NULL;
1502: PetscBool usecprow = a->compressedrow.use;
1504: PetscFunctionBegin;
1505: PetscCall(VecGetArrayRead(xx, &x));
1506: PetscCall(VecGetArrayWrite(zz, &zarray));
1508: v = a->a;
1509: if (usecprow) {
1510: mbs = a->compressedrow.nrows;
1511: ii = a->compressedrow.i;
1512: ridx = a->compressedrow.rindex;
1513: PetscCall(PetscArrayzero(zarray, 15 * a->mbs));
1514: } else {
1515: mbs = a->mbs;
1516: ii = a->i;
1517: z = zarray;
1518: }
1520: for (i = 0; i < mbs; i++) {
1521: n = ii[i + 1] - ii[i];
1522: idx = ij + ii[i];
1523: sum1 = 0.0;
1524: sum2 = 0.0;
1525: sum3 = 0.0;
1526: sum4 = 0.0;
1527: sum5 = 0.0;
1528: sum6 = 0.0;
1529: sum7 = 0.0;
1530: sum8 = 0.0;
1531: sum9 = 0.0;
1532: sum10 = 0.0;
1533: sum11 = 0.0;
1534: sum12 = 0.0;
1535: sum13 = 0.0;
1536: sum14 = 0.0;
1537: sum15 = 0.0;
1539: for (j = 0; j < n; j++) {
1540: xb = x + 15 * idx[j];
1542: for (k = 0; k < 15; k++) {
1543: xv = xb[k];
1544: sum1 += v[0] * xv;
1545: sum2 += v[1] * xv;
1546: sum3 += v[2] * xv;
1547: sum4 += v[3] * xv;
1548: sum5 += v[4] * xv;
1549: sum6 += v[5] * xv;
1550: sum7 += v[6] * xv;
1551: sum8 += v[7] * xv;
1552: sum9 += v[8] * xv;
1553: sum10 += v[9] * xv;
1554: sum11 += v[10] * xv;
1555: sum12 += v[11] * xv;
1556: sum13 += v[12] * xv;
1557: sum14 += v[13] * xv;
1558: sum15 += v[14] * xv;
1559: v += 15;
1560: }
1561: }
1562: if (usecprow) z = zarray + 15 * ridx[i];
1563: z[0] = sum1;
1564: z[1] = sum2;
1565: z[2] = sum3;
1566: z[3] = sum4;
1567: z[4] = sum5;
1568: z[5] = sum6;
1569: z[6] = sum7;
1570: z[7] = sum8;
1571: z[8] = sum9;
1572: z[9] = sum10;
1573: z[10] = sum11;
1574: z[11] = sum12;
1575: z[12] = sum13;
1576: z[13] = sum14;
1577: z[14] = sum15;
1579: if (!usecprow) z += 15;
1580: }
1582: PetscCall(VecRestoreArrayRead(xx, &x));
1583: PetscCall(VecRestoreArrayWrite(zz, &zarray));
1584: PetscCall(PetscLogFlops(450.0 * a->nz - 15.0 * a->nonzerorowcnt));
1585: PetscFunctionReturn(PETSC_SUCCESS);
1586: }
1588: /* MatMult_SeqBAIJ_15_ver2 : Columns in the block are accessed in sets of 4,4,4,3 */
1589: PetscErrorCode MatMult_SeqBAIJ_15_ver2(Mat A, Vec xx, Vec zz)
1590: {
1591: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
1592: PetscScalar *z = NULL, sum1, sum2, sum3, sum4, sum5, sum6, sum7, sum8, sum9, sum10, sum11, sum12, sum13, sum14, sum15;
1593: const PetscScalar *x, *xb;
1594: PetscScalar x1, x2, x3, x4, *zarray;
1595: const MatScalar *v;
1596: const PetscInt *ii, *ij = a->j, *idx;
1597: PetscInt mbs, i, j, n, *ridx = NULL;
1598: PetscBool usecprow = a->compressedrow.use;
1600: PetscFunctionBegin;
1601: PetscCall(VecGetArrayRead(xx, &x));
1602: PetscCall(VecGetArrayWrite(zz, &zarray));
1604: v = a->a;
1605: if (usecprow) {
1606: mbs = a->compressedrow.nrows;
1607: ii = a->compressedrow.i;
1608: ridx = a->compressedrow.rindex;
1609: PetscCall(PetscArrayzero(zarray, 15 * a->mbs));
1610: } else {
1611: mbs = a->mbs;
1612: ii = a->i;
1613: z = zarray;
1614: }
1616: for (i = 0; i < mbs; i++) {
1617: n = ii[i + 1] - ii[i];
1618: idx = ij + ii[i];
1619: sum1 = 0.0;
1620: sum2 = 0.0;
1621: sum3 = 0.0;
1622: sum4 = 0.0;
1623: sum5 = 0.0;
1624: sum6 = 0.0;
1625: sum7 = 0.0;
1626: sum8 = 0.0;
1627: sum9 = 0.0;
1628: sum10 = 0.0;
1629: sum11 = 0.0;
1630: sum12 = 0.0;
1631: sum13 = 0.0;
1632: sum14 = 0.0;
1633: sum15 = 0.0;
1635: for (j = 0; j < n; j++) {
1636: xb = x + 15 * idx[j];
1637: x1 = xb[0];
1638: x2 = xb[1];
1639: x3 = xb[2];
1640: x4 = xb[3];
1642: sum1 += v[0] * x1 + v[15] * x2 + v[30] * x3 + v[45] * x4;
1643: sum2 += v[1] * x1 + v[16] * x2 + v[31] * x3 + v[46] * x4;
1644: sum3 += v[2] * x1 + v[17] * x2 + v[32] * x3 + v[47] * x4;
1645: sum4 += v[3] * x1 + v[18] * x2 + v[33] * x3 + v[48] * x4;
1646: sum5 += v[4] * x1 + v[19] * x2 + v[34] * x3 + v[49] * x4;
1647: sum6 += v[5] * x1 + v[20] * x2 + v[35] * x3 + v[50] * x4;
1648: sum7 += v[6] * x1 + v[21] * x2 + v[36] * x3 + v[51] * x4;
1649: sum8 += v[7] * x1 + v[22] * x2 + v[37] * x3 + v[52] * x4;
1650: sum9 += v[8] * x1 + v[23] * x2 + v[38] * x3 + v[53] * x4;
1651: sum10 += v[9] * x1 + v[24] * x2 + v[39] * x3 + v[54] * x4;
1652: sum11 += v[10] * x1 + v[25] * x2 + v[40] * x3 + v[55] * x4;
1653: sum12 += v[11] * x1 + v[26] * x2 + v[41] * x3 + v[56] * x4;
1654: sum13 += v[12] * x1 + v[27] * x2 + v[42] * x3 + v[57] * x4;
1655: sum14 += v[13] * x1 + v[28] * x2 + v[43] * x3 + v[58] * x4;
1656: sum15 += v[14] * x1 + v[29] * x2 + v[44] * x3 + v[59] * x4;
1658: v += 60;
1660: x1 = xb[4];
1661: x2 = xb[5];
1662: x3 = xb[6];
1663: x4 = xb[7];
1665: sum1 += v[0] * x1 + v[15] * x2 + v[30] * x3 + v[45] * x4;
1666: sum2 += v[1] * x1 + v[16] * x2 + v[31] * x3 + v[46] * x4;
1667: sum3 += v[2] * x1 + v[17] * x2 + v[32] * x3 + v[47] * x4;
1668: sum4 += v[3] * x1 + v[18] * x2 + v[33] * x3 + v[48] * x4;
1669: sum5 += v[4] * x1 + v[19] * x2 + v[34] * x3 + v[49] * x4;
1670: sum6 += v[5] * x1 + v[20] * x2 + v[35] * x3 + v[50] * x4;
1671: sum7 += v[6] * x1 + v[21] * x2 + v[36] * x3 + v[51] * x4;
1672: sum8 += v[7] * x1 + v[22] * x2 + v[37] * x3 + v[52] * x4;
1673: sum9 += v[8] * x1 + v[23] * x2 + v[38] * x3 + v[53] * x4;
1674: sum10 += v[9] * x1 + v[24] * x2 + v[39] * x3 + v[54] * x4;
1675: sum11 += v[10] * x1 + v[25] * x2 + v[40] * x3 + v[55] * x4;
1676: sum12 += v[11] * x1 + v[26] * x2 + v[41] * x3 + v[56] * x4;
1677: sum13 += v[12] * x1 + v[27] * x2 + v[42] * x3 + v[57] * x4;
1678: sum14 += v[13] * x1 + v[28] * x2 + v[43] * x3 + v[58] * x4;
1679: sum15 += v[14] * x1 + v[29] * x2 + v[44] * x3 + v[59] * x4;
1680: v += 60;
1682: x1 = xb[8];
1683: x2 = xb[9];
1684: x3 = xb[10];
1685: x4 = xb[11];
1686: sum1 += v[0] * x1 + v[15] * x2 + v[30] * x3 + v[45] * x4;
1687: sum2 += v[1] * x1 + v[16] * x2 + v[31] * x3 + v[46] * x4;
1688: sum3 += v[2] * x1 + v[17] * x2 + v[32] * x3 + v[47] * x4;
1689: sum4 += v[3] * x1 + v[18] * x2 + v[33] * x3 + v[48] * x4;
1690: sum5 += v[4] * x1 + v[19] * x2 + v[34] * x3 + v[49] * x4;
1691: sum6 += v[5] * x1 + v[20] * x2 + v[35] * x3 + v[50] * x4;
1692: sum7 += v[6] * x1 + v[21] * x2 + v[36] * x3 + v[51] * x4;
1693: sum8 += v[7] * x1 + v[22] * x2 + v[37] * x3 + v[52] * x4;
1694: sum9 += v[8] * x1 + v[23] * x2 + v[38] * x3 + v[53] * x4;
1695: sum10 += v[9] * x1 + v[24] * x2 + v[39] * x3 + v[54] * x4;
1696: sum11 += v[10] * x1 + v[25] * x2 + v[40] * x3 + v[55] * x4;
1697: sum12 += v[11] * x1 + v[26] * x2 + v[41] * x3 + v[56] * x4;
1698: sum13 += v[12] * x1 + v[27] * x2 + v[42] * x3 + v[57] * x4;
1699: sum14 += v[13] * x1 + v[28] * x2 + v[43] * x3 + v[58] * x4;
1700: sum15 += v[14] * x1 + v[29] * x2 + v[44] * x3 + v[59] * x4;
1701: v += 60;
1703: x1 = xb[12];
1704: x2 = xb[13];
1705: x3 = xb[14];
1706: sum1 += v[0] * x1 + v[15] * x2 + v[30] * x3;
1707: sum2 += v[1] * x1 + v[16] * x2 + v[31] * x3;
1708: sum3 += v[2] * x1 + v[17] * x2 + v[32] * x3;
1709: sum4 += v[3] * x1 + v[18] * x2 + v[33] * x3;
1710: sum5 += v[4] * x1 + v[19] * x2 + v[34] * x3;
1711: sum6 += v[5] * x1 + v[20] * x2 + v[35] * x3;
1712: sum7 += v[6] * x1 + v[21] * x2 + v[36] * x3;
1713: sum8 += v[7] * x1 + v[22] * x2 + v[37] * x3;
1714: sum9 += v[8] * x1 + v[23] * x2 + v[38] * x3;
1715: sum10 += v[9] * x1 + v[24] * x2 + v[39] * x3;
1716: sum11 += v[10] * x1 + v[25] * x2 + v[40] * x3;
1717: sum12 += v[11] * x1 + v[26] * x2 + v[41] * x3;
1718: sum13 += v[12] * x1 + v[27] * x2 + v[42] * x3;
1719: sum14 += v[13] * x1 + v[28] * x2 + v[43] * x3;
1720: sum15 += v[14] * x1 + v[29] * x2 + v[44] * x3;
1721: v += 45;
1722: }
1723: if (usecprow) z = zarray + 15 * ridx[i];
1724: z[0] = sum1;
1725: z[1] = sum2;
1726: z[2] = sum3;
1727: z[3] = sum4;
1728: z[4] = sum5;
1729: z[5] = sum6;
1730: z[6] = sum7;
1731: z[7] = sum8;
1732: z[8] = sum9;
1733: z[9] = sum10;
1734: z[10] = sum11;
1735: z[11] = sum12;
1736: z[12] = sum13;
1737: z[13] = sum14;
1738: z[14] = sum15;
1740: if (!usecprow) z += 15;
1741: }
1743: PetscCall(VecRestoreArrayRead(xx, &x));
1744: PetscCall(VecRestoreArrayWrite(zz, &zarray));
1745: PetscCall(PetscLogFlops(450.0 * a->nz - 15.0 * a->nonzerorowcnt));
1746: PetscFunctionReturn(PETSC_SUCCESS);
1747: }
1749: /* MatMult_SeqBAIJ_15_ver3 : Columns in the block are accessed in sets of 8,7 */
1750: PetscErrorCode MatMult_SeqBAIJ_15_ver3(Mat A, Vec xx, Vec zz)
1751: {
1752: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
1753: PetscScalar *z = NULL, sum1, sum2, sum3, sum4, sum5, sum6, sum7, sum8, sum9, sum10, sum11, sum12, sum13, sum14, sum15;
1754: const PetscScalar *x, *xb;
1755: PetscScalar x1, x2, x3, x4, x5, x6, x7, x8, *zarray;
1756: const MatScalar *v;
1757: const PetscInt *ii, *ij = a->j, *idx;
1758: PetscInt mbs, i, j, n, *ridx = NULL;
1759: PetscBool usecprow = a->compressedrow.use;
1761: PetscFunctionBegin;
1762: PetscCall(VecGetArrayRead(xx, &x));
1763: PetscCall(VecGetArrayWrite(zz, &zarray));
1765: v = a->a;
1766: if (usecprow) {
1767: mbs = a->compressedrow.nrows;
1768: ii = a->compressedrow.i;
1769: ridx = a->compressedrow.rindex;
1770: PetscCall(PetscArrayzero(zarray, 15 * a->mbs));
1771: } else {
1772: mbs = a->mbs;
1773: ii = a->i;
1774: z = zarray;
1775: }
1777: for (i = 0; i < mbs; i++) {
1778: n = ii[i + 1] - ii[i];
1779: idx = ij + ii[i];
1780: sum1 = 0.0;
1781: sum2 = 0.0;
1782: sum3 = 0.0;
1783: sum4 = 0.0;
1784: sum5 = 0.0;
1785: sum6 = 0.0;
1786: sum7 = 0.0;
1787: sum8 = 0.0;
1788: sum9 = 0.0;
1789: sum10 = 0.0;
1790: sum11 = 0.0;
1791: sum12 = 0.0;
1792: sum13 = 0.0;
1793: sum14 = 0.0;
1794: sum15 = 0.0;
1796: for (j = 0; j < n; j++) {
1797: xb = x + 15 * idx[j];
1798: x1 = xb[0];
1799: x2 = xb[1];
1800: x3 = xb[2];
1801: x4 = xb[3];
1802: x5 = xb[4];
1803: x6 = xb[5];
1804: x7 = xb[6];
1805: x8 = xb[7];
1807: sum1 += v[0] * x1 + v[15] * x2 + v[30] * x3 + v[45] * x4 + v[60] * x5 + v[75] * x6 + v[90] * x7 + v[105] * x8;
1808: sum2 += v[1] * x1 + v[16] * x2 + v[31] * x3 + v[46] * x4 + v[61] * x5 + v[76] * x6 + v[91] * x7 + v[106] * x8;
1809: sum3 += v[2] * x1 + v[17] * x2 + v[32] * x3 + v[47] * x4 + v[62] * x5 + v[77] * x6 + v[92] * x7 + v[107] * x8;
1810: sum4 += v[3] * x1 + v[18] * x2 + v[33] * x3 + v[48] * x4 + v[63] * x5 + v[78] * x6 + v[93] * x7 + v[108] * x8;
1811: sum5 += v[4] * x1 + v[19] * x2 + v[34] * x3 + v[49] * x4 + v[64] * x5 + v[79] * x6 + v[94] * x7 + v[109] * x8;
1812: sum6 += v[5] * x1 + v[20] * x2 + v[35] * x3 + v[50] * x4 + v[65] * x5 + v[80] * x6 + v[95] * x7 + v[110] * x8;
1813: sum7 += v[6] * x1 + v[21] * x2 + v[36] * x3 + v[51] * x4 + v[66] * x5 + v[81] * x6 + v[96] * x7 + v[111] * x8;
1814: sum8 += v[7] * x1 + v[22] * x2 + v[37] * x3 + v[52] * x4 + v[67] * x5 + v[82] * x6 + v[97] * x7 + v[112] * x8;
1815: sum9 += v[8] * x1 + v[23] * x2 + v[38] * x3 + v[53] * x4 + v[68] * x5 + v[83] * x6 + v[98] * x7 + v[113] * x8;
1816: sum10 += v[9] * x1 + v[24] * x2 + v[39] * x3 + v[54] * x4 + v[69] * x5 + v[84] * x6 + v[99] * x7 + v[114] * x8;
1817: sum11 += v[10] * x1 + v[25] * x2 + v[40] * x3 + v[55] * x4 + v[70] * x5 + v[85] * x6 + v[100] * x7 + v[115] * x8;
1818: sum12 += v[11] * x1 + v[26] * x2 + v[41] * x3 + v[56] * x4 + v[71] * x5 + v[86] * x6 + v[101] * x7 + v[116] * x8;
1819: sum13 += v[12] * x1 + v[27] * x2 + v[42] * x3 + v[57] * x4 + v[72] * x5 + v[87] * x6 + v[102] * x7 + v[117] * x8;
1820: sum14 += v[13] * x1 + v[28] * x2 + v[43] * x3 + v[58] * x4 + v[73] * x5 + v[88] * x6 + v[103] * x7 + v[118] * x8;
1821: sum15 += v[14] * x1 + v[29] * x2 + v[44] * x3 + v[59] * x4 + v[74] * x5 + v[89] * x6 + v[104] * x7 + v[119] * x8;
1822: v += 120;
1824: x1 = xb[8];
1825: x2 = xb[9];
1826: x3 = xb[10];
1827: x4 = xb[11];
1828: x5 = xb[12];
1829: x6 = xb[13];
1830: x7 = xb[14];
1832: sum1 += v[0] * x1 + v[15] * x2 + v[30] * x3 + v[45] * x4 + v[60] * x5 + v[75] * x6 + v[90] * x7;
1833: sum2 += v[1] * x1 + v[16] * x2 + v[31] * x3 + v[46] * x4 + v[61] * x5 + v[76] * x6 + v[91] * x7;
1834: sum3 += v[2] * x1 + v[17] * x2 + v[32] * x3 + v[47] * x4 + v[62] * x5 + v[77] * x6 + v[92] * x7;
1835: sum4 += v[3] * x1 + v[18] * x2 + v[33] * x3 + v[48] * x4 + v[63] * x5 + v[78] * x6 + v[93] * x7;
1836: sum5 += v[4] * x1 + v[19] * x2 + v[34] * x3 + v[49] * x4 + v[64] * x5 + v[79] * x6 + v[94] * x7;
1837: sum6 += v[5] * x1 + v[20] * x2 + v[35] * x3 + v[50] * x4 + v[65] * x5 + v[80] * x6 + v[95] * x7;
1838: sum7 += v[6] * x1 + v[21] * x2 + v[36] * x3 + v[51] * x4 + v[66] * x5 + v[81] * x6 + v[96] * x7;
1839: sum8 += v[7] * x1 + v[22] * x2 + v[37] * x3 + v[52] * x4 + v[67] * x5 + v[82] * x6 + v[97] * x7;
1840: sum9 += v[8] * x1 + v[23] * x2 + v[38] * x3 + v[53] * x4 + v[68] * x5 + v[83] * x6 + v[98] * x7;
1841: sum10 += v[9] * x1 + v[24] * x2 + v[39] * x3 + v[54] * x4 + v[69] * x5 + v[84] * x6 + v[99] * x7;
1842: sum11 += v[10] * x1 + v[25] * x2 + v[40] * x3 + v[55] * x4 + v[70] * x5 + v[85] * x6 + v[100] * x7;
1843: sum12 += v[11] * x1 + v[26] * x2 + v[41] * x3 + v[56] * x4 + v[71] * x5 + v[86] * x6 + v[101] * x7;
1844: sum13 += v[12] * x1 + v[27] * x2 + v[42] * x3 + v[57] * x4 + v[72] * x5 + v[87] * x6 + v[102] * x7;
1845: sum14 += v[13] * x1 + v[28] * x2 + v[43] * x3 + v[58] * x4 + v[73] * x5 + v[88] * x6 + v[103] * x7;
1846: sum15 += v[14] * x1 + v[29] * x2 + v[44] * x3 + v[59] * x4 + v[74] * x5 + v[89] * x6 + v[104] * x7;
1847: v += 105;
1848: }
1849: if (usecprow) z = zarray + 15 * ridx[i];
1850: z[0] = sum1;
1851: z[1] = sum2;
1852: z[2] = sum3;
1853: z[3] = sum4;
1854: z[4] = sum5;
1855: z[5] = sum6;
1856: z[6] = sum7;
1857: z[7] = sum8;
1858: z[8] = sum9;
1859: z[9] = sum10;
1860: z[10] = sum11;
1861: z[11] = sum12;
1862: z[12] = sum13;
1863: z[13] = sum14;
1864: z[14] = sum15;
1866: if (!usecprow) z += 15;
1867: }
1869: PetscCall(VecRestoreArrayRead(xx, &x));
1870: PetscCall(VecRestoreArrayWrite(zz, &zarray));
1871: PetscCall(PetscLogFlops(450.0 * a->nz - 15.0 * a->nonzerorowcnt));
1872: PetscFunctionReturn(PETSC_SUCCESS);
1873: }
1875: /* MatMult_SeqBAIJ_15_ver4 : All columns in the block are accessed at once */
1876: PetscErrorCode MatMult_SeqBAIJ_15_ver4(Mat A, Vec xx, Vec zz)
1877: {
1878: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
1879: PetscScalar *z = NULL, sum1, sum2, sum3, sum4, sum5, sum6, sum7, sum8, sum9, sum10, sum11, sum12, sum13, sum14, sum15;
1880: const PetscScalar *x, *xb;
1881: PetscScalar x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15, *zarray;
1882: const MatScalar *v;
1883: const PetscInt *ii, *ij = a->j, *idx;
1884: PetscInt mbs, i, j, n, *ridx = NULL;
1885: PetscBool usecprow = a->compressedrow.use;
1887: PetscFunctionBegin;
1888: PetscCall(VecGetArrayRead(xx, &x));
1889: PetscCall(VecGetArrayWrite(zz, &zarray));
1891: v = a->a;
1892: if (usecprow) {
1893: mbs = a->compressedrow.nrows;
1894: ii = a->compressedrow.i;
1895: ridx = a->compressedrow.rindex;
1896: PetscCall(PetscArrayzero(zarray, 15 * a->mbs));
1897: } else {
1898: mbs = a->mbs;
1899: ii = a->i;
1900: z = zarray;
1901: }
1903: for (i = 0; i < mbs; i++) {
1904: n = ii[i + 1] - ii[i];
1905: idx = ij + ii[i];
1906: sum1 = 0.0;
1907: sum2 = 0.0;
1908: sum3 = 0.0;
1909: sum4 = 0.0;
1910: sum5 = 0.0;
1911: sum6 = 0.0;
1912: sum7 = 0.0;
1913: sum8 = 0.0;
1914: sum9 = 0.0;
1915: sum10 = 0.0;
1916: sum11 = 0.0;
1917: sum12 = 0.0;
1918: sum13 = 0.0;
1919: sum14 = 0.0;
1920: sum15 = 0.0;
1922: for (j = 0; j < n; j++) {
1923: xb = x + 15 * idx[j];
1924: x1 = xb[0];
1925: x2 = xb[1];
1926: x3 = xb[2];
1927: x4 = xb[3];
1928: x5 = xb[4];
1929: x6 = xb[5];
1930: x7 = xb[6];
1931: x8 = xb[7];
1932: x9 = xb[8];
1933: x10 = xb[9];
1934: x11 = xb[10];
1935: x12 = xb[11];
1936: x13 = xb[12];
1937: x14 = xb[13];
1938: x15 = xb[14];
1940: sum1 += v[0] * x1 + v[15] * x2 + v[30] * x3 + v[45] * x4 + v[60] * x5 + v[75] * x6 + v[90] * x7 + v[105] * x8 + v[120] * x9 + v[135] * x10 + v[150] * x11 + v[165] * x12 + v[180] * x13 + v[195] * x14 + v[210] * x15;
1941: sum2 += v[1] * x1 + v[16] * x2 + v[31] * x3 + v[46] * x4 + v[61] * x5 + v[76] * x6 + v[91] * x7 + v[106] * x8 + v[121] * x9 + v[136] * x10 + v[151] * x11 + v[166] * x12 + v[181] * x13 + v[196] * x14 + v[211] * x15;
1942: sum3 += v[2] * x1 + v[17] * x2 + v[32] * x3 + v[47] * x4 + v[62] * x5 + v[77] * x6 + v[92] * x7 + v[107] * x8 + v[122] * x9 + v[137] * x10 + v[152] * x11 + v[167] * x12 + v[182] * x13 + v[197] * x14 + v[212] * x15;
1943: sum4 += v[3] * x1 + v[18] * x2 + v[33] * x3 + v[48] * x4 + v[63] * x5 + v[78] * x6 + v[93] * x7 + v[108] * x8 + v[123] * x9 + v[138] * x10 + v[153] * x11 + v[168] * x12 + v[183] * x13 + v[198] * x14 + v[213] * x15;
1944: sum5 += v[4] * x1 + v[19] * x2 + v[34] * x3 + v[49] * x4 + v[64] * x5 + v[79] * x6 + v[94] * x7 + v[109] * x8 + v[124] * x9 + v[139] * x10 + v[154] * x11 + v[169] * x12 + v[184] * x13 + v[199] * x14 + v[214] * x15;
1945: sum6 += v[5] * x1 + v[20] * x2 + v[35] * x3 + v[50] * x4 + v[65] * x5 + v[80] * x6 + v[95] * x7 + v[110] * x8 + v[125] * x9 + v[140] * x10 + v[155] * x11 + v[170] * x12 + v[185] * x13 + v[200] * x14 + v[215] * x15;
1946: sum7 += v[6] * x1 + v[21] * x2 + v[36] * x3 + v[51] * x4 + v[66] * x5 + v[81] * x6 + v[96] * x7 + v[111] * x8 + v[126] * x9 + v[141] * x10 + v[156] * x11 + v[171] * x12 + v[186] * x13 + v[201] * x14 + v[216] * x15;
1947: sum8 += v[7] * x1 + v[22] * x2 + v[37] * x3 + v[52] * x4 + v[67] * x5 + v[82] * x6 + v[97] * x7 + v[112] * x8 + v[127] * x9 + v[142] * x10 + v[157] * x11 + v[172] * x12 + v[187] * x13 + v[202] * x14 + v[217] * x15;
1948: sum9 += v[8] * x1 + v[23] * x2 + v[38] * x3 + v[53] * x4 + v[68] * x5 + v[83] * x6 + v[98] * x7 + v[113] * x8 + v[128] * x9 + v[143] * x10 + v[158] * x11 + v[173] * x12 + v[188] * x13 + v[203] * x14 + v[218] * x15;
1949: sum10 += v[9] * x1 + v[24] * x2 + v[39] * x3 + v[54] * x4 + v[69] * x5 + v[84] * x6 + v[99] * x7 + v[114] * x8 + v[129] * x9 + v[144] * x10 + v[159] * x11 + v[174] * x12 + v[189] * x13 + v[204] * x14 + v[219] * x15;
1950: sum11 += v[10] * x1 + v[25] * x2 + v[40] * x3 + v[55] * x4 + v[70] * x5 + v[85] * x6 + v[100] * x7 + v[115] * x8 + v[130] * x9 + v[145] * x10 + v[160] * x11 + v[175] * x12 + v[190] * x13 + v[205] * x14 + v[220] * x15;
1951: sum12 += v[11] * x1 + v[26] * x2 + v[41] * x3 + v[56] * x4 + v[71] * x5 + v[86] * x6 + v[101] * x7 + v[116] * x8 + v[131] * x9 + v[146] * x10 + v[161] * x11 + v[176] * x12 + v[191] * x13 + v[206] * x14 + v[221] * x15;
1952: sum13 += v[12] * x1 + v[27] * x2 + v[42] * x3 + v[57] * x4 + v[72] * x5 + v[87] * x6 + v[102] * x7 + v[117] * x8 + v[132] * x9 + v[147] * x10 + v[162] * x11 + v[177] * x12 + v[192] * x13 + v[207] * x14 + v[222] * x15;
1953: sum14 += v[13] * x1 + v[28] * x2 + v[43] * x3 + v[58] * x4 + v[73] * x5 + v[88] * x6 + v[103] * x7 + v[118] * x8 + v[133] * x9 + v[148] * x10 + v[163] * x11 + v[178] * x12 + v[193] * x13 + v[208] * x14 + v[223] * x15;
1954: sum15 += v[14] * x1 + v[29] * x2 + v[44] * x3 + v[59] * x4 + v[74] * x5 + v[89] * x6 + v[104] * x7 + v[119] * x8 + v[134] * x9 + v[149] * x10 + v[164] * x11 + v[179] * x12 + v[194] * x13 + v[209] * x14 + v[224] * x15;
1955: v += 225;
1956: }
1957: if (usecprow) z = zarray + 15 * ridx[i];
1958: z[0] = sum1;
1959: z[1] = sum2;
1960: z[2] = sum3;
1961: z[3] = sum4;
1962: z[4] = sum5;
1963: z[5] = sum6;
1964: z[6] = sum7;
1965: z[7] = sum8;
1966: z[8] = sum9;
1967: z[9] = sum10;
1968: z[10] = sum11;
1969: z[11] = sum12;
1970: z[12] = sum13;
1971: z[13] = sum14;
1972: z[14] = sum15;
1974: if (!usecprow) z += 15;
1975: }
1977: PetscCall(VecRestoreArrayRead(xx, &x));
1978: PetscCall(VecRestoreArrayWrite(zz, &zarray));
1979: PetscCall(PetscLogFlops(450.0 * a->nz - 15.0 * a->nonzerorowcnt));
1980: PetscFunctionReturn(PETSC_SUCCESS);
1981: }
1983: /*
1984: This will not work with MatScalar == float because it calls the BLAS
1985: */
1986: PetscErrorCode MatMult_SeqBAIJ_N(Mat A, Vec xx, Vec zz)
1987: {
1988: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
1989: PetscScalar *z = NULL, *work, *workt, *zarray;
1990: const PetscScalar *x, *xb;
1991: const MatScalar *v;
1992: PetscInt mbs, i, bs = A->rmap->bs, j, n, bs2 = a->bs2;
1993: const PetscInt *idx, *ii, *ridx = NULL;
1994: PetscInt ncols, k;
1995: PetscBool usecprow = a->compressedrow.use;
1997: PetscFunctionBegin;
1998: PetscCall(VecGetArrayRead(xx, &x));
1999: PetscCall(VecGetArrayWrite(zz, &zarray));
2001: idx = a->j;
2002: v = a->a;
2003: if (usecprow) {
2004: mbs = a->compressedrow.nrows;
2005: ii = a->compressedrow.i;
2006: ridx = a->compressedrow.rindex;
2007: PetscCall(PetscArrayzero(zarray, bs * a->mbs));
2008: } else {
2009: mbs = a->mbs;
2010: ii = a->i;
2011: z = zarray;
2012: }
2014: if (!a->mult_work) {
2015: k = PetscMax(A->rmap->n, A->cmap->n);
2016: PetscCall(PetscMalloc1(k + 1, &a->mult_work));
2017: }
2018: work = a->mult_work;
2019: for (i = 0; i < mbs; i++) {
2020: n = ii[1] - ii[0];
2021: ii++;
2022: ncols = n * bs;
2023: workt = work;
2024: for (j = 0; j < n; j++) {
2025: xb = x + bs * (*idx++);
2026: for (k = 0; k < bs; k++) workt[k] = xb[k];
2027: workt += bs;
2028: }
2029: if (usecprow) z = zarray + bs * ridx[i];
2030: PetscKernel_w_gets_Ar_times_v(bs, ncols, work, v, z);
2031: v += n * bs2;
2032: if (!usecprow) z += bs;
2033: }
2034: PetscCall(VecRestoreArrayRead(xx, &x));
2035: PetscCall(VecRestoreArrayWrite(zz, &zarray));
2036: PetscCall(PetscLogFlops(2.0 * a->nz * bs2 - bs * a->nonzerorowcnt));
2037: PetscFunctionReturn(PETSC_SUCCESS);
2038: }
2040: PetscErrorCode MatMultAdd_SeqBAIJ_1(Mat A, Vec xx, Vec yy, Vec zz)
2041: {
2042: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
2043: const PetscScalar *x;
2044: PetscScalar *y, *z, sum;
2045: const MatScalar *v;
2046: PetscInt mbs = a->mbs, i, n, *ridx = NULL;
2047: const PetscInt *idx, *ii;
2048: PetscBool usecprow = a->compressedrow.use;
2050: PetscFunctionBegin;
2051: PetscCall(VecGetArrayRead(xx, &x));
2052: PetscCall(VecGetArrayPair(yy, zz, &y, &z));
2054: idx = a->j;
2055: v = a->a;
2056: if (usecprow) {
2057: if (zz != yy) PetscCall(PetscArraycpy(z, y, mbs));
2058: mbs = a->compressedrow.nrows;
2059: ii = a->compressedrow.i;
2060: ridx = a->compressedrow.rindex;
2061: } else {
2062: ii = a->i;
2063: }
2065: for (i = 0; i < mbs; i++) {
2066: n = ii[1] - ii[0];
2067: ii++;
2068: if (!usecprow) {
2069: sum = y[i];
2070: } else {
2071: sum = y[ridx[i]];
2072: }
2073: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
2074: PetscPrefetchBlock(v + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
2075: PetscSparseDensePlusDot(sum, x, v, idx, n);
2076: v += n;
2077: idx += n;
2078: if (usecprow) {
2079: z[ridx[i]] = sum;
2080: } else {
2081: z[i] = sum;
2082: }
2083: }
2084: PetscCall(VecRestoreArrayRead(xx, &x));
2085: PetscCall(VecRestoreArrayPair(yy, zz, &y, &z));
2086: PetscCall(PetscLogFlops(2.0 * a->nz));
2087: PetscFunctionReturn(PETSC_SUCCESS);
2088: }
2090: PetscErrorCode MatMultAdd_SeqBAIJ_2(Mat A, Vec xx, Vec yy, Vec zz)
2091: {
2092: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
2093: PetscScalar *y = NULL, *z = NULL, sum1, sum2;
2094: const PetscScalar *x, *xb;
2095: PetscScalar x1, x2, *yarray, *zarray;
2096: const MatScalar *v;
2097: PetscInt mbs = a->mbs, i, n, j;
2098: const PetscInt *idx, *ii, *ridx = NULL;
2099: PetscBool usecprow = a->compressedrow.use;
2101: PetscFunctionBegin;
2102: PetscCall(VecGetArrayRead(xx, &x));
2103: PetscCall(VecGetArrayPair(yy, zz, &yarray, &zarray));
2105: idx = a->j;
2106: v = a->a;
2107: if (usecprow) {
2108: if (zz != yy) PetscCall(PetscArraycpy(zarray, yarray, 2 * mbs));
2109: mbs = a->compressedrow.nrows;
2110: ii = a->compressedrow.i;
2111: ridx = a->compressedrow.rindex;
2112: } else {
2113: ii = a->i;
2114: y = yarray;
2115: z = zarray;
2116: }
2118: for (i = 0; i < mbs; i++) {
2119: n = ii[1] - ii[0];
2120: ii++;
2121: if (usecprow) {
2122: z = zarray + 2 * ridx[i];
2123: y = yarray + 2 * ridx[i];
2124: }
2125: sum1 = y[0];
2126: sum2 = y[1];
2127: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
2128: PetscPrefetchBlock(v + 4 * n, 4 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
2129: for (j = 0; j < n; j++) {
2130: xb = x + 2 * (*idx++);
2131: x1 = xb[0];
2132: x2 = xb[1];
2134: sum1 += v[0] * x1 + v[2] * x2;
2135: sum2 += v[1] * x1 + v[3] * x2;
2136: v += 4;
2137: }
2138: z[0] = sum1;
2139: z[1] = sum2;
2140: if (!usecprow) {
2141: z += 2;
2142: y += 2;
2143: }
2144: }
2145: PetscCall(VecRestoreArrayRead(xx, &x));
2146: PetscCall(VecRestoreArrayPair(yy, zz, &yarray, &zarray));
2147: PetscCall(PetscLogFlops(4.0 * a->nz));
2148: PetscFunctionReturn(PETSC_SUCCESS);
2149: }
2151: PetscErrorCode MatMultAdd_SeqBAIJ_3(Mat A, Vec xx, Vec yy, Vec zz)
2152: {
2153: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
2154: PetscScalar *y = NULL, *z = NULL, sum1, sum2, sum3, x1, x2, x3, *yarray, *zarray;
2155: const PetscScalar *x, *xb;
2156: const MatScalar *v;
2157: PetscInt mbs = a->mbs, i, j, n;
2158: const PetscInt *idx, *ii, *ridx = NULL;
2159: PetscBool usecprow = a->compressedrow.use;
2161: PetscFunctionBegin;
2162: PetscCall(VecGetArrayRead(xx, &x));
2163: PetscCall(VecGetArrayPair(yy, zz, &yarray, &zarray));
2165: idx = a->j;
2166: v = a->a;
2167: if (usecprow) {
2168: if (zz != yy) PetscCall(PetscArraycpy(zarray, yarray, 3 * mbs));
2169: mbs = a->compressedrow.nrows;
2170: ii = a->compressedrow.i;
2171: ridx = a->compressedrow.rindex;
2172: } else {
2173: ii = a->i;
2174: y = yarray;
2175: z = zarray;
2176: }
2178: for (i = 0; i < mbs; i++) {
2179: n = ii[1] - ii[0];
2180: ii++;
2181: if (usecprow) {
2182: z = zarray + 3 * ridx[i];
2183: y = yarray + 3 * ridx[i];
2184: }
2185: sum1 = y[0];
2186: sum2 = y[1];
2187: sum3 = y[2];
2188: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
2189: PetscPrefetchBlock(v + 9 * n, 9 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
2190: for (j = 0; j < n; j++) {
2191: xb = x + 3 * (*idx++);
2192: x1 = xb[0];
2193: x2 = xb[1];
2194: x3 = xb[2];
2195: sum1 += v[0] * x1 + v[3] * x2 + v[6] * x3;
2196: sum2 += v[1] * x1 + v[4] * x2 + v[7] * x3;
2197: sum3 += v[2] * x1 + v[5] * x2 + v[8] * x3;
2198: v += 9;
2199: }
2200: z[0] = sum1;
2201: z[1] = sum2;
2202: z[2] = sum3;
2203: if (!usecprow) {
2204: z += 3;
2205: y += 3;
2206: }
2207: }
2208: PetscCall(VecRestoreArrayRead(xx, &x));
2209: PetscCall(VecRestoreArrayPair(yy, zz, &yarray, &zarray));
2210: PetscCall(PetscLogFlops(18.0 * a->nz));
2211: PetscFunctionReturn(PETSC_SUCCESS);
2212: }
2214: PetscErrorCode MatMultAdd_SeqBAIJ_4(Mat A, Vec xx, Vec yy, Vec zz)
2215: {
2216: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
2217: PetscScalar *y = NULL, *z = NULL, sum1, sum2, sum3, sum4, x1, x2, x3, x4, *yarray, *zarray;
2218: const PetscScalar *x, *xb;
2219: const MatScalar *v;
2220: PetscInt mbs = a->mbs, i, j, n;
2221: const PetscInt *idx, *ii, *ridx = NULL;
2222: PetscBool usecprow = a->compressedrow.use;
2224: PetscFunctionBegin;
2225: PetscCall(VecGetArrayRead(xx, &x));
2226: PetscCall(VecGetArrayPair(yy, zz, &yarray, &zarray));
2228: idx = a->j;
2229: v = a->a;
2230: if (usecprow) {
2231: if (zz != yy) PetscCall(PetscArraycpy(zarray, yarray, 4 * mbs));
2232: mbs = a->compressedrow.nrows;
2233: ii = a->compressedrow.i;
2234: ridx = a->compressedrow.rindex;
2235: } else {
2236: ii = a->i;
2237: y = yarray;
2238: z = zarray;
2239: }
2241: for (i = 0; i < mbs; i++) {
2242: n = ii[1] - ii[0];
2243: ii++;
2244: if (usecprow) {
2245: z = zarray + 4 * ridx[i];
2246: y = yarray + 4 * ridx[i];
2247: }
2248: sum1 = y[0];
2249: sum2 = y[1];
2250: sum3 = y[2];
2251: sum4 = y[3];
2252: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
2253: PetscPrefetchBlock(v + 16 * n, 16 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
2254: for (j = 0; j < n; j++) {
2255: xb = x + 4 * (*idx++);
2256: x1 = xb[0];
2257: x2 = xb[1];
2258: x3 = xb[2];
2259: x4 = xb[3];
2260: sum1 += v[0] * x1 + v[4] * x2 + v[8] * x3 + v[12] * x4;
2261: sum2 += v[1] * x1 + v[5] * x2 + v[9] * x3 + v[13] * x4;
2262: sum3 += v[2] * x1 + v[6] * x2 + v[10] * x3 + v[14] * x4;
2263: sum4 += v[3] * x1 + v[7] * x2 + v[11] * x3 + v[15] * x4;
2264: v += 16;
2265: }
2266: z[0] = sum1;
2267: z[1] = sum2;
2268: z[2] = sum3;
2269: z[3] = sum4;
2270: if (!usecprow) {
2271: z += 4;
2272: y += 4;
2273: }
2274: }
2275: PetscCall(VecRestoreArrayRead(xx, &x));
2276: PetscCall(VecRestoreArrayPair(yy, zz, &yarray, &zarray));
2277: PetscCall(PetscLogFlops(32.0 * a->nz));
2278: PetscFunctionReturn(PETSC_SUCCESS);
2279: }
2281: PetscErrorCode MatMultAdd_SeqBAIJ_5(Mat A, Vec xx, Vec yy, Vec zz)
2282: {
2283: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
2284: PetscScalar *y = NULL, *z = NULL, sum1, sum2, sum3, sum4, sum5, x1, x2, x3, x4, x5;
2285: const PetscScalar *x, *xb;
2286: PetscScalar *yarray, *zarray;
2287: const MatScalar *v;
2288: PetscInt mbs = a->mbs, i, j, n;
2289: const PetscInt *idx, *ii, *ridx = NULL;
2290: PetscBool usecprow = a->compressedrow.use;
2292: PetscFunctionBegin;
2293: PetscCall(VecGetArrayRead(xx, &x));
2294: PetscCall(VecGetArrayPair(yy, zz, &yarray, &zarray));
2296: idx = a->j;
2297: v = a->a;
2298: if (usecprow) {
2299: if (zz != yy) PetscCall(PetscArraycpy(zarray, yarray, 5 * mbs));
2300: mbs = a->compressedrow.nrows;
2301: ii = a->compressedrow.i;
2302: ridx = a->compressedrow.rindex;
2303: } else {
2304: ii = a->i;
2305: y = yarray;
2306: z = zarray;
2307: }
2309: for (i = 0; i < mbs; i++) {
2310: n = ii[1] - ii[0];
2311: ii++;
2312: if (usecprow) {
2313: z = zarray + 5 * ridx[i];
2314: y = yarray + 5 * ridx[i];
2315: }
2316: sum1 = y[0];
2317: sum2 = y[1];
2318: sum3 = y[2];
2319: sum4 = y[3];
2320: sum5 = y[4];
2321: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
2322: PetscPrefetchBlock(v + 25 * n, 25 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
2323: for (j = 0; j < n; j++) {
2324: xb = x + 5 * (*idx++);
2325: x1 = xb[0];
2326: x2 = xb[1];
2327: x3 = xb[2];
2328: x4 = xb[3];
2329: x5 = xb[4];
2330: sum1 += v[0] * x1 + v[5] * x2 + v[10] * x3 + v[15] * x4 + v[20] * x5;
2331: sum2 += v[1] * x1 + v[6] * x2 + v[11] * x3 + v[16] * x4 + v[21] * x5;
2332: sum3 += v[2] * x1 + v[7] * x2 + v[12] * x3 + v[17] * x4 + v[22] * x5;
2333: sum4 += v[3] * x1 + v[8] * x2 + v[13] * x3 + v[18] * x4 + v[23] * x5;
2334: sum5 += v[4] * x1 + v[9] * x2 + v[14] * x3 + v[19] * x4 + v[24] * x5;
2335: v += 25;
2336: }
2337: z[0] = sum1;
2338: z[1] = sum2;
2339: z[2] = sum3;
2340: z[3] = sum4;
2341: z[4] = sum5;
2342: if (!usecprow) {
2343: z += 5;
2344: y += 5;
2345: }
2346: }
2347: PetscCall(VecRestoreArrayRead(xx, &x));
2348: PetscCall(VecRestoreArrayPair(yy, zz, &yarray, &zarray));
2349: PetscCall(PetscLogFlops(50.0 * a->nz));
2350: PetscFunctionReturn(PETSC_SUCCESS);
2351: }
2353: PetscErrorCode MatMultAdd_SeqBAIJ_6(Mat A, Vec xx, Vec yy, Vec zz)
2354: {
2355: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
2356: PetscScalar *y = NULL, *z = NULL, sum1, sum2, sum3, sum4, sum5, sum6;
2357: const PetscScalar *x, *xb;
2358: PetscScalar x1, x2, x3, x4, x5, x6, *yarray, *zarray;
2359: const MatScalar *v;
2360: PetscInt mbs = a->mbs, i, j, n;
2361: const PetscInt *idx, *ii, *ridx = NULL;
2362: PetscBool usecprow = a->compressedrow.use;
2364: PetscFunctionBegin;
2365: PetscCall(VecGetArrayRead(xx, &x));
2366: PetscCall(VecGetArrayPair(yy, zz, &yarray, &zarray));
2368: idx = a->j;
2369: v = a->a;
2370: if (usecprow) {
2371: if (zz != yy) PetscCall(PetscArraycpy(zarray, yarray, 6 * mbs));
2372: mbs = a->compressedrow.nrows;
2373: ii = a->compressedrow.i;
2374: ridx = a->compressedrow.rindex;
2375: } else {
2376: ii = a->i;
2377: y = yarray;
2378: z = zarray;
2379: }
2381: for (i = 0; i < mbs; i++) {
2382: n = ii[1] - ii[0];
2383: ii++;
2384: if (usecprow) {
2385: z = zarray + 6 * ridx[i];
2386: y = yarray + 6 * ridx[i];
2387: }
2388: sum1 = y[0];
2389: sum2 = y[1];
2390: sum3 = y[2];
2391: sum4 = y[3];
2392: sum5 = y[4];
2393: sum6 = y[5];
2394: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
2395: PetscPrefetchBlock(v + 36 * n, 36 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
2396: for (j = 0; j < n; j++) {
2397: xb = x + 6 * (*idx++);
2398: x1 = xb[0];
2399: x2 = xb[1];
2400: x3 = xb[2];
2401: x4 = xb[3];
2402: x5 = xb[4];
2403: x6 = xb[5];
2404: sum1 += v[0] * x1 + v[6] * x2 + v[12] * x3 + v[18] * x4 + v[24] * x5 + v[30] * x6;
2405: sum2 += v[1] * x1 + v[7] * x2 + v[13] * x3 + v[19] * x4 + v[25] * x5 + v[31] * x6;
2406: sum3 += v[2] * x1 + v[8] * x2 + v[14] * x3 + v[20] * x4 + v[26] * x5 + v[32] * x6;
2407: sum4 += v[3] * x1 + v[9] * x2 + v[15] * x3 + v[21] * x4 + v[27] * x5 + v[33] * x6;
2408: sum5 += v[4] * x1 + v[10] * x2 + v[16] * x3 + v[22] * x4 + v[28] * x5 + v[34] * x6;
2409: sum6 += v[5] * x1 + v[11] * x2 + v[17] * x3 + v[23] * x4 + v[29] * x5 + v[35] * x6;
2410: v += 36;
2411: }
2412: z[0] = sum1;
2413: z[1] = sum2;
2414: z[2] = sum3;
2415: z[3] = sum4;
2416: z[4] = sum5;
2417: z[5] = sum6;
2418: if (!usecprow) {
2419: z += 6;
2420: y += 6;
2421: }
2422: }
2423: PetscCall(VecRestoreArrayRead(xx, &x));
2424: PetscCall(VecRestoreArrayPair(yy, zz, &yarray, &zarray));
2425: PetscCall(PetscLogFlops(72.0 * a->nz));
2426: PetscFunctionReturn(PETSC_SUCCESS);
2427: }
2429: PetscErrorCode MatMultAdd_SeqBAIJ_7(Mat A, Vec xx, Vec yy, Vec zz)
2430: {
2431: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
2432: PetscScalar *y = NULL, *z = NULL, sum1, sum2, sum3, sum4, sum5, sum6, sum7;
2433: const PetscScalar *x, *xb;
2434: PetscScalar x1, x2, x3, x4, x5, x6, x7, *yarray, *zarray;
2435: const MatScalar *v;
2436: PetscInt mbs = a->mbs, i, j, n;
2437: const PetscInt *idx, *ii, *ridx = NULL;
2438: PetscBool usecprow = a->compressedrow.use;
2440: PetscFunctionBegin;
2441: PetscCall(VecGetArrayRead(xx, &x));
2442: PetscCall(VecGetArrayPair(yy, zz, &yarray, &zarray));
2444: idx = a->j;
2445: v = a->a;
2446: if (usecprow) {
2447: if (zz != yy) PetscCall(PetscArraycpy(zarray, yarray, 7 * mbs));
2448: mbs = a->compressedrow.nrows;
2449: ii = a->compressedrow.i;
2450: ridx = a->compressedrow.rindex;
2451: } else {
2452: ii = a->i;
2453: y = yarray;
2454: z = zarray;
2455: }
2457: for (i = 0; i < mbs; i++) {
2458: n = ii[1] - ii[0];
2459: ii++;
2460: if (usecprow) {
2461: z = zarray + 7 * ridx[i];
2462: y = yarray + 7 * ridx[i];
2463: }
2464: sum1 = y[0];
2465: sum2 = y[1];
2466: sum3 = y[2];
2467: sum4 = y[3];
2468: sum5 = y[4];
2469: sum6 = y[5];
2470: sum7 = y[6];
2471: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
2472: PetscPrefetchBlock(v + 49 * n, 49 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
2473: for (j = 0; j < n; j++) {
2474: xb = x + 7 * (*idx++);
2475: x1 = xb[0];
2476: x2 = xb[1];
2477: x3 = xb[2];
2478: x4 = xb[3];
2479: x5 = xb[4];
2480: x6 = xb[5];
2481: x7 = xb[6];
2482: sum1 += v[0] * x1 + v[7] * x2 + v[14] * x3 + v[21] * x4 + v[28] * x5 + v[35] * x6 + v[42] * x7;
2483: sum2 += v[1] * x1 + v[8] * x2 + v[15] * x3 + v[22] * x4 + v[29] * x5 + v[36] * x6 + v[43] * x7;
2484: sum3 += v[2] * x1 + v[9] * x2 + v[16] * x3 + v[23] * x4 + v[30] * x5 + v[37] * x6 + v[44] * x7;
2485: sum4 += v[3] * x1 + v[10] * x2 + v[17] * x3 + v[24] * x4 + v[31] * x5 + v[38] * x6 + v[45] * x7;
2486: sum5 += v[4] * x1 + v[11] * x2 + v[18] * x3 + v[25] * x4 + v[32] * x5 + v[39] * x6 + v[46] * x7;
2487: sum6 += v[5] * x1 + v[12] * x2 + v[19] * x3 + v[26] * x4 + v[33] * x5 + v[40] * x6 + v[47] * x7;
2488: sum7 += v[6] * x1 + v[13] * x2 + v[20] * x3 + v[27] * x4 + v[34] * x5 + v[41] * x6 + v[48] * x7;
2489: v += 49;
2490: }
2491: z[0] = sum1;
2492: z[1] = sum2;
2493: z[2] = sum3;
2494: z[3] = sum4;
2495: z[4] = sum5;
2496: z[5] = sum6;
2497: z[6] = sum7;
2498: if (!usecprow) {
2499: z += 7;
2500: y += 7;
2501: }
2502: }
2503: PetscCall(VecRestoreArrayRead(xx, &x));
2504: PetscCall(VecRestoreArrayPair(yy, zz, &yarray, &zarray));
2505: PetscCall(PetscLogFlops(98.0 * a->nz));
2506: PetscFunctionReturn(PETSC_SUCCESS);
2507: }
2509: #if PetscDefined(HAVE_IMMINTRIN_H) && defined(__AVX2__) && defined(__FMA__) && PetscDefined(USE_REAL_DOUBLE) && !PetscDefined(USE_COMPLEX) && !PetscDefined(USE_64BIT_INDICES)
2510: PetscErrorCode MatMultAdd_SeqBAIJ_9_AVX2(Mat A, Vec xx, Vec yy, Vec zz)
2511: {
2512: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
2513: PetscScalar *z = NULL, *work, *workt, *zarray;
2514: const PetscScalar *x, *xb;
2515: const MatScalar *v;
2516: PetscInt mbs, i, j, n;
2517: PetscInt k;
2518: PetscBool usecprow = a->compressedrow.use;
2519: const PetscInt *idx, *ii, *ridx = NULL, bs = 9, bs2 = 81;
2521: __m256d a0, a1, a2, a3, a4, a5;
2522: __m256d w0, w1, w2, w3;
2523: __m256d z0, z1, z2;
2524: __m256i mask1 = _mm256_set_epi64x(0LL, 0LL, 0LL, 1LL << 63);
2526: PetscFunctionBegin;
2527: PetscCall(VecCopy(yy, zz));
2528: PetscCall(VecGetArrayRead(xx, &x));
2529: PetscCall(VecGetArray(zz, &zarray));
2531: idx = a->j;
2532: v = a->a;
2533: if (usecprow) {
2534: mbs = a->compressedrow.nrows;
2535: ii = a->compressedrow.i;
2536: ridx = a->compressedrow.rindex;
2537: } else {
2538: mbs = a->mbs;
2539: ii = a->i;
2540: z = zarray;
2541: }
2543: if (!a->mult_work) {
2544: k = PetscMax(A->rmap->n, A->cmap->n);
2545: PetscCall(PetscMalloc1(k + 1, &a->mult_work));
2546: }
2548: work = a->mult_work;
2549: for (i = 0; i < mbs; i++) {
2550: n = ii[1] - ii[0];
2551: ii++;
2552: workt = work;
2553: for (j = 0; j < n; j++) {
2554: xb = x + bs * (*idx++);
2555: for (k = 0; k < bs; k++) workt[k] = xb[k];
2556: workt += bs;
2557: }
2558: if (usecprow) z = zarray + bs * ridx[i];
2560: z0 = _mm256_loadu_pd(&z[0]);
2561: z1 = _mm256_loadu_pd(&z[4]);
2562: z2 = _mm256_set1_pd(z[8]);
2564: for (j = 0; j < n; j++) {
2565: /* first column of a */
2566: w0 = _mm256_set1_pd(work[j * 9]);
2567: a0 = _mm256_loadu_pd(&v[j * 81]);
2568: z0 = _mm256_fmadd_pd(a0, w0, z0);
2569: a1 = _mm256_loadu_pd(&v[j * 81 + 4]);
2570: z1 = _mm256_fmadd_pd(a1, w0, z1);
2571: a2 = _mm256_loadu_pd(&v[j * 81 + 8]);
2572: z2 = _mm256_fmadd_pd(a2, w0, z2);
2574: /* second column of a */
2575: w1 = _mm256_set1_pd(work[j * 9 + 1]);
2576: a0 = _mm256_loadu_pd(&v[j * 81 + 9]);
2577: z0 = _mm256_fmadd_pd(a0, w1, z0);
2578: a1 = _mm256_loadu_pd(&v[j * 81 + 13]);
2579: z1 = _mm256_fmadd_pd(a1, w1, z1);
2580: a2 = _mm256_loadu_pd(&v[j * 81 + 17]);
2581: z2 = _mm256_fmadd_pd(a2, w1, z2);
2583: /* third column of a */
2584: w2 = _mm256_set1_pd(work[j * 9 + 2]);
2585: a3 = _mm256_loadu_pd(&v[j * 81 + 18]);
2586: z0 = _mm256_fmadd_pd(a3, w2, z0);
2587: a4 = _mm256_loadu_pd(&v[j * 81 + 22]);
2588: z1 = _mm256_fmadd_pd(a4, w2, z1);
2589: a5 = _mm256_loadu_pd(&v[j * 81 + 26]);
2590: z2 = _mm256_fmadd_pd(a5, w2, z2);
2592: /* fourth column of a */
2593: w3 = _mm256_set1_pd(work[j * 9 + 3]);
2594: a0 = _mm256_loadu_pd(&v[j * 81 + 27]);
2595: z0 = _mm256_fmadd_pd(a0, w3, z0);
2596: a1 = _mm256_loadu_pd(&v[j * 81 + 31]);
2597: z1 = _mm256_fmadd_pd(a1, w3, z1);
2598: a2 = _mm256_loadu_pd(&v[j * 81 + 35]);
2599: z2 = _mm256_fmadd_pd(a2, w3, z2);
2601: /* fifth column of a */
2602: w0 = _mm256_set1_pd(work[j * 9 + 4]);
2603: a3 = _mm256_loadu_pd(&v[j * 81 + 36]);
2604: z0 = _mm256_fmadd_pd(a3, w0, z0);
2605: a4 = _mm256_loadu_pd(&v[j * 81 + 40]);
2606: z1 = _mm256_fmadd_pd(a4, w0, z1);
2607: a5 = _mm256_loadu_pd(&v[j * 81 + 44]);
2608: z2 = _mm256_fmadd_pd(a5, w0, z2);
2610: /* sixth column of a */
2611: w1 = _mm256_set1_pd(work[j * 9 + 5]);
2612: a0 = _mm256_loadu_pd(&v[j * 81 + 45]);
2613: z0 = _mm256_fmadd_pd(a0, w1, z0);
2614: a1 = _mm256_loadu_pd(&v[j * 81 + 49]);
2615: z1 = _mm256_fmadd_pd(a1, w1, z1);
2616: a2 = _mm256_loadu_pd(&v[j * 81 + 53]);
2617: z2 = _mm256_fmadd_pd(a2, w1, z2);
2619: /* seventh column of a */
2620: w2 = _mm256_set1_pd(work[j * 9 + 6]);
2621: a0 = _mm256_loadu_pd(&v[j * 81 + 54]);
2622: z0 = _mm256_fmadd_pd(a0, w2, z0);
2623: a1 = _mm256_loadu_pd(&v[j * 81 + 58]);
2624: z1 = _mm256_fmadd_pd(a1, w2, z1);
2625: a2 = _mm256_loadu_pd(&v[j * 81 + 62]);
2626: z2 = _mm256_fmadd_pd(a2, w2, z2);
2628: /* eighth column of a */
2629: w3 = _mm256_set1_pd(work[j * 9 + 7]);
2630: a3 = _mm256_loadu_pd(&v[j * 81 + 63]);
2631: z0 = _mm256_fmadd_pd(a3, w3, z0);
2632: a4 = _mm256_loadu_pd(&v[j * 81 + 67]);
2633: z1 = _mm256_fmadd_pd(a4, w3, z1);
2634: a5 = _mm256_loadu_pd(&v[j * 81 + 71]);
2635: z2 = _mm256_fmadd_pd(a5, w3, z2);
2637: /* ninth column of a */
2638: w0 = _mm256_set1_pd(work[j * 9 + 8]);
2639: a0 = _mm256_loadu_pd(&v[j * 81 + 72]);
2640: z0 = _mm256_fmadd_pd(a0, w0, z0);
2641: a1 = _mm256_loadu_pd(&v[j * 81 + 76]);
2642: z1 = _mm256_fmadd_pd(a1, w0, z1);
2643: a2 = _mm256_maskload_pd(&v[j * 81 + 80], mask1);
2644: z2 = _mm256_fmadd_pd(a2, w0, z2);
2645: }
2647: _mm256_storeu_pd(&z[0], z0);
2648: _mm256_storeu_pd(&z[4], z1);
2649: _mm256_maskstore_pd(&z[8], mask1, z2);
2651: v += n * bs2;
2652: if (!usecprow) z += bs;
2653: }
2654: PetscCall(VecRestoreArrayRead(xx, &x));
2655: PetscCall(VecRestoreArray(zz, &zarray));
2656: PetscCall(PetscLogFlops(162.0 * a->nz));
2657: PetscFunctionReturn(PETSC_SUCCESS);
2658: }
2659: #endif
2661: PetscErrorCode MatMultAdd_SeqBAIJ_11(Mat A, Vec xx, Vec yy, Vec zz)
2662: {
2663: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
2664: PetscScalar *y = NULL, *z = NULL, sum1, sum2, sum3, sum4, sum5, sum6, sum7, sum8, sum9, sum10, sum11;
2665: const PetscScalar *x, *xb;
2666: PetscScalar x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, *yarray, *zarray;
2667: const MatScalar *v;
2668: PetscInt mbs = a->mbs, i, j, n;
2669: const PetscInt *idx, *ii, *ridx = NULL;
2670: PetscBool usecprow = a->compressedrow.use;
2672: PetscFunctionBegin;
2673: PetscCall(VecGetArrayRead(xx, &x));
2674: PetscCall(VecGetArrayPair(yy, zz, &yarray, &zarray));
2676: idx = a->j;
2677: v = a->a;
2678: if (usecprow) {
2679: if (zz != yy) PetscCall(PetscArraycpy(zarray, yarray, 7 * mbs));
2680: mbs = a->compressedrow.nrows;
2681: ii = a->compressedrow.i;
2682: ridx = a->compressedrow.rindex;
2683: } else {
2684: ii = a->i;
2685: y = yarray;
2686: z = zarray;
2687: }
2689: for (i = 0; i < mbs; i++) {
2690: n = ii[1] - ii[0];
2691: ii++;
2692: if (usecprow) {
2693: z = zarray + 11 * ridx[i];
2694: y = yarray + 11 * ridx[i];
2695: }
2696: sum1 = y[0];
2697: sum2 = y[1];
2698: sum3 = y[2];
2699: sum4 = y[3];
2700: sum5 = y[4];
2701: sum6 = y[5];
2702: sum7 = y[6];
2703: sum8 = y[7];
2704: sum9 = y[8];
2705: sum10 = y[9];
2706: sum11 = y[10];
2707: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
2708: PetscPrefetchBlock(v + 121 * n, 121 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
2709: for (j = 0; j < n; j++) {
2710: xb = x + 11 * (*idx++);
2711: x1 = xb[0];
2712: x2 = xb[1];
2713: x3 = xb[2];
2714: x4 = xb[3];
2715: x5 = xb[4];
2716: x6 = xb[5];
2717: x7 = xb[6];
2718: x8 = xb[7];
2719: x9 = xb[8];
2720: x10 = xb[9];
2721: x11 = xb[10];
2722: sum1 += v[0] * x1 + v[11] * x2 + v[2 * 11] * x3 + v[3 * 11] * x4 + v[4 * 11] * x5 + v[5 * 11] * x6 + v[6 * 11] * x7 + v[7 * 11] * x8 + v[8 * 11] * x9 + v[9 * 11] * x10 + v[10 * 11] * x11;
2723: sum2 += v[1 + 0] * x1 + v[1 + 11] * x2 + v[1 + 2 * 11] * x3 + v[1 + 3 * 11] * x4 + v[1 + 4 * 11] * x5 + v[1 + 5 * 11] * x6 + v[1 + 6 * 11] * x7 + v[1 + 7 * 11] * x8 + v[1 + 8 * 11] * x9 + v[1 + 9 * 11] * x10 + v[1 + 10 * 11] * x11;
2724: sum3 += v[2 + 0] * x1 + v[2 + 11] * x2 + v[2 + 2 * 11] * x3 + v[2 + 3 * 11] * x4 + v[2 + 4 * 11] * x5 + v[2 + 5 * 11] * x6 + v[2 + 6 * 11] * x7 + v[2 + 7 * 11] * x8 + v[2 + 8 * 11] * x9 + v[2 + 9 * 11] * x10 + v[2 + 10 * 11] * x11;
2725: sum4 += v[3 + 0] * x1 + v[3 + 11] * x2 + v[3 + 2 * 11] * x3 + v[3 + 3 * 11] * x4 + v[3 + 4 * 11] * x5 + v[3 + 5 * 11] * x6 + v[3 + 6 * 11] * x7 + v[3 + 7 * 11] * x8 + v[3 + 8 * 11] * x9 + v[3 + 9 * 11] * x10 + v[3 + 10 * 11] * x11;
2726: sum5 += v[4 + 0] * x1 + v[4 + 11] * x2 + v[4 + 2 * 11] * x3 + v[4 + 3 * 11] * x4 + v[4 + 4 * 11] * x5 + v[4 + 5 * 11] * x6 + v[4 + 6 * 11] * x7 + v[4 + 7 * 11] * x8 + v[4 + 8 * 11] * x9 + v[4 + 9 * 11] * x10 + v[4 + 10 * 11] * x11;
2727: sum6 += v[5 + 0] * x1 + v[5 + 11] * x2 + v[5 + 2 * 11] * x3 + v[5 + 3 * 11] * x4 + v[5 + 4 * 11] * x5 + v[5 + 5 * 11] * x6 + v[5 + 6 * 11] * x7 + v[5 + 7 * 11] * x8 + v[5 + 8 * 11] * x9 + v[5 + 9 * 11] * x10 + v[5 + 10 * 11] * x11;
2728: sum7 += v[6 + 0] * x1 + v[6 + 11] * x2 + v[6 + 2 * 11] * x3 + v[6 + 3 * 11] * x4 + v[6 + 4 * 11] * x5 + v[6 + 5 * 11] * x6 + v[6 + 6 * 11] * x7 + v[6 + 7 * 11] * x8 + v[6 + 8 * 11] * x9 + v[6 + 9 * 11] * x10 + v[6 + 10 * 11] * x11;
2729: sum8 += v[7 + 0] * x1 + v[7 + 11] * x2 + v[7 + 2 * 11] * x3 + v[7 + 3 * 11] * x4 + v[7 + 4 * 11] * x5 + v[7 + 5 * 11] * x6 + v[7 + 6 * 11] * x7 + v[7 + 7 * 11] * x8 + v[7 + 8 * 11] * x9 + v[7 + 9 * 11] * x10 + v[7 + 10 * 11] * x11;
2730: sum9 += v[8 + 0] * x1 + v[8 + 11] * x2 + v[8 + 2 * 11] * x3 + v[8 + 3 * 11] * x4 + v[8 + 4 * 11] * x5 + v[8 + 5 * 11] * x6 + v[8 + 6 * 11] * x7 + v[8 + 7 * 11] * x8 + v[8 + 8 * 11] * x9 + v[8 + 9 * 11] * x10 + v[8 + 10 * 11] * x11;
2731: sum10 += v[9 + 0] * x1 + v[9 + 11] * x2 + v[9 + 2 * 11] * x3 + v[9 + 3 * 11] * x4 + v[9 + 4 * 11] * x5 + v[9 + 5 * 11] * x6 + v[9 + 6 * 11] * x7 + v[9 + 7 * 11] * x8 + v[9 + 8 * 11] * x9 + v[9 + 9 * 11] * x10 + v[9 + 10 * 11] * x11;
2732: sum11 += v[10 + 0] * x1 + v[10 + 11] * x2 + v[10 + 2 * 11] * x3 + v[10 + 3 * 11] * x4 + v[10 + 4 * 11] * x5 + v[10 + 5 * 11] * x6 + v[10 + 6 * 11] * x7 + v[10 + 7 * 11] * x8 + v[10 + 8 * 11] * x9 + v[10 + 9 * 11] * x10 + v[10 + 10 * 11] * x11;
2733: v += 121;
2734: }
2735: z[0] = sum1;
2736: z[1] = sum2;
2737: z[2] = sum3;
2738: z[3] = sum4;
2739: z[4] = sum5;
2740: z[5] = sum6;
2741: z[6] = sum7;
2742: z[7] = sum8;
2743: z[8] = sum9;
2744: z[9] = sum10;
2745: z[10] = sum11;
2746: if (!usecprow) {
2747: z += 11;
2748: y += 11;
2749: }
2750: }
2751: PetscCall(VecRestoreArrayRead(xx, &x));
2752: PetscCall(VecRestoreArrayPair(yy, zz, &yarray, &zarray));
2753: PetscCall(PetscLogFlops(242.0 * a->nz));
2754: PetscFunctionReturn(PETSC_SUCCESS);
2755: }
2757: PetscErrorCode MatMultAdd_SeqBAIJ_N(Mat A, Vec xx, Vec yy, Vec zz)
2758: {
2759: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
2760: PetscScalar *z = NULL, *work, *workt, *zarray;
2761: const PetscScalar *x, *xb;
2762: const MatScalar *v;
2763: PetscInt mbs, i, bs = A->rmap->bs, j, n, bs2 = a->bs2;
2764: PetscInt ncols, k;
2765: const PetscInt *ridx = NULL, *idx, *ii;
2766: PetscBool usecprow = a->compressedrow.use;
2768: PetscFunctionBegin;
2769: PetscCall(VecCopy(yy, zz));
2770: PetscCall(VecGetArrayRead(xx, &x));
2771: PetscCall(VecGetArray(zz, &zarray));
2773: idx = a->j;
2774: v = a->a;
2775: if (usecprow) {
2776: mbs = a->compressedrow.nrows;
2777: ii = a->compressedrow.i;
2778: ridx = a->compressedrow.rindex;
2779: } else {
2780: mbs = a->mbs;
2781: ii = a->i;
2782: z = zarray;
2783: }
2785: if (!a->mult_work) {
2786: k = PetscMax(A->rmap->n, A->cmap->n);
2787: PetscCall(PetscMalloc1(k + 1, &a->mult_work));
2788: }
2789: work = a->mult_work;
2790: for (i = 0; i < mbs; i++) {
2791: n = ii[1] - ii[0];
2792: ii++;
2793: ncols = n * bs;
2794: workt = work;
2795: for (j = 0; j < n; j++) {
2796: xb = x + bs * (*idx++);
2797: for (k = 0; k < bs; k++) workt[k] = xb[k];
2798: workt += bs;
2799: }
2800: if (usecprow) z = zarray + bs * ridx[i];
2801: PetscKernel_w_gets_w_plus_Ar_times_v(bs, ncols, work, v, z);
2802: v += n * bs2;
2803: if (!usecprow) z += bs;
2804: }
2805: PetscCall(VecRestoreArrayRead(xx, &x));
2806: PetscCall(VecRestoreArray(zz, &zarray));
2807: PetscCall(PetscLogFlops(2.0 * a->nz * bs2));
2808: PetscFunctionReturn(PETSC_SUCCESS);
2809: }
2811: PetscErrorCode MatMultHermitianTranspose_SeqBAIJ(Mat A, Vec xx, Vec zz)
2812: {
2813: PetscScalar zero = 0.0;
2815: PetscFunctionBegin;
2816: PetscCall(VecSet(zz, zero));
2817: PetscCall(MatMultHermitianTransposeAdd_SeqBAIJ(A, xx, zz, zz));
2818: PetscFunctionReturn(PETSC_SUCCESS);
2819: }
2821: PetscErrorCode MatMultTranspose_SeqBAIJ(Mat A, Vec xx, Vec zz)
2822: {
2823: PetscScalar zero = 0.0;
2825: PetscFunctionBegin;
2826: PetscCall(VecSet(zz, zero));
2827: PetscCall(MatMultTransposeAdd_SeqBAIJ(A, xx, zz, zz));
2828: PetscFunctionReturn(PETSC_SUCCESS);
2829: }
2831: PetscErrorCode MatMultHermitianTransposeAdd_SeqBAIJ(Mat A, Vec xx, Vec yy, Vec zz)
2832: {
2833: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
2834: PetscScalar *z, x1, x2, x3, x4, x5;
2835: const PetscScalar *x, *xb = NULL;
2836: const MatScalar *v;
2837: PetscInt mbs, i, rval, bs = A->rmap->bs, j, n;
2838: const PetscInt *idx, *ii, *ib, *ridx = NULL;
2839: Mat_CompressedRow cprow = a->compressedrow;
2840: PetscBool usecprow = cprow.use;
2842: PetscFunctionBegin;
2843: if (yy != zz) PetscCall(VecCopy(yy, zz));
2844: PetscCall(VecGetArrayRead(xx, &x));
2845: PetscCall(VecGetArray(zz, &z));
2847: idx = a->j;
2848: v = a->a;
2849: if (usecprow) {
2850: mbs = cprow.nrows;
2851: ii = cprow.i;
2852: ridx = cprow.rindex;
2853: } else {
2854: mbs = a->mbs;
2855: ii = a->i;
2856: xb = x;
2857: }
2859: switch (bs) {
2860: case 1:
2861: for (i = 0; i < mbs; i++) {
2862: if (usecprow) xb = x + ridx[i];
2863: x1 = xb[0];
2864: ib = idx + ii[0];
2865: n = ii[1] - ii[0];
2866: ii++;
2867: for (j = 0; j < n; j++) {
2868: rval = ib[j];
2869: z[rval] += PetscConj(*v) * x1;
2870: v++;
2871: }
2872: if (!usecprow) xb++;
2873: }
2874: break;
2875: case 2:
2876: for (i = 0; i < mbs; i++) {
2877: if (usecprow) xb = x + 2 * ridx[i];
2878: x1 = xb[0];
2879: x2 = xb[1];
2880: ib = idx + ii[0];
2881: n = ii[1] - ii[0];
2882: ii++;
2883: for (j = 0; j < n; j++) {
2884: rval = ib[j] * 2;
2885: z[rval++] += PetscConj(v[0]) * x1 + PetscConj(v[1]) * x2;
2886: z[rval++] += PetscConj(v[2]) * x1 + PetscConj(v[3]) * x2;
2887: v += 4;
2888: }
2889: if (!usecprow) xb += 2;
2890: }
2891: break;
2892: case 3:
2893: for (i = 0; i < mbs; i++) {
2894: if (usecprow) xb = x + 3 * ridx[i];
2895: x1 = xb[0];
2896: x2 = xb[1];
2897: x3 = xb[2];
2898: ib = idx + ii[0];
2899: n = ii[1] - ii[0];
2900: ii++;
2901: for (j = 0; j < n; j++) {
2902: rval = ib[j] * 3;
2903: z[rval++] += PetscConj(v[0]) * x1 + PetscConj(v[1]) * x2 + PetscConj(v[2]) * x3;
2904: z[rval++] += PetscConj(v[3]) * x1 + PetscConj(v[4]) * x2 + PetscConj(v[5]) * x3;
2905: z[rval++] += PetscConj(v[6]) * x1 + PetscConj(v[7]) * x2 + PetscConj(v[8]) * x3;
2906: v += 9;
2907: }
2908: if (!usecprow) xb += 3;
2909: }
2910: break;
2911: case 4:
2912: for (i = 0; i < mbs; i++) {
2913: if (usecprow) xb = x + 4 * ridx[i];
2914: x1 = xb[0];
2915: x2 = xb[1];
2916: x3 = xb[2];
2917: x4 = xb[3];
2918: ib = idx + ii[0];
2919: n = ii[1] - ii[0];
2920: ii++;
2921: for (j = 0; j < n; j++) {
2922: rval = ib[j] * 4;
2923: z[rval++] += PetscConj(v[0]) * x1 + PetscConj(v[1]) * x2 + PetscConj(v[2]) * x3 + PetscConj(v[3]) * x4;
2924: z[rval++] += PetscConj(v[4]) * x1 + PetscConj(v[5]) * x2 + PetscConj(v[6]) * x3 + PetscConj(v[7]) * x4;
2925: z[rval++] += PetscConj(v[8]) * x1 + PetscConj(v[9]) * x2 + PetscConj(v[10]) * x3 + PetscConj(v[11]) * x4;
2926: z[rval++] += PetscConj(v[12]) * x1 + PetscConj(v[13]) * x2 + PetscConj(v[14]) * x3 + PetscConj(v[15]) * x4;
2927: v += 16;
2928: }
2929: if (!usecprow) xb += 4;
2930: }
2931: break;
2932: case 5:
2933: for (i = 0; i < mbs; i++) {
2934: if (usecprow) xb = x + 5 * ridx[i];
2935: x1 = xb[0];
2936: x2 = xb[1];
2937: x3 = xb[2];
2938: x4 = xb[3];
2939: x5 = xb[4];
2940: ib = idx + ii[0];
2941: n = ii[1] - ii[0];
2942: ii++;
2943: for (j = 0; j < n; j++) {
2944: rval = ib[j] * 5;
2945: z[rval++] += PetscConj(v[0]) * x1 + PetscConj(v[1]) * x2 + PetscConj(v[2]) * x3 + PetscConj(v[3]) * x4 + PetscConj(v[4]) * x5;
2946: z[rval++] += PetscConj(v[5]) * x1 + PetscConj(v[6]) * x2 + PetscConj(v[7]) * x3 + PetscConj(v[8]) * x4 + PetscConj(v[9]) * x5;
2947: z[rval++] += PetscConj(v[10]) * x1 + PetscConj(v[11]) * x2 + PetscConj(v[12]) * x3 + PetscConj(v[13]) * x4 + PetscConj(v[14]) * x5;
2948: z[rval++] += PetscConj(v[15]) * x1 + PetscConj(v[16]) * x2 + PetscConj(v[17]) * x3 + PetscConj(v[18]) * x4 + PetscConj(v[19]) * x5;
2949: z[rval++] += PetscConj(v[20]) * x1 + PetscConj(v[21]) * x2 + PetscConj(v[22]) * x3 + PetscConj(v[23]) * x4 + PetscConj(v[24]) * x5;
2950: v += 25;
2951: }
2952: if (!usecprow) xb += 5;
2953: }
2954: break;
2955: default: /* block sizes larger than 5 by 5 are handled by BLAS */
2956: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "block size larger than 5 is not supported yet");
2957: #if 0
2958: {
2959: PetscInt ncols,k,bs2=a->bs2;
2960: PetscScalar *work,*workt,zb;
2961: const PetscScalar *xtmp;
2962: if (!a->mult_work) {
2963: k = PetscMax(A->rmap->n,A->cmap->n);
2964: PetscCall(PetscMalloc1(k+1,&a->mult_work));
2965: }
2966: work = a->mult_work;
2967: xtmp = x;
2968: for (i=0; i<mbs; i++) {
2969: n = ii[1] - ii[0]; ii++;
2970: ncols = n*bs;
2971: PetscCall(PetscArrayzero(work,ncols));
2972: if (usecprow) xtmp = x + bs*ridx[i];
2973: PetscKernel_w_gets_w_plus_trans_Ar_times_v(bs,ncols,xtmp,v,work);
2974: v += n*bs2;
2975: if (!usecprow) xtmp += bs;
2976: workt = work;
2977: for (j=0; j<n; j++) {
2978: zb = z + bs*(*idx++);
2979: for (k=0; k<bs; k++) zb[k] += workt[k] ;
2980: workt += bs;
2981: }
2982: }
2983: }
2984: #endif
2985: }
2986: PetscCall(VecRestoreArrayRead(xx, &x));
2987: PetscCall(VecRestoreArray(zz, &z));
2988: PetscCall(PetscLogFlops(2.0 * a->nz * a->bs2));
2989: PetscFunctionReturn(PETSC_SUCCESS);
2990: }
2992: PetscErrorCode MatMultTransposeAdd_SeqBAIJ(Mat A, Vec xx, Vec yy, Vec zz)
2993: {
2994: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
2995: PetscScalar *zb, *z, x1, x2, x3, x4, x5;
2996: const PetscScalar *x, *xb = NULL;
2997: const MatScalar *v;
2998: PetscInt mbs, i, rval, bs = A->rmap->bs, j, n, bs2 = a->bs2;
2999: const PetscInt *idx, *ii, *ib, *ridx = NULL;
3000: Mat_CompressedRow cprow = a->compressedrow;
3001: PetscBool usecprow = cprow.use;
3003: PetscFunctionBegin;
3004: if (yy != zz) PetscCall(VecCopy(yy, zz));
3005: PetscCall(VecGetArrayRead(xx, &x));
3006: PetscCall(VecGetArray(zz, &z));
3008: idx = a->j;
3009: v = a->a;
3010: if (usecprow) {
3011: mbs = cprow.nrows;
3012: ii = cprow.i;
3013: ridx = cprow.rindex;
3014: } else {
3015: mbs = a->mbs;
3016: ii = a->i;
3017: xb = x;
3018: }
3020: switch (bs) {
3021: case 1:
3022: for (i = 0; i < mbs; i++) {
3023: if (usecprow) xb = x + ridx[i];
3024: x1 = xb[0];
3025: ib = idx + ii[0];
3026: n = ii[1] - ii[0];
3027: ii++;
3028: for (j = 0; j < n; j++) {
3029: rval = ib[j];
3030: z[rval] += *v * x1;
3031: v++;
3032: }
3033: if (!usecprow) xb++;
3034: }
3035: break;
3036: case 2:
3037: for (i = 0; i < mbs; i++) {
3038: if (usecprow) xb = x + 2 * ridx[i];
3039: x1 = xb[0];
3040: x2 = xb[1];
3041: ib = idx + ii[0];
3042: n = ii[1] - ii[0];
3043: ii++;
3044: for (j = 0; j < n; j++) {
3045: rval = ib[j] * 2;
3046: z[rval++] += v[0] * x1 + v[1] * x2;
3047: z[rval++] += v[2] * x1 + v[3] * x2;
3048: v += 4;
3049: }
3050: if (!usecprow) xb += 2;
3051: }
3052: break;
3053: case 3:
3054: for (i = 0; i < mbs; i++) {
3055: if (usecprow) xb = x + 3 * ridx[i];
3056: x1 = xb[0];
3057: x2 = xb[1];
3058: x3 = xb[2];
3059: ib = idx + ii[0];
3060: n = ii[1] - ii[0];
3061: ii++;
3062: for (j = 0; j < n; j++) {
3063: rval = ib[j] * 3;
3064: z[rval++] += v[0] * x1 + v[1] * x2 + v[2] * x3;
3065: z[rval++] += v[3] * x1 + v[4] * x2 + v[5] * x3;
3066: z[rval++] += v[6] * x1 + v[7] * x2 + v[8] * x3;
3067: v += 9;
3068: }
3069: if (!usecprow) xb += 3;
3070: }
3071: break;
3072: case 4:
3073: for (i = 0; i < mbs; i++) {
3074: if (usecprow) xb = x + 4 * ridx[i];
3075: x1 = xb[0];
3076: x2 = xb[1];
3077: x3 = xb[2];
3078: x4 = xb[3];
3079: ib = idx + ii[0];
3080: n = ii[1] - ii[0];
3081: ii++;
3082: for (j = 0; j < n; j++) {
3083: rval = ib[j] * 4;
3084: z[rval++] += v[0] * x1 + v[1] * x2 + v[2] * x3 + v[3] * x4;
3085: z[rval++] += v[4] * x1 + v[5] * x2 + v[6] * x3 + v[7] * x4;
3086: z[rval++] += v[8] * x1 + v[9] * x2 + v[10] * x3 + v[11] * x4;
3087: z[rval++] += v[12] * x1 + v[13] * x2 + v[14] * x3 + v[15] * x4;
3088: v += 16;
3089: }
3090: if (!usecprow) xb += 4;
3091: }
3092: break;
3093: case 5:
3094: for (i = 0; i < mbs; i++) {
3095: if (usecprow) xb = x + 5 * ridx[i];
3096: x1 = xb[0];
3097: x2 = xb[1];
3098: x3 = xb[2];
3099: x4 = xb[3];
3100: x5 = xb[4];
3101: ib = idx + ii[0];
3102: n = ii[1] - ii[0];
3103: ii++;
3104: for (j = 0; j < n; j++) {
3105: rval = ib[j] * 5;
3106: z[rval++] += v[0] * x1 + v[1] * x2 + v[2] * x3 + v[3] * x4 + v[4] * x5;
3107: z[rval++] += v[5] * x1 + v[6] * x2 + v[7] * x3 + v[8] * x4 + v[9] * x5;
3108: z[rval++] += v[10] * x1 + v[11] * x2 + v[12] * x3 + v[13] * x4 + v[14] * x5;
3109: z[rval++] += v[15] * x1 + v[16] * x2 + v[17] * x3 + v[18] * x4 + v[19] * x5;
3110: z[rval++] += v[20] * x1 + v[21] * x2 + v[22] * x3 + v[23] * x4 + v[24] * x5;
3111: v += 25;
3112: }
3113: if (!usecprow) xb += 5;
3114: }
3115: break;
3116: default: { /* block sizes larger than 5 by 5 are handled by BLAS */
3117: PetscInt ncols, k;
3118: PetscScalar *work, *workt;
3119: const PetscScalar *xtmp;
3120: if (!a->mult_work) {
3121: k = PetscMax(A->rmap->n, A->cmap->n);
3122: PetscCall(PetscMalloc1(k + 1, &a->mult_work));
3123: }
3124: work = a->mult_work;
3125: xtmp = x;
3126: for (i = 0; i < mbs; i++) {
3127: n = ii[1] - ii[0];
3128: ii++;
3129: ncols = n * bs;
3130: PetscCall(PetscArrayzero(work, ncols));
3131: if (usecprow) xtmp = x + bs * ridx[i];
3132: PetscKernel_w_gets_w_plus_trans_Ar_times_v(bs, ncols, xtmp, v, work);
3133: v += n * bs2;
3134: if (!usecprow) xtmp += bs;
3135: workt = work;
3136: for (j = 0; j < n; j++) {
3137: zb = z + bs * (*idx++);
3138: for (k = 0; k < bs; k++) zb[k] += workt[k];
3139: workt += bs;
3140: }
3141: }
3142: }
3143: }
3144: PetscCall(VecRestoreArrayRead(xx, &x));
3145: PetscCall(VecRestoreArray(zz, &z));
3146: PetscCall(PetscLogFlops(2.0 * a->nz * a->bs2));
3147: PetscFunctionReturn(PETSC_SUCCESS);
3148: }
3150: PetscErrorCode MatScale_SeqBAIJ(Mat inA, PetscScalar alpha)
3151: {
3152: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)inA->data;
3153: PetscInt totalnz = a->bs2 * a->nz;
3154: PetscScalar oalpha = alpha;
3155: PetscBLASInt one = 1, tnz;
3157: PetscFunctionBegin;
3158: PetscCall(PetscBLASIntCast(totalnz, &tnz));
3159: PetscCallBLAS("BLASscal", BLASscal_(&tnz, &oalpha, a->a, &one));
3160: PetscCall(PetscLogFlops(totalnz));
3161: PetscFunctionReturn(PETSC_SUCCESS);
3162: }
3164: PetscErrorCode MatNorm_SeqBAIJ(Mat A, NormType type, PetscReal *norm)
3165: {
3166: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
3167: MatScalar *v = a->a;
3168: PetscReal sum = 0.0;
3169: PetscInt i, j, k, bs = A->rmap->bs, nz = a->nz, bs2 = a->bs2, k1;
3171: PetscFunctionBegin;
3172: if (type == NORM_FROBENIUS) {
3173: #if PetscDefined(USE_REAL___FP16)
3174: PetscBLASInt one = 1, cnt = bs2 * nz;
3175: PetscCallBLAS("BLASnrm2", *norm = BLASnrm2_(&cnt, v, &one));
3176: #else
3177: for (i = 0; i < bs2 * nz; i++) {
3178: sum += PetscRealPart(PetscConj(*v) * (*v));
3179: v++;
3180: }
3181: #endif
3182: *norm = PetscSqrtReal(sum);
3183: PetscCall(PetscLogFlops(2.0 * bs2 * nz));
3184: } else if (type == NORM_1) { /* maximum column sum */
3185: PetscReal *tmp;
3186: PetscInt *bcol = a->j;
3187: PetscCall(PetscCalloc1(A->cmap->n + 1, &tmp));
3188: for (i = 0; i < nz; i++) {
3189: for (j = 0; j < bs; j++) {
3190: k1 = bs * (*bcol) + j; /* column index */
3191: for (k = 0; k < bs; k++) {
3192: tmp[k1] += PetscAbsScalar(*v);
3193: v++;
3194: }
3195: }
3196: bcol++;
3197: }
3198: *norm = 0.0;
3199: for (j = 0; j < A->cmap->n; j++) {
3200: if (tmp[j] > *norm) *norm = tmp[j];
3201: }
3202: PetscCall(PetscFree(tmp));
3203: PetscCall(PetscLogFlops(PetscMax(bs2 * nz - 1, 0)));
3204: } else if (type == NORM_INFINITY) { /* maximum row sum */
3205: *norm = 0.0;
3206: for (k = 0; k < bs; k++) {
3207: for (j = 0; j < a->mbs; j++) {
3208: v = a->a + bs2 * a->i[j] + k;
3209: sum = 0.0;
3210: for (i = 0; i < a->i[j + 1] - a->i[j]; i++) {
3211: for (k1 = 0; k1 < bs; k1++) {
3212: sum += PetscAbsScalar(*v);
3213: v += bs;
3214: }
3215: }
3216: if (sum > *norm) *norm = sum;
3217: }
3218: }
3219: PetscCall(PetscLogFlops(PetscMax(bs2 * nz - 1, 0)));
3220: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for this norm yet");
3221: PetscFunctionReturn(PETSC_SUCCESS);
3222: }
3224: PetscErrorCode MatGetDiagonal_SeqBAIJ(Mat A, Vec v)
3225: {
3226: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
3227: PetscInt n;
3228: const PetscInt bs = A->rmap->bs, ambs = a->mbs, bs2 = a->bs2;
3229: PetscScalar *x;
3230: const MatScalar *aa = a->a, *aa_j;
3231: const PetscInt *ai = a->i, *adiag;
3232: PetscBool diagDense;
3234: PetscFunctionBegin;
3235: PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
3237: PetscCall(VecGetLocalSize(v, &n));
3238: PetscCheck(n == A->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector");
3239: PetscCall(VecGetArrayWrite(v, &x));
3240: PetscCall(MatGetDiagonalMarkers_SeqBAIJ(A, &adiag, &diagDense));
3241: if (diagDense) {
3242: for (PetscInt i = 0, row = 0; i < ambs; i++) {
3243: aa_j = aa + adiag[i] * bs2;
3244: for (PetscInt k = 0; k < bs2; k += (bs + 1)) x[row++] = aa_j[k];
3245: }
3246: } else {
3247: for (PetscInt i = 0, row = 0; i < ambs; i++) {
3248: const PetscInt j = adiag[i];
3250: if (j != ai[i + 1]) {
3251: aa_j = aa + j * bs2;
3252: for (PetscInt k = 0; k < bs2; k += (bs + 1)) x[row++] = aa_j[k];
3253: } else {
3254: for (PetscInt k = 0; k < bs; k++) x[row++] = 0.0;
3255: }
3256: }
3257: }
3258: PetscCall(VecRestoreArrayWrite(v, &x));
3259: PetscFunctionReturn(PETSC_SUCCESS);
3260: }
3262: PetscErrorCode MatDiagonalScale_SeqBAIJ(Mat A, Vec ll, Vec rr)
3263: {
3264: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
3265: const PetscScalar *l, *r, *li, *ri;
3266: PetscScalar x;
3267: MatScalar *aa, *v;
3268: PetscInt i, j, k, lm, rn, M, m, n, mbs, tmp, bs, bs2, iai;
3269: const PetscInt *ai, *aj;
3271: PetscFunctionBegin;
3272: ai = a->i;
3273: aj = a->j;
3274: aa = a->a;
3275: m = A->rmap->n;
3276: n = A->cmap->n;
3277: bs = A->rmap->bs;
3278: mbs = a->mbs;
3279: bs2 = a->bs2;
3280: if (ll) {
3281: PetscCall(VecGetArrayRead(ll, &l));
3282: PetscCall(VecGetLocalSize(ll, &lm));
3283: PetscCheck(lm == m, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Left scaling vector wrong length");
3284: for (i = 0; i < mbs; i++) { /* for each block row */
3285: M = ai[i + 1] - ai[i];
3286: li = l + i * bs;
3287: v = PetscSafePointerPlusOffset(aa, bs2 * ai[i]);
3288: for (j = 0; j < M; j++) { /* for each block */
3289: for (k = 0; k < bs2; k++) (*v++) *= li[k % bs];
3290: }
3291: }
3292: PetscCall(VecRestoreArrayRead(ll, &l));
3293: PetscCall(PetscLogFlops(a->nz));
3294: }
3296: if (rr) {
3297: PetscCall(VecGetArrayRead(rr, &r));
3298: PetscCall(VecGetLocalSize(rr, &rn));
3299: PetscCheck(rn == n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Right scaling vector wrong length");
3300: for (i = 0; i < mbs; i++) { /* for each block row */
3301: iai = ai[i];
3302: M = ai[i + 1] - iai;
3303: v = PetscSafePointerPlusOffset(aa, bs2 * iai);
3304: for (j = 0; j < M; j++) { /* for each block */
3305: ri = r + bs * aj[iai + j];
3306: for (k = 0; k < bs; k++) {
3307: x = ri[k];
3308: for (tmp = 0; tmp < bs; tmp++) v[tmp] *= x;
3309: v += bs;
3310: }
3311: }
3312: }
3313: PetscCall(VecRestoreArrayRead(rr, &r));
3314: PetscCall(PetscLogFlops(a->nz));
3315: }
3316: PetscFunctionReturn(PETSC_SUCCESS);
3317: }
3319: PetscErrorCode MatGetInfo_SeqBAIJ(Mat A, MatInfoType flag, MatInfo *info)
3320: {
3321: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
3323: PetscFunctionBegin;
3324: info->block_size = a->bs2;
3325: info->nz_allocated = a->bs2 * a->maxnz;
3326: info->nz_used = a->bs2 * a->nz;
3327: info->nz_unneeded = info->nz_allocated - info->nz_used;
3328: info->assemblies = A->num_ass;
3329: info->mallocs = A->info.mallocs;
3330: info->memory = 0; /* REVIEW ME */
3331: if (A->factortype) {
3332: info->fill_ratio_given = A->info.fill_ratio_given;
3333: info->fill_ratio_needed = A->info.fill_ratio_needed;
3334: info->factor_mallocs = A->info.factor_mallocs;
3335: } else {
3336: info->fill_ratio_given = 0;
3337: info->fill_ratio_needed = 0;
3338: info->factor_mallocs = 0;
3339: }
3340: PetscFunctionReturn(PETSC_SUCCESS);
3341: }
3343: PetscErrorCode MatZeroEntries_SeqBAIJ(Mat A)
3344: {
3345: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
3347: PetscFunctionBegin;
3348: PetscCall(PetscArrayzero(a->a, a->bs2 * a->i[a->mbs]));
3349: PetscFunctionReturn(PETSC_SUCCESS);
3350: }
3352: PetscErrorCode MatMatMultSymbolic_SeqBAIJ_SeqDense(Mat A, Mat B, PetscReal fill, Mat C)
3353: {
3354: PetscFunctionBegin;
3355: PetscCall(MatMatMultSymbolic_SeqDense_SeqDense(A, B, 0.0, C));
3356: C->ops->matmultnumeric = MatMatMultNumeric_SeqBAIJ_SeqDense;
3357: PetscFunctionReturn(PETSC_SUCCESS);
3358: }
3360: PetscErrorCode MatTransposeMatMultSymbolic_SeqBAIJ_SeqDense(Mat A, Mat B, PetscReal fill, Mat C)
3361: {
3362: PetscFunctionBegin;
3363: PetscCall(MatTransposeMatMultSymbolic_SeqDense_SeqDense(A, B, 0.0, C));
3364: C->ops->transposematmultnumeric = MatTransposeMatMultNumeric_SeqBAIJ_SeqDense;
3365: PetscFunctionReturn(PETSC_SUCCESS);
3366: }
3368: static PetscErrorCode MatMatMult_SeqBAIJ_1_Private(Mat A, PetscScalar *b, PetscInt bm, PetscScalar *c, PetscInt cm, PetscInt cn)
3369: {
3370: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
3371: PetscScalar *z = NULL, sum1;
3372: const PetscScalar *xb;
3373: PetscScalar x1;
3374: const MatScalar *v, *vv;
3375: PetscInt mbs, i, *idx, *ii, j, *jj, n, k, *ridx = NULL;
3376: PetscBool usecprow = a->compressedrow.use;
3378: PetscFunctionBegin;
3379: idx = a->j;
3380: v = a->a;
3381: if (usecprow) {
3382: mbs = a->compressedrow.nrows;
3383: ii = a->compressedrow.i;
3384: ridx = a->compressedrow.rindex;
3385: } else {
3386: mbs = a->mbs;
3387: ii = a->i;
3388: z = c;
3389: }
3391: for (i = 0; i < mbs; i++) {
3392: n = ii[1] - ii[0];
3393: ii++;
3394: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
3395: PetscPrefetchBlock(v + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
3396: if (usecprow) z = c + ridx[i];
3397: jj = idx;
3398: vv = v;
3399: for (k = 0; k < cn; k++) {
3400: idx = jj;
3401: v = vv;
3402: sum1 = 0.0;
3403: for (j = 0; j < n; j++) {
3404: xb = b + (*idx++);
3405: x1 = xb[0 + k * bm];
3406: sum1 += v[0] * x1;
3407: v += 1;
3408: }
3409: z[0 + k * cm] = sum1;
3410: }
3411: if (!usecprow) z += 1;
3412: }
3413: PetscFunctionReturn(PETSC_SUCCESS);
3414: }
3416: static PetscErrorCode MatMatMult_SeqBAIJ_2_Private(Mat A, PetscScalar *b, PetscInt bm, PetscScalar *c, PetscInt cm, PetscInt cn)
3417: {
3418: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
3419: PetscScalar *z = NULL, sum1, sum2;
3420: const PetscScalar *xb;
3421: PetscScalar x1, x2;
3422: const MatScalar *v, *vv;
3423: PetscInt mbs, i, *idx, *ii, j, *jj, n, k, *ridx = NULL;
3424: PetscBool usecprow = a->compressedrow.use;
3426: PetscFunctionBegin;
3427: idx = a->j;
3428: v = a->a;
3429: if (usecprow) {
3430: mbs = a->compressedrow.nrows;
3431: ii = a->compressedrow.i;
3432: ridx = a->compressedrow.rindex;
3433: } else {
3434: mbs = a->mbs;
3435: ii = a->i;
3436: z = c;
3437: }
3439: for (i = 0; i < mbs; i++) {
3440: n = ii[1] - ii[0];
3441: ii++;
3442: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
3443: PetscPrefetchBlock(v + 4 * n, 4 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
3444: if (usecprow) z = c + 2 * ridx[i];
3445: jj = idx;
3446: vv = v;
3447: for (k = 0; k < cn; k++) {
3448: idx = jj;
3449: v = vv;
3450: sum1 = 0.0;
3451: sum2 = 0.0;
3452: for (j = 0; j < n; j++) {
3453: xb = b + 2 * (*idx++);
3454: x1 = xb[0 + k * bm];
3455: x2 = xb[1 + k * bm];
3456: sum1 += v[0] * x1 + v[2] * x2;
3457: sum2 += v[1] * x1 + v[3] * x2;
3458: v += 4;
3459: }
3460: z[0 + k * cm] = sum1;
3461: z[1 + k * cm] = sum2;
3462: }
3463: if (!usecprow) z += 2;
3464: }
3465: PetscFunctionReturn(PETSC_SUCCESS);
3466: }
3468: static PetscErrorCode MatMatMult_SeqBAIJ_3_Private(Mat A, PetscScalar *b, PetscInt bm, PetscScalar *c, PetscInt cm, PetscInt cn)
3469: {
3470: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
3471: PetscScalar *z = NULL, sum1, sum2, sum3;
3472: const PetscScalar *xb;
3473: PetscScalar x1, x2, x3;
3474: const MatScalar *v, *vv;
3475: PetscInt mbs, i, *idx, *ii, j, *jj, n, k, *ridx = NULL;
3476: PetscBool usecprow = a->compressedrow.use;
3478: PetscFunctionBegin;
3479: idx = a->j;
3480: v = a->a;
3481: if (usecprow) {
3482: mbs = a->compressedrow.nrows;
3483: ii = a->compressedrow.i;
3484: ridx = a->compressedrow.rindex;
3485: } else {
3486: mbs = a->mbs;
3487: ii = a->i;
3488: z = c;
3489: }
3491: for (i = 0; i < mbs; i++) {
3492: n = ii[1] - ii[0];
3493: ii++;
3494: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
3495: PetscPrefetchBlock(v + 9 * n, 9 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
3496: if (usecprow) z = c + 3 * ridx[i];
3497: jj = idx;
3498: vv = v;
3499: for (k = 0; k < cn; k++) {
3500: idx = jj;
3501: v = vv;
3502: sum1 = 0.0;
3503: sum2 = 0.0;
3504: sum3 = 0.0;
3505: for (j = 0; j < n; j++) {
3506: xb = b + 3 * (*idx++);
3507: x1 = xb[0 + k * bm];
3508: x2 = xb[1 + k * bm];
3509: x3 = xb[2 + k * bm];
3510: sum1 += v[0] * x1 + v[3] * x2 + v[6] * x3;
3511: sum2 += v[1] * x1 + v[4] * x2 + v[7] * x3;
3512: sum3 += v[2] * x1 + v[5] * x2 + v[8] * x3;
3513: v += 9;
3514: }
3515: z[0 + k * cm] = sum1;
3516: z[1 + k * cm] = sum2;
3517: z[2 + k * cm] = sum3;
3518: }
3519: if (!usecprow) z += 3;
3520: }
3521: PetscFunctionReturn(PETSC_SUCCESS);
3522: }
3524: static PetscErrorCode MatMatMult_SeqBAIJ_4_Private(Mat A, PetscScalar *b, PetscInt bm, PetscScalar *c, PetscInt cm, PetscInt cn)
3525: {
3526: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
3527: PetscScalar *z = NULL, sum1, sum2, sum3, sum4;
3528: const PetscScalar *xb;
3529: PetscScalar x1, x2, x3, x4;
3530: const MatScalar *v, *vv;
3531: PetscInt mbs, i, *idx, *ii, j, *jj, n, k, *ridx = NULL;
3532: PetscBool usecprow = a->compressedrow.use;
3534: PetscFunctionBegin;
3535: idx = a->j;
3536: v = a->a;
3537: if (usecprow) {
3538: mbs = a->compressedrow.nrows;
3539: ii = a->compressedrow.i;
3540: ridx = a->compressedrow.rindex;
3541: } else {
3542: mbs = a->mbs;
3543: ii = a->i;
3544: z = c;
3545: }
3547: for (i = 0; i < mbs; i++) {
3548: n = ii[1] - ii[0];
3549: ii++;
3550: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
3551: PetscPrefetchBlock(v + 16 * n, 16 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
3552: if (usecprow) z = c + 4 * ridx[i];
3553: jj = idx;
3554: vv = v;
3555: for (k = 0; k < cn; k++) {
3556: idx = jj;
3557: v = vv;
3558: sum1 = 0.0;
3559: sum2 = 0.0;
3560: sum3 = 0.0;
3561: sum4 = 0.0;
3562: for (j = 0; j < n; j++) {
3563: xb = b + 4 * (*idx++);
3564: x1 = xb[0 + k * bm];
3565: x2 = xb[1 + k * bm];
3566: x3 = xb[2 + k * bm];
3567: x4 = xb[3 + k * bm];
3568: sum1 += v[0] * x1 + v[4] * x2 + v[8] * x3 + v[12] * x4;
3569: sum2 += v[1] * x1 + v[5] * x2 + v[9] * x3 + v[13] * x4;
3570: sum3 += v[2] * x1 + v[6] * x2 + v[10] * x3 + v[14] * x4;
3571: sum4 += v[3] * x1 + v[7] * x2 + v[11] * x3 + v[15] * x4;
3572: v += 16;
3573: }
3574: z[0 + k * cm] = sum1;
3575: z[1 + k * cm] = sum2;
3576: z[2 + k * cm] = sum3;
3577: z[3 + k * cm] = sum4;
3578: }
3579: if (!usecprow) z += 4;
3580: }
3581: PetscFunctionReturn(PETSC_SUCCESS);
3582: }
3584: static PetscErrorCode MatMatMult_SeqBAIJ_5_Private(Mat A, PetscScalar *b, PetscInt bm, PetscScalar *c, PetscInt cm, PetscInt cn)
3585: {
3586: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
3587: PetscScalar *z = NULL, sum1, sum2, sum3, sum4, sum5;
3588: const PetscScalar *xb;
3589: PetscScalar x1, x2, x3, x4, x5;
3590: const MatScalar *v, *vv;
3591: PetscInt mbs, i, *idx, *ii, j, *jj, n, k, *ridx = NULL;
3592: PetscBool usecprow = a->compressedrow.use;
3594: PetscFunctionBegin;
3595: idx = a->j;
3596: v = a->a;
3597: if (usecprow) {
3598: mbs = a->compressedrow.nrows;
3599: ii = a->compressedrow.i;
3600: ridx = a->compressedrow.rindex;
3601: } else {
3602: mbs = a->mbs;
3603: ii = a->i;
3604: z = c;
3605: }
3607: for (i = 0; i < mbs; i++) {
3608: n = ii[1] - ii[0];
3609: ii++;
3610: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
3611: PetscPrefetchBlock(v + 25 * n, 25 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
3612: if (usecprow) z = c + 5 * ridx[i];
3613: jj = idx;
3614: vv = v;
3615: for (k = 0; k < cn; k++) {
3616: idx = jj;
3617: v = vv;
3618: sum1 = 0.0;
3619: sum2 = 0.0;
3620: sum3 = 0.0;
3621: sum4 = 0.0;
3622: sum5 = 0.0;
3623: for (j = 0; j < n; j++) {
3624: xb = b + 5 * (*idx++);
3625: x1 = xb[0 + k * bm];
3626: x2 = xb[1 + k * bm];
3627: x3 = xb[2 + k * bm];
3628: x4 = xb[3 + k * bm];
3629: x5 = xb[4 + k * bm];
3630: sum1 += v[0] * x1 + v[5] * x2 + v[10] * x3 + v[15] * x4 + v[20] * x5;
3631: sum2 += v[1] * x1 + v[6] * x2 + v[11] * x3 + v[16] * x4 + v[21] * x5;
3632: sum3 += v[2] * x1 + v[7] * x2 + v[12] * x3 + v[17] * x4 + v[22] * x5;
3633: sum4 += v[3] * x1 + v[8] * x2 + v[13] * x3 + v[18] * x4 + v[23] * x5;
3634: sum5 += v[4] * x1 + v[9] * x2 + v[14] * x3 + v[19] * x4 + v[24] * x5;
3635: v += 25;
3636: }
3637: z[0 + k * cm] = sum1;
3638: z[1 + k * cm] = sum2;
3639: z[2 + k * cm] = sum3;
3640: z[3 + k * cm] = sum4;
3641: z[4 + k * cm] = sum5;
3642: }
3643: if (!usecprow) z += 5;
3644: }
3645: PetscFunctionReturn(PETSC_SUCCESS);
3646: }
3648: PetscErrorCode MatMatMultNumeric_SeqBAIJ_SeqDense(Mat A, Mat B, Mat C)
3649: {
3650: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
3651: Mat_SeqDense *bd = (Mat_SeqDense *)B->data;
3652: Mat_SeqDense *cd = (Mat_SeqDense *)C->data;
3653: PetscInt cm = cd->lda, cn = B->cmap->n, bm = bd->lda;
3654: PetscInt mbs, i, bs = A->rmap->bs, j, n, bs2 = a->bs2;
3655: PetscBLASInt bbs, bcn, bbm, bcm;
3656: PetscScalar *z = NULL;
3657: PetscScalar *c, *b;
3658: const MatScalar *v;
3659: const PetscInt *idx, *ii, *ridx = NULL;
3660: PetscScalar _DZero = 0.0, _DOne = 1.0;
3661: PetscBool usecprow = a->compressedrow.use;
3663: PetscFunctionBegin;
3664: if (!cm || !cn) PetscFunctionReturn(PETSC_SUCCESS);
3665: 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);
3666: 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);
3667: 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);
3668: b = bd->v;
3669: if (a->nonzerorowcnt != A->rmap->n) PetscCall(MatZeroEntries(C));
3670: PetscCall(MatDenseGetArrayWrite(C, &c));
3671: switch (bs) {
3672: case 1:
3673: PetscCall(MatMatMult_SeqBAIJ_1_Private(A, b, bm, c, cm, cn));
3674: break;
3675: case 2:
3676: PetscCall(MatMatMult_SeqBAIJ_2_Private(A, b, bm, c, cm, cn));
3677: break;
3678: case 3:
3679: PetscCall(MatMatMult_SeqBAIJ_3_Private(A, b, bm, c, cm, cn));
3680: break;
3681: case 4:
3682: PetscCall(MatMatMult_SeqBAIJ_4_Private(A, b, bm, c, cm, cn));
3683: break;
3684: case 5:
3685: PetscCall(MatMatMult_SeqBAIJ_5_Private(A, b, bm, c, cm, cn));
3686: break;
3687: default: /* block sizes larger than 5 by 5 are handled by BLAS */
3688: PetscCall(PetscBLASIntCast(bs, &bbs));
3689: PetscCall(PetscBLASIntCast(cn, &bcn));
3690: PetscCall(PetscBLASIntCast(bm, &bbm));
3691: PetscCall(PetscBLASIntCast(cm, &bcm));
3692: idx = a->j;
3693: v = a->a;
3694: if (usecprow) {
3695: mbs = a->compressedrow.nrows;
3696: ii = a->compressedrow.i;
3697: ridx = a->compressedrow.rindex;
3698: } else {
3699: mbs = a->mbs;
3700: ii = a->i;
3701: z = c;
3702: }
3703: for (i = 0; i < mbs; i++) {
3704: n = ii[1] - ii[0];
3705: ii++;
3706: if (usecprow) z = c + bs * ridx[i];
3707: if (n) {
3708: PetscCallBLAS("BLASgemm", BLASgemm_("N", "N", &bbs, &bcn, &bbs, &_DOne, v, &bbs, b + bs * (*idx++), &bbm, &_DZero, z, &bcm));
3709: v += bs2;
3710: }
3711: for (j = 1; j < n; j++) {
3712: PetscCallBLAS("BLASgemm", BLASgemm_("N", "N", &bbs, &bcn, &bbs, &_DOne, v, &bbs, b + bs * (*idx++), &bbm, &_DOne, z, &bcm));
3713: v += bs2;
3714: }
3715: if (!usecprow) z += bs;
3716: }
3717: }
3718: PetscCall(MatDenseRestoreArrayWrite(C, &c));
3719: PetscCall(PetscLogFlops((2.0 * a->nz * bs2 - bs * a->nonzerorowcnt) * cn));
3720: PetscFunctionReturn(PETSC_SUCCESS);
3721: }
3723: static PetscErrorCode MatTransposeMatMult_SeqBAIJ_1_Private(Mat A, PetscScalar *b, PetscInt bm, PetscScalar *c, PetscInt cm, PetscInt cn)
3724: {
3725: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
3726: const MatScalar *v;
3727: PetscInt mbs, i, *idx, *ii, j, n, k, *ridx = NULL;
3728: PetscBool usecprow = a->compressedrow.use;
3729: const PetscScalar *bi;
3731: PetscFunctionBegin;
3732: idx = a->j;
3733: v = a->a;
3734: if (usecprow) {
3735: mbs = a->compressedrow.nrows;
3736: ii = a->compressedrow.i;
3737: ridx = a->compressedrow.rindex;
3738: } else {
3739: mbs = a->mbs;
3740: ii = a->i;
3741: }
3743: for (i = 0; i < mbs; i++) {
3744: PetscInt brow = usecprow ? ridx[i] : i;
3746: n = ii[1] - ii[0];
3747: ii++;
3748: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
3749: PetscPrefetchBlock(v + 1 * n, 1 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
3750: for (j = 0, bi = b + 1 * brow; j < n; j++) {
3751: PetscScalar *zcol = c + 1 * (*idx++);
3753: for (k = 0; k < cn; k++) zcol[0 + k * cm] += v[0] * bi[k * bm];
3754: ++v;
3755: }
3756: }
3757: PetscFunctionReturn(PETSC_SUCCESS);
3758: }
3760: static PetscErrorCode MatTransposeMatMult_SeqBAIJ_2_Private(Mat A, PetscScalar *b, PetscInt bm, PetscScalar *c, PetscInt cm, PetscInt cn)
3761: {
3762: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
3763: const MatScalar *v;
3764: PetscInt mbs, i, *idx, *ii, j, n, k, *ridx = NULL;
3765: PetscBool usecprow = a->compressedrow.use;
3766: const PetscScalar *bi;
3767: PetscScalar x1, x2;
3769: PetscFunctionBegin;
3770: idx = a->j;
3771: v = a->a;
3772: if (usecprow) {
3773: mbs = a->compressedrow.nrows;
3774: ii = a->compressedrow.i;
3775: ridx = a->compressedrow.rindex;
3776: } else {
3777: mbs = a->mbs;
3778: ii = a->i;
3779: }
3781: for (i = 0; i < mbs; i++) {
3782: PetscInt brow = usecprow ? ridx[i] : i;
3784: n = ii[1] - ii[0];
3785: ii++;
3786: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
3787: PetscPrefetchBlock(v + 4 * n, 4 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
3788: for (j = 0, bi = b + 2 * brow; j < n; j++) {
3789: PetscScalar *zcol = c + 2 * (*idx++);
3791: for (k = 0; k < cn; k++) {
3792: x1 = bi[0 + k * bm];
3793: x2 = bi[1 + k * bm];
3794: zcol[0 + k * cm] += v[0] * x1 + v[1] * x2;
3795: zcol[1 + k * cm] += v[2] * x1 + v[3] * x2;
3796: }
3797: v += 4;
3798: }
3799: }
3800: PetscFunctionReturn(PETSC_SUCCESS);
3801: }
3803: static PetscErrorCode MatTransposeMatMult_SeqBAIJ_3_Private(Mat A, PetscScalar *b, PetscInt bm, PetscScalar *c, PetscInt cm, PetscInt cn)
3804: {
3805: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
3806: const MatScalar *v;
3807: PetscInt mbs, i, *idx, *ii, j, n, k, *ridx = NULL;
3808: PetscBool usecprow = a->compressedrow.use;
3809: const PetscScalar *bi;
3810: PetscScalar x1, x2, x3;
3812: PetscFunctionBegin;
3813: idx = a->j;
3814: v = a->a;
3815: if (usecprow) {
3816: mbs = a->compressedrow.nrows;
3817: ii = a->compressedrow.i;
3818: ridx = a->compressedrow.rindex;
3819: } else {
3820: mbs = a->mbs;
3821: ii = a->i;
3822: }
3824: for (i = 0; i < mbs; i++) {
3825: PetscInt brow = usecprow ? ridx[i] : i;
3827: n = ii[1] - ii[0];
3828: ii++;
3829: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
3830: PetscPrefetchBlock(v + 9 * n, 9 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
3831: for (j = 0, bi = b + 3 * brow; j < n; j++) {
3832: PetscScalar *zcol = c + 3 * (*idx++);
3834: for (k = 0; k < cn; k++) {
3835: x1 = bi[0 + k * bm];
3836: x2 = bi[1 + k * bm];
3837: x3 = bi[2 + k * bm];
3838: zcol[0 + k * cm] += v[0] * x1 + v[1] * x2 + v[2] * x3;
3839: zcol[1 + k * cm] += v[3] * x1 + v[4] * x2 + v[5] * x3;
3840: zcol[2 + k * cm] += v[6] * x1 + v[7] * x2 + v[8] * x3;
3841: }
3842: v += 9;
3843: }
3844: }
3845: PetscFunctionReturn(PETSC_SUCCESS);
3846: }
3848: static PetscErrorCode MatTransposeMatMult_SeqBAIJ_4_Private(Mat A, PetscScalar *b, PetscInt bm, PetscScalar *c, PetscInt cm, PetscInt cn)
3849: {
3850: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
3851: const MatScalar *v;
3852: PetscInt mbs, i, *idx, *ii, j, n, k, *ridx = NULL;
3853: PetscBool usecprow = a->compressedrow.use;
3854: const PetscScalar *bi;
3855: PetscScalar x1, x2, x3, x4;
3857: PetscFunctionBegin;
3858: idx = a->j;
3859: v = a->a;
3860: if (usecprow) {
3861: mbs = a->compressedrow.nrows;
3862: ii = a->compressedrow.i;
3863: ridx = a->compressedrow.rindex;
3864: } else {
3865: mbs = a->mbs;
3866: ii = a->i;
3867: }
3869: for (i = 0; i < mbs; i++) {
3870: PetscInt brow = usecprow ? ridx[i] : i;
3872: n = ii[1] - ii[0];
3873: ii++;
3874: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
3875: PetscPrefetchBlock(v + 16 * n, 16 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
3876: for (j = 0, bi = b + 4 * brow; j < n; j++) {
3877: PetscScalar *zcol = c + 4 * (*idx++);
3879: for (k = 0; k < cn; k++) {
3880: x1 = bi[0 + k * bm];
3881: x2 = bi[1 + k * bm];
3882: x3 = bi[2 + k * bm];
3883: x4 = bi[3 + k * bm];
3884: zcol[0 + k * cm] += v[0] * x1 + v[1] * x2 + v[2] * x3 + v[3] * x4;
3885: zcol[1 + k * cm] += v[4] * x1 + v[5] * x2 + v[6] * x3 + v[7] * x4;
3886: zcol[2 + k * cm] += v[8] * x1 + v[9] * x2 + v[10] * x3 + v[11] * x4;
3887: zcol[3 + k * cm] += v[12] * x1 + v[13] * x2 + v[14] * x3 + v[15] * x4;
3888: }
3889: v += 16;
3890: }
3891: }
3892: PetscFunctionReturn(PETSC_SUCCESS);
3893: }
3895: static PetscErrorCode MatTransposeMatMult_SeqBAIJ_5_Private(Mat A, PetscScalar *b, PetscInt bm, PetscScalar *c, PetscInt cm, PetscInt cn)
3896: {
3897: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
3898: const MatScalar *v;
3899: PetscInt mbs, i, *idx, *ii, j, n, k, *ridx = NULL;
3900: PetscBool usecprow = a->compressedrow.use;
3901: const PetscScalar *bi;
3902: PetscScalar x1, x2, x3, x4, x5;
3904: PetscFunctionBegin;
3905: idx = a->j;
3906: v = a->a;
3907: if (usecprow) {
3908: mbs = a->compressedrow.nrows;
3909: ii = a->compressedrow.i;
3910: ridx = a->compressedrow.rindex;
3911: } else {
3912: mbs = a->mbs;
3913: ii = a->i;
3914: }
3916: for (i = 0; i < mbs; i++) {
3917: PetscInt brow = usecprow ? ridx[i] : i;
3919: n = ii[1] - ii[0];
3920: ii++;
3921: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
3922: PetscPrefetchBlock(v + 25 * n, 25 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
3923: for (j = 0, bi = b + 5 * brow; j < n; j++) {
3924: PetscScalar *zcol = c + 5 * (*idx++);
3926: for (k = 0; k < cn; k++) {
3927: x1 = bi[0 + k * bm];
3928: x2 = bi[1 + k * bm];
3929: x3 = bi[2 + k * bm];
3930: x4 = bi[3 + k * bm];
3931: x5 = bi[4 + k * bm];
3932: zcol[0 + k * cm] += v[0] * x1 + v[1] * x2 + v[2] * x3 + v[3] * x4 + v[4] * x5;
3933: zcol[1 + k * cm] += v[5] * x1 + v[6] * x2 + v[7] * x3 + v[8] * x4 + v[9] * x5;
3934: zcol[2 + k * cm] += v[10] * x1 + v[11] * x2 + v[12] * x3 + v[13] * x4 + v[14] * x5;
3935: zcol[3 + k * cm] += v[15] * x1 + v[16] * x2 + v[17] * x3 + v[18] * x4 + v[19] * x5;
3936: zcol[4 + k * cm] += v[20] * x1 + v[21] * x2 + v[22] * x3 + v[23] * x4 + v[24] * x5;
3937: }
3938: v += 25;
3939: }
3940: }
3941: PetscFunctionReturn(PETSC_SUCCESS);
3942: }
3944: PetscErrorCode MatTransposeMatMultNumeric_SeqBAIJ_SeqDense(Mat A, Mat B, Mat C)
3945: {
3946: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
3947: Mat_SeqDense *bd = (Mat_SeqDense *)B->data;
3948: Mat_SeqDense *cd = (Mat_SeqDense *)C->data;
3949: PetscInt cm = cd->lda, cn = B->cmap->n, bm = bd->lda;
3950: PetscInt mbs, i, bs = A->rmap->bs, j, n, bs2 = a->bs2;
3951: PetscBLASInt bbs, bcn, bbm, bcm;
3952: PetscScalar *c, *b;
3953: const MatScalar *v;
3954: const PetscInt *idx, *ii, *ridx = NULL;
3955: PetscScalar _DOne = 1.0;
3956: PetscBool usecprow = a->compressedrow.use;
3958: PetscFunctionBegin;
3959: if (!cm || !cn) PetscFunctionReturn(PETSC_SUCCESS);
3960: PetscCheck(B->rmap->n == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number rows in A %" PetscInt_FMT " not equal rows in B %" PetscInt_FMT, A->rmap->n, B->rmap->n);
3961: PetscCheck(A->cmap->n == C->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number rows in C %" PetscInt_FMT " not equal columns in A %" PetscInt_FMT, C->rmap->n, A->cmap->n);
3962: 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);
3963: b = bd->v;
3964: PetscCall(MatZeroEntries(C));
3965: PetscCall(MatDenseGetArrayWrite(C, &c));
3966: switch (bs) {
3967: case 1:
3968: PetscCall(MatTransposeMatMult_SeqBAIJ_1_Private(A, b, bm, c, cm, cn));
3969: break;
3970: case 2:
3971: PetscCall(MatTransposeMatMult_SeqBAIJ_2_Private(A, b, bm, c, cm, cn));
3972: break;
3973: case 3:
3974: PetscCall(MatTransposeMatMult_SeqBAIJ_3_Private(A, b, bm, c, cm, cn));
3975: break;
3976: case 4:
3977: PetscCall(MatTransposeMatMult_SeqBAIJ_4_Private(A, b, bm, c, cm, cn));
3978: break;
3979: case 5:
3980: PetscCall(MatTransposeMatMult_SeqBAIJ_5_Private(A, b, bm, c, cm, cn));
3981: break;
3982: default: /* block sizes larger than 5 by 5 are handled by BLAS */
3983: PetscCall(PetscBLASIntCast(bs, &bbs));
3984: PetscCall(PetscBLASIntCast(cn, &bcn));
3985: PetscCall(PetscBLASIntCast(bm, &bbm));
3986: PetscCall(PetscBLASIntCast(cm, &bcm));
3987: idx = a->j;
3988: v = a->a;
3989: if (usecprow) {
3990: mbs = a->compressedrow.nrows;
3991: ii = a->compressedrow.i;
3992: ridx = a->compressedrow.rindex;
3993: } else {
3994: mbs = a->mbs;
3995: ii = a->i;
3996: }
3997: for (i = 0; i < mbs; i++) {
3998: const PetscScalar *bi = b + bs * (usecprow ? ridx[i] : i);
4000: n = ii[1] - ii[0];
4001: ii++;
4002: for (j = 0; j < n; j++) {
4003: PetscCallBLAS("BLASgemm", BLASgemm_("T", "N", &bbs, &bcn, &bbs, &_DOne, v, &bbs, bi, &bbm, &_DOne, c + bs * (*idx++), &bcm));
4004: v += bs2;
4005: }
4006: }
4007: }
4008: PetscCall(MatDenseRestoreArrayWrite(C, &c));
4009: PetscCall(PetscLogFlops((2.0 * a->nz * bs2 - bs * a->nonzerorowcnt) * cn));
4010: PetscFunctionReturn(PETSC_SUCCESS);
4011: }