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], D;
230: IS is1, is2, intersect = NULL, sorted = NULL, perm = NULL, iperm = NULL, expanded = NULL;
231: PetscInt n1, n2, ni;
232: PetscBool implicit, sym, sameorder = PETSC_FALSE, issorted = PETSC_FALSE;
234: PetscFunctionBegin;
235: implicit = sym = (PetscBool)(A->rmap->N == A->cmap->N && (A->symmetric == PETSC_BOOL3_TRUE || A->hermitian == PETSC_BOOL3_TRUE));
236: PetscCall(ISCompressIndicesGeneral(A->rmap->N, A->rmap->n, A->rmap->bs, 1, &isrow, &is1));
237: if (isrow == iscol) {
238: is2 = is1;
239: PetscCall(PetscObjectReference((PetscObject)is2));
240: } else {
241: PetscCall(ISCompressIndicesGeneral(A->cmap->N, A->cmap->n, A->cmap->bs, 1, &iscol, &is2));
242: if (implicit == PETSC_TRUE) {
243: PetscCall(ISIntersect(is1, is2, &intersect));
244: PetscCall(ISGetLocalSize(intersect, &ni));
245: PetscCall(ISDestroy(&intersect));
246: if (ni == 0) sym = PETSC_FALSE;
247: else if (PetscDefined(USE_DEBUG)) {
248: PetscCall(ISGetLocalSize(is1, &n1));
249: PetscCall(ISGetLocalSize(is2, &n2));
250: PetscCheck(ni == n1 && ni == n2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot create such a submatrix");
251: }
252: }
253: }
254: // rectangular and nonsymmetric SeqSBAIJ matrices store their entries explicitly
255: if (sym == PETSC_TRUE) {
256: if (isrow == iscol) sameorder = PETSC_TRUE;
257: else PetscCall(ISEqual(isrow, iscol, &sameorder));
258: if (sameorder == PETSC_TRUE) PetscCall(ISSorted(is1, &issorted));
259: }
260: // keep the extracted matrix in upper-triangular storage before restoring the requested block order
261: if (sym == PETSC_TRUE && sameorder == PETSC_TRUE && issorted == PETSC_FALSE) {
262: PetscCheck(scall != MAT_INPLACE_MATRIX, PETSC_COMM_SELF, PETSC_ERR_SUP, "MAT_INPLACE_MATRIX not supported");
263: PetscCall(ISDuplicate(is1, &sorted));
264: PetscCall(ISSort(sorted));
265: PetscCall(MatCreateSubMatrix_SeqSBAIJ_Private(A, sorted, sorted, MAT_INITIAL_MATRIX, C, PETSC_TRUE));
266: PetscCall(MatPropagateSymmetryOptions(A, C[0]));
267: PetscCall(ISSortPermutation(is1, PETSC_TRUE, &perm));
268: PetscCall(ISInvertPermutation(perm, PETSC_DECIDE, &iperm));
269: PetscCall(ISExpandIndicesGeneral(A->rmap->N, A->rmap->n, A->rmap->bs, 1, &iperm, &expanded));
270: PetscCall(MatPermute(C[0], expanded, expanded, &D));
271: if (scall == MAT_REUSE_MATRIX) {
272: PetscCall(MatCopy(D, *B, DIFFERENT_NONZERO_PATTERN));
273: PetscCall(MatDestroy(&D));
274: } else *B = D;
275: PetscCall(MatDestroy(C));
276: PetscCall(ISDestroy(&expanded));
277: PetscCall(ISDestroy(&iperm));
278: PetscCall(ISDestroy(&perm));
279: PetscCall(ISDestroy(&sorted));
280: } else if (sym == PETSC_TRUE || implicit == PETSC_FALSE) PetscCall(MatCreateSubMatrix_SeqSBAIJ_Private(A, is1, is2, scall, B, !implicit ? PETSC_FALSE : sym));
281: else {
282: PetscCall(MatCreateSubMatrix_SeqSBAIJ_Private(A, is1, is2, MAT_INITIAL_MATRIX, C, sym));
283: PetscCall(MatCreateSubMatrix_SeqSBAIJ_Private(A, is2, is1, MAT_INITIAL_MATRIX, C + 1, sym));
284: PetscCall(MatTranspose(C[1], MAT_INPLACE_MATRIX, C + 1));
285: PetscCall(MatAXPY(C[0], 1.0, C[1], DIFFERENT_NONZERO_PATTERN));
286: PetscCheck(scall != MAT_INPLACE_MATRIX, PETSC_COMM_SELF, PETSC_ERR_SUP, "MAT_INPLACE_MATRIX not supported");
287: if (scall == MAT_REUSE_MATRIX) PetscCall(MatCopy(C[0], *B, SAME_NONZERO_PATTERN));
288: else if (A->rmap->bs == 1) PetscCall(MatConvert(C[0], MATAIJ, MAT_INITIAL_MATRIX, B));
289: else {
290: *B = C[0];
291: C[0] = NULL;
292: }
293: PetscCall(MatDestroy(C));
294: PetscCall(MatDestroy(C + 1));
295: }
296: PetscCall(ISDestroy(&is1));
297: PetscCall(ISDestroy(&is2));
299: if (implicit == PETSC_TRUE && sym == PETSC_TRUE && isrow != iscol) {
300: PetscBool isequal;
301: PetscCall(ISEqual(isrow, iscol, &isequal));
302: if (isequal == PETSC_FALSE) PetscCall(MatSeqSBAIJZeroOps_Private(*B));
303: }
304: PetscFunctionReturn(PETSC_SUCCESS);
305: }
307: PetscErrorCode MatCreateSubMatrices_SeqSBAIJ(Mat A, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *B[])
308: {
309: PetscInt i;
311: PetscFunctionBegin;
312: if (scall == MAT_INITIAL_MATRIX) PetscCall(PetscCalloc1(n + 1, B));
314: for (i = 0; i < n; i++) PetscCall(MatCreateSubMatrix_SeqSBAIJ(A, irow[i], icol[i], scall, &(*B)[i]));
315: PetscFunctionReturn(PETSC_SUCCESS);
316: }
318: /* Should check that shapes of vectors and matrices match */
319: PetscErrorCode MatMult_SeqSBAIJ_2(Mat A, Vec xx, Vec zz)
320: {
321: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
322: PetscScalar *z, x1, x2, zero = 0.0;
323: const PetscScalar *x, *xb;
324: const MatScalar *v;
325: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
326: const PetscInt *aj = a->j, *ai = a->i, *ib;
327: PetscInt nonzerorow = 0;
329: PetscFunctionBegin;
330: PetscCall(VecSet(zz, zero));
331: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
332: PetscCall(VecGetArrayRead(xx, &x));
333: PetscCall(VecGetArray(zz, &z));
335: v = a->a;
336: xb = x;
338: for (i = 0; i < mbs; i++, xb += 2, ai++) {
339: n = ai[1] - ai[0]; /* length of i_th block row of A */
340: if (!n) continue;
341: x1 = xb[0];
342: x2 = xb[1];
343: ib = aj + *ai;
344: jmin = 0;
345: nonzerorow++;
346: if (*ib == i) { /* (diag of A)*x */
347: z[2 * i] += v[0] * x1 + v[2] * x2;
348: z[2 * i + 1] += v[2] * x1 + v[3] * x2;
349: v += 4;
350: jmin++;
351: }
352: PetscPrefetchBlock(ib + jmin + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
353: PetscPrefetchBlock(v + 4 * n, 4 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
354: for (j = jmin; j < n; j++) {
355: /* (strict lower triangular part of A)*x */
356: cval = ib[j] * 2;
357: z[cval] += v[0] * x1 + v[1] * x2;
358: z[cval + 1] += v[2] * x1 + v[3] * x2;
359: /* (strict upper triangular part of A)*x */
360: z[2 * i] += v[0] * x[cval] + v[2] * x[cval + 1];
361: z[2 * i + 1] += v[1] * x[cval] + v[3] * x[cval + 1];
362: v += 4;
363: }
364: }
366: PetscCall(VecRestoreArrayRead(xx, &x));
367: PetscCall(VecRestoreArray(zz, &z));
368: PetscCall(PetscLogFlops(8.0 * (a->nz * 2.0 - nonzerorow) - nonzerorow));
369: PetscFunctionReturn(PETSC_SUCCESS);
370: }
372: PetscErrorCode MatMult_SeqSBAIJ_3(Mat A, Vec xx, Vec zz)
373: {
374: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
375: PetscScalar *z, x1, x2, x3, zero = 0.0;
376: const PetscScalar *x, *xb;
377: const MatScalar *v;
378: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
379: const PetscInt *aj = a->j, *ai = a->i, *ib;
380: PetscInt nonzerorow = 0;
382: PetscFunctionBegin;
383: PetscCall(VecSet(zz, zero));
384: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
385: PetscCall(VecGetArrayRead(xx, &x));
386: PetscCall(VecGetArray(zz, &z));
388: v = a->a;
389: xb = x;
391: for (i = 0; i < mbs; i++, xb += 3, ai++) {
392: n = ai[1] - ai[0]; /* length of i_th block row of A */
393: if (!n) continue;
394: x1 = xb[0];
395: x2 = xb[1];
396: x3 = xb[2];
397: ib = aj + *ai;
398: jmin = 0;
399: nonzerorow++;
400: if (*ib == i) { /* (diag of A)*x */
401: z[3 * i] += v[0] * x1 + v[3] * x2 + v[6] * x3;
402: z[3 * i + 1] += v[3] * x1 + v[4] * x2 + v[7] * x3;
403: z[3 * i + 2] += v[6] * x1 + v[7] * x2 + v[8] * x3;
404: v += 9;
405: jmin++;
406: }
407: PetscPrefetchBlock(ib + jmin + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
408: PetscPrefetchBlock(v + 9 * n, 9 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
409: for (j = jmin; j < n; j++) {
410: /* (strict lower triangular part of A)*x */
411: cval = ib[j] * 3;
412: z[cval] += v[0] * x1 + v[1] * x2 + v[2] * x3;
413: z[cval + 1] += v[3] * x1 + v[4] * x2 + v[5] * x3;
414: z[cval + 2] += v[6] * x1 + v[7] * x2 + v[8] * x3;
415: /* (strict upper triangular part of A)*x */
416: z[3 * i] += v[0] * x[cval] + v[3] * x[cval + 1] + v[6] * x[cval + 2];
417: z[3 * i + 1] += v[1] * x[cval] + v[4] * x[cval + 1] + v[7] * x[cval + 2];
418: z[3 * i + 2] += v[2] * x[cval] + v[5] * x[cval + 1] + v[8] * x[cval + 2];
419: v += 9;
420: }
421: }
423: PetscCall(VecRestoreArrayRead(xx, &x));
424: PetscCall(VecRestoreArray(zz, &z));
425: PetscCall(PetscLogFlops(18.0 * (a->nz * 2.0 - nonzerorow) - nonzerorow));
426: PetscFunctionReturn(PETSC_SUCCESS);
427: }
429: PetscErrorCode MatMult_SeqSBAIJ_4(Mat A, Vec xx, Vec zz)
430: {
431: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
432: PetscScalar *z, x1, x2, x3, x4, zero = 0.0;
433: const PetscScalar *x, *xb;
434: const MatScalar *v;
435: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
436: const PetscInt *aj = a->j, *ai = a->i, *ib;
437: PetscInt nonzerorow = 0;
439: PetscFunctionBegin;
440: PetscCall(VecSet(zz, zero));
441: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
442: PetscCall(VecGetArrayRead(xx, &x));
443: PetscCall(VecGetArray(zz, &z));
445: v = a->a;
446: xb = x;
448: for (i = 0; i < mbs; i++, xb += 4, ai++) {
449: n = ai[1] - ai[0]; /* length of i_th block row of A */
450: if (!n) continue;
451: x1 = xb[0];
452: x2 = xb[1];
453: x3 = xb[2];
454: x4 = xb[3];
455: ib = aj + *ai;
456: jmin = 0;
457: nonzerorow++;
458: if (*ib == i) { /* (diag of A)*x */
459: z[4 * i] += v[0] * x1 + v[4] * x2 + v[8] * x3 + v[12] * x4;
460: z[4 * i + 1] += v[4] * x1 + v[5] * x2 + v[9] * x3 + v[13] * x4;
461: z[4 * i + 2] += v[8] * x1 + v[9] * x2 + v[10] * x3 + v[14] * x4;
462: z[4 * i + 3] += v[12] * x1 + v[13] * x2 + v[14] * x3 + v[15] * x4;
463: v += 16;
464: jmin++;
465: }
466: PetscPrefetchBlock(ib + jmin + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
467: PetscPrefetchBlock(v + 16 * n, 16 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
468: for (j = jmin; j < n; j++) {
469: /* (strict lower triangular part of A)*x */
470: cval = ib[j] * 4;
471: z[cval] += v[0] * x1 + v[1] * x2 + v[2] * x3 + v[3] * x4;
472: z[cval + 1] += v[4] * x1 + v[5] * x2 + v[6] * x3 + v[7] * x4;
473: z[cval + 2] += v[8] * x1 + v[9] * x2 + v[10] * x3 + v[11] * x4;
474: z[cval + 3] += v[12] * x1 + v[13] * x2 + v[14] * x3 + v[15] * x4;
475: /* (strict upper triangular part of A)*x */
476: z[4 * i] += v[0] * x[cval] + v[4] * x[cval + 1] + v[8] * x[cval + 2] + v[12] * x[cval + 3];
477: z[4 * i + 1] += v[1] * x[cval] + v[5] * x[cval + 1] + v[9] * x[cval + 2] + v[13] * x[cval + 3];
478: z[4 * i + 2] += v[2] * x[cval] + v[6] * x[cval + 1] + v[10] * x[cval + 2] + v[14] * x[cval + 3];
479: z[4 * i + 3] += v[3] * x[cval] + v[7] * x[cval + 1] + v[11] * x[cval + 2] + v[15] * x[cval + 3];
480: v += 16;
481: }
482: }
484: PetscCall(VecRestoreArrayRead(xx, &x));
485: PetscCall(VecRestoreArray(zz, &z));
486: PetscCall(PetscLogFlops(32.0 * (a->nz * 2.0 - nonzerorow) - nonzerorow));
487: PetscFunctionReturn(PETSC_SUCCESS);
488: }
490: PetscErrorCode MatMult_SeqSBAIJ_5(Mat A, Vec xx, Vec zz)
491: {
492: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
493: PetscScalar *z, x1, x2, x3, x4, x5, zero = 0.0;
494: const PetscScalar *x, *xb;
495: const MatScalar *v;
496: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
497: const PetscInt *aj = a->j, *ai = a->i, *ib;
498: PetscInt nonzerorow = 0;
500: PetscFunctionBegin;
501: PetscCall(VecSet(zz, zero));
502: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
503: PetscCall(VecGetArrayRead(xx, &x));
504: PetscCall(VecGetArray(zz, &z));
506: v = a->a;
507: xb = x;
509: for (i = 0; i < mbs; i++, xb += 5, ai++) {
510: n = ai[1] - ai[0]; /* length of i_th block row of A */
511: if (!n) continue;
512: x1 = xb[0];
513: x2 = xb[1];
514: x3 = xb[2];
515: x4 = xb[3];
516: x5 = xb[4];
517: ib = aj + *ai;
518: jmin = 0;
519: nonzerorow++;
520: if (*ib == i) { /* (diag of A)*x */
521: z[5 * i] += v[0] * x1 + v[5] * x2 + v[10] * x3 + v[15] * x4 + v[20] * x5;
522: z[5 * i + 1] += v[5] * x1 + v[6] * x2 + v[11] * x3 + v[16] * x4 + v[21] * x5;
523: z[5 * i + 2] += v[10] * x1 + v[11] * x2 + v[12] * x3 + v[17] * x4 + v[22] * x5;
524: z[5 * i + 3] += v[15] * x1 + v[16] * x2 + v[17] * x3 + v[18] * x4 + v[23] * x5;
525: z[5 * i + 4] += v[20] * x1 + v[21] * x2 + v[22] * x3 + v[23] * x4 + v[24] * x5;
526: v += 25;
527: jmin++;
528: }
529: PetscPrefetchBlock(ib + jmin + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
530: PetscPrefetchBlock(v + 25 * n, 25 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
531: for (j = jmin; j < n; j++) {
532: /* (strict lower triangular part of A)*x */
533: cval = ib[j] * 5;
534: z[cval] += v[0] * x1 + v[1] * x2 + v[2] * x3 + v[3] * x4 + v[4] * x5;
535: z[cval + 1] += v[5] * x1 + v[6] * x2 + v[7] * x3 + v[8] * x4 + v[9] * x5;
536: z[cval + 2] += v[10] * x1 + v[11] * x2 + v[12] * x3 + v[13] * x4 + v[14] * x5;
537: z[cval + 3] += v[15] * x1 + v[16] * x2 + v[17] * x3 + v[18] * x4 + v[19] * x5;
538: z[cval + 4] += v[20] * x1 + v[21] * x2 + v[22] * x3 + v[23] * x4 + v[24] * x5;
539: /* (strict upper triangular part of A)*x */
540: 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];
541: 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];
542: 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];
543: 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];
544: 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];
545: v += 25;
546: }
547: }
549: PetscCall(VecRestoreArrayRead(xx, &x));
550: PetscCall(VecRestoreArray(zz, &z));
551: PetscCall(PetscLogFlops(50.0 * (a->nz * 2.0 - nonzerorow) - nonzerorow));
552: PetscFunctionReturn(PETSC_SUCCESS);
553: }
555: PetscErrorCode MatMult_SeqSBAIJ_6(Mat A, Vec xx, Vec zz)
556: {
557: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
558: PetscScalar *z, x1, x2, x3, x4, x5, x6, zero = 0.0;
559: const PetscScalar *x, *xb;
560: const MatScalar *v;
561: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
562: const PetscInt *aj = a->j, *ai = a->i, *ib;
563: PetscInt nonzerorow = 0;
565: PetscFunctionBegin;
566: PetscCall(VecSet(zz, zero));
567: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
568: PetscCall(VecGetArrayRead(xx, &x));
569: PetscCall(VecGetArray(zz, &z));
571: v = a->a;
572: xb = x;
574: for (i = 0; i < mbs; i++, xb += 6, ai++) {
575: n = ai[1] - ai[0]; /* length of i_th block row of A */
576: if (!n) continue;
577: x1 = xb[0];
578: x2 = xb[1];
579: x3 = xb[2];
580: x4 = xb[3];
581: x5 = xb[4];
582: x6 = xb[5];
583: ib = aj + *ai;
584: jmin = 0;
585: nonzerorow++;
586: if (*ib == i) { /* (diag of A)*x */
587: z[6 * i] += v[0] * x1 + v[6] * x2 + v[12] * x3 + v[18] * x4 + v[24] * x5 + v[30] * x6;
588: z[6 * i + 1] += v[6] * x1 + v[7] * x2 + v[13] * x3 + v[19] * x4 + v[25] * x5 + v[31] * x6;
589: z[6 * i + 2] += v[12] * x1 + v[13] * x2 + v[14] * x3 + v[20] * x4 + v[26] * x5 + v[32] * x6;
590: z[6 * i + 3] += v[18] * x1 + v[19] * x2 + v[20] * x3 + v[21] * x4 + v[27] * x5 + v[33] * x6;
591: z[6 * i + 4] += v[24] * x1 + v[25] * x2 + v[26] * x3 + v[27] * x4 + v[28] * x5 + v[34] * x6;
592: z[6 * i + 5] += v[30] * x1 + v[31] * x2 + v[32] * x3 + v[33] * x4 + v[34] * x5 + v[35] * x6;
593: v += 36;
594: jmin++;
595: }
596: PetscPrefetchBlock(ib + jmin + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
597: PetscPrefetchBlock(v + 36 * n, 36 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
598: for (j = jmin; j < n; j++) {
599: /* (strict lower triangular part of A)*x */
600: cval = ib[j] * 6;
601: z[cval] += v[0] * x1 + v[1] * x2 + v[2] * x3 + v[3] * x4 + v[4] * x5 + v[5] * x6;
602: z[cval + 1] += v[6] * x1 + v[7] * x2 + v[8] * x3 + v[9] * x4 + v[10] * x5 + v[11] * x6;
603: z[cval + 2] += v[12] * x1 + v[13] * x2 + v[14] * x3 + v[15] * x4 + v[16] * x5 + v[17] * x6;
604: z[cval + 3] += v[18] * x1 + v[19] * x2 + v[20] * x3 + v[21] * x4 + v[22] * x5 + v[23] * x6;
605: z[cval + 4] += v[24] * x1 + v[25] * x2 + v[26] * x3 + v[27] * x4 + v[28] * x5 + v[29] * x6;
606: z[cval + 5] += v[30] * x1 + v[31] * x2 + v[32] * x3 + v[33] * x4 + v[34] * x5 + v[35] * x6;
607: /* (strict upper triangular part of A)*x */
608: 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];
609: 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];
610: 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];
611: 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];
612: 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];
613: 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];
614: v += 36;
615: }
616: }
618: PetscCall(VecRestoreArrayRead(xx, &x));
619: PetscCall(VecRestoreArray(zz, &z));
620: PetscCall(PetscLogFlops(72.0 * (a->nz * 2.0 - nonzerorow) - nonzerorow));
621: PetscFunctionReturn(PETSC_SUCCESS);
622: }
624: PetscErrorCode MatMult_SeqSBAIJ_7(Mat A, Vec xx, Vec zz)
625: {
626: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
627: PetscScalar *z, x1, x2, x3, x4, x5, x6, x7, zero = 0.0;
628: const PetscScalar *x, *xb;
629: const MatScalar *v;
630: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
631: const PetscInt *aj = a->j, *ai = a->i, *ib;
632: PetscInt nonzerorow = 0;
634: PetscFunctionBegin;
635: PetscCall(VecSet(zz, zero));
636: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
637: PetscCall(VecGetArrayRead(xx, &x));
638: PetscCall(VecGetArray(zz, &z));
640: v = a->a;
641: xb = x;
643: for (i = 0; i < mbs; i++, xb += 7, ai++) {
644: n = ai[1] - ai[0]; /* length of i_th block row of A */
645: if (!n) continue;
646: x1 = xb[0];
647: x2 = xb[1];
648: x3 = xb[2];
649: x4 = xb[3];
650: x5 = xb[4];
651: x6 = xb[5];
652: x7 = xb[6];
653: ib = aj + *ai;
654: jmin = 0;
655: nonzerorow++;
656: if (*ib == i) { /* (diag of A)*x */
657: z[7 * i] += v[0] * x1 + v[7] * x2 + v[14] * x3 + v[21] * x4 + v[28] * x5 + v[35] * x6 + v[42] * x7;
658: 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;
659: 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;
660: 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;
661: 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;
662: 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;
663: 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;
664: v += 49;
665: jmin++;
666: }
667: PetscPrefetchBlock(ib + jmin + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
668: PetscPrefetchBlock(v + 49 * n, 49 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
669: for (j = jmin; j < n; j++) {
670: /* (strict lower triangular part of A)*x */
671: cval = ib[j] * 7;
672: z[cval] += v[0] * x1 + v[1] * x2 + v[2] * x3 + v[3] * x4 + v[4] * x5 + v[5] * x6 + v[6] * x7;
673: z[cval + 1] += v[7] * x1 + v[8] * x2 + v[9] * x3 + v[10] * x4 + v[11] * x5 + v[12] * x6 + v[13] * x7;
674: z[cval + 2] += v[14] * x1 + v[15] * x2 + v[16] * x3 + v[17] * x4 + v[18] * x5 + v[19] * x6 + v[20] * x7;
675: z[cval + 3] += v[21] * x1 + v[22] * x2 + v[23] * x3 + v[24] * x4 + v[25] * x5 + v[26] * x6 + v[27] * x7;
676: z[cval + 4] += v[28] * x1 + v[29] * x2 + v[30] * x3 + v[31] * x4 + v[32] * x5 + v[33] * x6 + v[34] * x7;
677: z[cval + 5] += v[35] * x1 + v[36] * x2 + v[37] * x3 + v[38] * x4 + v[39] * x5 + v[40] * x6 + v[41] * x7;
678: z[cval + 6] += v[42] * x1 + v[43] * x2 + v[44] * x3 + v[45] * x4 + v[46] * x5 + v[47] * x6 + v[48] * x7;
679: /* (strict upper triangular part of A)*x */
680: 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];
681: 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];
682: 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];
683: 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];
684: 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];
685: 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];
686: 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];
687: v += 49;
688: }
689: }
690: PetscCall(VecRestoreArrayRead(xx, &x));
691: PetscCall(VecRestoreArray(zz, &z));
692: PetscCall(PetscLogFlops(98.0 * (a->nz * 2.0 - nonzerorow) - nonzerorow));
693: PetscFunctionReturn(PETSC_SUCCESS);
694: }
696: /*
697: This will not work with MatScalar == float because it calls the BLAS
698: */
699: PetscErrorCode MatMult_SeqSBAIJ_N(Mat A, Vec xx, Vec zz)
700: {
701: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
702: PetscScalar *z, *z_ptr, *zb, *work, *workt, zero = 0.0;
703: const PetscScalar *x, *x_ptr, *xb;
704: const MatScalar *v;
705: PetscInt mbs = a->mbs, i, bs = A->rmap->bs, j, n, bs2 = a->bs2, ncols, k;
706: const PetscInt *idx, *aj, *ii;
707: PetscInt nonzerorow = 0;
709: PetscFunctionBegin;
710: PetscCall(VecSet(zz, zero));
711: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
712: PetscCall(VecGetArrayRead(xx, &x));
713: PetscCall(VecGetArray(zz, &z));
715: x_ptr = x;
716: z_ptr = z;
718: aj = a->j;
719: v = a->a;
720: ii = a->i;
722: if (!a->mult_work) PetscCall(PetscMalloc1(A->rmap->N + 1, &a->mult_work));
723: work = a->mult_work;
725: for (i = 0; i < mbs; i++) {
726: n = ii[1] - ii[0];
727: ncols = n * bs;
728: workt = work;
729: idx = aj + ii[0];
730: nonzerorow += (n > 0);
732: /* upper triangular part */
733: for (j = 0; j < n; j++) {
734: xb = x_ptr + bs * (*idx++);
735: for (k = 0; k < bs; k++) workt[k] = xb[k];
736: workt += bs;
737: }
738: /* z(i*bs:(i+1)*bs-1) += A(i,:)*x */
739: PetscKernel_w_gets_w_plus_Ar_times_v(bs, ncols, work, v, z);
741: /* strict lower triangular part */
742: idx = aj + ii[0];
743: if (n && *idx == i) {
744: ncols -= bs;
745: v += bs2;
746: idx++;
747: n--;
748: }
750: if (ncols > 0) {
751: workt = work;
752: PetscCall(PetscArrayzero(workt, ncols));
753: PetscKernel_w_gets_w_plus_trans_Ar_times_v(bs, ncols, x, v, workt);
754: for (j = 0; j < n; j++) {
755: zb = z_ptr + bs * (*idx++);
756: for (k = 0; k < bs; k++) zb[k] += workt[k];
757: workt += bs;
758: }
759: }
760: x += bs;
761: v += n * bs2;
762: z += bs;
763: ii++;
764: }
766: PetscCall(VecRestoreArrayRead(xx, &x));
767: PetscCall(VecRestoreArray(zz, &z));
768: PetscCall(PetscLogFlops(2.0 * (a->nz * 2.0 - nonzerorow) * bs2 - nonzerorow));
769: PetscFunctionReturn(PETSC_SUCCESS);
770: }
772: PetscErrorCode MatMultAdd_SeqSBAIJ_1(Mat A, Vec xx, Vec yy, Vec zz)
773: {
774: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
775: PetscScalar *z, x1;
776: const PetscScalar *x, *xb;
777: const MatScalar *v;
778: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
779: const PetscInt *aj = a->j, *ai = a->i, *ib;
780: PetscInt nonzerorow = 0;
781: const int aconj = PetscDefined(USE_COMPLEX) && A->hermitian == PETSC_BOOL3_TRUE ? 1 : 0;
783: PetscFunctionBegin;
784: PetscCall(VecCopy(yy, zz));
785: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
786: PetscCall(VecGetArrayRead(xx, &x));
787: PetscCall(VecGetArray(zz, &z));
788: v = a->a;
789: xb = x;
791: for (i = 0; i < mbs; i++, xb++, ai++) {
792: n = ai[1] - ai[0]; /* length of i_th row of A */
793: if (!n) continue;
794: x1 = xb[0];
795: ib = aj + *ai;
796: jmin = 0;
797: nonzerorow++;
798: if (*ib == i) { /* (diag of A)*x */
799: z[i] += *v++ * x[*ib++];
800: jmin++;
801: }
802: if (aconj) {
803: for (j = jmin; j < n; j++) {
804: cval = *ib;
805: z[cval] += PetscConj(*v) * x1; /* (strict lower triangular part of A)*x */
806: z[i] += *v++ * x[*ib++]; /* (strict upper triangular part of A)*x */
807: }
808: } else {
809: for (j = jmin; j < n; j++) {
810: cval = *ib;
811: z[cval] += *v * x1; /* (strict lower triangular part of A)*x */
812: z[i] += *v++ * x[*ib++]; /* (strict upper triangular part of A)*x */
813: }
814: }
815: }
817: PetscCall(VecRestoreArrayRead(xx, &x));
818: PetscCall(VecRestoreArray(zz, &z));
820: PetscCall(PetscLogFlops(2.0 * (a->nz * 2.0 - nonzerorow)));
821: PetscFunctionReturn(PETSC_SUCCESS);
822: }
824: PetscErrorCode MatMultAdd_SeqSBAIJ_2(Mat A, Vec xx, Vec yy, Vec zz)
825: {
826: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
827: PetscScalar *z, x1, x2;
828: const PetscScalar *x, *xb;
829: const MatScalar *v;
830: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
831: const PetscInt *aj = a->j, *ai = a->i, *ib;
832: PetscInt nonzerorow = 0;
834: PetscFunctionBegin;
835: PetscCall(VecCopy(yy, zz));
836: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
837: PetscCall(VecGetArrayRead(xx, &x));
838: PetscCall(VecGetArray(zz, &z));
840: v = a->a;
841: xb = x;
843: for (i = 0; i < mbs; i++, xb += 2, ai++) {
844: n = ai[1] - ai[0]; /* length of i_th block row of A */
845: if (!n) continue;
846: x1 = xb[0];
847: x2 = xb[1];
848: ib = aj + *ai;
849: jmin = 0;
850: nonzerorow++;
851: if (*ib == i) { /* (diag of A)*x */
852: z[2 * i] += v[0] * x1 + v[2] * x2;
853: z[2 * i + 1] += v[2] * x1 + v[3] * x2;
854: v += 4;
855: jmin++;
856: }
857: PetscPrefetchBlock(ib + jmin + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
858: PetscPrefetchBlock(v + 4 * n, 4 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
859: for (j = jmin; j < n; j++) {
860: /* (strict lower triangular part of A)*x */
861: cval = ib[j] * 2;
862: z[cval] += v[0] * x1 + v[1] * x2;
863: z[cval + 1] += v[2] * x1 + v[3] * x2;
864: /* (strict upper triangular part of A)*x */
865: z[2 * i] += v[0] * x[cval] + v[2] * x[cval + 1];
866: z[2 * i + 1] += v[1] * x[cval] + v[3] * x[cval + 1];
867: v += 4;
868: }
869: }
870: PetscCall(VecRestoreArrayRead(xx, &x));
871: PetscCall(VecRestoreArray(zz, &z));
873: PetscCall(PetscLogFlops(8.0 * (a->nz * 2.0 - nonzerorow)));
874: PetscFunctionReturn(PETSC_SUCCESS);
875: }
877: PetscErrorCode MatMultAdd_SeqSBAIJ_3(Mat A, Vec xx, Vec yy, Vec zz)
878: {
879: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
880: PetscScalar *z, x1, x2, x3;
881: const PetscScalar *x, *xb;
882: const MatScalar *v;
883: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
884: const PetscInt *aj = a->j, *ai = a->i, *ib;
885: PetscInt nonzerorow = 0;
887: PetscFunctionBegin;
888: PetscCall(VecCopy(yy, zz));
889: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
890: PetscCall(VecGetArrayRead(xx, &x));
891: PetscCall(VecGetArray(zz, &z));
893: v = a->a;
894: xb = x;
896: for (i = 0; i < mbs; i++, xb += 3, ai++) {
897: n = ai[1] - ai[0]; /* length of i_th block row of A */
898: if (!n) continue;
899: x1 = xb[0];
900: x2 = xb[1];
901: x3 = xb[2];
902: ib = aj + *ai;
903: jmin = 0;
904: nonzerorow++;
905: if (*ib == i) { /* (diag of A)*x */
906: z[3 * i] += v[0] * x1 + v[3] * x2 + v[6] * x3;
907: z[3 * i + 1] += v[3] * x1 + v[4] * x2 + v[7] * x3;
908: z[3 * i + 2] += v[6] * x1 + v[7] * x2 + v[8] * x3;
909: v += 9;
910: jmin++;
911: }
912: PetscPrefetchBlock(ib + jmin + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
913: PetscPrefetchBlock(v + 9 * n, 9 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
914: for (j = jmin; j < n; j++) {
915: /* (strict lower triangular part of A)*x */
916: cval = ib[j] * 3;
917: z[cval] += v[0] * x1 + v[1] * x2 + v[2] * x3;
918: z[cval + 1] += v[3] * x1 + v[4] * x2 + v[5] * x3;
919: z[cval + 2] += v[6] * x1 + v[7] * x2 + v[8] * x3;
920: /* (strict upper triangular part of A)*x */
921: z[3 * i] += v[0] * x[cval] + v[3] * x[cval + 1] + v[6] * x[cval + 2];
922: z[3 * i + 1] += v[1] * x[cval] + v[4] * x[cval + 1] + v[7] * x[cval + 2];
923: z[3 * i + 2] += v[2] * x[cval] + v[5] * x[cval + 1] + v[8] * x[cval + 2];
924: v += 9;
925: }
926: }
928: PetscCall(VecRestoreArrayRead(xx, &x));
929: PetscCall(VecRestoreArray(zz, &z));
931: PetscCall(PetscLogFlops(18.0 * (a->nz * 2.0 - nonzerorow)));
932: PetscFunctionReturn(PETSC_SUCCESS);
933: }
935: PetscErrorCode MatMultAdd_SeqSBAIJ_4(Mat A, Vec xx, Vec yy, Vec zz)
936: {
937: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
938: PetscScalar *z, x1, x2, x3, x4;
939: const PetscScalar *x, *xb;
940: const MatScalar *v;
941: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
942: const PetscInt *aj = a->j, *ai = a->i, *ib;
943: PetscInt nonzerorow = 0;
945: PetscFunctionBegin;
946: PetscCall(VecCopy(yy, zz));
947: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
948: PetscCall(VecGetArrayRead(xx, &x));
949: PetscCall(VecGetArray(zz, &z));
951: v = a->a;
952: xb = x;
954: for (i = 0; i < mbs; i++, xb += 4, ai++) {
955: n = ai[1] - ai[0]; /* length of i_th block row of A */
956: if (!n) continue;
957: x1 = xb[0];
958: x2 = xb[1];
959: x3 = xb[2];
960: x4 = xb[3];
961: ib = aj + *ai;
962: jmin = 0;
963: nonzerorow++;
964: if (*ib == i) { /* (diag of A)*x */
965: z[4 * i] += v[0] * x1 + v[4] * x2 + v[8] * x3 + v[12] * x4;
966: z[4 * i + 1] += v[4] * x1 + v[5] * x2 + v[9] * x3 + v[13] * x4;
967: z[4 * i + 2] += v[8] * x1 + v[9] * x2 + v[10] * x3 + v[14] * x4;
968: z[4 * i + 3] += v[12] * x1 + v[13] * x2 + v[14] * x3 + v[15] * x4;
969: v += 16;
970: jmin++;
971: }
972: PetscPrefetchBlock(ib + jmin + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
973: PetscPrefetchBlock(v + 16 * n, 16 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
974: for (j = jmin; j < n; j++) {
975: /* (strict lower triangular part of A)*x */
976: cval = ib[j] * 4;
977: z[cval] += v[0] * x1 + v[1] * x2 + v[2] * x3 + v[3] * x4;
978: z[cval + 1] += v[4] * x1 + v[5] * x2 + v[6] * x3 + v[7] * x4;
979: z[cval + 2] += v[8] * x1 + v[9] * x2 + v[10] * x3 + v[11] * x4;
980: z[cval + 3] += v[12] * x1 + v[13] * x2 + v[14] * x3 + v[15] * x4;
981: /* (strict upper triangular part of A)*x */
982: z[4 * i] += v[0] * x[cval] + v[4] * x[cval + 1] + v[8] * x[cval + 2] + v[12] * x[cval + 3];
983: z[4 * i + 1] += v[1] * x[cval] + v[5] * x[cval + 1] + v[9] * x[cval + 2] + v[13] * x[cval + 3];
984: z[4 * i + 2] += v[2] * x[cval] + v[6] * x[cval + 1] + v[10] * x[cval + 2] + v[14] * x[cval + 3];
985: z[4 * i + 3] += v[3] * x[cval] + v[7] * x[cval + 1] + v[11] * x[cval + 2] + v[15] * x[cval + 3];
986: v += 16;
987: }
988: }
990: PetscCall(VecRestoreArrayRead(xx, &x));
991: PetscCall(VecRestoreArray(zz, &z));
993: PetscCall(PetscLogFlops(32.0 * (a->nz * 2.0 - nonzerorow)));
994: PetscFunctionReturn(PETSC_SUCCESS);
995: }
997: PetscErrorCode MatMultAdd_SeqSBAIJ_5(Mat A, Vec xx, Vec yy, Vec zz)
998: {
999: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1000: PetscScalar *z, x1, x2, x3, x4, x5;
1001: const PetscScalar *x, *xb;
1002: const MatScalar *v;
1003: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
1004: const PetscInt *aj = a->j, *ai = a->i, *ib;
1005: PetscInt nonzerorow = 0;
1007: PetscFunctionBegin;
1008: PetscCall(VecCopy(yy, zz));
1009: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
1010: PetscCall(VecGetArrayRead(xx, &x));
1011: PetscCall(VecGetArray(zz, &z));
1013: v = a->a;
1014: xb = x;
1016: for (i = 0; i < mbs; i++, xb += 5, ai++) {
1017: n = ai[1] - ai[0]; /* length of i_th block row of A */
1018: if (!n) continue;
1019: x1 = xb[0];
1020: x2 = xb[1];
1021: x3 = xb[2];
1022: x4 = xb[3];
1023: x5 = xb[4];
1024: ib = aj + *ai;
1025: jmin = 0;
1026: nonzerorow++;
1027: if (*ib == i) { /* (diag of A)*x */
1028: z[5 * i] += v[0] * x1 + v[5] * x2 + v[10] * x3 + v[15] * x4 + v[20] * x5;
1029: z[5 * i + 1] += v[5] * x1 + v[6] * x2 + v[11] * x3 + v[16] * x4 + v[21] * x5;
1030: z[5 * i + 2] += v[10] * x1 + v[11] * x2 + v[12] * x3 + v[17] * x4 + v[22] * x5;
1031: z[5 * i + 3] += v[15] * x1 + v[16] * x2 + v[17] * x3 + v[18] * x4 + v[23] * x5;
1032: z[5 * i + 4] += v[20] * x1 + v[21] * x2 + v[22] * x3 + v[23] * x4 + v[24] * x5;
1033: v += 25;
1034: jmin++;
1035: }
1036: PetscPrefetchBlock(ib + jmin + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
1037: PetscPrefetchBlock(v + 25 * n, 25 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
1038: for (j = jmin; j < n; j++) {
1039: /* (strict lower triangular part of A)*x */
1040: cval = ib[j] * 5;
1041: z[cval] += v[0] * x1 + v[1] * x2 + v[2] * x3 + v[3] * x4 + v[4] * x5;
1042: z[cval + 1] += v[5] * x1 + v[6] * x2 + v[7] * x3 + v[8] * x4 + v[9] * x5;
1043: z[cval + 2] += v[10] * x1 + v[11] * x2 + v[12] * x3 + v[13] * x4 + v[14] * x5;
1044: z[cval + 3] += v[15] * x1 + v[16] * x2 + v[17] * x3 + v[18] * x4 + v[19] * x5;
1045: z[cval + 4] += v[20] * x1 + v[21] * x2 + v[22] * x3 + v[23] * x4 + v[24] * x5;
1046: /* (strict upper triangular part of A)*x */
1047: 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];
1048: 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];
1049: 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];
1050: 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];
1051: 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];
1052: v += 25;
1053: }
1054: }
1056: PetscCall(VecRestoreArrayRead(xx, &x));
1057: PetscCall(VecRestoreArray(zz, &z));
1059: PetscCall(PetscLogFlops(50.0 * (a->nz * 2.0 - nonzerorow)));
1060: PetscFunctionReturn(PETSC_SUCCESS);
1061: }
1063: PetscErrorCode MatMultAdd_SeqSBAIJ_6(Mat A, Vec xx, Vec yy, Vec zz)
1064: {
1065: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1066: PetscScalar *z, x1, x2, x3, x4, x5, x6;
1067: const PetscScalar *x, *xb;
1068: const MatScalar *v;
1069: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
1070: const PetscInt *aj = a->j, *ai = a->i, *ib;
1071: PetscInt nonzerorow = 0;
1073: PetscFunctionBegin;
1074: PetscCall(VecCopy(yy, zz));
1075: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
1076: PetscCall(VecGetArrayRead(xx, &x));
1077: PetscCall(VecGetArray(zz, &z));
1079: v = a->a;
1080: xb = x;
1082: for (i = 0; i < mbs; i++, xb += 6, ai++) {
1083: n = ai[1] - ai[0]; /* length of i_th block row of A */
1084: if (!n) continue;
1085: x1 = xb[0];
1086: x2 = xb[1];
1087: x3 = xb[2];
1088: x4 = xb[3];
1089: x5 = xb[4];
1090: x6 = xb[5];
1091: ib = aj + *ai;
1092: jmin = 0;
1093: nonzerorow++;
1094: if (*ib == i) { /* (diag of A)*x */
1095: z[6 * i] += v[0] * x1 + v[6] * x2 + v[12] * x3 + v[18] * x4 + v[24] * x5 + v[30] * x6;
1096: z[6 * i + 1] += v[6] * x1 + v[7] * x2 + v[13] * x3 + v[19] * x4 + v[25] * x5 + v[31] * x6;
1097: z[6 * i + 2] += v[12] * x1 + v[13] * x2 + v[14] * x3 + v[20] * x4 + v[26] * x5 + v[32] * x6;
1098: z[6 * i + 3] += v[18] * x1 + v[19] * x2 + v[20] * x3 + v[21] * x4 + v[27] * x5 + v[33] * x6;
1099: z[6 * i + 4] += v[24] * x1 + v[25] * x2 + v[26] * x3 + v[27] * x4 + v[28] * x5 + v[34] * x6;
1100: z[6 * i + 5] += v[30] * x1 + v[31] * x2 + v[32] * x3 + v[33] * x4 + v[34] * x5 + v[35] * x6;
1101: v += 36;
1102: jmin++;
1103: }
1104: PetscPrefetchBlock(ib + jmin + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
1105: PetscPrefetchBlock(v + 36 * n, 36 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
1106: for (j = jmin; j < n; j++) {
1107: /* (strict lower triangular part of A)*x */
1108: cval = ib[j] * 6;
1109: z[cval] += v[0] * x1 + v[1] * x2 + v[2] * x3 + v[3] * x4 + v[4] * x5 + v[5] * x6;
1110: z[cval + 1] += v[6] * x1 + v[7] * x2 + v[8] * x3 + v[9] * x4 + v[10] * x5 + v[11] * x6;
1111: z[cval + 2] += v[12] * x1 + v[13] * x2 + v[14] * x3 + v[15] * x4 + v[16] * x5 + v[17] * x6;
1112: z[cval + 3] += v[18] * x1 + v[19] * x2 + v[20] * x3 + v[21] * x4 + v[22] * x5 + v[23] * x6;
1113: z[cval + 4] += v[24] * x1 + v[25] * x2 + v[26] * x3 + v[27] * x4 + v[28] * x5 + v[29] * x6;
1114: z[cval + 5] += v[30] * x1 + v[31] * x2 + v[32] * x3 + v[33] * x4 + v[34] * x5 + v[35] * x6;
1115: /* (strict upper triangular part of A)*x */
1116: 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];
1117: 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];
1118: 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];
1119: 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];
1120: 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];
1121: 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];
1122: v += 36;
1123: }
1124: }
1126: PetscCall(VecRestoreArrayRead(xx, &x));
1127: PetscCall(VecRestoreArray(zz, &z));
1129: PetscCall(PetscLogFlops(72.0 * (a->nz * 2.0 - nonzerorow)));
1130: PetscFunctionReturn(PETSC_SUCCESS);
1131: }
1133: PetscErrorCode MatMultAdd_SeqSBAIJ_7(Mat A, Vec xx, Vec yy, Vec zz)
1134: {
1135: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1136: PetscScalar *z, x1, x2, x3, x4, x5, x6, x7;
1137: const PetscScalar *x, *xb;
1138: const MatScalar *v;
1139: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
1140: const PetscInt *aj = a->j, *ai = a->i, *ib;
1141: PetscInt nonzerorow = 0;
1143: PetscFunctionBegin;
1144: PetscCall(VecCopy(yy, zz));
1145: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
1146: PetscCall(VecGetArrayRead(xx, &x));
1147: PetscCall(VecGetArray(zz, &z));
1149: v = a->a;
1150: xb = x;
1152: for (i = 0; i < mbs; i++, xb += 7, ai++) {
1153: n = ai[1] - ai[0]; /* length of i_th block row of A */
1154: if (!n) continue;
1155: x1 = xb[0];
1156: x2 = xb[1];
1157: x3 = xb[2];
1158: x4 = xb[3];
1159: x5 = xb[4];
1160: x6 = xb[5];
1161: x7 = xb[6];
1162: ib = aj + *ai;
1163: jmin = 0;
1164: nonzerorow++;
1165: if (*ib == i) { /* (diag of A)*x */
1166: z[7 * i] += v[0] * x1 + v[7] * x2 + v[14] * x3 + v[21] * x4 + v[28] * x5 + v[35] * x6 + v[42] * x7;
1167: 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;
1168: 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;
1169: 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;
1170: 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;
1171: 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;
1172: 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;
1173: v += 49;
1174: jmin++;
1175: }
1176: PetscPrefetchBlock(ib + jmin + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
1177: PetscPrefetchBlock(v + 49 * n, 49 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
1178: for (j = jmin; j < n; j++) {
1179: /* (strict lower triangular part of A)*x */
1180: cval = ib[j] * 7;
1181: z[cval] += v[0] * x1 + v[1] * x2 + v[2] * x3 + v[3] * x4 + v[4] * x5 + v[5] * x6 + v[6] * x7;
1182: z[cval + 1] += v[7] * x1 + v[8] * x2 + v[9] * x3 + v[10] * x4 + v[11] * x5 + v[12] * x6 + v[13] * x7;
1183: z[cval + 2] += v[14] * x1 + v[15] * x2 + v[16] * x3 + v[17] * x4 + v[18] * x5 + v[19] * x6 + v[20] * x7;
1184: z[cval + 3] += v[21] * x1 + v[22] * x2 + v[23] * x3 + v[24] * x4 + v[25] * x5 + v[26] * x6 + v[27] * x7;
1185: z[cval + 4] += v[28] * x1 + v[29] * x2 + v[30] * x3 + v[31] * x4 + v[32] * x5 + v[33] * x6 + v[34] * x7;
1186: z[cval + 5] += v[35] * x1 + v[36] * x2 + v[37] * x3 + v[38] * x4 + v[39] * x5 + v[40] * x6 + v[41] * x7;
1187: z[cval + 6] += v[42] * x1 + v[43] * x2 + v[44] * x3 + v[45] * x4 + v[46] * x5 + v[47] * x6 + v[48] * x7;
1188: /* (strict upper triangular part of A)*x */
1189: 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];
1190: 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];
1191: 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];
1192: 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];
1193: 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];
1194: 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];
1195: 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];
1196: v += 49;
1197: }
1198: }
1200: PetscCall(VecRestoreArrayRead(xx, &x));
1201: PetscCall(VecRestoreArray(zz, &z));
1203: PetscCall(PetscLogFlops(98.0 * (a->nz * 2.0 - nonzerorow)));
1204: PetscFunctionReturn(PETSC_SUCCESS);
1205: }
1207: PetscErrorCode MatMultAdd_SeqSBAIJ_N(Mat A, Vec xx, Vec yy, Vec zz)
1208: {
1209: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1210: PetscScalar *z, *z_ptr = NULL, *zb, *work, *workt;
1211: const PetscScalar *x, *x_ptr, *xb;
1212: const MatScalar *v;
1213: PetscInt mbs = a->mbs, i, bs = A->rmap->bs, j, n, bs2 = a->bs2, ncols, k;
1214: const PetscInt *idx, *aj, *ii;
1215: PetscInt nonzerorow = 0;
1217: PetscFunctionBegin;
1218: PetscCall(VecCopy(yy, zz));
1219: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
1220: PetscCall(VecGetArrayRead(xx, &x));
1221: x_ptr = x;
1222: PetscCall(VecGetArray(zz, &z));
1223: z_ptr = z;
1225: aj = a->j;
1226: v = a->a;
1227: ii = a->i;
1229: if (!a->mult_work) PetscCall(PetscMalloc1(A->rmap->n + 1, &a->mult_work));
1230: work = a->mult_work;
1232: for (i = 0; i < mbs; i++) {
1233: n = ii[1] - ii[0];
1234: ncols = n * bs;
1235: workt = work;
1236: idx = aj + ii[0];
1237: nonzerorow += (n > 0);
1239: /* upper triangular part */
1240: for (j = 0; j < n; j++) {
1241: xb = x_ptr + bs * (*idx++);
1242: for (k = 0; k < bs; k++) workt[k] = xb[k];
1243: workt += bs;
1244: }
1245: /* z(i*bs:(i+1)*bs-1) += A(i,:)*x */
1246: PetscKernel_w_gets_w_plus_Ar_times_v(bs, ncols, work, v, z);
1248: /* strict lower triangular part */
1249: idx = aj + ii[0];
1250: if (n && *idx == i) {
1251: ncols -= bs;
1252: v += bs2;
1253: idx++;
1254: n--;
1255: }
1256: if (ncols > 0) {
1257: workt = work;
1258: PetscCall(PetscArrayzero(workt, ncols));
1259: PetscKernel_w_gets_w_plus_trans_Ar_times_v(bs, ncols, x, v, workt);
1260: for (j = 0; j < n; j++) {
1261: zb = z_ptr + bs * (*idx++);
1262: for (k = 0; k < bs; k++) zb[k] += workt[k];
1263: workt += bs;
1264: }
1265: }
1267: x += bs;
1268: v += n * bs2;
1269: z += bs;
1270: ii++;
1271: }
1273: PetscCall(VecRestoreArrayRead(xx, &x));
1274: PetscCall(VecRestoreArray(zz, &z));
1276: PetscCall(PetscLogFlops(2.0 * bs2 * (a->nz * 2.0 - nonzerorow)));
1277: PetscFunctionReturn(PETSC_SUCCESS);
1278: }
1280: PetscErrorCode MatScale_SeqSBAIJ(Mat inA, PetscScalar alpha)
1281: {
1282: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)inA->data;
1283: PetscScalar oalpha = alpha;
1284: PetscBLASInt one = 1, totalnz;
1286: PetscFunctionBegin;
1287: PetscCall(PetscBLASIntCast(a->bs2 * a->nz, &totalnz));
1288: PetscCallBLAS("BLASscal", BLASscal_(&totalnz, &oalpha, a->a, &one));
1289: PetscCall(PetscLogFlops(totalnz));
1290: PetscFunctionReturn(PETSC_SUCCESS);
1291: }
1293: PetscErrorCode MatNorm_SeqSBAIJ(Mat A, NormType type, PetscReal *norm)
1294: {
1295: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1296: const MatScalar *v = a->a;
1297: PetscReal sum_diag = 0.0, sum_off = 0.0, *sum;
1298: PetscInt i, j, k, bs = A->rmap->bs, bs2 = a->bs2, k1, mbs = a->mbs, jmin, jmax, nexti, ik, *jl, *il;
1299: const PetscInt *aj = a->j, *col;
1301: PetscFunctionBegin;
1302: if (!a->nz) {
1303: *norm = 0.0;
1304: PetscFunctionReturn(PETSC_SUCCESS);
1305: }
1306: if (type == NORM_FROBENIUS) {
1307: for (k = 0; k < mbs; k++) {
1308: jmin = a->i[k];
1309: jmax = a->i[k + 1];
1310: col = aj + jmin;
1311: if (jmax - jmin > 0 && *col == k) { /* diagonal block */
1312: for (i = 0; i < bs2; i++) {
1313: sum_diag += PetscRealPart(PetscConj(*v) * (*v));
1314: v++;
1315: }
1316: jmin++;
1317: }
1318: for (j = jmin; j < jmax; j++) { /* off-diagonal blocks */
1319: for (i = 0; i < bs2; i++) {
1320: sum_off += PetscRealPart(PetscConj(*v) * (*v));
1321: v++;
1322: }
1323: }
1324: }
1325: *norm = PetscSqrtReal(sum_diag + 2 * sum_off);
1326: PetscCall(PetscLogFlops(2.0 * bs2 * a->nz));
1327: } else if (type == NORM_INFINITY || type == NORM_1) { /* maximum row/column sum */
1328: PetscCall(PetscMalloc3(bs, &sum, mbs, &il, mbs, &jl));
1329: for (i = 0; i < mbs; i++) jl[i] = mbs;
1330: il[0] = 0;
1332: *norm = 0.0;
1333: for (k = 0; k < mbs; k++) { /* k_th block row */
1334: for (j = 0; j < bs; j++) sum[j] = 0.0;
1335: /*-- col sum --*/
1336: i = jl[k]; /* first |A(i,k)| to be added */
1337: /* jl[k]=i: first nonzero element in row i for submatrix A(1:k,k:n) (active window)
1338: at step k */
1339: while (i < mbs) {
1340: nexti = jl[i]; /* next block row to be added */
1341: ik = il[i]; /* block index of A(i,k) in the array a */
1342: for (j = 0; j < bs; j++) {
1343: v = a->a + ik * bs2 + j * bs;
1344: for (k1 = 0; k1 < bs; k1++) {
1345: sum[j] += PetscAbsScalar(*v);
1346: v++;
1347: }
1348: }
1349: /* update il, jl */
1350: jmin = ik + 1; /* block index of array a: points to the next nonzero of A in row i */
1351: jmax = a->i[i + 1];
1352: if (jmin < jmax) {
1353: il[i] = jmin;
1354: j = a->j[jmin];
1355: jl[i] = jl[j];
1356: jl[j] = i;
1357: }
1358: i = nexti;
1359: }
1360: /*-- row sum --*/
1361: jmin = a->i[k];
1362: jmax = a->i[k + 1];
1363: for (i = jmin; i < jmax; i++) {
1364: for (j = 0; j < bs; j++) {
1365: v = a->a + i * bs2 + j;
1366: for (k1 = 0; k1 < bs; k1++) {
1367: sum[j] += PetscAbsScalar(*v);
1368: v += bs;
1369: }
1370: }
1371: }
1372: /* add k_th block row to il, jl */
1373: col = aj + jmin;
1374: if (jmax - jmin > 0 && *col == k) jmin++;
1375: if (jmin < jmax) {
1376: il[k] = jmin;
1377: j = a->j[jmin];
1378: jl[k] = jl[j];
1379: jl[j] = k;
1380: }
1381: for (j = 0; j < bs; j++) {
1382: if (sum[j] > *norm) *norm = sum[j];
1383: }
1384: }
1385: PetscCall(PetscFree3(sum, il, jl));
1386: PetscCall(PetscLogFlops(PetscMax(mbs * a->nz - 1, 0)));
1387: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for this norm yet");
1388: PetscFunctionReturn(PETSC_SUCCESS);
1389: }
1391: PetscErrorCode MatEqual_SeqSBAIJ(Mat A, Mat B, PetscBool *flg)
1392: {
1393: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data, *b = (Mat_SeqSBAIJ *)B->data;
1395: PetscFunctionBegin;
1396: /* If the matrix/block dimensions are not equal, or no of nonzeros or shift */
1397: if ((A->rmap->N != B->rmap->N) || (A->cmap->n != B->cmap->n) || (A->rmap->bs != B->rmap->bs) || (a->nz != b->nz)) {
1398: *flg = PETSC_FALSE;
1399: PetscFunctionReturn(PETSC_SUCCESS);
1400: }
1402: /* if the a->i are the same */
1403: PetscCall(PetscArraycmp(a->i, b->i, a->mbs + 1, flg));
1404: if (!*flg) PetscFunctionReturn(PETSC_SUCCESS);
1406: /* if a->j are the same */
1407: PetscCall(PetscArraycmp(a->j, b->j, a->nz, flg));
1408: if (!*flg) PetscFunctionReturn(PETSC_SUCCESS);
1410: /* if a->a are the same */
1411: PetscCall(PetscArraycmp(a->a, b->a, a->nz * A->rmap->bs * A->rmap->bs, flg));
1412: PetscFunctionReturn(PETSC_SUCCESS);
1413: }
1415: PetscErrorCode MatGetDiagonal_SeqSBAIJ(Mat A, Vec v)
1416: {
1417: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1418: PetscInt n;
1419: const PetscInt bs = A->rmap->bs, ambs = a->mbs, bs2 = a->bs2;
1420: PetscScalar *x;
1421: const MatScalar *aa = a->a, *aa_j;
1422: const PetscInt *ai = a->i, *adiag;
1423: PetscBool diagDense;
1425: PetscFunctionBegin;
1426: PetscCheck(!A->factortype || bs <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix with bs>1");
1427: PetscCall(MatGetDiagonalMarkers_SeqSBAIJ(A, &adiag, &diagDense));
1428: if (A->factortype == MAT_FACTOR_CHOLESKY || A->factortype == MAT_FACTOR_ICC) {
1429: PetscCall(VecGetArrayWrite(v, &x));
1430: for (PetscInt i = 0; i < ambs; i++) x[i] = 1.0 / aa[adiag[i]];
1431: PetscCall(VecRestoreArrayWrite(v, &x));
1432: PetscFunctionReturn(PETSC_SUCCESS);
1433: }
1435: PetscCall(VecGetLocalSize(v, &n));
1436: PetscCheck(n == A->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector");
1437: PetscCall(VecGetArrayWrite(v, &x));
1439: if (diagDense) {
1440: for (PetscInt i = 0, row = 0; i < ambs; i++) {
1441: aa_j = aa + adiag[i] * bs2;
1442: for (PetscInt k = 0; k < bs2; k += (bs + 1)) x[row++] = aa_j[k];
1443: }
1444: } else {
1445: for (PetscInt i = 0, row = 0; i < ambs; i++) {
1446: const PetscInt j = adiag[i];
1448: if (j != ai[i + 1]) {
1449: aa_j = aa + j * bs2;
1450: for (PetscInt k = 0; k < bs2; k += (bs + 1)) x[row++] = aa_j[k];
1451: } else {
1452: for (PetscInt k = 0; k < bs; k++) x[row++] = 0.0;
1453: }
1454: }
1455: }
1456: PetscCall(VecRestoreArrayWrite(v, &x));
1457: PetscFunctionReturn(PETSC_SUCCESS);
1458: }
1460: PetscErrorCode MatDiagonalScale_SeqSBAIJ(Mat A, Vec ll, Vec rr)
1461: {
1462: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1463: const PetscScalar *l = NULL;
1464: MatScalar *aa = a->a;
1465: PetscInt lm, m = A->rmap->N, mbs = a->mbs, bs = A->rmap->bs, bs2 = a->bs2;
1466: const PetscInt *ai = a->i, *aj = a->j;
1468: PetscFunctionBegin;
1469: if (ll != rr) {
1470: Mat_SeqBAIJ *b;
1471: Mat B;
1472: const PetscScalar *r = NULL;
1473: PetscInt *browlengths, *browstart, *bj;
1474: MatScalar *ba;
1475: PetscInt n = A->cmap->N;
1476: PetscBool hermitian = (PetscBool)(PetscDefined(USE_COMPLEX) && A->hermitian == PETSC_BOOL3_TRUE);
1478: if (ll) {
1479: PetscCall(VecGetLocalSize(ll, &lm));
1480: PetscCheck(lm == m, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Left scaling vector wrong length");
1481: }
1482: if (rr) {
1483: PetscInt rn;
1485: PetscCall(VecGetLocalSize(rr, &rn));
1486: PetscCheck(rn == n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Right scaling vector wrong length");
1487: }
1488: if (ll) PetscCall(VecGetArrayRead(ll, &l));
1489: if (rr) PetscCall(VecGetArrayRead(rr, &r));
1490: PetscCall(PetscCalloc1(mbs, &browlengths));
1491: PetscCall(PetscMalloc1(mbs, &browstart));
1492: for (PetscInt i = 0; i < mbs; i++) {
1493: for (PetscInt k = ai[i]; k < ai[i + 1]; k++) {
1494: browlengths[i]++;
1495: if (aj[k] != i) browlengths[aj[k]]++;
1496: }
1497: }
1498: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
1499: PetscCall(MatSetSizes(B, m, n, m, n));
1500: PetscCall(MatSetType(B, MATSEQBAIJ));
1501: PetscCall(MatSeqBAIJSetPreallocation(B, bs, 0, browlengths));
1502: b = (Mat_SeqBAIJ *)B->data;
1503: ba = b->a;
1504: bj = b->j;
1505: for (PetscInt i = 0; i < mbs; i++) {
1506: b->ilen[i] = browlengths[i];
1507: browstart[i] = b->i[i];
1508: }
1509: PetscCall(PetscFree(browlengths));
1510: for (PetscInt i = 0; i < mbs; i++) {
1511: for (PetscInt k = ai[i]; k < ai[i + 1]; k++) {
1512: const PetscInt j = aj[k];
1513: const MatScalar *av = aa + k * bs2;
1514: MatScalar *v = ba + browstart[i] * bs2;
1515: const PetscScalar *li = PetscSafePointerPlusOffset(l, i * bs), *ri = PetscSafePointerPlusOffset(r, j * bs);
1517: bj[browstart[i]++] = j;
1518: for (PetscInt col = 0; col < bs; col++) {
1519: const PetscScalar x = r != NULL ? ri[col] : 1.0;
1521: for (PetscInt row = 0; row < bs; row++) v[col * bs + row] = av[col * bs + row] * (l != NULL ? li[row] : 1.0) * x;
1522: }
1523: if (j != i) {
1524: MatScalar *v = ba + browstart[j] * bs2;
1525: const PetscScalar *li = PetscSafePointerPlusOffset(l, j * bs), *ri = PetscSafePointerPlusOffset(r, i * bs);
1527: bj[browstart[j]++] = i;
1528: for (PetscInt col = 0; col < bs; col++) {
1529: const PetscScalar x = r != NULL ? ri[col] : 1.0;
1531: for (PetscInt row = 0; row < bs; row++) v[col * bs + row] = (hermitian == PETSC_TRUE ? PetscConj(av[row * bs + col]) : av[row * bs + col]) * (l != NULL ? li[row] : 1.0) * x;
1532: }
1533: }
1534: }
1535: }
1536: PetscCall(PetscFree(browstart));
1537: if (ll) PetscCall(VecRestoreArrayRead(ll, &l));
1538: if (rr) PetscCall(VecRestoreArrayRead(rr, &r));
1539: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
1540: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
1541: PetscCall(PetscLogFlops(((ll ? 1.0 : 0.0) + (rr ? 1.0 : 0.0)) * b->nz * bs2));
1542: B->symmetric = A->symmetric;
1543: B->structurally_symmetric = A->structurally_symmetric;
1544: B->hermitian = A->hermitian;
1545: PetscCall(MatHeaderReplace(A, &B));
1546: PetscFunctionReturn(PETSC_SUCCESS);
1547: }
1548: if (!ll) PetscFunctionReturn(PETSC_SUCCESS);
1549: PetscCall(VecGetArrayRead(ll, &l));
1550: PetscCall(VecGetLocalSize(ll, &lm));
1551: PetscCheck(lm == m, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Left scaling vector wrong length");
1552: for (PetscInt i = 0; i < mbs; i++) { /* for each block row */
1553: const PetscScalar *li = l + i * bs;
1554: MatScalar *v = aa + bs2 * ai[i];
1556: for (PetscInt j = 0; j < ai[i + 1] - ai[i]; j++) { /* for each block */
1557: const PetscScalar *ri = l + bs * aj[ai[i] + j];
1559: for (PetscInt k = 0; k < bs; k++) {
1560: for (PetscInt row = 0; row < bs; row++) (*v++) *= li[row] * ri[k];
1561: }
1562: }
1563: }
1564: PetscCall(VecRestoreArrayRead(ll, &l));
1565: PetscCall(PetscLogFlops(2.0 * a->nz * bs2));
1566: PetscFunctionReturn(PETSC_SUCCESS);
1567: }
1569: PetscErrorCode MatGetInfo_SeqSBAIJ(Mat A, MatInfoType flag, MatInfo *info)
1570: {
1571: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1573: PetscFunctionBegin;
1574: info->block_size = a->bs2;
1575: info->nz_allocated = a->bs2 * a->maxnz; /*num. of nonzeros in upper triangular part */
1576: info->nz_used = a->bs2 * a->nz; /*num. of nonzeros in upper triangular part */
1577: info->nz_unneeded = info->nz_allocated - info->nz_used;
1578: info->assemblies = A->num_ass;
1579: info->mallocs = A->info.mallocs;
1580: info->memory = 0; /* REVIEW ME */
1581: if (A->factortype) {
1582: info->fill_ratio_given = A->info.fill_ratio_given;
1583: info->fill_ratio_needed = A->info.fill_ratio_needed;
1584: info->factor_mallocs = A->info.factor_mallocs;
1585: } else {
1586: info->fill_ratio_given = 0;
1587: info->fill_ratio_needed = 0;
1588: info->factor_mallocs = 0;
1589: }
1590: PetscFunctionReturn(PETSC_SUCCESS);
1591: }
1593: PetscErrorCode MatZeroEntries_SeqSBAIJ(Mat A)
1594: {
1595: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1597: PetscFunctionBegin;
1598: PetscCall(PetscArrayzero(a->a, a->bs2 * a->i[a->mbs]));
1599: PetscFunctionReturn(PETSC_SUCCESS);
1600: }
1602: PetscErrorCode MatGetRowMaxAbs_SeqSBAIJ(Mat A, Vec v, PetscInt idx[])
1603: {
1604: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1605: PetscInt i, j, n, row, col, bs, mbs;
1606: const PetscInt *ai, *aj;
1607: PetscReal atmp;
1608: const MatScalar *aa;
1609: PetscScalar *x;
1610: PetscInt ncols, brow, bcol, krow, kcol;
1612: PetscFunctionBegin;
1613: PetscCheck(!idx, PETSC_COMM_SELF, PETSC_ERR_SUP, "Send email to petsc-maint@mcs.anl.gov");
1614: PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
1615: bs = A->rmap->bs;
1616: aa = a->a;
1617: ai = a->i;
1618: aj = a->j;
1619: mbs = a->mbs;
1621: PetscCall(VecSet(v, 0.0));
1622: PetscCall(VecGetArray(v, &x));
1623: PetscCall(VecGetLocalSize(v, &n));
1624: PetscCheck(n == A->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector");
1625: for (i = 0; i < mbs; i++) {
1626: ncols = ai[1] - ai[0];
1627: ai++;
1628: brow = bs * i;
1629: for (j = 0; j < ncols; j++) {
1630: bcol = bs * (*aj);
1631: for (kcol = 0; kcol < bs; kcol++) {
1632: col = bcol + kcol; /* col index */
1633: for (krow = 0; krow < bs; krow++) {
1634: atmp = PetscAbsScalar(*aa);
1635: aa++;
1636: row = brow + krow; /* row index */
1637: if (PetscRealPart(x[row]) < atmp) x[row] = atmp;
1638: if (*aj > i && PetscRealPart(x[col]) < atmp) x[col] = atmp;
1639: }
1640: }
1641: aj++;
1642: }
1643: }
1644: PetscCall(VecRestoreArray(v, &x));
1645: PetscFunctionReturn(PETSC_SUCCESS);
1646: }
1648: PetscErrorCode MatMatMultSymbolic_SeqSBAIJ_SeqDense(Mat A, Mat B, PetscReal fill, Mat C)
1649: {
1650: PetscFunctionBegin;
1651: PetscCall(MatMatMultSymbolic_SeqDense_SeqDense(A, B, 0.0, C));
1652: C->ops->matmultnumeric = MatMatMultNumeric_SeqSBAIJ_SeqDense;
1653: PetscFunctionReturn(PETSC_SUCCESS);
1654: }
1656: static PetscErrorCode MatMatMult_SeqSBAIJ_1_Private(Mat A, PetscScalar *b, PetscInt bm, PetscScalar *c, PetscInt cm, PetscInt cn)
1657: {
1658: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1659: PetscScalar *z = c;
1660: const PetscScalar *xb;
1661: PetscScalar x1;
1662: const MatScalar *v = a->a, *vv;
1663: PetscInt mbs = a->mbs, i, *idx = a->j, *ii = a->i, j, *jj, n, k;
1664: const int aconj = PetscDefined(USE_COMPLEX) && A->hermitian == PETSC_BOOL3_TRUE ? 1 : 0;
1666: PetscFunctionBegin;
1667: for (i = 0; i < mbs; i++) {
1668: n = ii[1] - ii[0];
1669: ii++;
1670: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
1671: PetscPrefetchBlock(v + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
1672: jj = idx;
1673: vv = v;
1674: for (k = 0; k < cn; k++) {
1675: idx = jj;
1676: v = vv;
1677: for (j = 0; j < n; j++) {
1678: xb = b + (*idx);
1679: x1 = xb[0 + k * bm];
1680: z[0 + k * cm] += v[0] * x1;
1681: if (*idx != i) c[(*idx) + k * cm] += (aconj ? PetscConj(v[0]) : v[0]) * b[i + k * bm];
1682: v += 1;
1683: ++idx;
1684: }
1685: }
1686: z += 1;
1687: }
1688: PetscFunctionReturn(PETSC_SUCCESS);
1689: }
1691: static PetscErrorCode MatMatMult_SeqSBAIJ_2_Private(Mat A, PetscScalar *b, PetscInt bm, PetscScalar *c, PetscInt cm, PetscInt cn)
1692: {
1693: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1694: PetscScalar *z = c;
1695: const PetscScalar *xb;
1696: PetscScalar x1, x2;
1697: const MatScalar *v = a->a, *vv;
1698: PetscInt mbs = a->mbs, i, *idx = a->j, *ii = a->i, j, *jj, n, k;
1700: PetscFunctionBegin;
1701: for (i = 0; i < mbs; i++) {
1702: n = ii[1] - ii[0];
1703: ii++;
1704: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
1705: PetscPrefetchBlock(v + 4 * n, 4 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
1706: jj = idx;
1707: vv = v;
1708: for (k = 0; k < cn; k++) {
1709: idx = jj;
1710: v = vv;
1711: for (j = 0; j < n; j++) {
1712: xb = b + 2 * (*idx);
1713: x1 = xb[0 + k * bm];
1714: x2 = xb[1 + k * bm];
1715: z[0 + k * cm] += v[0] * x1 + v[2] * x2;
1716: z[1 + k * cm] += v[1] * x1 + v[3] * x2;
1717: if (*idx != i) {
1718: c[2 * (*idx) + 0 + k * cm] += v[0] * b[2 * i + k * bm] + v[1] * b[2 * i + 1 + k * bm];
1719: c[2 * (*idx) + 1 + k * cm] += v[2] * b[2 * i + k * bm] + v[3] * b[2 * i + 1 + k * bm];
1720: }
1721: v += 4;
1722: ++idx;
1723: }
1724: }
1725: z += 2;
1726: }
1727: PetscFunctionReturn(PETSC_SUCCESS);
1728: }
1730: static PetscErrorCode MatMatMult_SeqSBAIJ_3_Private(Mat A, PetscScalar *b, PetscInt bm, PetscScalar *c, PetscInt cm, PetscInt cn)
1731: {
1732: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1733: PetscScalar *z = c;
1734: const PetscScalar *xb;
1735: PetscScalar x1, x2, x3;
1736: const MatScalar *v = a->a, *vv;
1737: PetscInt mbs = a->mbs, i, *idx = a->j, *ii = a->i, j, *jj, n, k;
1739: PetscFunctionBegin;
1740: for (i = 0; i < mbs; i++) {
1741: n = ii[1] - ii[0];
1742: ii++;
1743: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
1744: PetscPrefetchBlock(v + 9 * n, 9 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
1745: jj = idx;
1746: vv = v;
1747: for (k = 0; k < cn; k++) {
1748: idx = jj;
1749: v = vv;
1750: for (j = 0; j < n; j++) {
1751: xb = b + 3 * (*idx);
1752: x1 = xb[0 + k * bm];
1753: x2 = xb[1 + k * bm];
1754: x3 = xb[2 + k * bm];
1755: z[0 + k * cm] += v[0] * x1 + v[3] * x2 + v[6] * x3;
1756: z[1 + k * cm] += v[1] * x1 + v[4] * x2 + v[7] * x3;
1757: z[2 + k * cm] += v[2] * x1 + v[5] * x2 + v[8] * x3;
1758: if (*idx != i) {
1759: 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];
1760: 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];
1761: 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];
1762: }
1763: v += 9;
1764: ++idx;
1765: }
1766: }
1767: z += 3;
1768: }
1769: PetscFunctionReturn(PETSC_SUCCESS);
1770: }
1772: static PetscErrorCode MatMatMult_SeqSBAIJ_4_Private(Mat A, PetscScalar *b, PetscInt bm, PetscScalar *c, PetscInt cm, PetscInt cn)
1773: {
1774: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1775: PetscScalar *z = c;
1776: const PetscScalar *xb;
1777: PetscScalar x1, x2, x3, x4;
1778: const MatScalar *v = a->a, *vv;
1779: PetscInt mbs = a->mbs, i, *idx = a->j, *ii = a->i, j, *jj, n, k;
1781: PetscFunctionBegin;
1782: for (i = 0; i < mbs; i++) {
1783: n = ii[1] - ii[0];
1784: ii++;
1785: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
1786: PetscPrefetchBlock(v + 16 * n, 16 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
1787: jj = idx;
1788: vv = v;
1789: for (k = 0; k < cn; k++) {
1790: idx = jj;
1791: v = vv;
1792: for (j = 0; j < n; j++) {
1793: xb = b + 4 * (*idx);
1794: x1 = xb[0 + k * bm];
1795: x2 = xb[1 + k * bm];
1796: x3 = xb[2 + k * bm];
1797: x4 = xb[3 + k * bm];
1798: z[0 + k * cm] += v[0] * x1 + v[4] * x2 + v[8] * x3 + v[12] * x4;
1799: z[1 + k * cm] += v[1] * x1 + v[5] * x2 + v[9] * x3 + v[13] * x4;
1800: z[2 + k * cm] += v[2] * x1 + v[6] * x2 + v[10] * x3 + v[14] * x4;
1801: z[3 + k * cm] += v[3] * x1 + v[7] * x2 + v[11] * x3 + v[15] * x4;
1802: if (*idx != i) {
1803: 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];
1804: 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];
1805: 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];
1806: 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];
1807: }
1808: v += 16;
1809: ++idx;
1810: }
1811: }
1812: z += 4;
1813: }
1814: PetscFunctionReturn(PETSC_SUCCESS);
1815: }
1817: static PetscErrorCode MatMatMult_SeqSBAIJ_5_Private(Mat A, PetscScalar *b, PetscInt bm, PetscScalar *c, PetscInt cm, PetscInt cn)
1818: {
1819: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1820: PetscScalar *z = c;
1821: const PetscScalar *xb;
1822: PetscScalar x1, x2, x3, x4, x5;
1823: const MatScalar *v = a->a, *vv;
1824: PetscInt mbs = a->mbs, i, *idx = a->j, *ii = a->i, j, *jj, n, k;
1826: PetscFunctionBegin;
1827: for (i = 0; i < mbs; i++) {
1828: n = ii[1] - ii[0];
1829: ii++;
1830: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
1831: PetscPrefetchBlock(v + 25 * n, 25 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
1832: jj = idx;
1833: vv = v;
1834: for (k = 0; k < cn; k++) {
1835: idx = jj;
1836: v = vv;
1837: for (j = 0; j < n; j++) {
1838: xb = b + 5 * (*idx);
1839: x1 = xb[0 + k * bm];
1840: x2 = xb[1 + k * bm];
1841: x3 = xb[2 + k * bm];
1842: x4 = xb[3 + k * bm];
1843: x5 = xb[4 + k * cm];
1844: z[0 + k * cm] += v[0] * x1 + v[5] * x2 + v[10] * x3 + v[15] * x4 + v[20] * x5;
1845: z[1 + k * cm] += v[1] * x1 + v[6] * x2 + v[11] * x3 + v[16] * x4 + v[21] * x5;
1846: z[2 + k * cm] += v[2] * x1 + v[7] * x2 + v[12] * x3 + v[17] * x4 + v[22] * x5;
1847: z[3 + k * cm] += v[3] * x1 + v[8] * x2 + v[13] * x3 + v[18] * x4 + v[23] * x5;
1848: z[4 + k * cm] += v[4] * x1 + v[9] * x2 + v[14] * x3 + v[19] * x4 + v[24] * x5;
1849: if (*idx != i) {
1850: 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];
1851: 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];
1852: 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];
1853: 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];
1854: 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];
1855: }
1856: v += 25;
1857: ++idx;
1858: }
1859: }
1860: z += 5;
1861: }
1862: PetscFunctionReturn(PETSC_SUCCESS);
1863: }
1865: PetscErrorCode MatMatMultNumeric_SeqSBAIJ_SeqDense(Mat A, Mat B, Mat C)
1866: {
1867: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1868: Mat_SeqDense *bd = (Mat_SeqDense *)B->data;
1869: Mat_SeqDense *cd = (Mat_SeqDense *)C->data;
1870: PetscInt cm = cd->lda, cn = B->cmap->n, bm = bd->lda;
1871: PetscInt mbs, i, bs = A->rmap->bs, j, n, bs2 = a->bs2;
1872: PetscBLASInt bbs, bcn, bbm, bcm;
1873: PetscScalar *z = NULL;
1874: PetscScalar *c, *b;
1875: const MatScalar *v;
1876: const PetscInt *idx, *ii;
1877: PetscScalar _DOne = 1.0;
1879: PetscFunctionBegin;
1880: if (!cm || !cn) PetscFunctionReturn(PETSC_SUCCESS);
1881: 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);
1882: 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);
1883: 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);
1884: b = bd->v;
1885: PetscCall(MatZeroEntries(C));
1886: PetscCall(MatDenseGetArray(C, &c));
1887: switch (bs) {
1888: case 1:
1889: PetscCall(MatMatMult_SeqSBAIJ_1_Private(A, b, bm, c, cm, cn));
1890: break;
1891: case 2:
1892: PetscCall(MatMatMult_SeqSBAIJ_2_Private(A, b, bm, c, cm, cn));
1893: break;
1894: case 3:
1895: PetscCall(MatMatMult_SeqSBAIJ_3_Private(A, b, bm, c, cm, cn));
1896: break;
1897: case 4:
1898: PetscCall(MatMatMult_SeqSBAIJ_4_Private(A, b, bm, c, cm, cn));
1899: break;
1900: case 5:
1901: PetscCall(MatMatMult_SeqSBAIJ_5_Private(A, b, bm, c, cm, cn));
1902: break;
1903: default: /* block sizes larger than 5 by 5 are handled by BLAS */
1904: PetscCall(PetscBLASIntCast(bs, &bbs));
1905: PetscCall(PetscBLASIntCast(cn, &bcn));
1906: PetscCall(PetscBLASIntCast(bm, &bbm));
1907: PetscCall(PetscBLASIntCast(cm, &bcm));
1908: idx = a->j;
1909: v = a->a;
1910: mbs = a->mbs;
1911: ii = a->i;
1912: z = c;
1913: for (i = 0; i < mbs; i++) {
1914: n = ii[1] - ii[0];
1915: ii++;
1916: for (j = 0; j < n; j++) {
1917: if (*idx != i) PetscCallBLAS("BLASgemm", BLASgemm_("T", "N", &bbs, &bcn, &bbs, &_DOne, v, &bbs, b + bs * i, &bbm, &_DOne, c + bs * (*idx), &bcm));
1918: PetscCallBLAS("BLASgemm", BLASgemm_("N", "N", &bbs, &bcn, &bbs, &_DOne, v, &bbs, b + bs * (*idx++), &bbm, &_DOne, z, &bcm));
1919: v += bs2;
1920: }
1921: z += bs;
1922: }
1923: }
1924: PetscCall(MatDenseRestoreArray(C, &c));
1925: PetscCall(PetscLogFlops((2.0 * (a->nz * 2.0 - a->nonzerorowcnt) * bs2 - a->nonzerorowcnt) * cn));
1926: PetscFunctionReturn(PETSC_SUCCESS);
1927: }