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++) {
339: n = ai[1] - ai[0]; /* length of i_th block row of A */
340: x1 = xb[0];
341: x2 = xb[1];
342: ib = aj + *ai;
343: jmin = 0;
344: nonzerorow += (n > 0);
345: if (*ib == i) { /* (diag of A)*x */
346: z[2 * i] += v[0] * x1 + v[2] * x2;
347: z[2 * i + 1] += v[2] * x1 + v[3] * x2;
348: v += 4;
349: jmin++;
350: }
351: PetscPrefetchBlock(ib + jmin + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
352: PetscPrefetchBlock(v + 4 * n, 4 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
353: for (j = jmin; j < n; j++) {
354: /* (strict lower triangular part of A)*x */
355: cval = ib[j] * 2;
356: z[cval] += v[0] * x1 + v[1] * x2;
357: z[cval + 1] += v[2] * x1 + v[3] * x2;
358: /* (strict upper triangular part of A)*x */
359: z[2 * i] += v[0] * x[cval] + v[2] * x[cval + 1];
360: z[2 * i + 1] += v[1] * x[cval] + v[3] * x[cval + 1];
361: v += 4;
362: }
363: xb += 2;
364: ai++;
365: }
367: PetscCall(VecRestoreArrayRead(xx, &x));
368: PetscCall(VecRestoreArray(zz, &z));
369: PetscCall(PetscLogFlops(8.0 * (a->nz * 2.0 - nonzerorow) - nonzerorow));
370: PetscFunctionReturn(PETSC_SUCCESS);
371: }
373: PetscErrorCode MatMult_SeqSBAIJ_3(Mat A, Vec xx, Vec zz)
374: {
375: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
376: PetscScalar *z, x1, x2, x3, zero = 0.0;
377: const PetscScalar *x, *xb;
378: const MatScalar *v;
379: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
380: const PetscInt *aj = a->j, *ai = a->i, *ib;
381: PetscInt nonzerorow = 0;
383: PetscFunctionBegin;
384: PetscCall(VecSet(zz, zero));
385: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
386: PetscCall(VecGetArrayRead(xx, &x));
387: PetscCall(VecGetArray(zz, &z));
389: v = a->a;
390: xb = x;
392: for (i = 0; i < mbs; i++) {
393: n = ai[1] - ai[0]; /* length of i_th block row of A */
394: x1 = xb[0];
395: x2 = xb[1];
396: x3 = xb[2];
397: ib = aj + *ai;
398: jmin = 0;
399: nonzerorow += (n > 0);
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: xb += 3;
422: ai++;
423: }
425: PetscCall(VecRestoreArrayRead(xx, &x));
426: PetscCall(VecRestoreArray(zz, &z));
427: PetscCall(PetscLogFlops(18.0 * (a->nz * 2.0 - nonzerorow) - nonzerorow));
428: PetscFunctionReturn(PETSC_SUCCESS);
429: }
431: PetscErrorCode MatMult_SeqSBAIJ_4(Mat A, Vec xx, Vec zz)
432: {
433: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
434: PetscScalar *z, x1, x2, x3, x4, zero = 0.0;
435: const PetscScalar *x, *xb;
436: const MatScalar *v;
437: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
438: const PetscInt *aj = a->j, *ai = a->i, *ib;
439: PetscInt nonzerorow = 0;
441: PetscFunctionBegin;
442: PetscCall(VecSet(zz, zero));
443: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
444: PetscCall(VecGetArrayRead(xx, &x));
445: PetscCall(VecGetArray(zz, &z));
447: v = a->a;
448: xb = x;
450: for (i = 0; i < mbs; i++) {
451: n = ai[1] - ai[0]; /* length of i_th block row of A */
452: x1 = xb[0];
453: x2 = xb[1];
454: x3 = xb[2];
455: x4 = xb[3];
456: ib = aj + *ai;
457: jmin = 0;
458: nonzerorow += (n > 0);
459: if (*ib == i) { /* (diag of A)*x */
460: z[4 * i] += v[0] * x1 + v[4] * x2 + v[8] * x3 + v[12] * x4;
461: z[4 * i + 1] += v[4] * x1 + v[5] * x2 + v[9] * x3 + v[13] * x4;
462: z[4 * i + 2] += v[8] * x1 + v[9] * x2 + v[10] * x3 + v[14] * x4;
463: z[4 * i + 3] += v[12] * x1 + v[13] * x2 + v[14] * x3 + v[15] * x4;
464: v += 16;
465: jmin++;
466: }
467: PetscPrefetchBlock(ib + jmin + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
468: PetscPrefetchBlock(v + 16 * n, 16 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
469: for (j = jmin; j < n; j++) {
470: /* (strict lower triangular part of A)*x */
471: cval = ib[j] * 4;
472: z[cval] += v[0] * x1 + v[1] * x2 + v[2] * x3 + v[3] * x4;
473: z[cval + 1] += v[4] * x1 + v[5] * x2 + v[6] * x3 + v[7] * x4;
474: z[cval + 2] += v[8] * x1 + v[9] * x2 + v[10] * x3 + v[11] * x4;
475: z[cval + 3] += v[12] * x1 + v[13] * x2 + v[14] * x3 + v[15] * x4;
476: /* (strict upper triangular part of A)*x */
477: z[4 * i] += v[0] * x[cval] + v[4] * x[cval + 1] + v[8] * x[cval + 2] + v[12] * x[cval + 3];
478: z[4 * i + 1] += v[1] * x[cval] + v[5] * x[cval + 1] + v[9] * x[cval + 2] + v[13] * x[cval + 3];
479: z[4 * i + 2] += v[2] * x[cval] + v[6] * x[cval + 1] + v[10] * x[cval + 2] + v[14] * x[cval + 3];
480: z[4 * i + 3] += v[3] * x[cval] + v[7] * x[cval + 1] + v[11] * x[cval + 2] + v[15] * x[cval + 3];
481: v += 16;
482: }
483: xb += 4;
484: ai++;
485: }
487: PetscCall(VecRestoreArrayRead(xx, &x));
488: PetscCall(VecRestoreArray(zz, &z));
489: PetscCall(PetscLogFlops(32.0 * (a->nz * 2.0 - nonzerorow) - nonzerorow));
490: PetscFunctionReturn(PETSC_SUCCESS);
491: }
493: PetscErrorCode MatMult_SeqSBAIJ_5(Mat A, Vec xx, Vec zz)
494: {
495: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
496: PetscScalar *z, x1, x2, x3, x4, x5, zero = 0.0;
497: const PetscScalar *x, *xb;
498: const MatScalar *v;
499: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
500: const PetscInt *aj = a->j, *ai = a->i, *ib;
501: PetscInt nonzerorow = 0;
503: PetscFunctionBegin;
504: PetscCall(VecSet(zz, zero));
505: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
506: PetscCall(VecGetArrayRead(xx, &x));
507: PetscCall(VecGetArray(zz, &z));
509: v = a->a;
510: xb = x;
512: for (i = 0; i < mbs; i++) {
513: n = ai[1] - ai[0]; /* length of i_th block row of A */
514: x1 = xb[0];
515: x2 = xb[1];
516: x3 = xb[2];
517: x4 = xb[3];
518: x5 = xb[4];
519: ib = aj + *ai;
520: jmin = 0;
521: nonzerorow += (n > 0);
522: if (*ib == i) { /* (diag of A)*x */
523: z[5 * i] += v[0] * x1 + v[5] * x2 + v[10] * x3 + v[15] * x4 + v[20] * x5;
524: z[5 * i + 1] += v[5] * x1 + v[6] * x2 + v[11] * x3 + v[16] * x4 + v[21] * x5;
525: z[5 * i + 2] += v[10] * x1 + v[11] * x2 + v[12] * x3 + v[17] * x4 + v[22] * x5;
526: z[5 * i + 3] += v[15] * x1 + v[16] * x2 + v[17] * x3 + v[18] * x4 + v[23] * x5;
527: z[5 * i + 4] += v[20] * x1 + v[21] * x2 + v[22] * x3 + v[23] * x4 + v[24] * x5;
528: v += 25;
529: jmin++;
530: }
531: PetscPrefetchBlock(ib + jmin + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
532: PetscPrefetchBlock(v + 25 * n, 25 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
533: for (j = jmin; j < n; j++) {
534: /* (strict lower triangular part of A)*x */
535: cval = ib[j] * 5;
536: z[cval] += v[0] * x1 + v[1] * x2 + v[2] * x3 + v[3] * x4 + v[4] * x5;
537: z[cval + 1] += v[5] * x1 + v[6] * x2 + v[7] * x3 + v[8] * x4 + v[9] * x5;
538: z[cval + 2] += v[10] * x1 + v[11] * x2 + v[12] * x3 + v[13] * x4 + v[14] * x5;
539: z[cval + 3] += v[15] * x1 + v[16] * x2 + v[17] * x3 + v[18] * x4 + v[19] * x5;
540: z[cval + 4] += v[20] * x1 + v[21] * x2 + v[22] * x3 + v[23] * x4 + v[24] * x5;
541: /* (strict upper triangular part of A)*x */
542: 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];
543: 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];
544: 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];
545: 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];
546: 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];
547: v += 25;
548: }
549: xb += 5;
550: ai++;
551: }
553: PetscCall(VecRestoreArrayRead(xx, &x));
554: PetscCall(VecRestoreArray(zz, &z));
555: PetscCall(PetscLogFlops(50.0 * (a->nz * 2.0 - nonzerorow) - nonzerorow));
556: PetscFunctionReturn(PETSC_SUCCESS);
557: }
559: PetscErrorCode MatMult_SeqSBAIJ_6(Mat A, Vec xx, Vec zz)
560: {
561: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
562: PetscScalar *z, x1, x2, x3, x4, x5, x6, zero = 0.0;
563: const PetscScalar *x, *xb;
564: const MatScalar *v;
565: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
566: const PetscInt *aj = a->j, *ai = a->i, *ib;
567: PetscInt nonzerorow = 0;
569: PetscFunctionBegin;
570: PetscCall(VecSet(zz, zero));
571: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
572: PetscCall(VecGetArrayRead(xx, &x));
573: PetscCall(VecGetArray(zz, &z));
575: v = a->a;
576: xb = x;
578: for (i = 0; i < mbs; i++) {
579: n = ai[1] - ai[0]; /* length of i_th block row of A */
580: x1 = xb[0];
581: x2 = xb[1];
582: x3 = xb[2];
583: x4 = xb[3];
584: x5 = xb[4];
585: x6 = xb[5];
586: ib = aj + *ai;
587: jmin = 0;
588: nonzerorow += (n > 0);
589: if (*ib == i) { /* (diag of A)*x */
590: z[6 * i] += v[0] * x1 + v[6] * x2 + v[12] * x3 + v[18] * x4 + v[24] * x5 + v[30] * x6;
591: z[6 * i + 1] += v[6] * x1 + v[7] * x2 + v[13] * x3 + v[19] * x4 + v[25] * x5 + v[31] * x6;
592: z[6 * i + 2] += v[12] * x1 + v[13] * x2 + v[14] * x3 + v[20] * x4 + v[26] * x5 + v[32] * x6;
593: z[6 * i + 3] += v[18] * x1 + v[19] * x2 + v[20] * x3 + v[21] * x4 + v[27] * x5 + v[33] * x6;
594: z[6 * i + 4] += v[24] * x1 + v[25] * x2 + v[26] * x3 + v[27] * x4 + v[28] * x5 + v[34] * x6;
595: z[6 * i + 5] += v[30] * x1 + v[31] * x2 + v[32] * x3 + v[33] * x4 + v[34] * x5 + v[35] * x6;
596: v += 36;
597: jmin++;
598: }
599: PetscPrefetchBlock(ib + jmin + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
600: PetscPrefetchBlock(v + 36 * n, 36 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
601: for (j = jmin; j < n; j++) {
602: /* (strict lower triangular part of A)*x */
603: cval = ib[j] * 6;
604: z[cval] += v[0] * x1 + v[1] * x2 + v[2] * x3 + v[3] * x4 + v[4] * x5 + v[5] * x6;
605: z[cval + 1] += v[6] * x1 + v[7] * x2 + v[8] * x3 + v[9] * x4 + v[10] * x5 + v[11] * x6;
606: z[cval + 2] += v[12] * x1 + v[13] * x2 + v[14] * x3 + v[15] * x4 + v[16] * x5 + v[17] * x6;
607: z[cval + 3] += v[18] * x1 + v[19] * x2 + v[20] * x3 + v[21] * x4 + v[22] * x5 + v[23] * x6;
608: z[cval + 4] += v[24] * x1 + v[25] * x2 + v[26] * x3 + v[27] * x4 + v[28] * x5 + v[29] * x6;
609: z[cval + 5] += v[30] * x1 + v[31] * x2 + v[32] * x3 + v[33] * x4 + v[34] * x5 + v[35] * x6;
610: /* (strict upper triangular part of A)*x */
611: 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];
612: 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];
613: 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];
614: 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];
615: 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];
616: 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];
617: v += 36;
618: }
619: xb += 6;
620: ai++;
621: }
623: PetscCall(VecRestoreArrayRead(xx, &x));
624: PetscCall(VecRestoreArray(zz, &z));
625: PetscCall(PetscLogFlops(72.0 * (a->nz * 2.0 - nonzerorow) - nonzerorow));
626: PetscFunctionReturn(PETSC_SUCCESS);
627: }
629: PetscErrorCode MatMult_SeqSBAIJ_7(Mat A, Vec xx, Vec zz)
630: {
631: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
632: PetscScalar *z, x1, x2, x3, x4, x5, x6, x7, zero = 0.0;
633: const PetscScalar *x, *xb;
634: const MatScalar *v;
635: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
636: const PetscInt *aj = a->j, *ai = a->i, *ib;
637: PetscInt nonzerorow = 0;
639: PetscFunctionBegin;
640: PetscCall(VecSet(zz, zero));
641: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
642: PetscCall(VecGetArrayRead(xx, &x));
643: PetscCall(VecGetArray(zz, &z));
645: v = a->a;
646: xb = x;
648: for (i = 0; i < mbs; i++) {
649: n = ai[1] - ai[0]; /* length of i_th block row of A */
650: x1 = xb[0];
651: x2 = xb[1];
652: x3 = xb[2];
653: x4 = xb[3];
654: x5 = xb[4];
655: x6 = xb[5];
656: x7 = xb[6];
657: ib = aj + *ai;
658: jmin = 0;
659: nonzerorow += (n > 0);
660: if (*ib == i) { /* (diag of A)*x */
661: z[7 * i] += v[0] * x1 + v[7] * x2 + v[14] * x3 + v[21] * x4 + v[28] * x5 + v[35] * x6 + v[42] * x7;
662: 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;
663: 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;
664: 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;
665: 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;
666: 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;
667: 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;
668: v += 49;
669: jmin++;
670: }
671: PetscPrefetchBlock(ib + jmin + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
672: PetscPrefetchBlock(v + 49 * n, 49 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
673: for (j = jmin; j < n; j++) {
674: /* (strict lower triangular part of A)*x */
675: cval = ib[j] * 7;
676: z[cval] += v[0] * x1 + v[1] * x2 + v[2] * x3 + v[3] * x4 + v[4] * x5 + v[5] * x6 + v[6] * x7;
677: z[cval + 1] += v[7] * x1 + v[8] * x2 + v[9] * x3 + v[10] * x4 + v[11] * x5 + v[12] * x6 + v[13] * x7;
678: z[cval + 2] += v[14] * x1 + v[15] * x2 + v[16] * x3 + v[17] * x4 + v[18] * x5 + v[19] * x6 + v[20] * x7;
679: z[cval + 3] += v[21] * x1 + v[22] * x2 + v[23] * x3 + v[24] * x4 + v[25] * x5 + v[26] * x6 + v[27] * x7;
680: z[cval + 4] += v[28] * x1 + v[29] * x2 + v[30] * x3 + v[31] * x4 + v[32] * x5 + v[33] * x6 + v[34] * x7;
681: z[cval + 5] += v[35] * x1 + v[36] * x2 + v[37] * x3 + v[38] * x4 + v[39] * x5 + v[40] * x6 + v[41] * x7;
682: z[cval + 6] += v[42] * x1 + v[43] * x2 + v[44] * x3 + v[45] * x4 + v[46] * x5 + v[47] * x6 + v[48] * x7;
683: /* (strict upper triangular part of A)*x */
684: 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];
685: 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];
686: 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];
687: 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];
688: 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];
689: 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];
690: 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];
691: v += 49;
692: }
693: xb += 7;
694: ai++;
695: }
696: PetscCall(VecRestoreArrayRead(xx, &x));
697: PetscCall(VecRestoreArray(zz, &z));
698: PetscCall(PetscLogFlops(98.0 * (a->nz * 2.0 - nonzerorow) - nonzerorow));
699: PetscFunctionReturn(PETSC_SUCCESS);
700: }
702: /*
703: This will not work with MatScalar == float because it calls the BLAS
704: */
705: PetscErrorCode MatMult_SeqSBAIJ_N(Mat A, Vec xx, Vec zz)
706: {
707: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
708: PetscScalar *z, *z_ptr, *zb, *work, *workt, zero = 0.0;
709: const PetscScalar *x, *x_ptr, *xb;
710: const MatScalar *v;
711: PetscInt mbs = a->mbs, i, bs = A->rmap->bs, j, n, bs2 = a->bs2, ncols, k;
712: const PetscInt *idx, *aj, *ii;
713: PetscInt nonzerorow = 0;
715: PetscFunctionBegin;
716: PetscCall(VecSet(zz, zero));
717: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
718: PetscCall(VecGetArrayRead(xx, &x));
719: PetscCall(VecGetArray(zz, &z));
721: x_ptr = x;
722: z_ptr = z;
724: aj = a->j;
725: v = a->a;
726: ii = a->i;
728: if (!a->mult_work) PetscCall(PetscMalloc1(A->rmap->N + 1, &a->mult_work));
729: work = a->mult_work;
731: for (i = 0; i < mbs; i++) {
732: n = ii[1] - ii[0];
733: ncols = n * bs;
734: workt = work;
735: idx = aj + ii[0];
736: nonzerorow += (n > 0);
738: /* upper triangular part */
739: for (j = 0; j < n; j++) {
740: xb = x_ptr + bs * (*idx++);
741: for (k = 0; k < bs; k++) workt[k] = xb[k];
742: workt += bs;
743: }
744: /* z(i*bs:(i+1)*bs-1) += A(i,:)*x */
745: PetscKernel_w_gets_w_plus_Ar_times_v(bs, ncols, work, v, z);
747: /* strict lower triangular part */
748: idx = aj + ii[0];
749: if (n && *idx == i) {
750: ncols -= bs;
751: v += bs2;
752: idx++;
753: n--;
754: }
756: if (ncols > 0) {
757: workt = work;
758: PetscCall(PetscArrayzero(workt, ncols));
759: PetscKernel_w_gets_w_plus_trans_Ar_times_v(bs, ncols, x, v, workt);
760: for (j = 0; j < n; j++) {
761: zb = z_ptr + bs * (*idx++);
762: for (k = 0; k < bs; k++) zb[k] += workt[k];
763: workt += bs;
764: }
765: }
766: x += bs;
767: v += n * bs2;
768: z += bs;
769: ii++;
770: }
772: PetscCall(VecRestoreArrayRead(xx, &x));
773: PetscCall(VecRestoreArray(zz, &z));
774: PetscCall(PetscLogFlops(2.0 * (a->nz * 2.0 - nonzerorow) * bs2 - nonzerorow));
775: PetscFunctionReturn(PETSC_SUCCESS);
776: }
778: PetscErrorCode MatMultAdd_SeqSBAIJ_1(Mat A, Vec xx, Vec yy, Vec zz)
779: {
780: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
781: PetscScalar *z, x1;
782: const PetscScalar *x, *xb;
783: const MatScalar *v;
784: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
785: const PetscInt *aj = a->j, *ai = a->i, *ib;
786: PetscInt nonzerorow = 0;
787: const int aconj = PetscDefined(USE_COMPLEX) && A->hermitian == PETSC_BOOL3_TRUE ? 1 : 0;
789: PetscFunctionBegin;
790: PetscCall(VecCopy(yy, zz));
791: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
792: PetscCall(VecGetArrayRead(xx, &x));
793: PetscCall(VecGetArray(zz, &z));
794: v = a->a;
795: xb = x;
797: for (i = 0; i < mbs; i++) {
798: n = ai[1] - ai[0]; /* length of i_th row of A */
799: x1 = xb[0];
800: ib = aj + *ai;
801: jmin = 0;
802: nonzerorow += (n > 0);
803: if (n && *ib == i) { /* (diag of A)*x */
804: z[i] += *v++ * x[*ib++];
805: jmin++;
806: }
807: if (aconj) {
808: for (j = jmin; j < n; j++) {
809: cval = *ib;
810: z[cval] += PetscConj(*v) * x1; /* (strict lower triangular part of A)*x */
811: z[i] += *v++ * x[*ib++]; /* (strict upper triangular part of A)*x */
812: }
813: } else {
814: for (j = jmin; j < n; j++) {
815: cval = *ib;
816: z[cval] += *v * x1; /* (strict lower triangular part of A)*x */
817: z[i] += *v++ * x[*ib++]; /* (strict upper triangular part of A)*x */
818: }
819: }
820: xb++;
821: ai++;
822: }
824: PetscCall(VecRestoreArrayRead(xx, &x));
825: PetscCall(VecRestoreArray(zz, &z));
827: PetscCall(PetscLogFlops(2.0 * (a->nz * 2.0 - nonzerorow)));
828: PetscFunctionReturn(PETSC_SUCCESS);
829: }
831: PetscErrorCode MatMultAdd_SeqSBAIJ_2(Mat A, Vec xx, Vec yy, Vec zz)
832: {
833: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
834: PetscScalar *z, x1, x2;
835: const PetscScalar *x, *xb;
836: const MatScalar *v;
837: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
838: const PetscInt *aj = a->j, *ai = a->i, *ib;
839: PetscInt nonzerorow = 0;
841: PetscFunctionBegin;
842: PetscCall(VecCopy(yy, zz));
843: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
844: PetscCall(VecGetArrayRead(xx, &x));
845: PetscCall(VecGetArray(zz, &z));
847: v = a->a;
848: xb = x;
850: for (i = 0; i < mbs; i++) {
851: n = ai[1] - ai[0]; /* length of i_th block row of A */
852: x1 = xb[0];
853: x2 = xb[1];
854: ib = aj + *ai;
855: jmin = 0;
856: nonzerorow += (n > 0);
857: if (n && *ib == i) { /* (diag of A)*x */
858: z[2 * i] += v[0] * x1 + v[2] * x2;
859: z[2 * i + 1] += v[2] * x1 + v[3] * x2;
860: v += 4;
861: jmin++;
862: }
863: PetscPrefetchBlock(ib + jmin + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
864: PetscPrefetchBlock(v + 4 * n, 4 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
865: for (j = jmin; j < n; j++) {
866: /* (strict lower triangular part of A)*x */
867: cval = ib[j] * 2;
868: z[cval] += v[0] * x1 + v[1] * x2;
869: z[cval + 1] += v[2] * x1 + v[3] * x2;
870: /* (strict upper triangular part of A)*x */
871: z[2 * i] += v[0] * x[cval] + v[2] * x[cval + 1];
872: z[2 * i + 1] += v[1] * x[cval] + v[3] * x[cval + 1];
873: v += 4;
874: }
875: xb += 2;
876: ai++;
877: }
878: PetscCall(VecRestoreArrayRead(xx, &x));
879: PetscCall(VecRestoreArray(zz, &z));
881: PetscCall(PetscLogFlops(8.0 * (a->nz * 2.0 - nonzerorow)));
882: PetscFunctionReturn(PETSC_SUCCESS);
883: }
885: PetscErrorCode MatMultAdd_SeqSBAIJ_3(Mat A, Vec xx, Vec yy, Vec zz)
886: {
887: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
888: PetscScalar *z, x1, x2, x3;
889: const PetscScalar *x, *xb;
890: const MatScalar *v;
891: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
892: const PetscInt *aj = a->j, *ai = a->i, *ib;
893: PetscInt nonzerorow = 0;
895: PetscFunctionBegin;
896: PetscCall(VecCopy(yy, zz));
897: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
898: PetscCall(VecGetArrayRead(xx, &x));
899: PetscCall(VecGetArray(zz, &z));
901: v = a->a;
902: xb = x;
904: for (i = 0; i < mbs; i++) {
905: n = ai[1] - ai[0]; /* length of i_th block row of A */
906: x1 = xb[0];
907: x2 = xb[1];
908: x3 = xb[2];
909: ib = aj + *ai;
910: jmin = 0;
911: nonzerorow += (n > 0);
912: if (n && *ib == i) { /* (diag of A)*x */
913: z[3 * i] += v[0] * x1 + v[3] * x2 + v[6] * x3;
914: z[3 * i + 1] += v[3] * x1 + v[4] * x2 + v[7] * x3;
915: z[3 * i + 2] += v[6] * x1 + v[7] * x2 + v[8] * x3;
916: v += 9;
917: jmin++;
918: }
919: PetscPrefetchBlock(ib + jmin + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
920: PetscPrefetchBlock(v + 9 * n, 9 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
921: for (j = jmin; j < n; j++) {
922: /* (strict lower triangular part of A)*x */
923: cval = ib[j] * 3;
924: z[cval] += v[0] * x1 + v[1] * x2 + v[2] * x3;
925: z[cval + 1] += v[3] * x1 + v[4] * x2 + v[5] * x3;
926: z[cval + 2] += v[6] * x1 + v[7] * x2 + v[8] * x3;
927: /* (strict upper triangular part of A)*x */
928: z[3 * i] += v[0] * x[cval] + v[3] * x[cval + 1] + v[6] * x[cval + 2];
929: z[3 * i + 1] += v[1] * x[cval] + v[4] * x[cval + 1] + v[7] * x[cval + 2];
930: z[3 * i + 2] += v[2] * x[cval] + v[5] * x[cval + 1] + v[8] * x[cval + 2];
931: v += 9;
932: }
933: xb += 3;
934: ai++;
935: }
937: PetscCall(VecRestoreArrayRead(xx, &x));
938: PetscCall(VecRestoreArray(zz, &z));
940: PetscCall(PetscLogFlops(18.0 * (a->nz * 2.0 - nonzerorow)));
941: PetscFunctionReturn(PETSC_SUCCESS);
942: }
944: PetscErrorCode MatMultAdd_SeqSBAIJ_4(Mat A, Vec xx, Vec yy, Vec zz)
945: {
946: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
947: PetscScalar *z, x1, x2, x3, x4;
948: const PetscScalar *x, *xb;
949: const MatScalar *v;
950: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
951: const PetscInt *aj = a->j, *ai = a->i, *ib;
952: PetscInt nonzerorow = 0;
954: PetscFunctionBegin;
955: PetscCall(VecCopy(yy, zz));
956: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
957: PetscCall(VecGetArrayRead(xx, &x));
958: PetscCall(VecGetArray(zz, &z));
960: v = a->a;
961: xb = x;
963: for (i = 0; i < mbs; i++) {
964: n = ai[1] - ai[0]; /* length of i_th block row of A */
965: x1 = xb[0];
966: x2 = xb[1];
967: x3 = xb[2];
968: x4 = xb[3];
969: ib = aj + *ai;
970: jmin = 0;
971: nonzerorow += (n > 0);
972: if (n && *ib == i) { /* (diag of A)*x */
973: z[4 * i] += v[0] * x1 + v[4] * x2 + v[8] * x3 + v[12] * x4;
974: z[4 * i + 1] += v[4] * x1 + v[5] * x2 + v[9] * x3 + v[13] * x4;
975: z[4 * i + 2] += v[8] * x1 + v[9] * x2 + v[10] * x3 + v[14] * x4;
976: z[4 * i + 3] += v[12] * x1 + v[13] * x2 + v[14] * x3 + v[15] * x4;
977: v += 16;
978: jmin++;
979: }
980: PetscPrefetchBlock(ib + jmin + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
981: PetscPrefetchBlock(v + 16 * n, 16 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
982: for (j = jmin; j < n; j++) {
983: /* (strict lower triangular part of A)*x */
984: cval = ib[j] * 4;
985: z[cval] += v[0] * x1 + v[1] * x2 + v[2] * x3 + v[3] * x4;
986: z[cval + 1] += v[4] * x1 + v[5] * x2 + v[6] * x3 + v[7] * x4;
987: z[cval + 2] += v[8] * x1 + v[9] * x2 + v[10] * x3 + v[11] * x4;
988: z[cval + 3] += v[12] * x1 + v[13] * x2 + v[14] * x3 + v[15] * x4;
989: /* (strict upper triangular part of A)*x */
990: z[4 * i] += v[0] * x[cval] + v[4] * x[cval + 1] + v[8] * x[cval + 2] + v[12] * x[cval + 3];
991: z[4 * i + 1] += v[1] * x[cval] + v[5] * x[cval + 1] + v[9] * x[cval + 2] + v[13] * x[cval + 3];
992: z[4 * i + 2] += v[2] * x[cval] + v[6] * x[cval + 1] + v[10] * x[cval + 2] + v[14] * x[cval + 3];
993: z[4 * i + 3] += v[3] * x[cval] + v[7] * x[cval + 1] + v[11] * x[cval + 2] + v[15] * x[cval + 3];
994: v += 16;
995: }
996: xb += 4;
997: ai++;
998: }
1000: PetscCall(VecRestoreArrayRead(xx, &x));
1001: PetscCall(VecRestoreArray(zz, &z));
1003: PetscCall(PetscLogFlops(32.0 * (a->nz * 2.0 - nonzerorow)));
1004: PetscFunctionReturn(PETSC_SUCCESS);
1005: }
1007: PetscErrorCode MatMultAdd_SeqSBAIJ_5(Mat A, Vec xx, Vec yy, Vec zz)
1008: {
1009: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1010: PetscScalar *z, x1, x2, x3, x4, x5;
1011: const PetscScalar *x, *xb;
1012: const MatScalar *v;
1013: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
1014: const PetscInt *aj = a->j, *ai = a->i, *ib;
1015: PetscInt nonzerorow = 0;
1017: PetscFunctionBegin;
1018: PetscCall(VecCopy(yy, zz));
1019: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
1020: PetscCall(VecGetArrayRead(xx, &x));
1021: PetscCall(VecGetArray(zz, &z));
1023: v = a->a;
1024: xb = x;
1026: for (i = 0; i < mbs; i++) {
1027: n = ai[1] - ai[0]; /* length of i_th block row of A */
1028: x1 = xb[0];
1029: x2 = xb[1];
1030: x3 = xb[2];
1031: x4 = xb[3];
1032: x5 = xb[4];
1033: ib = aj + *ai;
1034: jmin = 0;
1035: nonzerorow += (n > 0);
1036: if (n && *ib == i) { /* (diag of A)*x */
1037: z[5 * i] += v[0] * x1 + v[5] * x2 + v[10] * x3 + v[15] * x4 + v[20] * x5;
1038: z[5 * i + 1] += v[5] * x1 + v[6] * x2 + v[11] * x3 + v[16] * x4 + v[21] * x5;
1039: z[5 * i + 2] += v[10] * x1 + v[11] * x2 + v[12] * x3 + v[17] * x4 + v[22] * x5;
1040: z[5 * i + 3] += v[15] * x1 + v[16] * x2 + v[17] * x3 + v[18] * x4 + v[23] * x5;
1041: z[5 * i + 4] += v[20] * x1 + v[21] * x2 + v[22] * x3 + v[23] * x4 + v[24] * x5;
1042: v += 25;
1043: jmin++;
1044: }
1045: PetscPrefetchBlock(ib + jmin + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
1046: PetscPrefetchBlock(v + 25 * n, 25 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
1047: for (j = jmin; j < n; j++) {
1048: /* (strict lower triangular part of A)*x */
1049: cval = ib[j] * 5;
1050: z[cval] += v[0] * x1 + v[1] * x2 + v[2] * x3 + v[3] * x4 + v[4] * x5;
1051: z[cval + 1] += v[5] * x1 + v[6] * x2 + v[7] * x3 + v[8] * x4 + v[9] * x5;
1052: z[cval + 2] += v[10] * x1 + v[11] * x2 + v[12] * x3 + v[13] * x4 + v[14] * x5;
1053: z[cval + 3] += v[15] * x1 + v[16] * x2 + v[17] * x3 + v[18] * x4 + v[19] * x5;
1054: z[cval + 4] += v[20] * x1 + v[21] * x2 + v[22] * x3 + v[23] * x4 + v[24] * x5;
1055: /* (strict upper triangular part of A)*x */
1056: 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];
1057: 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];
1058: 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];
1059: 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];
1060: 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];
1061: v += 25;
1062: }
1063: xb += 5;
1064: ai++;
1065: }
1067: PetscCall(VecRestoreArrayRead(xx, &x));
1068: PetscCall(VecRestoreArray(zz, &z));
1070: PetscCall(PetscLogFlops(50.0 * (a->nz * 2.0 - nonzerorow)));
1071: PetscFunctionReturn(PETSC_SUCCESS);
1072: }
1074: PetscErrorCode MatMultAdd_SeqSBAIJ_6(Mat A, Vec xx, Vec yy, Vec zz)
1075: {
1076: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1077: PetscScalar *z, x1, x2, x3, x4, x5, x6;
1078: const PetscScalar *x, *xb;
1079: const MatScalar *v;
1080: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
1081: const PetscInt *aj = a->j, *ai = a->i, *ib;
1082: PetscInt nonzerorow = 0;
1084: PetscFunctionBegin;
1085: PetscCall(VecCopy(yy, zz));
1086: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
1087: PetscCall(VecGetArrayRead(xx, &x));
1088: PetscCall(VecGetArray(zz, &z));
1090: v = a->a;
1091: xb = x;
1093: for (i = 0; i < mbs; i++) {
1094: n = ai[1] - ai[0]; /* length of i_th block row of A */
1095: x1 = xb[0];
1096: x2 = xb[1];
1097: x3 = xb[2];
1098: x4 = xb[3];
1099: x5 = xb[4];
1100: x6 = xb[5];
1101: ib = aj + *ai;
1102: jmin = 0;
1103: nonzerorow += (n > 0);
1104: if (n && *ib == i) { /* (diag of A)*x */
1105: z[6 * i] += v[0] * x1 + v[6] * x2 + v[12] * x3 + v[18] * x4 + v[24] * x5 + v[30] * x6;
1106: z[6 * i + 1] += v[6] * x1 + v[7] * x2 + v[13] * x3 + v[19] * x4 + v[25] * x5 + v[31] * x6;
1107: z[6 * i + 2] += v[12] * x1 + v[13] * x2 + v[14] * x3 + v[20] * x4 + v[26] * x5 + v[32] * x6;
1108: z[6 * i + 3] += v[18] * x1 + v[19] * x2 + v[20] * x3 + v[21] * x4 + v[27] * x5 + v[33] * x6;
1109: z[6 * i + 4] += v[24] * x1 + v[25] * x2 + v[26] * x3 + v[27] * x4 + v[28] * x5 + v[34] * x6;
1110: z[6 * i + 5] += v[30] * x1 + v[31] * x2 + v[32] * x3 + v[33] * x4 + v[34] * x5 + v[35] * x6;
1111: v += 36;
1112: jmin++;
1113: }
1114: PetscPrefetchBlock(ib + jmin + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
1115: PetscPrefetchBlock(v + 36 * n, 36 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
1116: for (j = jmin; j < n; j++) {
1117: /* (strict lower triangular part of A)*x */
1118: cval = ib[j] * 6;
1119: z[cval] += v[0] * x1 + v[1] * x2 + v[2] * x3 + v[3] * x4 + v[4] * x5 + v[5] * x6;
1120: z[cval + 1] += v[6] * x1 + v[7] * x2 + v[8] * x3 + v[9] * x4 + v[10] * x5 + v[11] * x6;
1121: z[cval + 2] += v[12] * x1 + v[13] * x2 + v[14] * x3 + v[15] * x4 + v[16] * x5 + v[17] * x6;
1122: z[cval + 3] += v[18] * x1 + v[19] * x2 + v[20] * x3 + v[21] * x4 + v[22] * x5 + v[23] * x6;
1123: z[cval + 4] += v[24] * x1 + v[25] * x2 + v[26] * x3 + v[27] * x4 + v[28] * x5 + v[29] * x6;
1124: z[cval + 5] += v[30] * x1 + v[31] * x2 + v[32] * x3 + v[33] * x4 + v[34] * x5 + v[35] * x6;
1125: /* (strict upper triangular part of A)*x */
1126: 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];
1127: 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];
1128: 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];
1129: 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];
1130: 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];
1131: 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];
1132: v += 36;
1133: }
1134: xb += 6;
1135: ai++;
1136: }
1138: PetscCall(VecRestoreArrayRead(xx, &x));
1139: PetscCall(VecRestoreArray(zz, &z));
1141: PetscCall(PetscLogFlops(72.0 * (a->nz * 2.0 - nonzerorow)));
1142: PetscFunctionReturn(PETSC_SUCCESS);
1143: }
1145: PetscErrorCode MatMultAdd_SeqSBAIJ_7(Mat A, Vec xx, Vec yy, Vec zz)
1146: {
1147: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1148: PetscScalar *z, x1, x2, x3, x4, x5, x6, x7;
1149: const PetscScalar *x, *xb;
1150: const MatScalar *v;
1151: PetscInt mbs = a->mbs, i, n, cval, j, jmin;
1152: const PetscInt *aj = a->j, *ai = a->i, *ib;
1153: PetscInt nonzerorow = 0;
1155: PetscFunctionBegin;
1156: PetscCall(VecCopy(yy, zz));
1157: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
1158: PetscCall(VecGetArrayRead(xx, &x));
1159: PetscCall(VecGetArray(zz, &z));
1161: v = a->a;
1162: xb = x;
1164: for (i = 0; i < mbs; i++) {
1165: n = ai[1] - ai[0]; /* length of i_th block row of A */
1166: x1 = xb[0];
1167: x2 = xb[1];
1168: x3 = xb[2];
1169: x4 = xb[3];
1170: x5 = xb[4];
1171: x6 = xb[5];
1172: x7 = xb[6];
1173: ib = aj + *ai;
1174: jmin = 0;
1175: nonzerorow += (n > 0);
1176: if (n && *ib == i) { /* (diag of A)*x */
1177: z[7 * i] += v[0] * x1 + v[7] * x2 + v[14] * x3 + v[21] * x4 + v[28] * x5 + v[35] * x6 + v[42] * x7;
1178: 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;
1179: 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;
1180: 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;
1181: 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;
1182: 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;
1183: 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;
1184: v += 49;
1185: jmin++;
1186: }
1187: PetscPrefetchBlock(ib + jmin + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
1188: PetscPrefetchBlock(v + 49 * n, 49 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
1189: for (j = jmin; j < n; j++) {
1190: /* (strict lower triangular part of A)*x */
1191: cval = ib[j] * 7;
1192: z[cval] += v[0] * x1 + v[1] * x2 + v[2] * x3 + v[3] * x4 + v[4] * x5 + v[5] * x6 + v[6] * x7;
1193: z[cval + 1] += v[7] * x1 + v[8] * x2 + v[9] * x3 + v[10] * x4 + v[11] * x5 + v[12] * x6 + v[13] * x7;
1194: z[cval + 2] += v[14] * x1 + v[15] * x2 + v[16] * x3 + v[17] * x4 + v[18] * x5 + v[19] * x6 + v[20] * x7;
1195: z[cval + 3] += v[21] * x1 + v[22] * x2 + v[23] * x3 + v[24] * x4 + v[25] * x5 + v[26] * x6 + v[27] * x7;
1196: z[cval + 4] += v[28] * x1 + v[29] * x2 + v[30] * x3 + v[31] * x4 + v[32] * x5 + v[33] * x6 + v[34] * x7;
1197: z[cval + 5] += v[35] * x1 + v[36] * x2 + v[37] * x3 + v[38] * x4 + v[39] * x5 + v[40] * x6 + v[41] * x7;
1198: z[cval + 6] += v[42] * x1 + v[43] * x2 + v[44] * x3 + v[45] * x4 + v[46] * x5 + v[47] * x6 + v[48] * x7;
1199: /* (strict upper triangular part of A)*x */
1200: 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];
1201: 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];
1202: 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];
1203: 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];
1204: 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];
1205: 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];
1206: 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];
1207: v += 49;
1208: }
1209: xb += 7;
1210: ai++;
1211: }
1213: PetscCall(VecRestoreArrayRead(xx, &x));
1214: PetscCall(VecRestoreArray(zz, &z));
1216: PetscCall(PetscLogFlops(98.0 * (a->nz * 2.0 - nonzerorow)));
1217: PetscFunctionReturn(PETSC_SUCCESS);
1218: }
1220: PetscErrorCode MatMultAdd_SeqSBAIJ_N(Mat A, Vec xx, Vec yy, Vec zz)
1221: {
1222: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1223: PetscScalar *z, *z_ptr = NULL, *zb, *work, *workt;
1224: const PetscScalar *x, *x_ptr, *xb;
1225: const MatScalar *v;
1226: PetscInt mbs = a->mbs, i, bs = A->rmap->bs, j, n, bs2 = a->bs2, ncols, k;
1227: const PetscInt *idx, *aj, *ii;
1228: PetscInt nonzerorow = 0;
1230: PetscFunctionBegin;
1231: PetscCall(VecCopy(yy, zz));
1232: if (!a->nz) PetscFunctionReturn(PETSC_SUCCESS);
1233: PetscCall(VecGetArrayRead(xx, &x));
1234: x_ptr = x;
1235: PetscCall(VecGetArray(zz, &z));
1236: z_ptr = z;
1238: aj = a->j;
1239: v = a->a;
1240: ii = a->i;
1242: if (!a->mult_work) PetscCall(PetscMalloc1(A->rmap->n + 1, &a->mult_work));
1243: work = a->mult_work;
1245: for (i = 0; i < mbs; i++) {
1246: n = ii[1] - ii[0];
1247: ncols = n * bs;
1248: workt = work;
1249: idx = aj + ii[0];
1250: nonzerorow += (n > 0);
1252: /* upper triangular part */
1253: for (j = 0; j < n; j++) {
1254: xb = x_ptr + bs * (*idx++);
1255: for (k = 0; k < bs; k++) workt[k] = xb[k];
1256: workt += bs;
1257: }
1258: /* z(i*bs:(i+1)*bs-1) += A(i,:)*x */
1259: PetscKernel_w_gets_w_plus_Ar_times_v(bs, ncols, work, v, z);
1261: /* strict lower triangular part */
1262: idx = aj + ii[0];
1263: if (n && *idx == i) {
1264: ncols -= bs;
1265: v += bs2;
1266: idx++;
1267: n--;
1268: }
1269: if (ncols > 0) {
1270: workt = work;
1271: PetscCall(PetscArrayzero(workt, ncols));
1272: PetscKernel_w_gets_w_plus_trans_Ar_times_v(bs, ncols, x, v, workt);
1273: for (j = 0; j < n; j++) {
1274: zb = z_ptr + bs * (*idx++);
1275: for (k = 0; k < bs; k++) zb[k] += workt[k];
1276: workt += bs;
1277: }
1278: }
1280: x += bs;
1281: v += n * bs2;
1282: z += bs;
1283: ii++;
1284: }
1286: PetscCall(VecRestoreArrayRead(xx, &x));
1287: PetscCall(VecRestoreArray(zz, &z));
1289: PetscCall(PetscLogFlops(2.0 * bs2 * (a->nz * 2.0 - nonzerorow)));
1290: PetscFunctionReturn(PETSC_SUCCESS);
1291: }
1293: PetscErrorCode MatScale_SeqSBAIJ(Mat inA, PetscScalar alpha)
1294: {
1295: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)inA->data;
1296: PetscScalar oalpha = alpha;
1297: PetscBLASInt one = 1, totalnz;
1299: PetscFunctionBegin;
1300: PetscCall(PetscBLASIntCast(a->bs2 * a->nz, &totalnz));
1301: PetscCallBLAS("BLASscal", BLASscal_(&totalnz, &oalpha, a->a, &one));
1302: PetscCall(PetscLogFlops(totalnz));
1303: PetscFunctionReturn(PETSC_SUCCESS);
1304: }
1306: PetscErrorCode MatNorm_SeqSBAIJ(Mat A, NormType type, PetscReal *norm)
1307: {
1308: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1309: const MatScalar *v = a->a;
1310: PetscReal sum_diag = 0.0, sum_off = 0.0, *sum;
1311: PetscInt i, j, k, bs = A->rmap->bs, bs2 = a->bs2, k1, mbs = a->mbs, jmin, jmax, nexti, ik, *jl, *il;
1312: const PetscInt *aj = a->j, *col;
1314: PetscFunctionBegin;
1315: if (!a->nz) {
1316: *norm = 0.0;
1317: PetscFunctionReturn(PETSC_SUCCESS);
1318: }
1319: if (type == NORM_FROBENIUS) {
1320: for (k = 0; k < mbs; k++) {
1321: jmin = a->i[k];
1322: jmax = a->i[k + 1];
1323: col = aj + jmin;
1324: if (jmax - jmin > 0 && *col == k) { /* diagonal block */
1325: for (i = 0; i < bs2; i++) {
1326: sum_diag += PetscRealPart(PetscConj(*v) * (*v));
1327: v++;
1328: }
1329: jmin++;
1330: }
1331: for (j = jmin; j < jmax; j++) { /* off-diagonal blocks */
1332: for (i = 0; i < bs2; i++) {
1333: sum_off += PetscRealPart(PetscConj(*v) * (*v));
1334: v++;
1335: }
1336: }
1337: }
1338: *norm = PetscSqrtReal(sum_diag + 2 * sum_off);
1339: PetscCall(PetscLogFlops(2.0 * bs2 * a->nz));
1340: } else if (type == NORM_INFINITY || type == NORM_1) { /* maximum row/column sum */
1341: PetscCall(PetscMalloc3(bs, &sum, mbs, &il, mbs, &jl));
1342: for (i = 0; i < mbs; i++) jl[i] = mbs;
1343: il[0] = 0;
1345: *norm = 0.0;
1346: for (k = 0; k < mbs; k++) { /* k_th block row */
1347: for (j = 0; j < bs; j++) sum[j] = 0.0;
1348: /*-- col sum --*/
1349: i = jl[k]; /* first |A(i,k)| to be added */
1350: /* jl[k]=i: first nonzero element in row i for submatrix A(1:k,k:n) (active window)
1351: at step k */
1352: while (i < mbs) {
1353: nexti = jl[i]; /* next block row to be added */
1354: ik = il[i]; /* block index of A(i,k) in the array a */
1355: for (j = 0; j < bs; j++) {
1356: v = a->a + ik * bs2 + j * bs;
1357: for (k1 = 0; k1 < bs; k1++) {
1358: sum[j] += PetscAbsScalar(*v);
1359: v++;
1360: }
1361: }
1362: /* update il, jl */
1363: jmin = ik + 1; /* block index of array a: points to the next nonzero of A in row i */
1364: jmax = a->i[i + 1];
1365: if (jmin < jmax) {
1366: il[i] = jmin;
1367: j = a->j[jmin];
1368: jl[i] = jl[j];
1369: jl[j] = i;
1370: }
1371: i = nexti;
1372: }
1373: /*-- row sum --*/
1374: jmin = a->i[k];
1375: jmax = a->i[k + 1];
1376: for (i = jmin; i < jmax; i++) {
1377: for (j = 0; j < bs; j++) {
1378: v = a->a + i * bs2 + j;
1379: for (k1 = 0; k1 < bs; k1++) {
1380: sum[j] += PetscAbsScalar(*v);
1381: v += bs;
1382: }
1383: }
1384: }
1385: /* add k_th block row to il, jl */
1386: col = aj + jmin;
1387: if (jmax - jmin > 0 && *col == k) jmin++;
1388: if (jmin < jmax) {
1389: il[k] = jmin;
1390: j = a->j[jmin];
1391: jl[k] = jl[j];
1392: jl[j] = k;
1393: }
1394: for (j = 0; j < bs; j++) {
1395: if (sum[j] > *norm) *norm = sum[j];
1396: }
1397: }
1398: PetscCall(PetscFree3(sum, il, jl));
1399: PetscCall(PetscLogFlops(PetscMax(mbs * a->nz - 1, 0)));
1400: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for this norm yet");
1401: PetscFunctionReturn(PETSC_SUCCESS);
1402: }
1404: PetscErrorCode MatEqual_SeqSBAIJ(Mat A, Mat B, PetscBool *flg)
1405: {
1406: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data, *b = (Mat_SeqSBAIJ *)B->data;
1408: PetscFunctionBegin;
1409: /* If the matrix/block dimensions are not equal, or no of nonzeros or shift */
1410: if ((A->rmap->N != B->rmap->N) || (A->cmap->n != B->cmap->n) || (A->rmap->bs != B->rmap->bs) || (a->nz != b->nz)) {
1411: *flg = PETSC_FALSE;
1412: PetscFunctionReturn(PETSC_SUCCESS);
1413: }
1415: /* if the a->i are the same */
1416: PetscCall(PetscArraycmp(a->i, b->i, a->mbs + 1, flg));
1417: if (!*flg) PetscFunctionReturn(PETSC_SUCCESS);
1419: /* if a->j are the same */
1420: PetscCall(PetscArraycmp(a->j, b->j, a->nz, flg));
1421: if (!*flg) PetscFunctionReturn(PETSC_SUCCESS);
1423: /* if a->a are the same */
1424: PetscCall(PetscArraycmp(a->a, b->a, a->nz * A->rmap->bs * A->rmap->bs, flg));
1425: PetscFunctionReturn(PETSC_SUCCESS);
1426: }
1428: PetscErrorCode MatGetDiagonal_SeqSBAIJ(Mat A, Vec v)
1429: {
1430: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1431: PetscInt n;
1432: const PetscInt bs = A->rmap->bs, ambs = a->mbs, bs2 = a->bs2;
1433: PetscScalar *x;
1434: const MatScalar *aa = a->a, *aa_j;
1435: const PetscInt *ai = a->i, *adiag;
1436: PetscBool diagDense;
1438: PetscFunctionBegin;
1439: PetscCheck(!A->factortype || bs <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix with bs>1");
1440: PetscCall(MatGetDiagonalMarkers_SeqSBAIJ(A, &adiag, &diagDense));
1441: if (A->factortype == MAT_FACTOR_CHOLESKY || A->factortype == MAT_FACTOR_ICC) {
1442: PetscCall(VecGetArrayWrite(v, &x));
1443: for (PetscInt i = 0; i < ambs; i++) x[i] = 1.0 / aa[adiag[i]];
1444: PetscCall(VecRestoreArrayWrite(v, &x));
1445: PetscFunctionReturn(PETSC_SUCCESS);
1446: }
1448: PetscCall(VecGetLocalSize(v, &n));
1449: PetscCheck(n == A->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector");
1450: PetscCall(VecGetArrayWrite(v, &x));
1452: if (diagDense) {
1453: for (PetscInt i = 0, row = 0; i < ambs; i++) {
1454: aa_j = aa + adiag[i] * bs2;
1455: for (PetscInt k = 0; k < bs2; k += (bs + 1)) x[row++] = aa_j[k];
1456: }
1457: } else {
1458: for (PetscInt i = 0, row = 0; i < ambs; i++) {
1459: const PetscInt j = adiag[i];
1461: if (j != ai[i + 1]) {
1462: aa_j = aa + j * bs2;
1463: for (PetscInt k = 0; k < bs2; k += (bs + 1)) x[row++] = aa_j[k];
1464: } else {
1465: for (PetscInt k = 0; k < bs; k++) x[row++] = 0.0;
1466: }
1467: }
1468: }
1469: PetscCall(VecRestoreArrayWrite(v, &x));
1470: PetscFunctionReturn(PETSC_SUCCESS);
1471: }
1473: PetscErrorCode MatDiagonalScale_SeqSBAIJ(Mat A, Vec ll, Vec rr)
1474: {
1475: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1476: const PetscScalar *l = NULL;
1477: MatScalar *aa = a->a;
1478: PetscInt lm, m = A->rmap->N, mbs = a->mbs, bs = A->rmap->bs, bs2 = a->bs2;
1479: const PetscInt *ai = a->i, *aj = a->j;
1481: PetscFunctionBegin;
1482: if (ll != rr) {
1483: Mat_SeqBAIJ *b;
1484: Mat B;
1485: const PetscScalar *r = NULL;
1486: PetscInt *browlengths, *browstart, *bj;
1487: MatScalar *ba;
1488: PetscInt n = A->cmap->N;
1489: PetscBool hermitian = (PetscBool)(PetscDefined(USE_COMPLEX) && A->hermitian == PETSC_BOOL3_TRUE);
1491: if (ll) {
1492: PetscCall(VecGetLocalSize(ll, &lm));
1493: PetscCheck(lm == m, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Left scaling vector wrong length");
1494: }
1495: if (rr) {
1496: PetscInt rn;
1498: PetscCall(VecGetLocalSize(rr, &rn));
1499: PetscCheck(rn == n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Right scaling vector wrong length");
1500: }
1501: if (ll) PetscCall(VecGetArrayRead(ll, &l));
1502: if (rr) PetscCall(VecGetArrayRead(rr, &r));
1503: PetscCall(PetscCalloc1(mbs, &browlengths));
1504: PetscCall(PetscMalloc1(mbs, &browstart));
1505: for (PetscInt i = 0; i < mbs; i++) {
1506: for (PetscInt k = ai[i]; k < ai[i + 1]; k++) {
1507: browlengths[i]++;
1508: if (aj[k] != i) browlengths[aj[k]]++;
1509: }
1510: }
1511: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
1512: PetscCall(MatSetSizes(B, m, n, m, n));
1513: PetscCall(MatSetType(B, MATSEQBAIJ));
1514: PetscCall(MatSeqBAIJSetPreallocation(B, bs, 0, browlengths));
1515: b = (Mat_SeqBAIJ *)B->data;
1516: ba = b->a;
1517: bj = b->j;
1518: for (PetscInt i = 0; i < mbs; i++) {
1519: b->ilen[i] = browlengths[i];
1520: browstart[i] = b->i[i];
1521: }
1522: PetscCall(PetscFree(browlengths));
1523: for (PetscInt i = 0; i < mbs; i++) {
1524: for (PetscInt k = ai[i]; k < ai[i + 1]; k++) {
1525: const PetscInt j = aj[k];
1526: const MatScalar *av = aa + k * bs2;
1527: MatScalar *v = ba + browstart[i] * bs2;
1528: const PetscScalar *li = PetscSafePointerPlusOffset(l, i * bs), *ri = PetscSafePointerPlusOffset(r, j * bs);
1530: bj[browstart[i]++] = j;
1531: for (PetscInt col = 0; col < bs; col++) {
1532: const PetscScalar x = r != NULL ? ri[col] : 1.0;
1534: for (PetscInt row = 0; row < bs; row++) v[col * bs + row] = av[col * bs + row] * (l != NULL ? li[row] : 1.0) * x;
1535: }
1536: if (j != i) {
1537: MatScalar *v = ba + browstart[j] * bs2;
1538: const PetscScalar *li = PetscSafePointerPlusOffset(l, j * bs), *ri = PetscSafePointerPlusOffset(r, i * bs);
1540: bj[browstart[j]++] = i;
1541: for (PetscInt col = 0; col < bs; col++) {
1542: const PetscScalar x = r != NULL ? ri[col] : 1.0;
1544: 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;
1545: }
1546: }
1547: }
1548: }
1549: PetscCall(PetscFree(browstart));
1550: if (ll) PetscCall(VecRestoreArrayRead(ll, &l));
1551: if (rr) PetscCall(VecRestoreArrayRead(rr, &r));
1552: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
1553: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
1554: PetscCall(PetscLogFlops(((ll ? 1.0 : 0.0) + (rr ? 1.0 : 0.0)) * b->nz * bs2));
1555: B->symmetric = A->symmetric;
1556: B->structurally_symmetric = A->structurally_symmetric;
1557: B->hermitian = A->hermitian;
1558: PetscCall(MatHeaderReplace(A, &B));
1559: PetscFunctionReturn(PETSC_SUCCESS);
1560: }
1561: if (!ll) PetscFunctionReturn(PETSC_SUCCESS);
1562: PetscCall(VecGetArrayRead(ll, &l));
1563: PetscCall(VecGetLocalSize(ll, &lm));
1564: PetscCheck(lm == m, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Left scaling vector wrong length");
1565: for (PetscInt i = 0; i < mbs; i++) { /* for each block row */
1566: const PetscScalar *li = l + i * bs;
1567: MatScalar *v = aa + bs2 * ai[i];
1569: for (PetscInt j = 0; j < ai[i + 1] - ai[i]; j++) { /* for each block */
1570: const PetscScalar *ri = l + bs * aj[ai[i] + j];
1572: for (PetscInt k = 0; k < bs; k++) {
1573: for (PetscInt row = 0; row < bs; row++) (*v++) *= li[row] * ri[k];
1574: }
1575: }
1576: }
1577: PetscCall(VecRestoreArrayRead(ll, &l));
1578: PetscCall(PetscLogFlops(2.0 * a->nz * bs2));
1579: PetscFunctionReturn(PETSC_SUCCESS);
1580: }
1582: PetscErrorCode MatGetInfo_SeqSBAIJ(Mat A, MatInfoType flag, MatInfo *info)
1583: {
1584: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1586: PetscFunctionBegin;
1587: info->block_size = a->bs2;
1588: info->nz_allocated = a->bs2 * a->maxnz; /*num. of nonzeros in upper triangular part */
1589: info->nz_used = a->bs2 * a->nz; /*num. of nonzeros in upper triangular part */
1590: info->nz_unneeded = info->nz_allocated - info->nz_used;
1591: info->assemblies = A->num_ass;
1592: info->mallocs = A->info.mallocs;
1593: info->memory = 0; /* REVIEW ME */
1594: if (A->factortype) {
1595: info->fill_ratio_given = A->info.fill_ratio_given;
1596: info->fill_ratio_needed = A->info.fill_ratio_needed;
1597: info->factor_mallocs = A->info.factor_mallocs;
1598: } else {
1599: info->fill_ratio_given = 0;
1600: info->fill_ratio_needed = 0;
1601: info->factor_mallocs = 0;
1602: }
1603: PetscFunctionReturn(PETSC_SUCCESS);
1604: }
1606: PetscErrorCode MatZeroEntries_SeqSBAIJ(Mat A)
1607: {
1608: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1610: PetscFunctionBegin;
1611: PetscCall(PetscArrayzero(a->a, a->bs2 * a->i[a->mbs]));
1612: PetscFunctionReturn(PETSC_SUCCESS);
1613: }
1615: PetscErrorCode MatGetRowMaxAbs_SeqSBAIJ(Mat A, Vec v, PetscInt idx[])
1616: {
1617: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1618: PetscInt i, j, n, row, col, bs, mbs;
1619: const PetscInt *ai, *aj;
1620: PetscReal atmp;
1621: const MatScalar *aa;
1622: PetscScalar *x;
1623: PetscInt ncols, brow, bcol, krow, kcol;
1625: PetscFunctionBegin;
1626: PetscCheck(!idx, PETSC_COMM_SELF, PETSC_ERR_SUP, "Send email to petsc-maint@mcs.anl.gov");
1627: PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
1628: bs = A->rmap->bs;
1629: aa = a->a;
1630: ai = a->i;
1631: aj = a->j;
1632: mbs = a->mbs;
1634: PetscCall(VecSet(v, 0.0));
1635: PetscCall(VecGetArray(v, &x));
1636: PetscCall(VecGetLocalSize(v, &n));
1637: PetscCheck(n == A->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector");
1638: for (i = 0; i < mbs; i++) {
1639: ncols = ai[1] - ai[0];
1640: ai++;
1641: brow = bs * i;
1642: for (j = 0; j < ncols; j++) {
1643: bcol = bs * (*aj);
1644: for (kcol = 0; kcol < bs; kcol++) {
1645: col = bcol + kcol; /* col index */
1646: for (krow = 0; krow < bs; krow++) {
1647: atmp = PetscAbsScalar(*aa);
1648: aa++;
1649: row = brow + krow; /* row index */
1650: if (PetscRealPart(x[row]) < atmp) x[row] = atmp;
1651: if (*aj > i && PetscRealPart(x[col]) < atmp) x[col] = atmp;
1652: }
1653: }
1654: aj++;
1655: }
1656: }
1657: PetscCall(VecRestoreArray(v, &x));
1658: PetscFunctionReturn(PETSC_SUCCESS);
1659: }
1661: PetscErrorCode MatMatMultSymbolic_SeqSBAIJ_SeqDense(Mat A, Mat B, PetscReal fill, Mat C)
1662: {
1663: PetscFunctionBegin;
1664: PetscCall(MatMatMultSymbolic_SeqDense_SeqDense(A, B, 0.0, C));
1665: C->ops->matmultnumeric = MatMatMultNumeric_SeqSBAIJ_SeqDense;
1666: PetscFunctionReturn(PETSC_SUCCESS);
1667: }
1669: static PetscErrorCode MatMatMult_SeqSBAIJ_1_Private(Mat A, PetscScalar *b, PetscInt bm, PetscScalar *c, PetscInt cm, PetscInt cn)
1670: {
1671: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1672: PetscScalar *z = c;
1673: const PetscScalar *xb;
1674: PetscScalar x1;
1675: const MatScalar *v = a->a, *vv;
1676: PetscInt mbs = a->mbs, i, *idx = a->j, *ii = a->i, j, *jj, n, k;
1677: const int aconj = PetscDefined(USE_COMPLEX) && A->hermitian == PETSC_BOOL3_TRUE ? 1 : 0;
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 + n, 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 + (*idx);
1692: x1 = xb[0 + k * bm];
1693: z[0 + k * cm] += v[0] * x1;
1694: if (*idx != i) c[(*idx) + k * cm] += (aconj ? PetscConj(v[0]) : v[0]) * b[i + k * bm];
1695: v += 1;
1696: ++idx;
1697: }
1698: }
1699: z += 1;
1700: }
1701: PetscFunctionReturn(PETSC_SUCCESS);
1702: }
1704: static PetscErrorCode MatMatMult_SeqSBAIJ_2_Private(Mat A, PetscScalar *b, PetscInt bm, PetscScalar *c, PetscInt cm, PetscInt cn)
1705: {
1706: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1707: PetscScalar *z = c;
1708: const PetscScalar *xb;
1709: PetscScalar x1, x2;
1710: const MatScalar *v = a->a, *vv;
1711: PetscInt mbs = a->mbs, i, *idx = a->j, *ii = a->i, j, *jj, n, k;
1713: PetscFunctionBegin;
1714: for (i = 0; i < mbs; i++) {
1715: n = ii[1] - ii[0];
1716: ii++;
1717: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
1718: PetscPrefetchBlock(v + 4 * n, 4 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
1719: jj = idx;
1720: vv = v;
1721: for (k = 0; k < cn; k++) {
1722: idx = jj;
1723: v = vv;
1724: for (j = 0; j < n; j++) {
1725: xb = b + 2 * (*idx);
1726: x1 = xb[0 + k * bm];
1727: x2 = xb[1 + k * bm];
1728: z[0 + k * cm] += v[0] * x1 + v[2] * x2;
1729: z[1 + k * cm] += v[1] * x1 + v[3] * x2;
1730: if (*idx != i) {
1731: c[2 * (*idx) + 0 + k * cm] += v[0] * b[2 * i + k * bm] + v[1] * b[2 * i + 1 + k * bm];
1732: c[2 * (*idx) + 1 + k * cm] += v[2] * b[2 * i + k * bm] + v[3] * b[2 * i + 1 + k * bm];
1733: }
1734: v += 4;
1735: ++idx;
1736: }
1737: }
1738: z += 2;
1739: }
1740: PetscFunctionReturn(PETSC_SUCCESS);
1741: }
1743: static PetscErrorCode MatMatMult_SeqSBAIJ_3_Private(Mat A, PetscScalar *b, PetscInt bm, PetscScalar *c, PetscInt cm, PetscInt cn)
1744: {
1745: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1746: PetscScalar *z = c;
1747: const PetscScalar *xb;
1748: PetscScalar x1, x2, x3;
1749: const MatScalar *v = a->a, *vv;
1750: PetscInt mbs = a->mbs, i, *idx = a->j, *ii = a->i, j, *jj, n, k;
1752: PetscFunctionBegin;
1753: for (i = 0; i < mbs; i++) {
1754: n = ii[1] - ii[0];
1755: ii++;
1756: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
1757: PetscPrefetchBlock(v + 9 * n, 9 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
1758: jj = idx;
1759: vv = v;
1760: for (k = 0; k < cn; k++) {
1761: idx = jj;
1762: v = vv;
1763: for (j = 0; j < n; j++) {
1764: xb = b + 3 * (*idx);
1765: x1 = xb[0 + k * bm];
1766: x2 = xb[1 + k * bm];
1767: x3 = xb[2 + k * bm];
1768: z[0 + k * cm] += v[0] * x1 + v[3] * x2 + v[6] * x3;
1769: z[1 + k * cm] += v[1] * x1 + v[4] * x2 + v[7] * x3;
1770: z[2 + k * cm] += v[2] * x1 + v[5] * x2 + v[8] * x3;
1771: if (*idx != i) {
1772: 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];
1773: 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];
1774: 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];
1775: }
1776: v += 9;
1777: ++idx;
1778: }
1779: }
1780: z += 3;
1781: }
1782: PetscFunctionReturn(PETSC_SUCCESS);
1783: }
1785: static PetscErrorCode MatMatMult_SeqSBAIJ_4_Private(Mat A, PetscScalar *b, PetscInt bm, PetscScalar *c, PetscInt cm, PetscInt cn)
1786: {
1787: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1788: PetscScalar *z = c;
1789: const PetscScalar *xb;
1790: PetscScalar x1, x2, x3, x4;
1791: const MatScalar *v = a->a, *vv;
1792: PetscInt mbs = a->mbs, i, *idx = a->j, *ii = a->i, j, *jj, n, k;
1794: PetscFunctionBegin;
1795: for (i = 0; i < mbs; i++) {
1796: n = ii[1] - ii[0];
1797: ii++;
1798: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
1799: PetscPrefetchBlock(v + 16 * n, 16 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
1800: jj = idx;
1801: vv = v;
1802: for (k = 0; k < cn; k++) {
1803: idx = jj;
1804: v = vv;
1805: for (j = 0; j < n; j++) {
1806: xb = b + 4 * (*idx);
1807: x1 = xb[0 + k * bm];
1808: x2 = xb[1 + k * bm];
1809: x3 = xb[2 + k * bm];
1810: x4 = xb[3 + k * bm];
1811: z[0 + k * cm] += v[0] * x1 + v[4] * x2 + v[8] * x3 + v[12] * x4;
1812: z[1 + k * cm] += v[1] * x1 + v[5] * x2 + v[9] * x3 + v[13] * x4;
1813: z[2 + k * cm] += v[2] * x1 + v[6] * x2 + v[10] * x3 + v[14] * x4;
1814: z[3 + k * cm] += v[3] * x1 + v[7] * x2 + v[11] * x3 + v[15] * x4;
1815: if (*idx != i) {
1816: 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];
1817: 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];
1818: 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];
1819: 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];
1820: }
1821: v += 16;
1822: ++idx;
1823: }
1824: }
1825: z += 4;
1826: }
1827: PetscFunctionReturn(PETSC_SUCCESS);
1828: }
1830: static PetscErrorCode MatMatMult_SeqSBAIJ_5_Private(Mat A, PetscScalar *b, PetscInt bm, PetscScalar *c, PetscInt cm, PetscInt cn)
1831: {
1832: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1833: PetscScalar *z = c;
1834: const PetscScalar *xb;
1835: PetscScalar x1, x2, x3, x4, x5;
1836: const MatScalar *v = a->a, *vv;
1837: PetscInt mbs = a->mbs, i, *idx = a->j, *ii = a->i, j, *jj, n, k;
1839: PetscFunctionBegin;
1840: for (i = 0; i < mbs; i++) {
1841: n = ii[1] - ii[0];
1842: ii++;
1843: PetscPrefetchBlock(idx + n, n, 0, PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
1844: PetscPrefetchBlock(v + 25 * n, 25 * n, 0, PETSC_PREFETCH_HINT_NTA); /* Entries for the next row */
1845: jj = idx;
1846: vv = v;
1847: for (k = 0; k < cn; k++) {
1848: idx = jj;
1849: v = vv;
1850: for (j = 0; j < n; j++) {
1851: xb = b + 5 * (*idx);
1852: x1 = xb[0 + k * bm];
1853: x2 = xb[1 + k * bm];
1854: x3 = xb[2 + k * bm];
1855: x4 = xb[3 + k * bm];
1856: x5 = xb[4 + k * cm];
1857: z[0 + k * cm] += v[0] * x1 + v[5] * x2 + v[10] * x3 + v[15] * x4 + v[20] * x5;
1858: z[1 + k * cm] += v[1] * x1 + v[6] * x2 + v[11] * x3 + v[16] * x4 + v[21] * x5;
1859: z[2 + k * cm] += v[2] * x1 + v[7] * x2 + v[12] * x3 + v[17] * x4 + v[22] * x5;
1860: z[3 + k * cm] += v[3] * x1 + v[8] * x2 + v[13] * x3 + v[18] * x4 + v[23] * x5;
1861: z[4 + k * cm] += v[4] * x1 + v[9] * x2 + v[14] * x3 + v[19] * x4 + v[24] * x5;
1862: if (*idx != i) {
1863: 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];
1864: 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];
1865: 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];
1866: 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];
1867: 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];
1868: }
1869: v += 25;
1870: ++idx;
1871: }
1872: }
1873: z += 5;
1874: }
1875: PetscFunctionReturn(PETSC_SUCCESS);
1876: }
1878: PetscErrorCode MatMatMultNumeric_SeqSBAIJ_SeqDense(Mat A, Mat B, Mat C)
1879: {
1880: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
1881: Mat_SeqDense *bd = (Mat_SeqDense *)B->data;
1882: Mat_SeqDense *cd = (Mat_SeqDense *)C->data;
1883: PetscInt cm = cd->lda, cn = B->cmap->n, bm = bd->lda;
1884: PetscInt mbs, i, bs = A->rmap->bs, j, n, bs2 = a->bs2;
1885: PetscBLASInt bbs, bcn, bbm, bcm;
1886: PetscScalar *z = NULL;
1887: PetscScalar *c, *b;
1888: const MatScalar *v;
1889: const PetscInt *idx, *ii;
1890: PetscScalar _DOne = 1.0;
1892: PetscFunctionBegin;
1893: if (!cm || !cn) PetscFunctionReturn(PETSC_SUCCESS);
1894: 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);
1895: 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);
1896: 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);
1897: b = bd->v;
1898: PetscCall(MatZeroEntries(C));
1899: PetscCall(MatDenseGetArray(C, &c));
1900: switch (bs) {
1901: case 1:
1902: PetscCall(MatMatMult_SeqSBAIJ_1_Private(A, b, bm, c, cm, cn));
1903: break;
1904: case 2:
1905: PetscCall(MatMatMult_SeqSBAIJ_2_Private(A, b, bm, c, cm, cn));
1906: break;
1907: case 3:
1908: PetscCall(MatMatMult_SeqSBAIJ_3_Private(A, b, bm, c, cm, cn));
1909: break;
1910: case 4:
1911: PetscCall(MatMatMult_SeqSBAIJ_4_Private(A, b, bm, c, cm, cn));
1912: break;
1913: case 5:
1914: PetscCall(MatMatMult_SeqSBAIJ_5_Private(A, b, bm, c, cm, cn));
1915: break;
1916: default: /* block sizes larger than 5 by 5 are handled by BLAS */
1917: PetscCall(PetscBLASIntCast(bs, &bbs));
1918: PetscCall(PetscBLASIntCast(cn, &bcn));
1919: PetscCall(PetscBLASIntCast(bm, &bbm));
1920: PetscCall(PetscBLASIntCast(cm, &bcm));
1921: idx = a->j;
1922: v = a->a;
1923: mbs = a->mbs;
1924: ii = a->i;
1925: z = c;
1926: for (i = 0; i < mbs; i++) {
1927: n = ii[1] - ii[0];
1928: ii++;
1929: for (j = 0; j < n; j++) {
1930: if (*idx != i) PetscCallBLAS("BLASgemm", BLASgemm_("T", "N", &bbs, &bcn, &bbs, &_DOne, v, &bbs, b + bs * i, &bbm, &_DOne, c + bs * (*idx), &bcm));
1931: PetscCallBLAS("BLASgemm", BLASgemm_("N", "N", &bbs, &bcn, &bbs, &_DOne, v, &bbs, b + bs * (*idx++), &bbm, &_DOne, z, &bcm));
1932: v += bs2;
1933: }
1934: z += bs;
1935: }
1936: }
1937: PetscCall(MatDenseRestoreArray(C, &c));
1938: PetscCall(PetscLogFlops((2.0 * (a->nz * 2.0 - a->nonzerorowcnt) * bs2 - a->nonzerorowcnt) * cn));
1939: PetscFunctionReturn(PETSC_SUCCESS);
1940: }