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