Actual source code: sbaij2.c
1: #include <../src/mat/impls/baij/seq/baij.h>
2: #include <../src/mat/impls/dense/seq/dense.h>
3: #include <../src/mat/impls/sbaij/seq/sbaij.h>
4: #include <petsc/private/kernels/blockinvert.h>
5: #include <petscbt.h>
6: #include <petscblaslapack.h>
8: PetscErrorCode MatIncreaseOverlap_SeqSBAIJ(Mat A, PetscInt is_max, IS is[], PetscInt ov)
9: {
10: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
11: PetscInt brow, i, j, k, l, mbs, n, *nidx, isz, bcol, bcol_max, start, end, *ai, *aj, bs;
12: const PetscInt *idx;
13: PetscBT table_out, table_in;
15: PetscFunctionBegin;
16: PetscCheck(ov >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative overlap specified");
17: mbs = a->mbs;
18: ai = a->i;
19: aj = a->j;
20: bs = A->rmap->bs;
21: PetscCall(PetscBTCreate(mbs, &table_out));
22: PetscCall(PetscMalloc1(mbs + 1, &nidx));
23: PetscCall(PetscBTCreate(mbs, &table_in));
25: for (i = 0; i < is_max; i++) { /* for each is */
26: isz = 0;
27: PetscCall(PetscBTMemzero(mbs, table_out));
29: /* Extract the indices, assume there can be duplicate entries */
30: PetscCall(ISGetIndices(is[i], &idx));
31: PetscCall(ISGetLocalSize(is[i], &n));
33: /* Enter these into the temp arrays i.e mark table_out[brow], enter brow into new index */
34: bcol_max = 0;
35: for (j = 0; j < n; ++j) {
36: brow = idx[j] / bs; /* convert the indices into block indices */
37: PetscCheck(brow < mbs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "index greater than mat-dim");
38: if (!PetscBTLookupSet(table_out, brow)) {
39: nidx[isz++] = brow;
40: if (bcol_max < brow) bcol_max = brow;
41: }
42: }
43: PetscCall(ISRestoreIndices(is[i], &idx));
44: PetscCall(ISDestroy(&is[i]));
46: k = 0;
47: for (j = 0; j < ov; j++) { /* for each overlap */
48: /* set table_in for lookup - only mark entries that are added onto nidx in (j-1)-th overlap */
49: PetscCall(PetscBTMemzero(mbs, table_in));
50: for (l = k; l < isz; l++) PetscCall(PetscBTSet(table_in, nidx[l]));
52: n = isz; /* length of the updated is[i] */
53: for (brow = 0; brow < mbs; brow++) {
54: start = ai[brow];
55: end = ai[brow + 1];
56: if (PetscBTLookup(table_in, brow)) { /* brow is on nidx - row search: collect all bcol in this brow */
57: for (l = start; l < end; l++) {
58: bcol = aj[l];
59: if (!PetscBTLookupSet(table_out, bcol)) {
60: nidx[isz++] = bcol;
61: if (bcol_max < bcol) bcol_max = bcol;
62: }
63: }
64: k++;
65: if (k >= n) break; /* for (brow=0; brow<mbs; brow++) */
66: } else { /* brow is not on nidx - col search: add brow onto nidx if there is a bcol in nidx */
67: for (l = start; l < end; l++) {
68: bcol = aj[l];
69: if (bcol > bcol_max) break;
70: if (PetscBTLookup(table_in, bcol)) {
71: if (!PetscBTLookupSet(table_out, brow)) nidx[isz++] = brow;
72: break; /* for l = start; l<end ; l++) */
73: }
74: }
75: }
76: }
77: } /* for each overlap */
78: PetscCall(ISCreateBlock(PETSC_COMM_SELF, bs, isz, nidx, PETSC_COPY_VALUES, is + i));
79: } /* for each is */
80: PetscCall(PetscBTDestroy(&table_out));
81: PetscCall(PetscFree(nidx));
82: PetscCall(PetscBTDestroy(&table_in));
83: PetscFunctionReturn(PETSC_SUCCESS);
84: }
86: /* Bseq is non-symmetric SBAIJ matrix, only used internally by PETSc.
87: Zero some ops' to avoid invalid use */
88: PetscErrorCode MatSeqSBAIJZeroOps_Private(Mat Bseq)
89: {
90: PetscFunctionBegin;
91: PetscCall(MatSetOption(Bseq, MAT_SYMMETRIC, PETSC_FALSE));
92: Bseq->ops->mult = NULL;
93: Bseq->ops->multadd = NULL;
94: Bseq->ops->multtranspose = NULL;
95: Bseq->ops->multtransposeadd = NULL;
96: Bseq->ops->lufactor = NULL;
97: Bseq->ops->choleskyfactor = NULL;
98: Bseq->ops->lufactorsymbolic = NULL;
99: Bseq->ops->choleskyfactorsymbolic = NULL;
100: Bseq->ops->getinertia = NULL;
101: PetscFunctionReturn(PETSC_SUCCESS);
102: }
104: /* same as MatCreateSubMatrices_SeqBAIJ(), except cast Mat_SeqSBAIJ */
105: static PetscErrorCode MatCreateSubMatrix_SeqSBAIJ_Private(Mat A, IS isrow, IS iscol, MatReuse scall, Mat *B, PetscBool sym)
106: {
107: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data, *c = NULL;
108: Mat_SeqBAIJ *d = NULL;
109: PetscInt *smap, i, k, kstart, kend, oldcols = a->nbs, *lens;
110: PetscInt row, mat_i, *mat_j, tcol, *mat_ilen;
111: const PetscInt *irow, *icol;
112: PetscInt nrows, ncols, *ssmap, bs = A->rmap->bs, bs2 = a->bs2;
113: PetscInt *aj = a->j, *ai = a->i;
114: MatScalar *mat_a;
115: Mat C;
116: PetscBool flag;
118: PetscFunctionBegin;
119: PetscCall(ISGetIndices(isrow, &irow));
120: PetscCall(ISGetIndices(iscol, &icol));
121: PetscCall(ISGetLocalSize(isrow, &nrows));
122: PetscCall(ISGetLocalSize(iscol, &ncols));
124: PetscCall(PetscCalloc1(1 + oldcols, &smap));
125: ssmap = smap;
126: PetscCall(PetscMalloc1(1 + nrows, &lens));
127: for (i = 0; i < ncols; i++) smap[icol[i]] = i + 1;
128: /* determine lens of each row */
129: for (i = 0; i < nrows; i++) {
130: kstart = ai[irow[i]];
131: kend = kstart + a->ilen[irow[i]];
132: lens[i] = 0;
133: for (k = kstart; k < kend; k++) {
134: if (ssmap[aj[k]]) lens[i]++;
135: }
136: }
137: /* Create and fill new matrix */
138: if (scall == MAT_REUSE_MATRIX) {
139: if (sym) {
140: c = (Mat_SeqSBAIJ *)((*B)->data);
142: PetscCheck(c->mbs == nrows && c->nbs == ncols && (*B)->rmap->bs == bs, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Submatrix wrong size");
143: PetscCall(PetscArraycmp(c->ilen, lens, c->mbs, &flag));
144: PetscCheck(flag, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Cannot reuse matrix. wrong number of nonzeros");
145: PetscCall(PetscArrayzero(c->ilen, c->mbs));
146: } else {
147: d = (Mat_SeqBAIJ *)((*B)->data);
149: PetscCheck(d->mbs == nrows && d->nbs == ncols && (*B)->rmap->bs == bs, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Submatrix wrong size");
150: PetscCall(PetscArraycmp(d->ilen, lens, d->mbs, &flag));
151: PetscCheck(flag, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Cannot reuse matrix. wrong number of nonzeros");
152: PetscCall(PetscArrayzero(d->ilen, d->mbs));
153: }
154: C = *B;
155: } else {
156: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C));
157: PetscCall(MatSetSizes(C, nrows * bs, ncols * bs, PETSC_DETERMINE, PETSC_DETERMINE));
158: if (sym) {
159: PetscCall(MatSetType(C, ((PetscObject)A)->type_name));
160: PetscCall(MatSeqSBAIJSetPreallocation(C, bs, 0, lens));
161: } else {
162: PetscCall(MatSetType(C, MATSEQBAIJ));
163: PetscCall(MatSeqBAIJSetPreallocation(C, bs, 0, lens));
164: }
165: }
166: if (sym) c = (Mat_SeqSBAIJ *)C->data;
167: else d = (Mat_SeqBAIJ *)C->data;
168: for (i = 0; i < nrows; i++) {
169: row = irow[i];
170: kstart = ai[row];
171: kend = kstart + a->ilen[row];
172: if (sym) {
173: mat_i = c->i[i];
174: mat_j = PetscSafePointerPlusOffset(c->j, mat_i);
175: mat_a = PetscSafePointerPlusOffset(c->a, mat_i * bs2);
176: mat_ilen = c->ilen + i;
177: } else {
178: mat_i = d->i[i];
179: mat_j = PetscSafePointerPlusOffset(d->j, mat_i);
180: mat_a = PetscSafePointerPlusOffset(d->a, mat_i * bs2);
181: mat_ilen = d->ilen + i;
182: }
183: for (k = kstart; k < kend; k++) {
184: if ((tcol = ssmap[a->j[k]])) {
185: *mat_j++ = tcol - 1;
186: PetscCall(PetscArraycpy(mat_a, a->a + k * bs2, bs2));
187: mat_a += bs2;
188: (*mat_ilen)++;
189: }
190: }
191: }
192: /* sort */
193: {
194: MatScalar *work;
196: PetscCall(PetscMalloc1(bs2, &work));
197: for (i = 0; i < nrows; i++) {
198: PetscInt ilen;
199: if (sym) {
200: mat_i = c->i[i];
201: mat_j = PetscSafePointerPlusOffset(c->j, mat_i);
202: mat_a = PetscSafePointerPlusOffset(c->a, mat_i * bs2);
203: ilen = c->ilen[i];
204: } else {
205: mat_i = d->i[i];
206: mat_j = PetscSafePointerPlusOffset(d->j, mat_i);
207: mat_a = PetscSafePointerPlusOffset(d->a, mat_i * bs2);
208: ilen = d->ilen[i];
209: }
210: PetscCall(PetscSortIntWithDataArray(ilen, mat_j, mat_a, bs2 * sizeof(MatScalar), work));
211: }
212: PetscCall(PetscFree(work));
213: }
215: /* Free work space */
216: PetscCall(ISRestoreIndices(iscol, &icol));
217: PetscCall(PetscFree(smap));
218: PetscCall(PetscFree(lens));
219: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
220: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
222: PetscCall(ISRestoreIndices(isrow, &irow));
223: *B = C;
224: PetscFunctionReturn(PETSC_SUCCESS);
225: }
227: PetscErrorCode MatCreateSubMatrix_SeqSBAIJ(Mat A, IS isrow, IS iscol, MatReuse scall, Mat *B)
228: {
229: Mat C[2];
230: IS is1, is2, intersect = NULL;
231: PetscInt n1, n2, ni;
232: PetscBool sym = PETSC_TRUE;
234: PetscFunctionBegin;
235: PetscCall(ISCompressIndicesGeneral(A->rmap->N, A->rmap->n, A->rmap->bs, 1, &isrow, &is1));
236: if (isrow == iscol) {
237: is2 = is1;
238: PetscCall(PetscObjectReference((PetscObject)is2));
239: } else {
240: PetscCall(ISCompressIndicesGeneral(A->cmap->N, A->cmap->n, A->cmap->bs, 1, &iscol, &is2));
241: PetscCall(ISIntersect(is1, is2, &intersect));
242: PetscCall(ISGetLocalSize(intersect, &ni));
243: PetscCall(ISDestroy(&intersect));
244: if (ni == 0) sym = PETSC_FALSE;
245: else if (PetscDefined(USE_DEBUG)) {
246: PetscCall(ISGetLocalSize(is1, &n1));
247: PetscCall(ISGetLocalSize(is2, &n2));
248: PetscCheck(ni == n1 && ni == n2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot create such a submatrix");
249: }
250: }
251: if (sym) PetscCall(MatCreateSubMatrix_SeqSBAIJ_Private(A, is1, is2, scall, B, sym));
252: else {
253: PetscCall(MatCreateSubMatrix_SeqSBAIJ_Private(A, is1, is2, MAT_INITIAL_MATRIX, C, sym));
254: PetscCall(MatCreateSubMatrix_SeqSBAIJ_Private(A, is2, is1, MAT_INITIAL_MATRIX, C + 1, sym));
255: PetscCall(MatTranspose(C[1], MAT_INPLACE_MATRIX, C + 1));
256: PetscCall(MatAXPY(C[0], 1.0, C[1], DIFFERENT_NONZERO_PATTERN));
257: PetscCheck(scall != MAT_INPLACE_MATRIX, PETSC_COMM_SELF, PETSC_ERR_SUP, "MAT_INPLACE_MATRIX not supported");
258: if (scall == MAT_REUSE_MATRIX) PetscCall(MatCopy(C[0], *B, SAME_NONZERO_PATTERN));
259: else if (A->rmap->bs == 1) PetscCall(MatConvert(C[0], MATAIJ, MAT_INITIAL_MATRIX, B));
260: else PetscCall(MatCopy(C[0], *B, SAME_NONZERO_PATTERN));
261: PetscCall(MatDestroy(C));
262: PetscCall(MatDestroy(C + 1));
263: }
264: PetscCall(ISDestroy(&is1));
265: PetscCall(ISDestroy(&is2));
267: if (sym && isrow != iscol) {
268: PetscBool isequal;
269: PetscCall(ISEqual(isrow, iscol, &isequal));
270: if (!isequal) PetscCall(MatSeqSBAIJZeroOps_Private(*B));
271: }
272: PetscFunctionReturn(PETSC_SUCCESS);
273: }
275: PetscErrorCode MatCreateSubMatrices_SeqSBAIJ(Mat A, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *B[])
276: {
277: PetscInt i;
279: PetscFunctionBegin;
280: if (scall == MAT_INITIAL_MATRIX) PetscCall(PetscCalloc1(n, B));
282: for (i = 0; i < n; i++) PetscCall(MatCreateSubMatrix_SeqSBAIJ(A, irow[i], icol[i], scall, &(*B)[i]));
283: PetscFunctionReturn(PETSC_SUCCESS);
284: }
286: /* Should check that shapes of vectors and matrices match */
287: PetscErrorCode MatMult_SeqSBAIJ_2(Mat A, Vec xx, Vec zz)
288: {
289: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
290: PetscScalar *z, x1, x2, zero = 0.0;
291: const PetscScalar *x, *xb;
292: const MatScalar *v;
293: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
294: const PetscInt *aj = a->j, *ai = a->i, *ib;
295: PetscInt nonzerorow = 0;
297: PetscFunctionBegin;
298: PetscCall(VecSet(zz, zero));
299: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
300: PetscCall(VecGetArrayRead(xx, &x));
301: PetscCall(VecGetArray(zz, &z));
303: v = a->a;
304: xb = x;
306: for (i = 0; i < mbs; i++, xb += 2, ai++) {
307: n = ai[1] - ai[0]; /* length of i_th block row of A */
308: if (!n) continue;
309: x1 = xb[0];
310: x2 = xb[1];
311: ib = aj + *ai;
312: jmin = 0;
313: nonzerorow++;
314: if (*ib == i) { /* (diag of A)*x */
315: z[2 * i] += v[0] * x1 + v[2] * x2;
316: z[2 * i + 1] += v[2] * x1 + v[3] * x2;
317: v += 4;
318: jmin++;
319: }
320: PetscPrefetchBlock(ib + jmin + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
321: PetscPrefetchBlock(v + 4 * n, 4 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
322: for (j = jmin; j < n; j++) {
323: /* (strict lower triangular part of A)*x */
324: cval = ib[j] * 2;
325: z[cval] += v[0] * x1 + v[1] * x2;
326: z[cval + 1] += v[2] * x1 + v[3] * x2;
327: /* (strict upper triangular part of A)*x */
328: z[2 * i] += v[0] * x[cval] + v[2] * x[cval + 1];
329: z[2 * i + 1] += v[1] * x[cval] + v[3] * x[cval + 1];
330: v += 4;
331: }
332: }
334: PetscCall(VecRestoreArrayRead(xx, &x));
335: PetscCall(VecRestoreArray(zz, &z));
336: PetscCall(PetscLogFlops(8.0 * (a->nz * 2.0 - nonzerorow) - nonzerorow));
337: PetscFunctionReturn(PETSC_SUCCESS);
338: }
340: PetscErrorCode MatMult_SeqSBAIJ_3(Mat A, Vec xx, Vec zz)
341: {
342: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
343: PetscScalar *z, x1, x2, x3, zero = 0.0;
344: const PetscScalar *x, *xb;
345: const MatScalar *v;
346: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
347: const PetscInt *aj = a->j, *ai = a->i, *ib;
348: PetscInt nonzerorow = 0;
350: PetscFunctionBegin;
351: PetscCall(VecSet(zz, zero));
352: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
353: PetscCall(VecGetArrayRead(xx, &x));
354: PetscCall(VecGetArray(zz, &z));
356: v = a->a;
357: xb = x;
359: for (i = 0; i < mbs; i++, xb += 3, ai++) {
360: n = ai[1] - ai[0]; /* length of i_th block row of A */
361: if (!n) continue;
362: x1 = xb[0];
363: x2 = xb[1];
364: x3 = xb[2];
365: ib = aj + *ai;
366: jmin = 0;
367: nonzerorow++;
368: if (*ib == i) { /* (diag of A)*x */
369: z[3 * i] += v[0] * x1 + v[3] * x2 + v[6] * x3;
370: z[3 * i + 1] += v[3] * x1 + v[4] * x2 + v[7] * x3;
371: z[3 * i + 2] += v[6] * x1 + v[7] * x2 + v[8] * x3;
372: v += 9;
373: jmin++;
374: }
375: PetscPrefetchBlock(ib + jmin + 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 = jmin; j < n; j++) {
378: /* (strict lower triangular part of A)*x */
379: cval = ib[j] * 3;
380: z[cval] += v[0] * x1 + v[1] * x2 + v[2] * x3;
381: z[cval + 1] += v[3] * x1 + v[4] * x2 + v[5] * x3;
382: z[cval + 2] += v[6] * x1 + v[7] * x2 + v[8] * x3;
383: /* (strict upper triangular part of A)*x */
384: z[3 * i] += v[0] * x[cval] + v[3] * x[cval + 1] + v[6] * x[cval + 2];
385: z[3 * i + 1] += v[1] * x[cval] + v[4] * x[cval + 1] + v[7] * x[cval + 2];
386: z[3 * i + 2] += v[2] * x[cval] + v[5] * x[cval + 1] + v[8] * x[cval + 2];
387: v += 9;
388: }
389: }
391: PetscCall(VecRestoreArrayRead(xx, &x));
392: PetscCall(VecRestoreArray(zz, &z));
393: PetscCall(PetscLogFlops(18.0 * (a->nz * 2.0 - nonzerorow) - nonzerorow));
394: PetscFunctionReturn(PETSC_SUCCESS);
395: }
397: PetscErrorCode MatMult_SeqSBAIJ_4(Mat A, Vec xx, Vec zz)
398: {
399: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
400: PetscScalar *z, x1, x2, x3, x4, zero = 0.0;
401: const PetscScalar *x, *xb;
402: const MatScalar *v;
403: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
404: const PetscInt *aj = a->j, *ai = a->i, *ib;
405: PetscInt nonzerorow = 0;
407: PetscFunctionBegin;
408: PetscCall(VecSet(zz, zero));
409: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
410: PetscCall(VecGetArrayRead(xx, &x));
411: PetscCall(VecGetArray(zz, &z));
413: v = a->a;
414: xb = x;
416: for (i = 0; i < mbs; i++, xb += 4, ai++) {
417: n = ai[1] - ai[0]; /* length of i_th block row of A */
418: if (!n) continue;
419: x1 = xb[0];
420: x2 = xb[1];
421: x3 = xb[2];
422: x4 = xb[3];
423: ib = aj + *ai;
424: jmin = 0;
425: nonzerorow++;
426: if (*ib == i) { /* (diag of A)*x */
427: z[4 * i] += v[0] * x1 + v[4] * x2 + v[8] * x3 + v[12] * x4;
428: z[4 * i + 1] += v[4] * x1 + v[5] * x2 + v[9] * x3 + v[13] * x4;
429: z[4 * i + 2] += v[8] * x1 + v[9] * x2 + v[10] * x3 + v[14] * x4;
430: z[4 * i + 3] += v[12] * x1 + v[13] * x2 + v[14] * x3 + v[15] * x4;
431: v += 16;
432: jmin++;
433: }
434: PetscPrefetchBlock(ib + jmin + 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 = jmin; j < n; j++) {
437: /* (strict lower triangular part of A)*x */
438: cval = ib[j] * 4;
439: z[cval] += v[0] * x1 + v[1] * x2 + v[2] * x3 + v[3] * x4;
440: z[cval + 1] += v[4] * x1 + v[5] * x2 + v[6] * x3 + v[7] * x4;
441: z[cval + 2] += v[8] * x1 + v[9] * x2 + v[10] * x3 + v[11] * x4;
442: z[cval + 3] += v[12] * x1 + v[13] * x2 + v[14] * x3 + v[15] * x4;
443: /* (strict upper triangular part of A)*x */
444: z[4 * i] += v[0] * x[cval] + v[4] * x[cval + 1] + v[8] * x[cval + 2] + v[12] * x[cval + 3];
445: z[4 * i + 1] += v[1] * x[cval] + v[5] * x[cval + 1] + v[9] * x[cval + 2] + v[13] * x[cval + 3];
446: z[4 * i + 2] += v[2] * x[cval] + v[6] * x[cval + 1] + v[10] * x[cval + 2] + v[14] * x[cval + 3];
447: z[4 * i + 3] += v[3] * x[cval] + v[7] * x[cval + 1] + v[11] * x[cval + 2] + v[15] * x[cval + 3];
448: v += 16;
449: }
450: }
452: PetscCall(VecRestoreArrayRead(xx, &x));
453: PetscCall(VecRestoreArray(zz, &z));
454: PetscCall(PetscLogFlops(32.0 * (a->nz * 2.0 - nonzerorow) - nonzerorow));
455: PetscFunctionReturn(PETSC_SUCCESS);
456: }
458: PetscErrorCode MatMult_SeqSBAIJ_5(Mat A, Vec xx, Vec zz)
459: {
460: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
461: PetscScalar *z, x1, x2, x3, x4, x5, zero = 0.0;
462: const PetscScalar *x, *xb;
463: const MatScalar *v;
464: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
465: const PetscInt *aj = a->j, *ai = a->i, *ib;
466: PetscInt nonzerorow = 0;
468: PetscFunctionBegin;
469: PetscCall(VecSet(zz, zero));
470: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
471: PetscCall(VecGetArrayRead(xx, &x));
472: PetscCall(VecGetArray(zz, &z));
474: v = a->a;
475: xb = x;
477: for (i = 0; i < mbs; i++, xb += 5, ai++) {
478: n = ai[1] - ai[0]; /* length of i_th block row of A */
479: if (!n) continue;
480: x1 = xb[0];
481: x2 = xb[1];
482: x3 = xb[2];
483: x4 = xb[3];
484: x5 = xb[4];
485: ib = aj + *ai;
486: jmin = 0;
487: nonzerorow++;
488: if (*ib == i) { /* (diag of A)*x */
489: z[5 * i] += v[0] * x1 + v[5] * x2 + v[10] * x3 + v[15] * x4 + v[20] * x5;
490: z[5 * i + 1] += v[5] * x1 + v[6] * x2 + v[11] * x3 + v[16] * x4 + v[21] * x5;
491: z[5 * i + 2] += v[10] * x1 + v[11] * x2 + v[12] * x3 + v[17] * x4 + v[22] * x5;
492: z[5 * i + 3] += v[15] * x1 + v[16] * x2 + v[17] * x3 + v[18] * x4 + v[23] * x5;
493: z[5 * i + 4] += v[20] * x1 + v[21] * x2 + v[22] * x3 + v[23] * x4 + v[24] * x5;
494: v += 25;
495: jmin++;
496: }
497: PetscPrefetchBlock(ib + jmin + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
498: PetscPrefetchBlock(v + 25 * n, 25 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
499: for (j = jmin; j < n; j++) {
500: /* (strict lower triangular part of A)*x */
501: cval = ib[j] * 5;
502: z[cval] += v[0] * x1 + v[1] * x2 + v[2] * x3 + v[3] * x4 + v[4] * x5;
503: z[cval + 1] += v[5] * x1 + v[6] * x2 + v[7] * x3 + v[8] * x4 + v[9] * x5;
504: z[cval + 2] += v[10] * x1 + v[11] * x2 + v[12] * x3 + v[13] * x4 + v[14] * x5;
505: z[cval + 3] += v[15] * x1 + v[16] * x2 + v[17] * x3 + v[18] * x4 + v[19] * x5;
506: z[cval + 4] += v[20] * x1 + v[21] * x2 + v[22] * x3 + v[23] * x4 + v[24] * x5;
507: /* (strict upper triangular part of A)*x */
508: z[5 * i] += v[0] * x[cval] + v[5] * x[cval + 1] + v[10] * x[cval + 2] + v[15] * x[cval + 3] + v[20] * x[cval + 4];
509: z[5 * i + 1] += v[1] * x[cval] + v[6] * x[cval + 1] + v[11] * x[cval + 2] + v[16] * x[cval + 3] + v[21] * x[cval + 4];
510: z[5 * i + 2] += v[2] * x[cval] + v[7] * x[cval + 1] + v[12] * x[cval + 2] + v[17] * x[cval + 3] + v[22] * x[cval + 4];
511: z[5 * i + 3] += v[3] * x[cval] + v[8] * x[cval + 1] + v[13] * x[cval + 2] + v[18] * x[cval + 3] + v[23] * x[cval + 4];
512: z[5 * i + 4] += v[4] * x[cval] + v[9] * x[cval + 1] + v[14] * x[cval + 2] + v[19] * x[cval + 3] + v[24] * x[cval + 4];
513: v += 25;
514: }
515: }
517: PetscCall(VecRestoreArrayRead(xx, &x));
518: PetscCall(VecRestoreArray(zz, &z));
519: PetscCall(PetscLogFlops(50.0 * (a->nz * 2.0 - nonzerorow) - nonzerorow));
520: PetscFunctionReturn(PETSC_SUCCESS);
521: }
523: PetscErrorCode MatMult_SeqSBAIJ_6(Mat A, Vec xx, Vec zz)
524: {
525: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
526: PetscScalar *z, x1, x2, x3, x4, x5, x6, zero = 0.0;
527: const PetscScalar *x, *xb;
528: const MatScalar *v;
529: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
530: const PetscInt *aj = a->j, *ai = a->i, *ib;
531: PetscInt nonzerorow = 0;
533: PetscFunctionBegin;
534: PetscCall(VecSet(zz, zero));
535: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
536: PetscCall(VecGetArrayRead(xx, &x));
537: PetscCall(VecGetArray(zz, &z));
539: v = a->a;
540: xb = x;
542: for (i = 0; i < mbs; i++, xb += 6, ai++) {
543: n = ai[1] - ai[0]; /* length of i_th block row of A */
544: if (!n) continue;
545: x1 = xb[0];
546: x2 = xb[1];
547: x3 = xb[2];
548: x4 = xb[3];
549: x5 = xb[4];
550: x6 = xb[5];
551: ib = aj + *ai;
552: jmin = 0;
553: nonzerorow++;
554: if (*ib == i) { /* (diag of A)*x */
555: z[6 * i] += v[0] * x1 + v[6] * x2 + v[12] * x3 + v[18] * x4 + v[24] * x5 + v[30] * x6;
556: z[6 * i + 1] += v[6] * x1 + v[7] * x2 + v[13] * x3 + v[19] * x4 + v[25] * x5 + v[31] * x6;
557: z[6 * i + 2] += v[12] * x1 + v[13] * x2 + v[14] * x3 + v[20] * x4 + v[26] * x5 + v[32] * x6;
558: z[6 * i + 3] += v[18] * x1 + v[19] * x2 + v[20] * x3 + v[21] * x4 + v[27] * x5 + v[33] * x6;
559: z[6 * i + 4] += v[24] * x1 + v[25] * x2 + v[26] * x3 + v[27] * x4 + v[28] * x5 + v[34] * x6;
560: z[6 * i + 5] += v[30] * x1 + v[31] * x2 + v[32] * x3 + v[33] * x4 + v[34] * x5 + v[35] * x6;
561: v += 36;
562: jmin++;
563: }
564: PetscPrefetchBlock(ib + jmin + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
565: PetscPrefetchBlock(v + 36 * n, 36 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
566: for (j = jmin; j < n; j++) {
567: /* (strict lower triangular part of A)*x */
568: cval = ib[j] * 6;
569: z[cval] += v[0] * x1 + v[1] * x2 + v[2] * x3 + v[3] * x4 + v[4] * x5 + v[5] * x6;
570: z[cval + 1] += v[6] * x1 + v[7] * x2 + v[8] * x3 + v[9] * x4 + v[10] * x5 + v[11] * x6;
571: z[cval + 2] += v[12] * x1 + v[13] * x2 + v[14] * x3 + v[15] * x4 + v[16] * x5 + v[17] * x6;
572: z[cval + 3] += v[18] * x1 + v[19] * x2 + v[20] * x3 + v[21] * x4 + v[22] * x5 + v[23] * x6;
573: z[cval + 4] += v[24] * x1 + v[25] * x2 + v[26] * x3 + v[27] * x4 + v[28] * x5 + v[29] * x6;
574: z[cval + 5] += v[30] * x1 + v[31] * x2 + v[32] * x3 + v[33] * x4 + v[34] * x5 + v[35] * x6;
575: /* (strict upper triangular part of A)*x */
576: z[6 * i] += v[0] * x[cval] + v[6] * x[cval + 1] + v[12] * x[cval + 2] + v[18] * x[cval + 3] + v[24] * x[cval + 4] + v[30] * x[cval + 5];
577: z[6 * i + 1] += v[1] * x[cval] + v[7] * x[cval + 1] + v[13] * x[cval + 2] + v[19] * x[cval + 3] + v[25] * x[cval + 4] + v[31] * x[cval + 5];
578: z[6 * i + 2] += v[2] * x[cval] + v[8] * x[cval + 1] + v[14] * x[cval + 2] + v[20] * x[cval + 3] + v[26] * x[cval + 4] + v[32] * x[cval + 5];
579: z[6 * i + 3] += v[3] * x[cval] + v[9] * x[cval + 1] + v[15] * x[cval + 2] + v[21] * x[cval + 3] + v[27] * x[cval + 4] + v[33] * x[cval + 5];
580: z[6 * i + 4] += v[4] * x[cval] + v[10] * x[cval + 1] + v[16] * x[cval + 2] + v[22] * x[cval + 3] + v[28] * x[cval + 4] + v[34] * x[cval + 5];
581: z[6 * i + 5] += v[5] * x[cval] + v[11] * x[cval + 1] + v[17] * x[cval + 2] + v[23] * x[cval + 3] + v[29] * x[cval + 4] + v[35] * x[cval + 5];
582: v += 36;
583: }
584: }
586: PetscCall(VecRestoreArrayRead(xx, &x));
587: PetscCall(VecRestoreArray(zz, &z));
588: PetscCall(PetscLogFlops(72.0 * (a->nz * 2.0 - nonzerorow) - nonzerorow));
589: PetscFunctionReturn(PETSC_SUCCESS);
590: }
592: PetscErrorCode MatMult_SeqSBAIJ_7(Mat A, Vec xx, Vec zz)
593: {
594: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
595: PetscScalar *z, x1, x2, x3, x4, x5, x6, x7, zero = 0.0;
596: const PetscScalar *x, *xb;
597: const MatScalar *v;
598: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
599: const PetscInt *aj = a->j, *ai = a->i, *ib;
600: PetscInt nonzerorow = 0;
602: PetscFunctionBegin;
603: PetscCall(VecSet(zz, zero));
604: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
605: PetscCall(VecGetArrayRead(xx, &x));
606: PetscCall(VecGetArray(zz, &z));
608: v = a->a;
609: xb = x;
611: for (i = 0; i < mbs; i++, xb += 7, ai++) {
612: n = ai[1] - ai[0]; /* length of i_th block row of A */
613: if (!n) continue;
614: x1 = xb[0];
615: x2 = xb[1];
616: x3 = xb[2];
617: x4 = xb[3];
618: x5 = xb[4];
619: x6 = xb[5];
620: x7 = xb[6];
621: ib = aj + *ai;
622: jmin = 0;
623: nonzerorow++;
624: if (*ib == i) { /* (diag of A)*x */
625: z[7 * i] += v[0] * x1 + v[7] * x2 + v[14] * x3 + v[21] * x4 + v[28] * x5 + v[35] * x6 + v[42] * x7;
626: z[7 * i + 1] += v[7] * x1 + v[8] * x2 + v[15] * x3 + v[22] * x4 + v[29] * x5 + v[36] * x6 + v[43] * x7;
627: z[7 * i + 2] += v[14] * x1 + v[15] * x2 + v[16] * x3 + v[23] * x4 + v[30] * x5 + v[37] * x6 + v[44] * x7;
628: z[7 * i + 3] += v[21] * x1 + v[22] * x2 + v[23] * x3 + v[24] * x4 + v[31] * x5 + v[38] * x6 + v[45] * x7;
629: z[7 * i + 4] += v[28] * x1 + v[29] * x2 + v[30] * x3 + v[31] * x4 + v[32] * x5 + v[39] * x6 + v[46] * x7;
630: z[7 * i + 5] += v[35] * x1 + v[36] * x2 + v[37] * x3 + v[38] * x4 + v[39] * x5 + v[40] * x6 + v[47] * x7;
631: z[7 * i + 6] += v[42] * x1 + v[43] * x2 + v[44] * x3 + v[45] * x4 + v[46] * x5 + v[47] * x6 + v[48] * x7;
632: v += 49;
633: jmin++;
634: }
635: PetscPrefetchBlock(ib + jmin + 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 = jmin; j < n; j++) {
638: /* (strict lower triangular part of A)*x */
639: cval = ib[j] * 7;
640: z[cval] += v[0] * x1 + v[1] * x2 + v[2] * x3 + v[3] * x4 + v[4] * x5 + v[5] * x6 + v[6] * x7;
641: z[cval + 1] += v[7] * x1 + v[8] * x2 + v[9] * x3 + v[10] * x4 + v[11] * x5 + v[12] * x6 + v[13] * x7;
642: z[cval + 2] += v[14] * x1 + v[15] * x2 + v[16] * x3 + v[17] * x4 + v[18] * x5 + v[19] * x6 + v[20] * x7;
643: z[cval + 3] += v[21] * x1 + v[22] * x2 + v[23] * x3 + v[24] * x4 + v[25] * x5 + v[26] * x6 + v[27] * x7;
644: z[cval + 4] += v[28] * x1 + v[29] * x2 + v[30] * x3 + v[31] * x4 + v[32] * x5 + v[33] * x6 + v[34] * x7;
645: z[cval + 5] += v[35] * x1 + v[36] * x2 + v[37] * x3 + v[38] * x4 + v[39] * x5 + v[40] * x6 + v[41] * x7;
646: z[cval + 6] += v[42] * x1 + v[43] * x2 + v[44] * x3 + v[45] * x4 + v[46] * x5 + v[47] * x6 + v[48] * x7;
647: /* (strict upper triangular part of A)*x */
648: z[7 * i] += v[0] * x[cval] + v[7] * x[cval + 1] + v[14] * x[cval + 2] + v[21] * x[cval + 3] + v[28] * x[cval + 4] + v[35] * x[cval + 5] + v[42] * x[cval + 6];
649: z[7 * i + 1] += v[1] * x[cval] + v[8] * x[cval + 1] + v[15] * x[cval + 2] + v[22] * x[cval + 3] + v[29] * x[cval + 4] + v[36] * x[cval + 5] + v[43] * x[cval + 6];
650: z[7 * i + 2] += v[2] * x[cval] + v[9] * x[cval + 1] + v[16] * x[cval + 2] + v[23] * x[cval + 3] + v[30] * x[cval + 4] + v[37] * x[cval + 5] + v[44] * x[cval + 6];
651: z[7 * i + 3] += v[3] * x[cval] + v[10] * x[cval + 1] + v[17] * x[cval + 2] + v[24] * x[cval + 3] + v[31] * x[cval + 4] + v[38] * x[cval + 5] + v[45] * x[cval + 6];
652: z[7 * i + 4] += v[4] * x[cval] + v[11] * x[cval + 1] + v[18] * x[cval + 2] + v[25] * x[cval + 3] + v[32] * x[cval + 4] + v[39] * x[cval + 5] + v[46] * x[cval + 6];
653: z[7 * i + 5] += v[5] * x[cval] + v[12] * x[cval + 1] + v[19] * x[cval + 2] + v[26] * x[cval + 3] + v[33] * x[cval + 4] + v[40] * x[cval + 5] + v[47] * x[cval + 6];
654: z[7 * i + 6] += v[6] * x[cval] + v[13] * x[cval + 1] + v[20] * x[cval + 2] + v[27] * x[cval + 3] + v[34] * x[cval + 4] + v[41] * x[cval + 5] + v[48] * x[cval + 6];
655: v += 49;
656: }
657: }
658: PetscCall(VecRestoreArrayRead(xx, &x));
659: PetscCall(VecRestoreArray(zz, &z));
660: PetscCall(PetscLogFlops(98.0 * (a->nz * 2.0 - nonzerorow) - nonzerorow));
661: PetscFunctionReturn(PETSC_SUCCESS);
662: }
664: /*
665: This will not work with MatScalar == float because it calls the BLAS
666: */
667: PetscErrorCode MatMult_SeqSBAIJ_N(Mat A, Vec xx, Vec zz)
668: {
669: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
670: PetscScalar *z, *z_ptr, *zb, *work, *workt, zero = 0.0;
671: const PetscScalar *x, *x_ptr, *xb;
672: const MatScalar *v;
673: PetscInt mbs = a->mbs, i, bs = A->rmap->bs, j, n, bs2 = a->bs2, ncols, k;
674: const PetscInt *idx, *aj, *ii;
675: PetscInt nonzerorow = 0;
677: PetscFunctionBegin;
678: PetscCall(VecSet(zz, zero));
679: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
680: PetscCall(VecGetArrayRead(xx, &x));
681: PetscCall(VecGetArray(zz, &z));
683: x_ptr = x;
684: z_ptr = z;
686: aj = a->j;
687: v = a->a;
688: ii = a->i;
690: if (!a->mult_work) PetscCall(PetscMalloc1(A->rmap->N + 1, &a->mult_work));
691: work = a->mult_work;
693: for (i = 0; i < mbs; i++) {
694: n = ii[1] - ii[0];
695: ncols = n * bs;
696: workt = work;
697: idx = aj + ii[0];
698: nonzerorow += (n > 0);
700: /* upper triangular part */
701: for (j = 0; j < n; j++) {
702: xb = x_ptr + bs * (*idx++);
703: for (k = 0; k < bs; k++) workt[k] = xb[k];
704: workt += bs;
705: }
706: /* z(i*bs:(i+1)*bs-1) += A(i,:)*x */
707: PetscKernel_w_gets_w_plus_Ar_times_v(bs, ncols, work, v, z);
709: /* strict lower triangular part */
710: idx = aj + ii[0];
711: if (n && *idx == i) {
712: ncols -= bs;
713: v += bs2;
714: idx++;
715: n--;
716: }
718: if (ncols > 0) {
719: workt = work;
720: PetscCall(PetscArrayzero(workt, ncols));
721: PetscKernel_w_gets_w_plus_trans_Ar_times_v(bs, ncols, x, v, workt);
722: for (j = 0; j < n; j++) {
723: zb = z_ptr + bs * (*idx++);
724: for (k = 0; k < bs; k++) zb[k] += workt[k];
725: workt += bs;
726: }
727: }
728: x += bs;
729: v += n * bs2;
730: z += bs;
731: ii++;
732: }
734: PetscCall(VecRestoreArrayRead(xx, &x));
735: PetscCall(VecRestoreArray(zz, &z));
736: PetscCall(PetscLogFlops(2.0 * (a->nz * 2.0 - nonzerorow) * bs2 - nonzerorow));
737: PetscFunctionReturn(PETSC_SUCCESS);
738: }
740: PetscErrorCode MatMultAdd_SeqSBAIJ_1(Mat A, Vec xx, Vec yy, Vec zz)
741: {
742: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
743: PetscScalar *z, x1;
744: const PetscScalar *x, *xb;
745: const MatScalar *v;
746: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
747: const PetscInt *aj = a->j, *ai = a->i, *ib;
748: PetscInt nonzerorow = 0;
749: const int aconj = PetscDefined(USE_COMPLEX) && A->hermitian == PETSC_BOOL3_TRUE ? 1 : 0;
751: PetscFunctionBegin;
752: PetscCall(VecCopy(yy, zz));
753: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
754: PetscCall(VecGetArrayRead(xx, &x));
755: PetscCall(VecGetArray(zz, &z));
756: v = a->a;
757: xb = x;
759: for (i = 0; i < mbs; i++, xb++, ai++) {
760: n = ai[1] - ai[0]; /* length of i_th row of A */
761: if (!n) continue;
762: x1 = xb[0];
763: ib = aj + *ai;
764: jmin = 0;
765: nonzerorow++;
766: if (*ib == i) { /* (diag of A)*x */
767: z[i] += *v++ * x[*ib++];
768: jmin++;
769: }
770: if (aconj) {
771: for (j = jmin; j < n; j++) {
772: cval = *ib;
773: z[cval] += PetscConj(*v) * x1; /* (strict lower triangular part of A)*x */
774: z[i] += *v++ * x[*ib++]; /* (strict upper triangular part of A)*x */
775: }
776: } else {
777: for (j = jmin; j < n; j++) {
778: cval = *ib;
779: z[cval] += *v * x1; /* (strict lower triangular part of A)*x */
780: z[i] += *v++ * x[*ib++]; /* (strict upper triangular part of A)*x */
781: }
782: }
783: }
785: PetscCall(VecRestoreArrayRead(xx, &x));
786: PetscCall(VecRestoreArray(zz, &z));
788: PetscCall(PetscLogFlops(2.0 * (a->nz * 2.0 - nonzerorow)));
789: PetscFunctionReturn(PETSC_SUCCESS);
790: }
792: PetscErrorCode MatMultAdd_SeqSBAIJ_2(Mat A, Vec xx, Vec yy, Vec zz)
793: {
794: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
795: PetscScalar *z, x1, x2;
796: const PetscScalar *x, *xb;
797: const MatScalar *v;
798: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
799: const PetscInt *aj = a->j, *ai = a->i, *ib;
800: PetscInt nonzerorow = 0;
802: PetscFunctionBegin;
803: PetscCall(VecCopy(yy, zz));
804: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
805: PetscCall(VecGetArrayRead(xx, &x));
806: PetscCall(VecGetArray(zz, &z));
808: v = a->a;
809: xb = x;
811: for (i = 0; i < mbs; i++, xb += 2, ai++) {
812: n = ai[1] - ai[0]; /* length of i_th block row of A */
813: if (!n) continue;
814: x1 = xb[0];
815: x2 = xb[1];
816: ib = aj + *ai;
817: jmin = 0;
818: nonzerorow++;
819: if (*ib == i) { /* (diag of A)*x */
820: z[2 * i] += v[0] * x1 + v[2] * x2;
821: z[2 * i + 1] += v[2] * x1 + v[3] * x2;
822: v += 4;
823: jmin++;
824: }
825: PetscPrefetchBlock(ib + jmin + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
826: PetscPrefetchBlock(v + 4 * n, 4 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
827: for (j = jmin; j < n; j++) {
828: /* (strict lower triangular part of A)*x */
829: cval = ib[j] * 2;
830: z[cval] += v[0] * x1 + v[1] * x2;
831: z[cval + 1] += v[2] * x1 + v[3] * x2;
832: /* (strict upper triangular part of A)*x */
833: z[2 * i] += v[0] * x[cval] + v[2] * x[cval + 1];
834: z[2 * i + 1] += v[1] * x[cval] + v[3] * x[cval + 1];
835: v += 4;
836: }
837: }
838: PetscCall(VecRestoreArrayRead(xx, &x));
839: PetscCall(VecRestoreArray(zz, &z));
841: PetscCall(PetscLogFlops(4.0 * (a->nz * 2.0 - nonzerorow)));
842: PetscFunctionReturn(PETSC_SUCCESS);
843: }
845: PetscErrorCode MatMultAdd_SeqSBAIJ_3(Mat A, Vec xx, Vec yy, Vec zz)
846: {
847: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
848: PetscScalar *z, x1, x2, x3;
849: const PetscScalar *x, *xb;
850: const MatScalar *v;
851: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
852: const PetscInt *aj = a->j, *ai = a->i, *ib;
853: PetscInt nonzerorow = 0;
855: PetscFunctionBegin;
856: PetscCall(VecCopy(yy, zz));
857: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
858: PetscCall(VecGetArrayRead(xx, &x));
859: PetscCall(VecGetArray(zz, &z));
861: v = a->a;
862: xb = x;
864: for (i = 0; i < mbs; i++, xb += 3, ai++) {
865: n = ai[1] - ai[0]; /* length of i_th block row of A */
866: if (!n) continue;
867: x1 = xb[0];
868: x2 = xb[1];
869: x3 = xb[2];
870: ib = aj + *ai;
871: jmin = 0;
872: nonzerorow++;
873: if (*ib == i) { /* (diag of A)*x */
874: z[3 * i] += v[0] * x1 + v[3] * x2 + v[6] * x3;
875: z[3 * i + 1] += v[3] * x1 + v[4] * x2 + v[7] * x3;
876: z[3 * i + 2] += v[6] * x1 + v[7] * x2 + v[8] * x3;
877: v += 9;
878: jmin++;
879: }
880: PetscPrefetchBlock(ib + jmin + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
881: PetscPrefetchBlock(v + 9 * n, 9 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
882: for (j = jmin; j < n; j++) {
883: /* (strict lower triangular part of A)*x */
884: cval = ib[j] * 3;
885: z[cval] += v[0] * x1 + v[1] * x2 + v[2] * x3;
886: z[cval + 1] += v[3] * x1 + v[4] * x2 + v[5] * x3;
887: z[cval + 2] += v[6] * x1 + v[7] * x2 + v[8] * x3;
888: /* (strict upper triangular part of A)*x */
889: z[3 * i] += v[0] * x[cval] + v[3] * x[cval + 1] + v[6] * x[cval + 2];
890: z[3 * i + 1] += v[1] * x[cval] + v[4] * x[cval + 1] + v[7] * x[cval + 2];
891: z[3 * i + 2] += v[2] * x[cval] + v[5] * x[cval + 1] + v[8] * x[cval + 2];
892: v += 9;
893: }
894: }
896: PetscCall(VecRestoreArrayRead(xx, &x));
897: PetscCall(VecRestoreArray(zz, &z));
899: PetscCall(PetscLogFlops(18.0 * (a->nz * 2.0 - nonzerorow)));
900: PetscFunctionReturn(PETSC_SUCCESS);
901: }
903: PetscErrorCode MatMultAdd_SeqSBAIJ_4(Mat A, Vec xx, Vec yy, Vec zz)
904: {
905: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
906: PetscScalar *z, x1, x2, x3, x4;
907: const PetscScalar *x, *xb;
908: const MatScalar *v;
909: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
910: const PetscInt *aj = a->j, *ai = a->i, *ib;
911: PetscInt nonzerorow = 0;
913: PetscFunctionBegin;
914: PetscCall(VecCopy(yy, zz));
915: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
916: PetscCall(VecGetArrayRead(xx, &x));
917: PetscCall(VecGetArray(zz, &z));
919: v = a->a;
920: xb = x;
922: for (i = 0; i < mbs; i++, xb += 4, ai++) {
923: n = ai[1] - ai[0]; /* length of i_th block row of A */
924: if (!n) continue;
925: x1 = xb[0];
926: x2 = xb[1];
927: x3 = xb[2];
928: x4 = xb[3];
929: ib = aj + *ai;
930: jmin = 0;
931: nonzerorow++;
932: if (*ib == i) { /* (diag of A)*x */
933: z[4 * i] += v[0] * x1 + v[4] * x2 + v[8] * x3 + v[12] * x4;
934: z[4 * i + 1] += v[4] * x1 + v[5] * x2 + v[9] * x3 + v[13] * x4;
935: z[4 * i + 2] += v[8] * x1 + v[9] * x2 + v[10] * x3 + v[14] * x4;
936: z[4 * i + 3] += v[12] * x1 + v[13] * x2 + v[14] * x3 + v[15] * x4;
937: v += 16;
938: jmin++;
939: }
940: PetscPrefetchBlock(ib + jmin + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
941: PetscPrefetchBlock(v + 16 * n, 16 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
942: for (j = jmin; j < n; j++) {
943: /* (strict lower triangular part of A)*x */
944: cval = ib[j] * 4;
945: z[cval] += v[0] * x1 + v[1] * x2 + v[2] * x3 + v[3] * x4;
946: z[cval + 1] += v[4] * x1 + v[5] * x2 + v[6] * x3 + v[7] * x4;
947: z[cval + 2] += v[8] * x1 + v[9] * x2 + v[10] * x3 + v[11] * x4;
948: z[cval + 3] += v[12] * x1 + v[13] * x2 + v[14] * x3 + v[15] * x4;
949: /* (strict upper triangular part of A)*x */
950: z[4 * i] += v[0] * x[cval] + v[4] * x[cval + 1] + v[8] * x[cval + 2] + v[12] * x[cval + 3];
951: z[4 * i + 1] += v[1] * x[cval] + v[5] * x[cval + 1] + v[9] * x[cval + 2] + v[13] * x[cval + 3];
952: z[4 * i + 2] += v[2] * x[cval] + v[6] * x[cval + 1] + v[10] * x[cval + 2] + v[14] * x[cval + 3];
953: z[4 * i + 3] += v[3] * x[cval] + v[7] * x[cval + 1] + v[11] * x[cval + 2] + v[15] * x[cval + 3];
954: v += 16;
955: }
956: }
958: PetscCall(VecRestoreArrayRead(xx, &x));
959: PetscCall(VecRestoreArray(zz, &z));
961: PetscCall(PetscLogFlops(32.0 * (a->nz * 2.0 - nonzerorow)));
962: PetscFunctionReturn(PETSC_SUCCESS);
963: }
965: PetscErrorCode MatMultAdd_SeqSBAIJ_5(Mat A, Vec xx, Vec yy, Vec zz)
966: {
967: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
968: PetscScalar *z, x1, x2, x3, x4, x5;
969: const PetscScalar *x, *xb;
970: const MatScalar *v;
971: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
972: const PetscInt *aj = a->j, *ai = a->i, *ib;
973: PetscInt nonzerorow = 0;
975: PetscFunctionBegin;
976: PetscCall(VecCopy(yy, zz));
977: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
978: PetscCall(VecGetArrayRead(xx, &x));
979: PetscCall(VecGetArray(zz, &z));
981: v = a->a;
982: xb = x;
984: for (i = 0; i < mbs; i++, xb += 5, ai++) {
985: n = ai[1] - ai[0]; /* length of i_th block row of A */
986: if (!n) continue;
987: x1 = xb[0];
988: x2 = xb[1];
989: x3 = xb[2];
990: x4 = xb[3];
991: x5 = xb[4];
992: ib = aj + *ai;
993: jmin = 0;
994: nonzerorow++;
995: if (*ib == i) { /* (diag of A)*x */
996: z[5 * i] += v[0] * x1 + v[5] * x2 + v[10] * x3 + v[15] * x4 + v[20] * x5;
997: z[5 * i + 1] += v[5] * x1 + v[6] * x2 + v[11] * x3 + v[16] * x4 + v[21] * x5;
998: z[5 * i + 2] += v[10] * x1 + v[11] * x2 + v[12] * x3 + v[17] * x4 + v[22] * x5;
999: z[5 * i + 3] += v[15] * x1 + v[16] * x2 + v[17] * x3 + v[18] * x4 + v[23] * x5;
1000: z[5 * i + 4] += v[20] * x1 + v[21] * x2 + v[22] * x3 + v[23] * x4 + v[24] * x5;
1001: v += 25;
1002: jmin++;
1003: }
1004: PetscPrefetchBlock(ib + jmin + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
1005: PetscPrefetchBlock(v + 25 * n, 25 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
1006: for (j = jmin; j < n; j++) {
1007: /* (strict lower triangular part of A)*x */
1008: cval = ib[j] * 5;
1009: z[cval] += v[0] * x1 + v[1] * x2 + v[2] * x3 + v[3] * x4 + v[4] * x5;
1010: z[cval + 1] += v[5] * x1 + v[6] * x2 + v[7] * x3 + v[8] * x4 + v[9] * x5;
1011: z[cval + 2] += v[10] * x1 + v[11] * x2 + v[12] * x3 + v[13] * x4 + v[14] * x5;
1012: z[cval + 3] += v[15] * x1 + v[16] * x2 + v[17] * x3 + v[18] * x4 + v[19] * x5;
1013: z[cval + 4] += v[20] * x1 + v[21] * x2 + v[22] * x3 + v[23] * x4 + v[24] * x5;
1014: /* (strict upper triangular part of A)*x */
1015: z[5 * i] += v[0] * x[cval] + v[5] * x[cval + 1] + v[10] * x[cval + 2] + v[15] * x[cval + 3] + v[20] * x[cval + 4];
1016: z[5 * i + 1] += v[1] * x[cval] + v[6] * x[cval + 1] + v[11] * x[cval + 2] + v[16] * x[cval + 3] + v[21] * x[cval + 4];
1017: z[5 * i + 2] += v[2] * x[cval] + v[7] * x[cval + 1] + v[12] * x[cval + 2] + v[17] * x[cval + 3] + v[22] * x[cval + 4];
1018: z[5 * i + 3] += v[3] * x[cval] + v[8] * x[cval + 1] + v[13] * x[cval + 2] + v[18] * x[cval + 3] + v[23] * x[cval + 4];
1019: z[5 * i + 4] += v[4] * x[cval] + v[9] * x[cval + 1] + v[14] * x[cval + 2] + v[19] * x[cval + 3] + v[24] * x[cval + 4];
1020: v += 25;
1021: }
1022: }
1024: PetscCall(VecRestoreArrayRead(xx, &x));
1025: PetscCall(VecRestoreArray(zz, &z));
1027: PetscCall(PetscLogFlops(50.0 * (a->nz * 2.0 - nonzerorow)));
1028: PetscFunctionReturn(PETSC_SUCCESS);
1029: }
1031: PetscErrorCode MatMultAdd_SeqSBAIJ_6(Mat A, Vec xx, Vec yy, Vec zz)
1032: {
1033: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1034: PetscScalar *z, x1, x2, x3, x4, x5, x6;
1035: const PetscScalar *x, *xb;
1036: const MatScalar *v;
1037: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
1038: const PetscInt *aj = a->j, *ai = a->i, *ib;
1039: PetscInt nonzerorow = 0;
1041: PetscFunctionBegin;
1042: PetscCall(VecCopy(yy, zz));
1043: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
1044: PetscCall(VecGetArrayRead(xx, &x));
1045: PetscCall(VecGetArray(zz, &z));
1047: v = a->a;
1048: xb = x;
1050: for (i = 0; i < mbs; i++, xb += 6, ai++) {
1051: n = ai[1] - ai[0]; /* length of i_th block row of A */
1052: if (!n) continue;
1053: x1 = xb[0];
1054: x2 = xb[1];
1055: x3 = xb[2];
1056: x4 = xb[3];
1057: x5 = xb[4];
1058: x6 = xb[5];
1059: ib = aj + *ai;
1060: jmin = 0;
1061: nonzerorow++;
1062: if (*ib == i) { /* (diag of A)*x */
1063: z[6 * i] += v[0] * x1 + v[6] * x2 + v[12] * x3 + v[18] * x4 + v[24] * x5 + v[30] * x6;
1064: z[6 * i + 1] += v[6] * x1 + v[7] * x2 + v[13] * x3 + v[19] * x4 + v[25] * x5 + v[31] * x6;
1065: z[6 * i + 2] += v[12] * x1 + v[13] * x2 + v[14] * x3 + v[20] * x4 + v[26] * x5 + v[32] * x6;
1066: z[6 * i + 3] += v[18] * x1 + v[19] * x2 + v[20] * x3 + v[21] * x4 + v[27] * x5 + v[33] * x6;
1067: z[6 * i + 4] += v[24] * x1 + v[25] * x2 + v[26] * x3 + v[27] * x4 + v[28] * x5 + v[34] * x6;
1068: z[6 * i + 5] += v[30] * x1 + v[31] * x2 + v[32] * x3 + v[33] * x4 + v[34] * x5 + v[35] * x6;
1069: v += 36;
1070: jmin++;
1071: }
1072: PetscPrefetchBlock(ib + jmin + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
1073: PetscPrefetchBlock(v + 36 * n, 36 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
1074: for (j = jmin; j < n; j++) {
1075: /* (strict lower triangular part of A)*x */
1076: cval = ib[j] * 6;
1077: z[cval] += v[0] * x1 + v[1] * x2 + v[2] * x3 + v[3] * x4 + v[4] * x5 + v[5] * x6;
1078: z[cval + 1] += v[6] * x1 + v[7] * x2 + v[8] * x3 + v[9] * x4 + v[10] * x5 + v[11] * x6;
1079: z[cval + 2] += v[12] * x1 + v[13] * x2 + v[14] * x3 + v[15] * x4 + v[16] * x5 + v[17] * x6;
1080: z[cval + 3] += v[18] * x1 + v[19] * x2 + v[20] * x3 + v[21] * x4 + v[22] * x5 + v[23] * x6;
1081: z[cval + 4] += v[24] * x1 + v[25] * x2 + v[26] * x3 + v[27] * x4 + v[28] * x5 + v[29] * x6;
1082: z[cval + 5] += v[30] * x1 + v[31] * x2 + v[32] * x3 + v[33] * x4 + v[34] * x5 + v[35] * x6;
1083: /* (strict upper triangular part of A)*x */
1084: z[6 * i] += v[0] * x[cval] + v[6] * x[cval + 1] + v[12] * x[cval + 2] + v[18] * x[cval + 3] + v[24] * x[cval + 4] + v[30] * x[cval + 5];
1085: z[6 * i + 1] += v[1] * x[cval] + v[7] * x[cval + 1] + v[13] * x[cval + 2] + v[19] * x[cval + 3] + v[25] * x[cval + 4] + v[31] * x[cval + 5];
1086: z[6 * i + 2] += v[2] * x[cval] + v[8] * x[cval + 1] + v[14] * x[cval + 2] + v[20] * x[cval + 3] + v[26] * x[cval + 4] + v[32] * x[cval + 5];
1087: z[6 * i + 3] += v[3] * x[cval] + v[9] * x[cval + 1] + v[15] * x[cval + 2] + v[21] * x[cval + 3] + v[27] * x[cval + 4] + v[33] * x[cval + 5];
1088: z[6 * i + 4] += v[4] * x[cval] + v[10] * x[cval + 1] + v[16] * x[cval + 2] + v[22] * x[cval + 3] + v[28] * x[cval + 4] + v[34] * x[cval + 5];
1089: z[6 * i + 5] += v[5] * x[cval] + v[11] * x[cval + 1] + v[17] * x[cval + 2] + v[23] * x[cval + 3] + v[29] * x[cval + 4] + v[35] * x[cval + 5];
1090: v += 36;
1091: }
1092: }
1094: PetscCall(VecRestoreArrayRead(xx, &x));
1095: PetscCall(VecRestoreArray(zz, &z));
1097: PetscCall(PetscLogFlops(72.0 * (a->nz * 2.0 - nonzerorow)));
1098: PetscFunctionReturn(PETSC_SUCCESS);
1099: }
1101: PetscErrorCode MatMultAdd_SeqSBAIJ_7(Mat A, Vec xx, Vec yy, Vec zz)
1102: {
1103: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1104: PetscScalar *z, x1, x2, x3, x4, x5, x6, x7;
1105: const PetscScalar *x, *xb;
1106: const MatScalar *v;
1107: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
1108: const PetscInt *aj = a->j, *ai = a->i, *ib;
1109: PetscInt nonzerorow = 0;
1111: PetscFunctionBegin;
1112: PetscCall(VecCopy(yy, zz));
1113: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
1114: PetscCall(VecGetArrayRead(xx, &x));
1115: PetscCall(VecGetArray(zz, &z));
1117: v = a->a;
1118: xb = x;
1120: for (i = 0; i < mbs; i++, xb += 7, ai++) {
1121: n = ai[1] - ai[0]; /* length of i_th block row of A */
1122: if (!n) continue;
1123: x1 = xb[0];
1124: x2 = xb[1];
1125: x3 = xb[2];
1126: x4 = xb[3];
1127: x5 = xb[4];
1128: x6 = xb[5];
1129: x7 = xb[6];
1130: ib = aj + *ai;
1131: jmin = 0;
1132: nonzerorow++;
1133: if (*ib == i) { /* (diag of A)*x */
1134: z[7 * i] += v[0] * x1 + v[7] * x2 + v[14] * x3 + v[21] * x4 + v[28] * x5 + v[35] * x6 + v[42] * x7;
1135: z[7 * i + 1] += v[7] * x1 + v[8] * x2 + v[15] * x3 + v[22] * x4 + v[29] * x5 + v[36] * x6 + v[43] * x7;
1136: z[7 * i + 2] += v[14] * x1 + v[15] * x2 + v[16] * x3 + v[23] * x4 + v[30] * x5 + v[37] * x6 + v[44] * x7;
1137: z[7 * i + 3] += v[21] * x1 + v[22] * x2 + v[23] * x3 + v[24] * x4 + v[31] * x5 + v[38] * x6 + v[45] * x7;
1138: z[7 * i + 4] += v[28] * x1 + v[29] * x2 + v[30] * x3 + v[31] * x4 + v[32] * x5 + v[39] * x6 + v[46] * x7;
1139: z[7 * i + 5] += v[35] * x1 + v[36] * x2 + v[37] * x3 + v[38] * x4 + v[39] * x5 + v[40] * x6 + v[47] * x7;
1140: z[7 * i + 6] += v[42] * x1 + v[43] * x2 + v[44] * x3 + v[45] * x4 + v[46] * x5 + v[47] * x6 + v[48] * x7;
1141: v += 49;
1142: jmin++;
1143: }
1144: PetscPrefetchBlock(ib + jmin + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
1145: PetscPrefetchBlock(v + 49 * n, 49 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
1146: for (j = jmin; j < n; j++) {
1147: /* (strict lower triangular part of A)*x */
1148: cval = ib[j] * 7;
1149: z[cval] += v[0] * x1 + v[1] * x2 + v[2] * x3 + v[3] * x4 + v[4] * x5 + v[5] * x6 + v[6] * x7;
1150: z[cval + 1] += v[7] * x1 + v[8] * x2 + v[9] * x3 + v[10] * x4 + v[11] * x5 + v[12] * x6 + v[13] * x7;
1151: z[cval + 2] += v[14] * x1 + v[15] * x2 + v[16] * x3 + v[17] * x4 + v[18] * x5 + v[19] * x6 + v[20] * x7;
1152: z[cval + 3] += v[21] * x1 + v[22] * x2 + v[23] * x3 + v[24] * x4 + v[25] * x5 + v[26] * x6 + v[27] * x7;
1153: z[cval + 4] += v[28] * x1 + v[29] * x2 + v[30] * x3 + v[31] * x4 + v[32] * x5 + v[33] * x6 + v[34] * x7;
1154: z[cval + 5] += v[35] * x1 + v[36] * x2 + v[37] * x3 + v[38] * x4 + v[39] * x5 + v[40] * x6 + v[41] * x7;
1155: z[cval + 6] += v[42] * x1 + v[43] * x2 + v[44] * x3 + v[45] * x4 + v[46] * x5 + v[47] * x6 + v[48] * x7;
1156: /* (strict upper triangular part of A)*x */
1157: z[7 * i] += v[0] * x[cval] + v[7] * x[cval + 1] + v[14] * x[cval + 2] + v[21] * x[cval + 3] + v[28] * x[cval + 4] + v[35] * x[cval + 5] + v[42] * x[cval + 6];
1158: z[7 * i + 1] += v[1] * x[cval] + v[8] * x[cval + 1] + v[15] * x[cval + 2] + v[22] * x[cval + 3] + v[29] * x[cval + 4] + v[36] * x[cval + 5] + v[43] * x[cval + 6];
1159: z[7 * i + 2] += v[2] * x[cval] + v[9] * x[cval + 1] + v[16] * x[cval + 2] + v[23] * x[cval + 3] + v[30] * x[cval + 4] + v[37] * x[cval + 5] + v[44] * x[cval + 6];
1160: z[7 * i + 3] += v[3] * x[cval] + v[10] * x[cval + 1] + v[17] * x[cval + 2] + v[24] * x[cval + 3] + v[31] * x[cval + 4] + v[38] * x[cval + 5] + v[45] * x[cval + 6];
1161: z[7 * i + 4] += v[4] * x[cval] + v[11] * x[cval + 1] + v[18] * x[cval + 2] + v[25] * x[cval + 3] + v[32] * x[cval + 4] + v[39] * x[cval + 5] + v[46] * x[cval + 6];
1162: z[7 * i + 5] += v[5] * x[cval] + v[12] * x[cval + 1] + v[19] * x[cval + 2] + v[26] * x[cval + 3] + v[33] * x[cval + 4] + v[40] * x[cval + 5] + v[47] * x[cval + 6];
1163: z[7 * i + 6] += v[6] * x[cval] + v[13] * x[cval + 1] + v[20] * x[cval + 2] + v[27] * x[cval + 3] + v[34] * x[cval + 4] + v[41] * x[cval + 5] + v[48] * x[cval + 6];
1164: v += 49;
1165: }
1166: }
1168: PetscCall(VecRestoreArrayRead(xx, &x));
1169: PetscCall(VecRestoreArray(zz, &z));
1171: PetscCall(PetscLogFlops(98.0 * (a->nz * 2.0 - nonzerorow)));
1172: PetscFunctionReturn(PETSC_SUCCESS);
1173: }
1175: PetscErrorCode MatMultAdd_SeqSBAIJ_N(Mat A, Vec xx, Vec yy, Vec zz)
1176: {
1177: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1178: PetscScalar *z, *z_ptr = NULL, *zb, *work, *workt;
1179: const PetscScalar *x, *x_ptr, *xb;
1180: const MatScalar *v;
1181: PetscInt mbs = a->mbs, i, bs = A->rmap->bs, j, n, bs2 = a->bs2, ncols, k;
1182: const PetscInt *idx, *aj, *ii;
1183: PetscInt nonzerorow = 0;
1185: PetscFunctionBegin;
1186: PetscCall(VecCopy(yy, zz));
1187: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
1188: PetscCall(VecGetArrayRead(xx, &x));
1189: x_ptr = x;
1190: PetscCall(VecGetArray(zz, &z));
1191: z_ptr = z;
1193: aj = a->j;
1194: v = a->a;
1195: ii = a->i;
1197: if (!a->mult_work) PetscCall(PetscMalloc1(A->rmap->n + 1, &a->mult_work));
1198: work = a->mult_work;
1200: for (i = 0; i < mbs; i++) {
1201: n = ii[1] - ii[0];
1202: ncols = n * bs;
1203: workt = work;
1204: idx = aj + ii[0];
1205: nonzerorow += (n > 0);
1207: /* upper triangular part */
1208: for (j = 0; j < n; j++) {
1209: xb = x_ptr + bs * (*idx++);
1210: for (k = 0; k < bs; k++) workt[k] = xb[k];
1211: workt += bs;
1212: }
1213: /* z(i*bs:(i+1)*bs-1) += A(i,:)*x */
1214: PetscKernel_w_gets_w_plus_Ar_times_v(bs, ncols, work, v, z);
1216: /* strict lower triangular part */
1217: idx = aj + ii[0];
1218: if (n && *idx == i) {
1219: ncols -= bs;
1220: v += bs2;
1221: idx++;
1222: n--;
1223: }
1224: if (ncols > 0) {
1225: workt = work;
1226: PetscCall(PetscArrayzero(workt, ncols));
1227: PetscKernel_w_gets_w_plus_trans_Ar_times_v(bs, ncols, x, v, workt);
1228: for (j = 0; j < n; j++) {
1229: zb = z_ptr + bs * (*idx++);
1230: for (k = 0; k < bs; k++) zb[k] += workt[k];
1231: workt += bs;
1232: }
1233: }
1235: x += bs;
1236: v += n * bs2;
1237: z += bs;
1238: ii++;
1239: }
1241: PetscCall(VecRestoreArrayRead(xx, &x));
1242: PetscCall(VecRestoreArray(zz, &z));
1244: PetscCall(PetscLogFlops(2.0 * (a->nz * 2.0 - nonzerorow)));
1245: PetscFunctionReturn(PETSC_SUCCESS);
1246: }
1248: PetscErrorCode MatScale_SeqSBAIJ(Mat inA, PetscScalar alpha)
1249: {
1250: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)inA->data;
1251: PetscScalar oalpha = alpha;
1252: PetscBLASInt one = 1, totalnz;
1254: PetscFunctionBegin;
1255: PetscCall(PetscBLASIntCast(a->bs2 * a->nz, &totalnz));
1256: PetscCallBLAS("BLASscal", BLASscal_(&totalnz, &oalpha, a->a, &one));
1257: PetscCall(PetscLogFlops(totalnz));
1258: PetscFunctionReturn(PETSC_SUCCESS);
1259: }
1261: PetscErrorCode MatNorm_SeqSBAIJ(Mat A, NormType type, PetscReal *norm)
1262: {
1263: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1264: const MatScalar *v = a->a;
1265: PetscReal sum_diag = 0.0, sum_off = 0.0, *sum;
1266: PetscInt i, j, k, bs = A->rmap->bs, bs2 = a->bs2, k1, mbs = a->mbs, jmin, jmax, nexti, ik, *jl, *il;
1267: const PetscInt *aj = a->j, *col;
1269: PetscFunctionBegin;
1270: if (!a->nz) {
1271: *norm = 0.0;
1272: PetscFunctionReturn(PETSC_SUCCESS);
1273: }
1274: if (type == NORM_FROBENIUS) {
1275: for (k = 0; k < mbs; k++) {
1276: jmin = a->i[k];
1277: jmax = a->i[k + 1];
1278: col = aj + jmin;
1279: if (jmax - jmin > 0 && *col == k) { /* diagonal block */
1280: for (i = 0; i < bs2; i++) {
1281: sum_diag += PetscRealPart(PetscConj(*v) * (*v));
1282: v++;
1283: }
1284: jmin++;
1285: }
1286: for (j = jmin; j < jmax; j++) { /* off-diagonal blocks */
1287: for (i = 0; i < bs2; i++) {
1288: sum_off += PetscRealPart(PetscConj(*v) * (*v));
1289: v++;
1290: }
1291: }
1292: }
1293: *norm = PetscSqrtReal(sum_diag + 2 * sum_off);
1294: PetscCall(PetscLogFlops(2.0 * bs2 * a->nz));
1295: } else if (type == NORM_INFINITY || type == NORM_1) { /* maximum row/column sum */
1296: PetscCall(PetscMalloc3(bs, &sum, mbs, &il, mbs, &jl));
1297: for (i = 0; i < mbs; i++) jl[i] = mbs;
1298: il[0] = 0;
1300: *norm = 0.0;
1301: for (k = 0; k < mbs; k++) { /* k_th block row */
1302: for (j = 0; j < bs; j++) sum[j] = 0.0;
1303: /*-- col sum --*/
1304: i = jl[k]; /* first |A(i,k)| to be added */
1305: /* jl[k]=i: first nonzero element in row i for submatrix A(1:k,k:n) (active window)
1306: at step k */
1307: while (i < mbs) {
1308: nexti = jl[i]; /* next block row to be added */
1309: ik = il[i]; /* block index of A(i,k) in the array a */
1310: for (j = 0; j < bs; j++) {
1311: v = a->a + ik * bs2 + j * bs;
1312: for (k1 = 0; k1 < bs; k1++) {
1313: sum[j] += PetscAbsScalar(*v);
1314: v++;
1315: }
1316: }
1317: /* update il, jl */
1318: jmin = ik + 1; /* block index of array a: points to the next nonzero of A in row i */
1319: jmax = a->i[i + 1];
1320: if (jmin < jmax) {
1321: il[i] = jmin;
1322: j = a->j[jmin];
1323: jl[i] = jl[j];
1324: jl[j] = i;
1325: }
1326: i = nexti;
1327: }
1328: /*-- row sum --*/
1329: jmin = a->i[k];
1330: jmax = a->i[k + 1];
1331: for (i = jmin; i < jmax; i++) {
1332: for (j = 0; j < bs; j++) {
1333: v = a->a + i * bs2 + j;
1334: for (k1 = 0; k1 < bs; k1++) {
1335: sum[j] += PetscAbsScalar(*v);
1336: v += bs;
1337: }
1338: }
1339: }
1340: /* add k_th block row to il, jl */
1341: col = aj + jmin;
1342: if (jmax - jmin > 0 && *col == k) jmin++;
1343: if (jmin < jmax) {
1344: il[k] = jmin;
1345: j = a->j[jmin];
1346: jl[k] = jl[j];
1347: jl[j] = k;
1348: }
1349: for (j = 0; j < bs; j++) {
1350: if (sum[j] > *norm) *norm = sum[j];
1351: }
1352: }
1353: PetscCall(PetscFree3(sum, il, jl));
1354: PetscCall(PetscLogFlops(PetscMax(mbs * a->nz - 1, 0)));
1355: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for this norm yet");
1356: PetscFunctionReturn(PETSC_SUCCESS);
1357: }
1359: PetscErrorCode MatEqual_SeqSBAIJ(Mat A, Mat B, PetscBool *flg)
1360: {
1361: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data, *b = (Mat_SeqSBAIJ *)B->data;
1363: PetscFunctionBegin;
1364: /* If the matrix/block dimensions are not equal, or no of nonzeros or shift */
1365: if ((A->rmap->N != B->rmap->N) || (A->cmap->n != B->cmap->n) || (A->rmap->bs != B->rmap->bs) || (a->nz != b->nz)) {
1366: *flg = PETSC_FALSE;
1367: PetscFunctionReturn(PETSC_SUCCESS);
1368: }
1370: /* if the a->i are the same */
1371: PetscCall(PetscArraycmp(a->i, b->i, a->mbs + 1, flg));
1372: if (!*flg) PetscFunctionReturn(PETSC_SUCCESS);
1374: /* if a->j are the same */
1375: PetscCall(PetscArraycmp(a->j, b->j, a->nz, flg));
1376: if (!*flg) PetscFunctionReturn(PETSC_SUCCESS);
1378: /* if a->a are the same */
1379: PetscCall(PetscArraycmp(a->a, b->a, a->nz * A->rmap->bs * A->rmap->bs, flg));
1380: PetscFunctionReturn(PETSC_SUCCESS);
1381: }
1383: PetscErrorCode MatGetDiagonal_SeqSBAIJ(Mat A, Vec v)
1384: {
1385: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1386: PetscInt n;
1387: const PetscInt bs = A->rmap->bs, ambs = a->mbs, bs2 = a->bs2;
1388: PetscScalar *x;
1389: const MatScalar *aa = a->a, *aa_j;
1390: const PetscInt *ai = a->i, *adiag;
1391: PetscBool diagDense;
1393: PetscFunctionBegin;
1394: PetscCheck(!A->factortype || bs <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix with bs>1");
1395: PetscCall(MatGetDiagonalMarkers_SeqSBAIJ(A, &adiag, &diagDense));
1396: if (A->factortype == MAT_FACTOR_CHOLESKY || A->factortype == MAT_FACTOR_ICC) {
1397: PetscCall(VecGetArrayWrite(v, &x));
1398: for (PetscInt i = 0; i < ambs; i++) x[i] = 1.0 / aa[adiag[i]];
1399: PetscCall(VecRestoreArrayWrite(v, &x));
1400: PetscFunctionReturn(PETSC_SUCCESS);
1401: }
1403: PetscCall(VecGetLocalSize(v, &n));
1404: PetscCheck(n == A->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector");
1405: PetscCall(VecGetArrayWrite(v, &x));
1407: if (diagDense) {
1408: for (PetscInt i = 0, row = 0; i < ambs; i++) {
1409: aa_j = aa + adiag[i] * bs2;
1410: for (PetscInt k = 0; k < bs2; k += (bs + 1)) x[row++] = aa_j[k];
1411: }
1412: } else {
1413: for (PetscInt i = 0, row = 0; i < ambs; i++) {
1414: const PetscInt j = adiag[i];
1416: if (j != ai[i + 1]) {
1417: aa_j = aa + j * bs2;
1418: for (PetscInt k = 0; k < bs2; k += (bs + 1)) x[row++] = aa_j[k];
1419: } else {
1420: for (PetscInt k = 0; k < bs; k++) x[row++] = 0.0;
1421: }
1422: }
1423: }
1424: PetscCall(VecRestoreArrayWrite(v, &x));
1425: PetscFunctionReturn(PETSC_SUCCESS);
1426: }
1428: PetscErrorCode MatDiagonalScale_SeqSBAIJ(Mat A, Vec ll, Vec rr)
1429: {
1430: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1431: PetscScalar x;
1432: const PetscScalar *l, *li, *ri;
1433: MatScalar *aa, *v;
1434: PetscInt i, j, k, lm, M, m, mbs, tmp, bs, bs2;
1435: const PetscInt *ai, *aj;
1437: PetscFunctionBegin;
1438: if (!ll) PetscFunctionReturn(PETSC_SUCCESS);
1439: ai = a->i;
1440: aj = a->j;
1441: aa = a->a;
1442: m = A->rmap->N;
1443: bs = A->rmap->bs;
1444: mbs = a->mbs;
1445: bs2 = a->bs2;
1447: PetscCall(VecGetArrayRead(ll, &l));
1448: PetscCall(VecGetLocalSize(ll, &lm));
1449: PetscCheck(lm == m, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Left scaling vector wrong length");
1450: for (i = 0; i < mbs; i++) { /* for each block row */
1451: M = ai[i + 1] - ai[i];
1452: li = l + i * bs;
1453: v = aa + bs2 * ai[i];
1454: for (j = 0; j < M; j++) { /* for each block */
1455: ri = l + bs * aj[ai[i] + j];
1456: for (k = 0; k < bs; k++) {
1457: x = ri[k];
1458: for (tmp = 0; tmp < bs; tmp++) (*v++) *= li[tmp] * x;
1459: }
1460: }
1461: }
1462: PetscCall(VecRestoreArrayRead(ll, &l));
1463: PetscCall(PetscLogFlops(2.0 * a->nz));
1464: PetscFunctionReturn(PETSC_SUCCESS);
1465: }
1467: PetscErrorCode MatGetInfo_SeqSBAIJ(Mat A, MatInfoType flag, MatInfo *info)
1468: {
1469: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1471: PetscFunctionBegin;
1472: info->block_size = a->bs2;
1473: info->nz_allocated = a->bs2 * a->maxnz; /*num. of nonzeros in upper triangular part */
1474: info->nz_used = a->bs2 * a->nz; /*num. of nonzeros in upper triangular part */
1475: info->nz_unneeded = info->nz_allocated - info->nz_used;
1476: info->assemblies = A->num_ass;
1477: info->mallocs = A->info.mallocs;
1478: info->memory = 0; /* REVIEW ME */
1479: if (A->factortype) {
1480: info->fill_ratio_given = A->info.fill_ratio_given;
1481: info->fill_ratio_needed = A->info.fill_ratio_needed;
1482: info->factor_mallocs = A->info.factor_mallocs;
1483: } else {
1484: info->fill_ratio_given = 0;
1485: info->fill_ratio_needed = 0;
1486: info->factor_mallocs = 0;
1487: }
1488: PetscFunctionReturn(PETSC_SUCCESS);
1489: }
1491: PetscErrorCode MatZeroEntries_SeqSBAIJ(Mat A)
1492: {
1493: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1495: PetscFunctionBegin;
1496: PetscCall(PetscArrayzero(a->a, a->bs2 * a->i[a->mbs]));
1497: PetscFunctionReturn(PETSC_SUCCESS);
1498: }
1500: PetscErrorCode MatGetRowMaxAbs_SeqSBAIJ(Mat A, Vec v, PetscInt idx[])
1501: {
1502: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1503: PetscInt i, j, n, row, col, bs, mbs;
1504: const PetscInt *ai, *aj;
1505: PetscReal atmp;
1506: const MatScalar *aa;
1507: PetscScalar *x;
1508: PetscInt ncols, brow, bcol, krow, kcol;
1510: PetscFunctionBegin;
1511: PetscCheck(!idx, PETSC_COMM_SELF, PETSC_ERR_SUP, "Send email to petsc-maint@mcs.anl.gov");
1512: PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
1513: bs = A->rmap->bs;
1514: aa = a->a;
1515: ai = a->i;
1516: aj = a->j;
1517: mbs = a->mbs;
1519: PetscCall(VecSet(v, 0.0));
1520: PetscCall(VecGetArray(v, &x));
1521: PetscCall(VecGetLocalSize(v, &n));
1522: PetscCheck(n == A->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector");
1523: for (i = 0; i < mbs; i++) {
1524: ncols = ai[1] - ai[0];
1525: ai++;
1526: brow = bs * i;
1527: for (j = 0; j < ncols; j++) {
1528: bcol = bs * (*aj);
1529: for (kcol = 0; kcol < bs; kcol++) {
1530: col = bcol + kcol; /* col index */
1531: for (krow = 0; krow < bs; krow++) {
1532: atmp = PetscAbsScalar(*aa);
1533: aa++;
1534: row = brow + krow; /* row index */
1535: if (PetscRealPart(x[row]) < atmp) x[row] = atmp;
1536: if (*aj > i && PetscRealPart(x[col]) < atmp) x[col] = atmp;
1537: }
1538: }
1539: aj++;
1540: }
1541: }
1542: PetscCall(VecRestoreArray(v, &x));
1543: PetscFunctionReturn(PETSC_SUCCESS);
1544: }
1546: PetscErrorCode MatMatMultSymbolic_SeqSBAIJ_SeqDense(Mat A, Mat B, PetscReal fill, Mat C)
1547: {
1548: PetscFunctionBegin;
1549: PetscCall(MatMatMultSymbolic_SeqDense_SeqDense(A, B, 0.0, C));
1550: C->ops->matmultnumeric = MatMatMultNumeric_SeqSBAIJ_SeqDense;
1551: PetscFunctionReturn(PETSC_SUCCESS);
1552: }
1554: static PetscErrorCode MatMatMult_SeqSBAIJ_1_Private(Mat A, PetscScalar *b, PetscInt bm, PetscScalar *c, PetscInt cm, PetscInt cn)
1555: {
1556: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1557: PetscScalar *z = c;
1558: const PetscScalar *xb;
1559: PetscScalar x1;
1560: const MatScalar *v = a->a, *vv;
1561: PetscInt mbs = a->mbs, i, *idx = a->j, *ii = a->i, j, *jj, n, k;
1562: const int aconj = PetscDefined(USE_COMPLEX) && A->hermitian == PETSC_BOOL3_TRUE ? 1 : 0;
1564: PetscFunctionBegin;
1565: for (i = 0; i < mbs; i++) {
1566: n = ii[1] - ii[0];
1567: ii++;
1568: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
1569: PetscPrefetchBlock(v + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
1570: jj = idx;
1571: vv = v;
1572: for (k = 0; k < cn; k++) {
1573: idx = jj;
1574: v = vv;
1575: for (j = 0; j < n; j++) {
1576: xb = b + (*idx);
1577: x1 = xb[0 + k * bm];
1578: z[0 + k * cm] += v[0] * x1;
1579: if (*idx != i) c[(*idx) + k * cm] += (aconj ? PetscConj(v[0]) : v[0]) * b[i + k * bm];
1580: v += 1;
1581: ++idx;
1582: }
1583: }
1584: z += 1;
1585: }
1586: PetscFunctionReturn(PETSC_SUCCESS);
1587: }
1589: static PetscErrorCode MatMatMult_SeqSBAIJ_2_Private(Mat A, PetscScalar *b, PetscInt bm, PetscScalar *c, PetscInt cm, PetscInt cn)
1590: {
1591: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1592: PetscScalar *z = c;
1593: const PetscScalar *xb;
1594: PetscScalar x1, x2;
1595: const MatScalar *v = a->a, *vv;
1596: PetscInt mbs = a->mbs, i, *idx = a->j, *ii = a->i, j, *jj, n, k;
1598: PetscFunctionBegin;
1599: for (i = 0; i < mbs; i++) {
1600: n = ii[1] - ii[0];
1601: ii++;
1602: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
1603: PetscPrefetchBlock(v + 4 * n, 4 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
1604: jj = idx;
1605: vv = v;
1606: for (k = 0; k < cn; k++) {
1607: idx = jj;
1608: v = vv;
1609: for (j = 0; j < n; j++) {
1610: xb = b + 2 * (*idx);
1611: x1 = xb[0 + k * bm];
1612: x2 = xb[1 + k * bm];
1613: z[0 + k * cm] += v[0] * x1 + v[2] * x2;
1614: z[1 + k * cm] += v[1] * x1 + v[3] * x2;
1615: if (*idx != i) {
1616: c[2 * (*idx) + 0 + k * cm] += v[0] * b[2 * i + k * bm] + v[1] * b[2 * i + 1 + k * bm];
1617: c[2 * (*idx) + 1 + k * cm] += v[2] * b[2 * i + k * bm] + v[3] * b[2 * i + 1 + k * bm];
1618: }
1619: v += 4;
1620: ++idx;
1621: }
1622: }
1623: z += 2;
1624: }
1625: PetscFunctionReturn(PETSC_SUCCESS);
1626: }
1628: static PetscErrorCode MatMatMult_SeqSBAIJ_3_Private(Mat A, PetscScalar *b, PetscInt bm, PetscScalar *c, PetscInt cm, PetscInt cn)
1629: {
1630: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1631: PetscScalar *z = c;
1632: const PetscScalar *xb;
1633: PetscScalar x1, x2, x3;
1634: const MatScalar *v = a->a, *vv;
1635: PetscInt mbs = a->mbs, i, *idx = a->j, *ii = a->i, j, *jj, n, k;
1637: PetscFunctionBegin;
1638: for (i = 0; i < mbs; i++) {
1639: n = ii[1] - ii[0];
1640: ii++;
1641: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
1642: PetscPrefetchBlock(v + 9 * n, 9 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
1643: jj = idx;
1644: vv = v;
1645: for (k = 0; k < cn; k++) {
1646: idx = jj;
1647: v = vv;
1648: for (j = 0; j < n; j++) {
1649: xb = b + 3 * (*idx);
1650: x1 = xb[0 + k * bm];
1651: x2 = xb[1 + k * bm];
1652: x3 = xb[2 + k * bm];
1653: z[0 + k * cm] += v[0] * x1 + v[3] * x2 + v[6] * x3;
1654: z[1 + k * cm] += v[1] * x1 + v[4] * x2 + v[7] * x3;
1655: z[2 + k * cm] += v[2] * x1 + v[5] * x2 + v[8] * x3;
1656: if (*idx != i) {
1657: c[3 * (*idx) + 0 + k * cm] += v[0] * b[3 * i + k * bm] + v[3] * b[3 * i + 1 + k * bm] + v[6] * b[3 * i + 2 + k * bm];
1658: c[3 * (*idx) + 1 + k * cm] += v[1] * b[3 * i + k * bm] + v[4] * b[3 * i + 1 + k * bm] + v[7] * b[3 * i + 2 + k * bm];
1659: c[3 * (*idx) + 2 + k * cm] += v[2] * b[3 * i + k * bm] + v[5] * b[3 * i + 1 + k * bm] + v[8] * b[3 * i + 2 + k * bm];
1660: }
1661: v += 9;
1662: ++idx;
1663: }
1664: }
1665: z += 3;
1666: }
1667: PetscFunctionReturn(PETSC_SUCCESS);
1668: }
1670: static PetscErrorCode MatMatMult_SeqSBAIJ_4_Private(Mat A, PetscScalar *b, PetscInt bm, PetscScalar *c, PetscInt cm, PetscInt cn)
1671: {
1672: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1673: PetscScalar *z = c;
1674: const PetscScalar *xb;
1675: PetscScalar x1, x2, x3, x4;
1676: const MatScalar *v = a->a, *vv;
1677: PetscInt mbs = a->mbs, i, *idx = a->j, *ii = a->i, j, *jj, n, k;
1679: PetscFunctionBegin;
1680: for (i = 0; i < mbs; i++) {
1681: n = ii[1] - ii[0];
1682: ii++;
1683: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
1684: PetscPrefetchBlock(v + 16 * n, 16 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
1685: jj = idx;
1686: vv = v;
1687: for (k = 0; k < cn; k++) {
1688: idx = jj;
1689: v = vv;
1690: for (j = 0; j < n; j++) {
1691: xb = b + 4 * (*idx);
1692: x1 = xb[0 + k * bm];
1693: x2 = xb[1 + k * bm];
1694: x3 = xb[2 + k * bm];
1695: x4 = xb[3 + k * bm];
1696: z[0 + k * cm] += v[0] * x1 + v[4] * x2 + v[8] * x3 + v[12] * x4;
1697: z[1 + k * cm] += v[1] * x1 + v[5] * x2 + v[9] * x3 + v[13] * x4;
1698: z[2 + k * cm] += v[2] * x1 + v[6] * x2 + v[10] * x3 + v[14] * x4;
1699: z[3 + k * cm] += v[3] * x1 + v[7] * x2 + v[11] * x3 + v[15] * x4;
1700: if (*idx != i) {
1701: c[4 * (*idx) + 0 + k * cm] += v[0] * b[4 * i + k * bm] + v[4] * b[4 * i + 1 + k * bm] + v[8] * b[4 * i + 2 + k * bm] + v[12] * b[4 * i + 3 + k * bm];
1702: c[4 * (*idx) + 1 + k * cm] += v[1] * b[4 * i + k * bm] + v[5] * b[4 * i + 1 + k * bm] + v[9] * b[4 * i + 2 + k * bm] + v[13] * b[4 * i + 3 + k * bm];
1703: c[4 * (*idx) + 2 + k * cm] += v[2] * b[4 * i + k * bm] + v[6] * b[4 * i + 1 + k * bm] + v[10] * b[4 * i + 2 + k * bm] + v[14] * b[4 * i + 3 + k * bm];
1704: c[4 * (*idx) + 3 + k * cm] += v[3] * b[4 * i + k * bm] + v[7] * b[4 * i + 1 + k * bm] + v[11] * b[4 * i + 2 + k * bm] + v[15] * b[4 * i + 3 + k * bm];
1705: }
1706: v += 16;
1707: ++idx;
1708: }
1709: }
1710: z += 4;
1711: }
1712: PetscFunctionReturn(PETSC_SUCCESS);
1713: }
1715: static PetscErrorCode MatMatMult_SeqSBAIJ_5_Private(Mat A, PetscScalar *b, PetscInt bm, PetscScalar *c, PetscInt cm, PetscInt cn)
1716: {
1717: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1718: PetscScalar *z = c;
1719: const PetscScalar *xb;
1720: PetscScalar x1, x2, x3, x4, x5;
1721: const MatScalar *v = a->a, *vv;
1722: PetscInt mbs = a->mbs, i, *idx = a->j, *ii = a->i, j, *jj, n, k;
1724: PetscFunctionBegin;
1725: for (i = 0; i < mbs; i++) {
1726: n = ii[1] - ii[0];
1727: ii++;
1728: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
1729: PetscPrefetchBlock(v + 25 * n, 25 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
1730: jj = idx;
1731: vv = v;
1732: for (k = 0; k < cn; k++) {
1733: idx = jj;
1734: v = vv;
1735: for (j = 0; j < n; j++) {
1736: xb = b + 5 * (*idx);
1737: x1 = xb[0 + k * bm];
1738: x2 = xb[1 + k * bm];
1739: x3 = xb[2 + k * bm];
1740: x4 = xb[3 + k * bm];
1741: x5 = xb[4 + k * cm];
1742: z[0 + k * cm] += v[0] * x1 + v[5] * x2 + v[10] * x3 + v[15] * x4 + v[20] * x5;
1743: z[1 + k * cm] += v[1] * x1 + v[6] * x2 + v[11] * x3 + v[16] * x4 + v[21] * x5;
1744: z[2 + k * cm] += v[2] * x1 + v[7] * x2 + v[12] * x3 + v[17] * x4 + v[22] * x5;
1745: z[3 + k * cm] += v[3] * x1 + v[8] * x2 + v[13] * x3 + v[18] * x4 + v[23] * x5;
1746: z[4 + k * cm] += v[4] * x1 + v[9] * x2 + v[14] * x3 + v[19] * x4 + v[24] * x5;
1747: if (*idx != i) {
1748: c[5 * (*idx) + 0 + k * cm] += v[0] * b[5 * i + k * bm] + v[5] * b[5 * i + 1 + k * bm] + v[10] * b[5 * i + 2 + k * bm] + v[15] * b[5 * i + 3 + k * bm] + v[20] * b[5 * i + 4 + k * bm];
1749: c[5 * (*idx) + 1 + k * cm] += v[1] * b[5 * i + k * bm] + v[6] * b[5 * i + 1 + k * bm] + v[11] * b[5 * i + 2 + k * bm] + v[16] * b[5 * i + 3 + k * bm] + v[21] * b[5 * i + 4 + k * bm];
1750: c[5 * (*idx) + 2 + k * cm] += v[2] * b[5 * i + k * bm] + v[7] * b[5 * i + 1 + k * bm] + v[12] * b[5 * i + 2 + k * bm] + v[17] * b[5 * i + 3 + k * bm] + v[22] * b[5 * i + 4 + k * bm];
1751: c[5 * (*idx) + 3 + k * cm] += v[3] * b[5 * i + k * bm] + v[8] * b[5 * i + 1 + k * bm] + v[13] * b[5 * i + 2 + k * bm] + v[18] * b[5 * i + 3 + k * bm] + v[23] * b[5 * i + 4 + k * bm];
1752: c[5 * (*idx) + 4 + k * cm] += v[4] * b[5 * i + k * bm] + v[9] * b[5 * i + 1 + k * bm] + v[14] * b[5 * i + 2 + k * bm] + v[19] * b[5 * i + 3 + k * bm] + v[24] * b[5 * i + 4 + k * bm];
1753: }
1754: v += 25;
1755: ++idx;
1756: }
1757: }
1758: z += 5;
1759: }
1760: PetscFunctionReturn(PETSC_SUCCESS);
1761: }
1763: PetscErrorCode MatMatMultNumeric_SeqSBAIJ_SeqDense(Mat A, Mat B, Mat C)
1764: {
1765: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1766: Mat_SeqDense *bd = (Mat_SeqDense *)B->data;
1767: Mat_SeqDense *cd = (Mat_SeqDense *)C->data;
1768: PetscInt cm = cd->lda, cn = B->cmap->n, bm = bd->lda;
1769: PetscInt mbs, i, bs = A->rmap->bs, j, n, bs2 = a->bs2;
1770: PetscBLASInt bbs, bcn, bbm, bcm;
1771: PetscScalar *z = NULL;
1772: PetscScalar *c, *b;
1773: const MatScalar *v;
1774: const PetscInt *idx, *ii;
1775: PetscScalar _DOne = 1.0;
1777: PetscFunctionBegin;
1778: if (!cm || !cn) PetscFunctionReturn(PETSC_SUCCESS);
1779: 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);
1780: 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);
1781: 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);
1782: b = bd->v;
1783: PetscCall(MatZeroEntries(C));
1784: PetscCall(MatDenseGetArray(C, &c));
1785: switch (bs) {
1786: case 1:
1787: PetscCall(MatMatMult_SeqSBAIJ_1_Private(A, b, bm, c, cm, cn));
1788: break;
1789: case 2:
1790: PetscCall(MatMatMult_SeqSBAIJ_2_Private(A, b, bm, c, cm, cn));
1791: break;
1792: case 3:
1793: PetscCall(MatMatMult_SeqSBAIJ_3_Private(A, b, bm, c, cm, cn));
1794: break;
1795: case 4:
1796: PetscCall(MatMatMult_SeqSBAIJ_4_Private(A, b, bm, c, cm, cn));
1797: break;
1798: case 5:
1799: PetscCall(MatMatMult_SeqSBAIJ_5_Private(A, b, bm, c, cm, cn));
1800: break;
1801: default: /* block sizes larger than 5 by 5 are handled by BLAS */
1802: PetscCall(PetscBLASIntCast(bs, &bbs));
1803: PetscCall(PetscBLASIntCast(cn, &bcn));
1804: PetscCall(PetscBLASIntCast(bm, &bbm));
1805: PetscCall(PetscBLASIntCast(cm, &bcm));
1806: idx = a->j;
1807: v = a->a;
1808: mbs = a->mbs;
1809: ii = a->i;
1810: z = c;
1811: for (i = 0; i < mbs; i++) {
1812: n = ii[1] - ii[0];
1813: ii++;
1814: for (j = 0; j < n; j++) {
1815: if (*idx != i) PetscCallBLAS("BLASgemm", BLASgemm_("T", "N", &bbs, &bcn, &bbs, &_DOne, v, &bbs, b + bs * i, &bbm, &_DOne, c + bs * (*idx), &bcm));
1816: PetscCallBLAS("BLASgemm", BLASgemm_("N", "N", &bbs, &bcn, &bbs, &_DOne, v, &bbs, b + bs * (*idx++), &bbm, &_DOne, z, &bcm));
1817: v += bs2;
1818: }
1819: z += bs;
1820: }
1821: }
1822: PetscCall(MatDenseRestoreArray(C, &c));
1823: PetscCall(PetscLogFlops((2.0 * (a->nz * 2.0 - a->nonzerorowcnt) * bs2 - a->nonzerorowcnt) * cn));
1824: PetscFunctionReturn(PETSC_SUCCESS);
1825: }