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