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