Actual source code: baij.c
1: /*
2: Defines the basic matrix operations for the BAIJ (compressed row)
3: matrix storage format.
4: */
5: #include <../src/mat/impls/baij/seq/baij.h>
6: #include <petscblaslapack.h>
7: #include <petsc/private/kernels/blockinvert.h>
8: #include <petsc/private/kernels/blockmatmult.h>
10: /* defines MatSetValues_Seq_Hash(), MatAssemblyEnd_Seq_Hash(), MatSetUp_Seq_Hash() */
11: #define TYPE BAIJ
12: #define TYPE_BS
13: #include "../src/mat/impls/aij/seq/seqhashmatsetvalues.h"
14: #undef TYPE_BS
15: #define TYPE_BS _BS
16: #define TYPE_BS_ON
17: #include "../src/mat/impls/aij/seq/seqhashmatsetvalues.h"
18: #undef TYPE_BS
19: #include "../src/mat/impls/aij/seq/seqhashmat.h"
20: #undef TYPE
21: #undef TYPE_BS_ON
23: #if PetscDefined(HAVE_HYPRE)
24: PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat, MatType, MatReuse, Mat *);
25: #endif
27: #if PetscDefined(HAVE_MKL_SPARSE_OPTIMIZE)
28: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqBAIJMKL(Mat, MatType, MatReuse, Mat *);
29: #endif
30: #if PetscDefined(HAVE_LIBXSMM)
31: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqBAIJLIBXSMM(Mat, MatType, MatReuse, Mat *);
32: #endif
33: PETSC_INTERN PetscErrorCode MatConvert_XAIJ_IS(Mat, MatType, MatReuse, Mat *);
35: MatGetDiagonalMarkers(SeqBAIJ, A->rmap->bs)
37: static PetscErrorCode MatGetColumnReductions_SeqBAIJ(Mat A, PetscInt type, PetscReal *reductions)
38: {
39: Mat_SeqBAIJ *a_aij = (Mat_SeqBAIJ *)A->data;
40: PetscInt m, n, ib, jb, bs = A->rmap->bs;
41: MatScalar *a_val = a_aij->a;
43: PetscFunctionBegin;
44: PetscCall(MatGetSize(A, &m, &n));
45: PetscCall(PetscArrayzero(reductions, n));
46: if (type == NORM_2) {
47: for (PetscInt i = a_aij->i[0]; i < a_aij->i[A->rmap->n / bs]; i++) {
48: for (jb = 0; jb < bs; jb++) {
49: for (ib = 0; ib < bs; ib++) {
50: reductions[A->cmap->rstart + a_aij->j[i] * bs + jb] += PetscAbsScalar(*a_val * *a_val);
51: a_val++;
52: }
53: }
54: }
55: } else if (type == NORM_1) {
56: for (PetscInt i = a_aij->i[0]; i < a_aij->i[A->rmap->n / bs]; i++) {
57: for (jb = 0; jb < bs; jb++) {
58: for (ib = 0; ib < bs; ib++) {
59: reductions[A->cmap->rstart + a_aij->j[i] * bs + jb] += PetscAbsScalar(*a_val);
60: a_val++;
61: }
62: }
63: }
64: } else if (type == NORM_INFINITY) {
65: for (PetscInt i = a_aij->i[0]; i < a_aij->i[A->rmap->n / bs]; i++) {
66: for (jb = 0; jb < bs; jb++) {
67: for (ib = 0; ib < bs; ib++) {
68: PetscInt col = A->cmap->rstart + a_aij->j[i] * bs + jb;
69: reductions[col] = PetscMax(PetscAbsScalar(*a_val), reductions[col]);
70: a_val++;
71: }
72: }
73: }
74: } else if (type == REDUCTION_SUM_REALPART || type == REDUCTION_MEAN_REALPART) {
75: for (PetscInt i = a_aij->i[0]; i < a_aij->i[A->rmap->n / bs]; i++) {
76: for (jb = 0; jb < bs; jb++) {
77: for (ib = 0; ib < bs; ib++) {
78: reductions[A->cmap->rstart + a_aij->j[i] * bs + jb] += PetscRealPart(*a_val);
79: a_val++;
80: }
81: }
82: }
83: } else {
84: PetscCheck(type == REDUCTION_SUM_IMAGINARYPART || type == REDUCTION_MEAN_IMAGINARYPART, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Unknown reduction type");
85: for (PetscInt i = a_aij->i[0]; i < a_aij->i[A->rmap->n / bs]; i++) {
86: for (jb = 0; jb < bs; jb++) {
87: for (ib = 0; ib < bs; ib++) {
88: reductions[A->cmap->rstart + a_aij->j[i] * bs + jb] += PetscImaginaryPart(*a_val);
89: a_val++;
90: }
91: }
92: }
93: }
94: if (type == NORM_2) {
95: for (PetscInt i = 0; i < n; i++) reductions[i] = PetscSqrtReal(reductions[i]);
96: } else if (type == REDUCTION_MEAN_REALPART || type == REDUCTION_MEAN_IMAGINARYPART) {
97: for (PetscInt i = 0; i < n; i++) reductions[i] /= m;
98: }
99: PetscFunctionReturn(PETSC_SUCCESS);
100: }
102: static PetscErrorCode MatInvertBlockDiagonal_SeqBAIJ(Mat A, const PetscScalar **values)
103: {
104: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
105: PetscInt i, bs = A->rmap->bs, mbs = a->mbs, ipvt[5], bs2 = bs * bs, *v_pivots;
106: MatScalar *v = a->a, *odiag, *diag, work[25], *v_work;
107: PetscReal shift = 0.0;
108: PetscBool allowzeropivot, zeropivotdetected = PETSC_FALSE;
109: const PetscInt *adiag;
111: PetscFunctionBegin;
112: allowzeropivot = PetscNot(A->erroriffailure);
114: if (a->idiag && a->idiagState == ((PetscObject)A)->state) {
115: if (values) *values = a->idiag;
116: PetscFunctionReturn(PETSC_SUCCESS);
117: }
118: PetscCall(MatGetDiagonalMarkers_SeqBAIJ(A, &adiag, NULL));
119: if (!a->idiag) PetscCall(PetscMalloc1(bs2 * mbs, &a->idiag));
120: diag = a->idiag;
121: if (values) *values = a->idiag;
122: /* factor and invert each block */
123: switch (bs) {
124: case 1:
125: for (i = 0; i < mbs; i++) {
126: odiag = v + 1 * adiag[i];
127: diag[0] = odiag[0];
129: if (PetscAbsScalar(diag[0] + shift) < PETSC_MACHINE_EPSILON) {
130: PetscCheck(allowzeropivot, PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Zero pivot, row %" PetscInt_FMT " pivot value %g tolerance %g", i, (double)PetscAbsScalar(diag[0]), (double)PETSC_MACHINE_EPSILON);
131: A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
132: A->factorerror_zeropivot_value = PetscAbsScalar(diag[0]);
133: A->factorerror_zeropivot_row = i;
134: PetscCall(PetscInfo(A, "Zero pivot, row %" PetscInt_FMT "\n", i));
135: }
137: diag[0] = (PetscScalar)1.0 / (diag[0] + shift);
138: diag += 1;
139: }
140: break;
141: case 2:
142: for (i = 0; i < mbs; i++) {
143: odiag = v + 4 * adiag[i];
144: diag[0] = odiag[0];
145: diag[1] = odiag[1];
146: diag[2] = odiag[2];
147: diag[3] = odiag[3];
148: PetscCall(PetscKernel_A_gets_inverse_A_2(diag, shift, allowzeropivot, &zeropivotdetected));
149: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
150: diag += 4;
151: }
152: break;
153: case 3:
154: for (i = 0; i < mbs; i++) {
155: odiag = v + 9 * adiag[i];
156: diag[0] = odiag[0];
157: diag[1] = odiag[1];
158: diag[2] = odiag[2];
159: diag[3] = odiag[3];
160: diag[4] = odiag[4];
161: diag[5] = odiag[5];
162: diag[6] = odiag[6];
163: diag[7] = odiag[7];
164: diag[8] = odiag[8];
165: PetscCall(PetscKernel_A_gets_inverse_A_3(diag, shift, allowzeropivot, &zeropivotdetected));
166: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
167: diag += 9;
168: }
169: break;
170: case 4:
171: for (i = 0; i < mbs; i++) {
172: odiag = v + 16 * adiag[i];
173: PetscCall(PetscArraycpy(diag, odiag, 16));
174: PetscCall(PetscKernel_A_gets_inverse_A_4(diag, shift, allowzeropivot, &zeropivotdetected));
175: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
176: diag += 16;
177: }
178: break;
179: case 5:
180: for (i = 0; i < mbs; i++) {
181: odiag = v + 25 * adiag[i];
182: PetscCall(PetscArraycpy(diag, odiag, 25));
183: PetscCall(PetscKernel_A_gets_inverse_A_5(diag, ipvt, work, shift, allowzeropivot, &zeropivotdetected));
184: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
185: diag += 25;
186: }
187: break;
188: case 6:
189: for (i = 0; i < mbs; i++) {
190: odiag = v + 36 * adiag[i];
191: PetscCall(PetscArraycpy(diag, odiag, 36));
192: PetscCall(PetscKernel_A_gets_inverse_A_6(diag, shift, allowzeropivot, &zeropivotdetected));
193: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
194: diag += 36;
195: }
196: break;
197: case 7:
198: for (i = 0; i < mbs; i++) {
199: odiag = v + 49 * adiag[i];
200: PetscCall(PetscArraycpy(diag, odiag, 49));
201: PetscCall(PetscKernel_A_gets_inverse_A_7(diag, shift, allowzeropivot, &zeropivotdetected));
202: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
203: diag += 49;
204: }
205: break;
206: default:
207: PetscCall(PetscMalloc2(bs, &v_work, bs, &v_pivots));
208: for (i = 0; i < mbs; i++) {
209: odiag = v + bs2 * adiag[i];
210: PetscCall(PetscArraycpy(diag, odiag, bs2));
211: PetscCall(PetscKernel_A_gets_inverse_A(bs, diag, v_pivots, v_work, allowzeropivot, &zeropivotdetected));
212: if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
213: diag += bs2;
214: }
215: PetscCall(PetscFree2(v_work, v_pivots));
216: }
217: a->idiagState = ((PetscObject)A)->state;
218: PetscFunctionReturn(PETSC_SUCCESS);
219: }
221: static PetscErrorCode MatSOR_SeqBAIJ(Mat A, Vec bb, PetscReal omega, MatSORType flag, PetscReal fshift, PetscInt its, PetscInt lits, Vec xx)
222: {
223: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
224: PetscScalar *x, *work, *w, *workt, *t;
225: const MatScalar *v, *aa = a->a, *idiag;
226: const PetscScalar *b, *xb;
227: PetscScalar s[7], xw[7] = {0}; /* avoid some compilers thinking xw is uninitialized */
228: PetscInt m = a->mbs, i, i2, nz, bs = A->rmap->bs, bs2 = bs * bs, k, j, idx, it;
229: const PetscInt *diag, *ai = a->i, *aj = a->j, *vi;
231: PetscFunctionBegin;
232: its = its * lits;
233: PetscCheck(!(flag & SOR_EISENSTAT), PETSC_COMM_SELF, PETSC_ERR_SUP, "No support yet for Eisenstat");
234: PetscCheck(its > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Relaxation requires global its %" PetscInt_FMT " and local its %" PetscInt_FMT " both positive", its, lits);
235: PetscCheck(!fshift, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for diagonal shift");
236: PetscCheck(omega == 1.0, PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for non-trivial relaxation factor");
237: PetscCheck(!(flag & SOR_APPLY_UPPER) && !(flag & SOR_APPLY_LOWER), PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for applying upper or lower triangular parts");
239: PetscCall(MatInvertBlockDiagonal(A, NULL)); /* a no-op if the cached inverse is still current */
241: if (!m) PetscFunctionReturn(PETSC_SUCCESS);
242: diag = a->diag;
243: idiag = a->idiag;
244: k = PetscMax(A->rmap->n, A->cmap->n);
245: if (!a->mult_work) PetscCall(PetscMalloc1(k + 1, &a->mult_work));
246: if (!a->sor_workt) PetscCall(PetscMalloc1(k, &a->sor_workt));
247: if (!a->sor_work) PetscCall(PetscMalloc1(bs, &a->sor_work));
248: work = a->mult_work;
249: t = a->sor_workt;
250: w = a->sor_work;
252: PetscCall(VecGetArray(xx, &x));
253: PetscCall(VecGetArrayRead(bb, &b));
255: if (flag & SOR_ZERO_INITIAL_GUESS) {
256: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
257: switch (bs) {
258: case 1:
259: PetscKernel_v_gets_A_times_w_1(x, idiag, b);
260: t[0] = b[0];
261: i2 = 1;
262: idiag += 1;
263: for (i = 1; i < m; i++) {
264: v = aa + ai[i];
265: vi = aj + ai[i];
266: nz = diag[i] - ai[i];
267: s[0] = b[i2];
268: for (j = 0; j < nz; j++) {
269: xw[0] = x[vi[j]];
270: PetscKernel_v_gets_v_minus_A_times_w_1(s, (v + j), xw);
271: }
272: t[i2] = s[0];
273: PetscKernel_v_gets_A_times_w_1(xw, idiag, s);
274: x[i2] = xw[0];
275: idiag += 1;
276: i2 += 1;
277: }
278: break;
279: case 2:
280: PetscKernel_v_gets_A_times_w_2(x, idiag, b);
281: t[0] = b[0];
282: t[1] = b[1];
283: i2 = 2;
284: idiag += 4;
285: for (i = 1; i < m; i++) {
286: v = aa + 4 * ai[i];
287: vi = aj + ai[i];
288: nz = diag[i] - ai[i];
289: s[0] = b[i2];
290: s[1] = b[i2 + 1];
291: for (j = 0; j < nz; j++) {
292: idx = 2 * vi[j];
293: it = 4 * j;
294: xw[0] = x[idx];
295: xw[1] = x[1 + idx];
296: PetscKernel_v_gets_v_minus_A_times_w_2(s, (v + it), xw);
297: }
298: t[i2] = s[0];
299: t[i2 + 1] = s[1];
300: PetscKernel_v_gets_A_times_w_2(xw, idiag, s);
301: x[i2] = xw[0];
302: x[i2 + 1] = xw[1];
303: idiag += 4;
304: i2 += 2;
305: }
306: break;
307: case 3:
308: PetscKernel_v_gets_A_times_w_3(x, idiag, b);
309: t[0] = b[0];
310: t[1] = b[1];
311: t[2] = b[2];
312: i2 = 3;
313: idiag += 9;
314: for (i = 1; i < m; i++) {
315: v = aa + 9 * ai[i];
316: vi = aj + ai[i];
317: nz = diag[i] - ai[i];
318: s[0] = b[i2];
319: s[1] = b[i2 + 1];
320: s[2] = b[i2 + 2];
321: while (nz--) {
322: idx = 3 * (*vi++);
323: xw[0] = x[idx];
324: xw[1] = x[1 + idx];
325: xw[2] = x[2 + idx];
326: PetscKernel_v_gets_v_minus_A_times_w_3(s, v, xw);
327: v += 9;
328: }
329: t[i2] = s[0];
330: t[i2 + 1] = s[1];
331: t[i2 + 2] = s[2];
332: PetscKernel_v_gets_A_times_w_3(xw, idiag, s);
333: x[i2] = xw[0];
334: x[i2 + 1] = xw[1];
335: x[i2 + 2] = xw[2];
336: idiag += 9;
337: i2 += 3;
338: }
339: break;
340: case 4:
341: PetscKernel_v_gets_A_times_w_4(x, idiag, b);
342: t[0] = b[0];
343: t[1] = b[1];
344: t[2] = b[2];
345: t[3] = b[3];
346: i2 = 4;
347: idiag += 16;
348: for (i = 1; i < m; i++) {
349: v = aa + 16 * ai[i];
350: vi = aj + ai[i];
351: nz = diag[i] - ai[i];
352: s[0] = b[i2];
353: s[1] = b[i2 + 1];
354: s[2] = b[i2 + 2];
355: s[3] = b[i2 + 3];
356: while (nz--) {
357: idx = 4 * (*vi++);
358: xw[0] = x[idx];
359: xw[1] = x[1 + idx];
360: xw[2] = x[2 + idx];
361: xw[3] = x[3 + idx];
362: PetscKernel_v_gets_v_minus_A_times_w_4(s, v, xw);
363: v += 16;
364: }
365: t[i2] = s[0];
366: t[i2 + 1] = s[1];
367: t[i2 + 2] = s[2];
368: t[i2 + 3] = s[3];
369: PetscKernel_v_gets_A_times_w_4(xw, idiag, s);
370: x[i2] = xw[0];
371: x[i2 + 1] = xw[1];
372: x[i2 + 2] = xw[2];
373: x[i2 + 3] = xw[3];
374: idiag += 16;
375: i2 += 4;
376: }
377: break;
378: case 5:
379: PetscKernel_v_gets_A_times_w_5(x, idiag, b);
380: t[0] = b[0];
381: t[1] = b[1];
382: t[2] = b[2];
383: t[3] = b[3];
384: t[4] = b[4];
385: i2 = 5;
386: idiag += 25;
387: for (i = 1; i < m; i++) {
388: v = aa + 25 * ai[i];
389: vi = aj + ai[i];
390: nz = diag[i] - ai[i];
391: s[0] = b[i2];
392: s[1] = b[i2 + 1];
393: s[2] = b[i2 + 2];
394: s[3] = b[i2 + 3];
395: s[4] = b[i2 + 4];
396: while (nz--) {
397: idx = 5 * (*vi++);
398: xw[0] = x[idx];
399: xw[1] = x[1 + idx];
400: xw[2] = x[2 + idx];
401: xw[3] = x[3 + idx];
402: xw[4] = x[4 + idx];
403: PetscKernel_v_gets_v_minus_A_times_w_5(s, v, xw);
404: v += 25;
405: }
406: t[i2] = s[0];
407: t[i2 + 1] = s[1];
408: t[i2 + 2] = s[2];
409: t[i2 + 3] = s[3];
410: t[i2 + 4] = s[4];
411: PetscKernel_v_gets_A_times_w_5(xw, idiag, s);
412: x[i2] = xw[0];
413: x[i2 + 1] = xw[1];
414: x[i2 + 2] = xw[2];
415: x[i2 + 3] = xw[3];
416: x[i2 + 4] = xw[4];
417: idiag += 25;
418: i2 += 5;
419: }
420: break;
421: case 6:
422: PetscKernel_v_gets_A_times_w_6(x, idiag, b);
423: t[0] = b[0];
424: t[1] = b[1];
425: t[2] = b[2];
426: t[3] = b[3];
427: t[4] = b[4];
428: t[5] = b[5];
429: i2 = 6;
430: idiag += 36;
431: for (i = 1; i < m; i++) {
432: v = aa + 36 * ai[i];
433: vi = aj + ai[i];
434: nz = diag[i] - ai[i];
435: s[0] = b[i2];
436: s[1] = b[i2 + 1];
437: s[2] = b[i2 + 2];
438: s[3] = b[i2 + 3];
439: s[4] = b[i2 + 4];
440: s[5] = b[i2 + 5];
441: while (nz--) {
442: idx = 6 * (*vi++);
443: xw[0] = x[idx];
444: xw[1] = x[1 + idx];
445: xw[2] = x[2 + idx];
446: xw[3] = x[3 + idx];
447: xw[4] = x[4 + idx];
448: xw[5] = x[5 + idx];
449: PetscKernel_v_gets_v_minus_A_times_w_6(s, v, xw);
450: v += 36;
451: }
452: t[i2] = s[0];
453: t[i2 + 1] = s[1];
454: t[i2 + 2] = s[2];
455: t[i2 + 3] = s[3];
456: t[i2 + 4] = s[4];
457: t[i2 + 5] = s[5];
458: PetscKernel_v_gets_A_times_w_6(xw, idiag, s);
459: x[i2] = xw[0];
460: x[i2 + 1] = xw[1];
461: x[i2 + 2] = xw[2];
462: x[i2 + 3] = xw[3];
463: x[i2 + 4] = xw[4];
464: x[i2 + 5] = xw[5];
465: idiag += 36;
466: i2 += 6;
467: }
468: break;
469: case 7:
470: PetscKernel_v_gets_A_times_w_7(x, idiag, b);
471: t[0] = b[0];
472: t[1] = b[1];
473: t[2] = b[2];
474: t[3] = b[3];
475: t[4] = b[4];
476: t[5] = b[5];
477: t[6] = b[6];
478: i2 = 7;
479: idiag += 49;
480: for (i = 1; i < m; i++) {
481: v = aa + 49 * ai[i];
482: vi = aj + ai[i];
483: nz = diag[i] - ai[i];
484: s[0] = b[i2];
485: s[1] = b[i2 + 1];
486: s[2] = b[i2 + 2];
487: s[3] = b[i2 + 3];
488: s[4] = b[i2 + 4];
489: s[5] = b[i2 + 5];
490: s[6] = b[i2 + 6];
491: while (nz--) {
492: idx = 7 * (*vi++);
493: xw[0] = x[idx];
494: xw[1] = x[1 + idx];
495: xw[2] = x[2 + idx];
496: xw[3] = x[3 + idx];
497: xw[4] = x[4 + idx];
498: xw[5] = x[5 + idx];
499: xw[6] = x[6 + idx];
500: PetscKernel_v_gets_v_minus_A_times_w_7(s, v, xw);
501: v += 49;
502: }
503: t[i2] = s[0];
504: t[i2 + 1] = s[1];
505: t[i2 + 2] = s[2];
506: t[i2 + 3] = s[3];
507: t[i2 + 4] = s[4];
508: t[i2 + 5] = s[5];
509: t[i2 + 6] = s[6];
510: PetscKernel_v_gets_A_times_w_7(xw, idiag, s);
511: x[i2] = xw[0];
512: x[i2 + 1] = xw[1];
513: x[i2 + 2] = xw[2];
514: x[i2 + 3] = xw[3];
515: x[i2 + 4] = xw[4];
516: x[i2 + 5] = xw[5];
517: x[i2 + 6] = xw[6];
518: idiag += 49;
519: i2 += 7;
520: }
521: break;
522: default:
523: PetscKernel_w_gets_Ar_times_v(bs, bs, b, idiag, x);
524: PetscCall(PetscArraycpy(t, b, bs));
525: i2 = bs;
526: idiag += bs2;
527: for (i = 1; i < m; i++) {
528: v = aa + bs2 * ai[i];
529: vi = aj + ai[i];
530: nz = diag[i] - ai[i];
532: PetscCall(PetscArraycpy(w, b + i2, bs));
533: /* copy all rows of x that are needed into contiguous space */
534: workt = work;
535: for (j = 0; j < nz; j++) {
536: PetscCall(PetscArraycpy(workt, x + bs * (*vi++), bs));
537: workt += bs;
538: }
539: PetscKernel_w_gets_w_minus_Ar_times_v(bs, bs * nz, w, v, work);
540: PetscCall(PetscArraycpy(t + i2, w, bs));
541: PetscKernel_w_gets_Ar_times_v(bs, bs, w, idiag, x + i2);
543: idiag += bs2;
544: i2 += bs;
545: }
546: break;
547: }
548: /* for logging purposes assume number of nonzero in lower half is 1/2 of total */
549: PetscCall(PetscLogFlops(1.0 * bs2 * a->nz));
550: xb = t;
551: } else xb = b;
552: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
553: idiag = a->idiag + bs2 * (a->mbs - 1);
554: i2 = bs * (m - 1);
555: switch (bs) {
556: case 1:
557: s[0] = xb[i2];
558: PetscKernel_v_gets_A_times_w_1(xw, idiag, s);
559: x[i2] = xw[0];
560: i2 -= 1;
561: for (i = m - 2; i >= 0; i--) {
562: v = aa + (diag[i] + 1);
563: vi = aj + diag[i] + 1;
564: nz = ai[i + 1] - diag[i] - 1;
565: s[0] = xb[i2];
566: for (j = 0; j < nz; j++) {
567: xw[0] = x[vi[j]];
568: PetscKernel_v_gets_v_minus_A_times_w_1(s, (v + j), xw);
569: }
570: PetscKernel_v_gets_A_times_w_1(xw, idiag, s);
571: x[i2] = xw[0];
572: idiag -= 1;
573: i2 -= 1;
574: }
575: break;
576: case 2:
577: s[0] = xb[i2];
578: s[1] = xb[i2 + 1];
579: PetscKernel_v_gets_A_times_w_2(xw, idiag, s);
580: x[i2] = xw[0];
581: x[i2 + 1] = xw[1];
582: i2 -= 2;
583: idiag -= 4;
584: for (i = m - 2; i >= 0; i--) {
585: v = aa + 4 * (diag[i] + 1);
586: vi = aj + diag[i] + 1;
587: nz = ai[i + 1] - diag[i] - 1;
588: s[0] = xb[i2];
589: s[1] = xb[i2 + 1];
590: for (j = 0; j < nz; j++) {
591: idx = 2 * vi[j];
592: it = 4 * j;
593: xw[0] = x[idx];
594: xw[1] = x[1 + idx];
595: PetscKernel_v_gets_v_minus_A_times_w_2(s, (v + it), xw);
596: }
597: PetscKernel_v_gets_A_times_w_2(xw, idiag, s);
598: x[i2] = xw[0];
599: x[i2 + 1] = xw[1];
600: idiag -= 4;
601: i2 -= 2;
602: }
603: break;
604: case 3:
605: s[0] = xb[i2];
606: s[1] = xb[i2 + 1];
607: s[2] = xb[i2 + 2];
608: PetscKernel_v_gets_A_times_w_3(xw, idiag, s);
609: x[i2] = xw[0];
610: x[i2 + 1] = xw[1];
611: x[i2 + 2] = xw[2];
612: i2 -= 3;
613: idiag -= 9;
614: for (i = m - 2; i >= 0; i--) {
615: v = aa + 9 * (diag[i] + 1);
616: vi = aj + diag[i] + 1;
617: nz = ai[i + 1] - diag[i] - 1;
618: s[0] = xb[i2];
619: s[1] = xb[i2 + 1];
620: s[2] = xb[i2 + 2];
621: while (nz--) {
622: idx = 3 * (*vi++);
623: xw[0] = x[idx];
624: xw[1] = x[1 + idx];
625: xw[2] = x[2 + idx];
626: PetscKernel_v_gets_v_minus_A_times_w_3(s, v, xw);
627: v += 9;
628: }
629: PetscKernel_v_gets_A_times_w_3(xw, idiag, s);
630: x[i2] = xw[0];
631: x[i2 + 1] = xw[1];
632: x[i2 + 2] = xw[2];
633: idiag -= 9;
634: i2 -= 3;
635: }
636: break;
637: case 4:
638: s[0] = xb[i2];
639: s[1] = xb[i2 + 1];
640: s[2] = xb[i2 + 2];
641: s[3] = xb[i2 + 3];
642: PetscKernel_v_gets_A_times_w_4(xw, idiag, s);
643: x[i2] = xw[0];
644: x[i2 + 1] = xw[1];
645: x[i2 + 2] = xw[2];
646: x[i2 + 3] = xw[3];
647: i2 -= 4;
648: idiag -= 16;
649: for (i = m - 2; i >= 0; i--) {
650: v = aa + 16 * (diag[i] + 1);
651: vi = aj + diag[i] + 1;
652: nz = ai[i + 1] - diag[i] - 1;
653: s[0] = xb[i2];
654: s[1] = xb[i2 + 1];
655: s[2] = xb[i2 + 2];
656: s[3] = xb[i2 + 3];
657: while (nz--) {
658: idx = 4 * (*vi++);
659: xw[0] = x[idx];
660: xw[1] = x[1 + idx];
661: xw[2] = x[2 + idx];
662: xw[3] = x[3 + idx];
663: PetscKernel_v_gets_v_minus_A_times_w_4(s, v, xw);
664: v += 16;
665: }
666: PetscKernel_v_gets_A_times_w_4(xw, idiag, s);
667: x[i2] = xw[0];
668: x[i2 + 1] = xw[1];
669: x[i2 + 2] = xw[2];
670: x[i2 + 3] = xw[3];
671: idiag -= 16;
672: i2 -= 4;
673: }
674: break;
675: case 5:
676: s[0] = xb[i2];
677: s[1] = xb[i2 + 1];
678: s[2] = xb[i2 + 2];
679: s[3] = xb[i2 + 3];
680: s[4] = xb[i2 + 4];
681: PetscKernel_v_gets_A_times_w_5(xw, idiag, s);
682: x[i2] = xw[0];
683: x[i2 + 1] = xw[1];
684: x[i2 + 2] = xw[2];
685: x[i2 + 3] = xw[3];
686: x[i2 + 4] = xw[4];
687: i2 -= 5;
688: idiag -= 25;
689: for (i = m - 2; i >= 0; i--) {
690: v = aa + 25 * (diag[i] + 1);
691: vi = aj + diag[i] + 1;
692: nz = ai[i + 1] - diag[i] - 1;
693: s[0] = xb[i2];
694: s[1] = xb[i2 + 1];
695: s[2] = xb[i2 + 2];
696: s[3] = xb[i2 + 3];
697: s[4] = xb[i2 + 4];
698: while (nz--) {
699: idx = 5 * (*vi++);
700: xw[0] = x[idx];
701: xw[1] = x[1 + idx];
702: xw[2] = x[2 + idx];
703: xw[3] = x[3 + idx];
704: xw[4] = x[4 + idx];
705: PetscKernel_v_gets_v_minus_A_times_w_5(s, v, xw);
706: v += 25;
707: }
708: PetscKernel_v_gets_A_times_w_5(xw, idiag, s);
709: x[i2] = xw[0];
710: x[i2 + 1] = xw[1];
711: x[i2 + 2] = xw[2];
712: x[i2 + 3] = xw[3];
713: x[i2 + 4] = xw[4];
714: idiag -= 25;
715: i2 -= 5;
716: }
717: break;
718: case 6:
719: s[0] = xb[i2];
720: s[1] = xb[i2 + 1];
721: s[2] = xb[i2 + 2];
722: s[3] = xb[i2 + 3];
723: s[4] = xb[i2 + 4];
724: s[5] = xb[i2 + 5];
725: PetscKernel_v_gets_A_times_w_6(xw, idiag, s);
726: x[i2] = xw[0];
727: x[i2 + 1] = xw[1];
728: x[i2 + 2] = xw[2];
729: x[i2 + 3] = xw[3];
730: x[i2 + 4] = xw[4];
731: x[i2 + 5] = xw[5];
732: i2 -= 6;
733: idiag -= 36;
734: for (i = m - 2; i >= 0; i--) {
735: v = aa + 36 * (diag[i] + 1);
736: vi = aj + diag[i] + 1;
737: nz = ai[i + 1] - diag[i] - 1;
738: s[0] = xb[i2];
739: s[1] = xb[i2 + 1];
740: s[2] = xb[i2 + 2];
741: s[3] = xb[i2 + 3];
742: s[4] = xb[i2 + 4];
743: s[5] = xb[i2 + 5];
744: while (nz--) {
745: idx = 6 * (*vi++);
746: xw[0] = x[idx];
747: xw[1] = x[1 + idx];
748: xw[2] = x[2 + idx];
749: xw[3] = x[3 + idx];
750: xw[4] = x[4 + idx];
751: xw[5] = x[5 + idx];
752: PetscKernel_v_gets_v_minus_A_times_w_6(s, v, xw);
753: v += 36;
754: }
755: PetscKernel_v_gets_A_times_w_6(xw, idiag, s);
756: x[i2] = xw[0];
757: x[i2 + 1] = xw[1];
758: x[i2 + 2] = xw[2];
759: x[i2 + 3] = xw[3];
760: x[i2 + 4] = xw[4];
761: x[i2 + 5] = xw[5];
762: idiag -= 36;
763: i2 -= 6;
764: }
765: break;
766: case 7:
767: s[0] = xb[i2];
768: s[1] = xb[i2 + 1];
769: s[2] = xb[i2 + 2];
770: s[3] = xb[i2 + 3];
771: s[4] = xb[i2 + 4];
772: s[5] = xb[i2 + 5];
773: s[6] = xb[i2 + 6];
774: PetscKernel_v_gets_A_times_w_7(x, idiag, b);
775: x[i2] = xw[0];
776: x[i2 + 1] = xw[1];
777: x[i2 + 2] = xw[2];
778: x[i2 + 3] = xw[3];
779: x[i2 + 4] = xw[4];
780: x[i2 + 5] = xw[5];
781: x[i2 + 6] = xw[6];
782: i2 -= 7;
783: idiag -= 49;
784: for (i = m - 2; i >= 0; i--) {
785: v = aa + 49 * (diag[i] + 1);
786: vi = aj + diag[i] + 1;
787: nz = ai[i + 1] - diag[i] - 1;
788: s[0] = xb[i2];
789: s[1] = xb[i2 + 1];
790: s[2] = xb[i2 + 2];
791: s[3] = xb[i2 + 3];
792: s[4] = xb[i2 + 4];
793: s[5] = xb[i2 + 5];
794: s[6] = xb[i2 + 6];
795: while (nz--) {
796: idx = 7 * (*vi++);
797: xw[0] = x[idx];
798: xw[1] = x[1 + idx];
799: xw[2] = x[2 + idx];
800: xw[3] = x[3 + idx];
801: xw[4] = x[4 + idx];
802: xw[5] = x[5 + idx];
803: xw[6] = x[6 + idx];
804: PetscKernel_v_gets_v_minus_A_times_w_7(s, v, xw);
805: v += 49;
806: }
807: PetscKernel_v_gets_A_times_w_7(xw, idiag, s);
808: x[i2] = xw[0];
809: x[i2 + 1] = xw[1];
810: x[i2 + 2] = xw[2];
811: x[i2 + 3] = xw[3];
812: x[i2 + 4] = xw[4];
813: x[i2 + 5] = xw[5];
814: x[i2 + 6] = xw[6];
815: idiag -= 49;
816: i2 -= 7;
817: }
818: break;
819: default:
820: PetscCall(PetscArraycpy(w, xb + i2, bs));
821: PetscKernel_w_gets_Ar_times_v(bs, bs, w, idiag, x + i2);
822: i2 -= bs;
823: idiag -= bs2;
824: for (i = m - 2; i >= 0; i--) {
825: v = aa + bs2 * (diag[i] + 1);
826: vi = aj + diag[i] + 1;
827: nz = ai[i + 1] - diag[i] - 1;
829: PetscCall(PetscArraycpy(w, xb + i2, bs));
830: /* copy all rows of x that are needed into contiguous space */
831: workt = work;
832: for (j = 0; j < nz; j++) {
833: PetscCall(PetscArraycpy(workt, x + bs * (*vi++), bs));
834: workt += bs;
835: }
836: PetscKernel_w_gets_w_minus_Ar_times_v(bs, bs * nz, w, v, work);
837: PetscKernel_w_gets_Ar_times_v(bs, bs, w, idiag, x + i2);
839: idiag -= bs2;
840: i2 -= bs;
841: }
842: break;
843: }
844: PetscCall(PetscLogFlops(1.0 * bs2 * (a->nz)));
845: }
846: its--;
847: }
848: while (its--) {
849: if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {
850: idiag = a->idiag;
851: i2 = 0;
852: switch (bs) {
853: case 1:
854: for (i = 0; i < m; i++) {
855: v = aa + ai[i];
856: vi = aj + ai[i];
857: nz = ai[i + 1] - ai[i];
858: s[0] = b[i2];
859: for (j = 0; j < nz; j++) {
860: xw[0] = x[vi[j]];
861: PetscKernel_v_gets_v_minus_A_times_w_1(s, (v + j), xw);
862: }
863: PetscKernel_v_gets_A_times_w_1(xw, idiag, s);
864: x[i2] += xw[0];
865: idiag += 1;
866: i2 += 1;
867: }
868: break;
869: case 2:
870: for (i = 0; i < m; i++) {
871: v = aa + 4 * ai[i];
872: vi = aj + ai[i];
873: nz = ai[i + 1] - ai[i];
874: s[0] = b[i2];
875: s[1] = b[i2 + 1];
876: for (j = 0; j < nz; j++) {
877: idx = 2 * vi[j];
878: it = 4 * j;
879: xw[0] = x[idx];
880: xw[1] = x[1 + idx];
881: PetscKernel_v_gets_v_minus_A_times_w_2(s, (v + it), xw);
882: }
883: PetscKernel_v_gets_A_times_w_2(xw, idiag, s);
884: x[i2] += xw[0];
885: x[i2 + 1] += xw[1];
886: idiag += 4;
887: i2 += 2;
888: }
889: break;
890: case 3:
891: for (i = 0; i < m; i++) {
892: v = aa + 9 * ai[i];
893: vi = aj + ai[i];
894: nz = ai[i + 1] - ai[i];
895: s[0] = b[i2];
896: s[1] = b[i2 + 1];
897: s[2] = b[i2 + 2];
898: while (nz--) {
899: idx = 3 * (*vi++);
900: xw[0] = x[idx];
901: xw[1] = x[1 + idx];
902: xw[2] = x[2 + idx];
903: PetscKernel_v_gets_v_minus_A_times_w_3(s, v, xw);
904: v += 9;
905: }
906: PetscKernel_v_gets_A_times_w_3(xw, idiag, s);
907: x[i2] += xw[0];
908: x[i2 + 1] += xw[1];
909: x[i2 + 2] += xw[2];
910: idiag += 9;
911: i2 += 3;
912: }
913: break;
914: case 4:
915: for (i = 0; i < m; i++) {
916: v = aa + 16 * ai[i];
917: vi = aj + ai[i];
918: nz = ai[i + 1] - ai[i];
919: s[0] = b[i2];
920: s[1] = b[i2 + 1];
921: s[2] = b[i2 + 2];
922: s[3] = b[i2 + 3];
923: while (nz--) {
924: idx = 4 * (*vi++);
925: xw[0] = x[idx];
926: xw[1] = x[1 + idx];
927: xw[2] = x[2 + idx];
928: xw[3] = x[3 + idx];
929: PetscKernel_v_gets_v_minus_A_times_w_4(s, v, xw);
930: v += 16;
931: }
932: PetscKernel_v_gets_A_times_w_4(xw, idiag, s);
933: x[i2] += xw[0];
934: x[i2 + 1] += xw[1];
935: x[i2 + 2] += xw[2];
936: x[i2 + 3] += xw[3];
937: idiag += 16;
938: i2 += 4;
939: }
940: break;
941: case 5:
942: for (i = 0; i < m; i++) {
943: v = aa + 25 * ai[i];
944: vi = aj + ai[i];
945: nz = ai[i + 1] - ai[i];
946: s[0] = b[i2];
947: s[1] = b[i2 + 1];
948: s[2] = b[i2 + 2];
949: s[3] = b[i2 + 3];
950: s[4] = b[i2 + 4];
951: while (nz--) {
952: idx = 5 * (*vi++);
953: xw[0] = x[idx];
954: xw[1] = x[1 + idx];
955: xw[2] = x[2 + idx];
956: xw[3] = x[3 + idx];
957: xw[4] = x[4 + idx];
958: PetscKernel_v_gets_v_minus_A_times_w_5(s, v, xw);
959: v += 25;
960: }
961: PetscKernel_v_gets_A_times_w_5(xw, idiag, s);
962: x[i2] += xw[0];
963: x[i2 + 1] += xw[1];
964: x[i2 + 2] += xw[2];
965: x[i2 + 3] += xw[3];
966: x[i2 + 4] += xw[4];
967: idiag += 25;
968: i2 += 5;
969: }
970: break;
971: case 6:
972: for (i = 0; i < m; i++) {
973: v = aa + 36 * ai[i];
974: vi = aj + ai[i];
975: nz = ai[i + 1] - ai[i];
976: s[0] = b[i2];
977: s[1] = b[i2 + 1];
978: s[2] = b[i2 + 2];
979: s[3] = b[i2 + 3];
980: s[4] = b[i2 + 4];
981: s[5] = b[i2 + 5];
982: while (nz--) {
983: idx = 6 * (*vi++);
984: xw[0] = x[idx];
985: xw[1] = x[1 + idx];
986: xw[2] = x[2 + idx];
987: xw[3] = x[3 + idx];
988: xw[4] = x[4 + idx];
989: xw[5] = x[5 + idx];
990: PetscKernel_v_gets_v_minus_A_times_w_6(s, v, xw);
991: v += 36;
992: }
993: PetscKernel_v_gets_A_times_w_6(xw, idiag, s);
994: x[i2] += xw[0];
995: x[i2 + 1] += xw[1];
996: x[i2 + 2] += xw[2];
997: x[i2 + 3] += xw[3];
998: x[i2 + 4] += xw[4];
999: x[i2 + 5] += xw[5];
1000: idiag += 36;
1001: i2 += 6;
1002: }
1003: break;
1004: case 7:
1005: for (i = 0; i < m; i++) {
1006: v = aa + 49 * ai[i];
1007: vi = aj + ai[i];
1008: nz = ai[i + 1] - ai[i];
1009: s[0] = b[i2];
1010: s[1] = b[i2 + 1];
1011: s[2] = b[i2 + 2];
1012: s[3] = b[i2 + 3];
1013: s[4] = b[i2 + 4];
1014: s[5] = b[i2 + 5];
1015: s[6] = b[i2 + 6];
1016: while (nz--) {
1017: idx = 7 * (*vi++);
1018: xw[0] = x[idx];
1019: xw[1] = x[1 + idx];
1020: xw[2] = x[2 + idx];
1021: xw[3] = x[3 + idx];
1022: xw[4] = x[4 + idx];
1023: xw[5] = x[5 + idx];
1024: xw[6] = x[6 + idx];
1025: PetscKernel_v_gets_v_minus_A_times_w_7(s, v, xw);
1026: v += 49;
1027: }
1028: PetscKernel_v_gets_A_times_w_7(xw, idiag, s);
1029: x[i2] += xw[0];
1030: x[i2 + 1] += xw[1];
1031: x[i2 + 2] += xw[2];
1032: x[i2 + 3] += xw[3];
1033: x[i2 + 4] += xw[4];
1034: x[i2 + 5] += xw[5];
1035: x[i2 + 6] += xw[6];
1036: idiag += 49;
1037: i2 += 7;
1038: }
1039: break;
1040: default:
1041: for (i = 0; i < m; i++) {
1042: v = aa + bs2 * ai[i];
1043: vi = aj + ai[i];
1044: nz = ai[i + 1] - ai[i];
1046: PetscCall(PetscArraycpy(w, b + i2, bs));
1047: /* copy all rows of x that are needed into contiguous space */
1048: workt = work;
1049: for (j = 0; j < nz; j++) {
1050: PetscCall(PetscArraycpy(workt, x + bs * (*vi++), bs));
1051: workt += bs;
1052: }
1053: PetscKernel_w_gets_w_minus_Ar_times_v(bs, bs * nz, w, v, work);
1054: PetscKernel_w_gets_w_plus_Ar_times_v(bs, bs, w, idiag, x + i2);
1056: idiag += bs2;
1057: i2 += bs;
1058: }
1059: break;
1060: }
1061: PetscCall(PetscLogFlops(2.0 * bs2 * a->nz));
1062: }
1063: if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {
1064: idiag = a->idiag + bs2 * (a->mbs - 1);
1065: i2 = bs * (m - 1);
1066: switch (bs) {
1067: case 1:
1068: for (i = m - 1; i >= 0; i--) {
1069: v = aa + ai[i];
1070: vi = aj + ai[i];
1071: nz = ai[i + 1] - ai[i];
1072: s[0] = b[i2];
1073: for (j = 0; j < nz; j++) {
1074: xw[0] = x[vi[j]];
1075: PetscKernel_v_gets_v_minus_A_times_w_1(s, (v + j), xw);
1076: }
1077: PetscKernel_v_gets_A_times_w_1(xw, idiag, s);
1078: x[i2] += xw[0];
1079: idiag -= 1;
1080: i2 -= 1;
1081: }
1082: break;
1083: case 2:
1084: for (i = m - 1; i >= 0; i--) {
1085: v = aa + 4 * ai[i];
1086: vi = aj + ai[i];
1087: nz = ai[i + 1] - ai[i];
1088: s[0] = b[i2];
1089: s[1] = b[i2 + 1];
1090: for (j = 0; j < nz; j++) {
1091: idx = 2 * vi[j];
1092: it = 4 * j;
1093: xw[0] = x[idx];
1094: xw[1] = x[1 + idx];
1095: PetscKernel_v_gets_v_minus_A_times_w_2(s, (v + it), xw);
1096: }
1097: PetscKernel_v_gets_A_times_w_2(xw, idiag, s);
1098: x[i2] += xw[0];
1099: x[i2 + 1] += xw[1];
1100: idiag -= 4;
1101: i2 -= 2;
1102: }
1103: break;
1104: case 3:
1105: for (i = m - 1; i >= 0; i--) {
1106: v = aa + 9 * ai[i];
1107: vi = aj + ai[i];
1108: nz = ai[i + 1] - ai[i];
1109: s[0] = b[i2];
1110: s[1] = b[i2 + 1];
1111: s[2] = b[i2 + 2];
1112: while (nz--) {
1113: idx = 3 * (*vi++);
1114: xw[0] = x[idx];
1115: xw[1] = x[1 + idx];
1116: xw[2] = x[2 + idx];
1117: PetscKernel_v_gets_v_minus_A_times_w_3(s, v, xw);
1118: v += 9;
1119: }
1120: PetscKernel_v_gets_A_times_w_3(xw, idiag, s);
1121: x[i2] += xw[0];
1122: x[i2 + 1] += xw[1];
1123: x[i2 + 2] += xw[2];
1124: idiag -= 9;
1125: i2 -= 3;
1126: }
1127: break;
1128: case 4:
1129: for (i = m - 1; i >= 0; i--) {
1130: v = aa + 16 * ai[i];
1131: vi = aj + ai[i];
1132: nz = ai[i + 1] - ai[i];
1133: s[0] = b[i2];
1134: s[1] = b[i2 + 1];
1135: s[2] = b[i2 + 2];
1136: s[3] = b[i2 + 3];
1137: while (nz--) {
1138: idx = 4 * (*vi++);
1139: xw[0] = x[idx];
1140: xw[1] = x[1 + idx];
1141: xw[2] = x[2 + idx];
1142: xw[3] = x[3 + idx];
1143: PetscKernel_v_gets_v_minus_A_times_w_4(s, v, xw);
1144: v += 16;
1145: }
1146: PetscKernel_v_gets_A_times_w_4(xw, idiag, s);
1147: x[i2] += xw[0];
1148: x[i2 + 1] += xw[1];
1149: x[i2 + 2] += xw[2];
1150: x[i2 + 3] += xw[3];
1151: idiag -= 16;
1152: i2 -= 4;
1153: }
1154: break;
1155: case 5:
1156: for (i = m - 1; i >= 0; i--) {
1157: v = aa + 25 * ai[i];
1158: vi = aj + ai[i];
1159: nz = ai[i + 1] - ai[i];
1160: s[0] = b[i2];
1161: s[1] = b[i2 + 1];
1162: s[2] = b[i2 + 2];
1163: s[3] = b[i2 + 3];
1164: s[4] = b[i2 + 4];
1165: while (nz--) {
1166: idx = 5 * (*vi++);
1167: xw[0] = x[idx];
1168: xw[1] = x[1 + idx];
1169: xw[2] = x[2 + idx];
1170: xw[3] = x[3 + idx];
1171: xw[4] = x[4 + idx];
1172: PetscKernel_v_gets_v_minus_A_times_w_5(s, v, xw);
1173: v += 25;
1174: }
1175: PetscKernel_v_gets_A_times_w_5(xw, idiag, s);
1176: x[i2] += xw[0];
1177: x[i2 + 1] += xw[1];
1178: x[i2 + 2] += xw[2];
1179: x[i2 + 3] += xw[3];
1180: x[i2 + 4] += xw[4];
1181: idiag -= 25;
1182: i2 -= 5;
1183: }
1184: break;
1185: case 6:
1186: for (i = m - 1; i >= 0; i--) {
1187: v = aa + 36 * ai[i];
1188: vi = aj + ai[i];
1189: nz = ai[i + 1] - ai[i];
1190: s[0] = b[i2];
1191: s[1] = b[i2 + 1];
1192: s[2] = b[i2 + 2];
1193: s[3] = b[i2 + 3];
1194: s[4] = b[i2 + 4];
1195: s[5] = b[i2 + 5];
1196: while (nz--) {
1197: idx = 6 * (*vi++);
1198: xw[0] = x[idx];
1199: xw[1] = x[1 + idx];
1200: xw[2] = x[2 + idx];
1201: xw[3] = x[3 + idx];
1202: xw[4] = x[4 + idx];
1203: xw[5] = x[5 + idx];
1204: PetscKernel_v_gets_v_minus_A_times_w_6(s, v, xw);
1205: v += 36;
1206: }
1207: PetscKernel_v_gets_A_times_w_6(xw, idiag, s);
1208: x[i2] += xw[0];
1209: x[i2 + 1] += xw[1];
1210: x[i2 + 2] += xw[2];
1211: x[i2 + 3] += xw[3];
1212: x[i2 + 4] += xw[4];
1213: x[i2 + 5] += xw[5];
1214: idiag -= 36;
1215: i2 -= 6;
1216: }
1217: break;
1218: case 7:
1219: for (i = m - 1; i >= 0; i--) {
1220: v = aa + 49 * ai[i];
1221: vi = aj + ai[i];
1222: nz = ai[i + 1] - ai[i];
1223: s[0] = b[i2];
1224: s[1] = b[i2 + 1];
1225: s[2] = b[i2 + 2];
1226: s[3] = b[i2 + 3];
1227: s[4] = b[i2 + 4];
1228: s[5] = b[i2 + 5];
1229: s[6] = b[i2 + 6];
1230: while (nz--) {
1231: idx = 7 * (*vi++);
1232: xw[0] = x[idx];
1233: xw[1] = x[1 + idx];
1234: xw[2] = x[2 + idx];
1235: xw[3] = x[3 + idx];
1236: xw[4] = x[4 + idx];
1237: xw[5] = x[5 + idx];
1238: xw[6] = x[6 + idx];
1239: PetscKernel_v_gets_v_minus_A_times_w_7(s, v, xw);
1240: v += 49;
1241: }
1242: PetscKernel_v_gets_A_times_w_7(xw, idiag, s);
1243: x[i2] += xw[0];
1244: x[i2 + 1] += xw[1];
1245: x[i2 + 2] += xw[2];
1246: x[i2 + 3] += xw[3];
1247: x[i2 + 4] += xw[4];
1248: x[i2 + 5] += xw[5];
1249: x[i2 + 6] += xw[6];
1250: idiag -= 49;
1251: i2 -= 7;
1252: }
1253: break;
1254: default:
1255: for (i = m - 1; i >= 0; i--) {
1256: v = aa + bs2 * ai[i];
1257: vi = aj + ai[i];
1258: nz = ai[i + 1] - ai[i];
1260: PetscCall(PetscArraycpy(w, b + i2, bs));
1261: /* copy all rows of x that are needed into contiguous space */
1262: workt = work;
1263: for (j = 0; j < nz; j++) {
1264: PetscCall(PetscArraycpy(workt, x + bs * (*vi++), bs));
1265: workt += bs;
1266: }
1267: PetscKernel_w_gets_w_minus_Ar_times_v(bs, bs * nz, w, v, work);
1268: PetscKernel_w_gets_w_plus_Ar_times_v(bs, bs, w, idiag, x + i2);
1270: idiag -= bs2;
1271: i2 -= bs;
1272: }
1273: break;
1274: }
1275: PetscCall(PetscLogFlops(2.0 * bs2 * (a->nz)));
1276: }
1277: }
1278: PetscCall(VecRestoreArray(xx, &x));
1279: PetscCall(VecRestoreArrayRead(bb, &b));
1280: PetscFunctionReturn(PETSC_SUCCESS);
1281: }
1283: /*
1284: Special version for direct calls from Fortran (Used in PETSc-fun3d)
1285: */
1286: #if PetscDefined(HAVE_FORTRAN_CAPS)
1287: #define matsetvaluesblocked4_ MATSETVALUESBLOCKED4
1288: #elif !PetscDefined(HAVE_FORTRAN_UNDERSCORE)
1289: #define matsetvaluesblocked4_ matsetvaluesblocked4
1290: #endif
1292: PETSC_EXTERN void matsetvaluesblocked4_(Mat *AA, PetscInt *mm, const PetscInt im[], PetscInt *nn, const PetscInt in[], const PetscScalar v[])
1293: {
1294: Mat A = *AA;
1295: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
1296: PetscInt *rp, k, low, high, t, ii, jj, row, nrow, i, col, l, N, m = *mm, n = *nn;
1297: PetscInt *ai = a->i, *ailen = a->ilen;
1298: PetscInt *aj = a->j, stepval, lastcol = -1;
1299: const PetscScalar *value = v;
1300: MatScalar *ap, *aa = a->a, *bap;
1302: PetscFunctionBegin;
1303: if (A->rmap->bs != 4) SETERRABORT(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Can only be called with a block size of 4");
1304: stepval = (n - 1) * 4;
1305: for (k = 0; k < m; k++) { /* loop over added rows */
1306: row = im[k];
1307: rp = aj + ai[row];
1308: ap = aa + 16 * ai[row];
1309: nrow = ailen[row];
1310: low = 0;
1311: high = nrow;
1312: for (l = 0; l < n; l++) { /* loop over added columns */
1313: col = in[l];
1314: if (col <= lastcol) low = 0;
1315: else high = nrow;
1316: lastcol = col;
1317: value = v + k * (stepval + 4 + l) * 4;
1318: while (high - low > 7) {
1319: t = (low + high) / 2;
1320: if (rp[t] > col) high = t;
1321: else low = t;
1322: }
1323: for (i = low; i < high; i++) {
1324: if (rp[i] > col) break;
1325: if (rp[i] == col) {
1326: bap = ap + 16 * i;
1327: for (ii = 0; ii < 4; ii++, value += stepval) {
1328: for (jj = ii; jj < 16; jj += 4) bap[jj] += *value++;
1329: }
1330: goto noinsert2;
1331: }
1332: }
1333: N = nrow++ - 1;
1334: high++; /* added new column index thus must search to one higher than before */
1335: /* shift up all the later entries in this row */
1336: for (ii = N; ii >= i; ii--) {
1337: rp[ii + 1] = rp[ii];
1338: PetscCallVoid(PetscArraycpy(ap + 16 * (ii + 1), ap + 16 * (ii), 16));
1339: }
1340: if (N >= i) PetscCallVoid(PetscArrayzero(ap + 16 * i, 16));
1341: rp[i] = col;
1342: bap = ap + 16 * i;
1343: for (ii = 0; ii < 4; ii++, value += stepval) {
1344: for (jj = ii; jj < 16; jj += 4) bap[jj] = *value++;
1345: }
1346: noinsert2:;
1347: low = i;
1348: }
1349: ailen[row] = nrow;
1350: }
1351: PetscFunctionReturnVoid();
1352: }
1354: #if PetscDefined(HAVE_FORTRAN_CAPS)
1355: #define matsetvalues4_ MATSETVALUES4
1356: #elif !PetscDefined(HAVE_FORTRAN_UNDERSCORE)
1357: #define matsetvalues4_ matsetvalues4
1358: #endif
1360: PETSC_EXTERN void matsetvalues4_(Mat *AA, PetscInt *mm, PetscInt *im, PetscInt *nn, PetscInt *in, PetscScalar *v)
1361: {
1362: Mat A = *AA;
1363: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
1364: PetscInt *rp, k, low, high, t, row, nrow, i, col, l, N, n = *nn, m = *mm;
1365: PetscInt *ai = a->i, *ailen = a->ilen;
1366: PetscInt *aj = a->j, brow, bcol;
1367: PetscInt ridx, cidx, lastcol = -1;
1368: MatScalar *ap, value, *aa = a->a, *bap;
1370: PetscFunctionBegin;
1371: for (k = 0; k < m; k++) { /* loop over added rows */
1372: row = im[k];
1373: brow = row / 4;
1374: rp = aj + ai[brow];
1375: ap = aa + 16 * ai[brow];
1376: nrow = ailen[brow];
1377: low = 0;
1378: high = nrow;
1379: for (l = 0; l < n; l++) { /* loop over added columns */
1380: col = in[l];
1381: bcol = col / 4;
1382: ridx = row % 4;
1383: cidx = col % 4;
1384: value = v[l + k * n];
1385: if (col <= lastcol) low = 0;
1386: else high = nrow;
1387: lastcol = col;
1388: while (high - low > 7) {
1389: t = (low + high) / 2;
1390: if (rp[t] > bcol) high = t;
1391: else low = t;
1392: }
1393: for (i = low; i < high; i++) {
1394: if (rp[i] > bcol) break;
1395: if (rp[i] == bcol) {
1396: bap = ap + 16 * i + 4 * cidx + ridx;
1397: *bap += value;
1398: goto noinsert1;
1399: }
1400: }
1401: N = nrow++ - 1;
1402: high++; /* added new column thus must search to one higher than before */
1403: /* shift up all the later entries in this row */
1404: PetscCallVoid(PetscArraymove(rp + i + 1, rp + i, N - i + 1));
1405: PetscCallVoid(PetscArraymove(ap + 16 * i + 16, ap + 16 * i, 16 * (N - i + 1)));
1406: PetscCallVoid(PetscArrayzero(ap + 16 * i, 16));
1407: rp[i] = bcol;
1408: ap[16 * i + 4 * cidx + ridx] = value;
1409: noinsert1:;
1410: low = i;
1411: }
1412: ailen[brow] = nrow;
1413: }
1414: PetscFunctionReturnVoid();
1415: }
1417: static PetscErrorCode MatGetRowIJ_SeqBAIJ(Mat A, PetscInt oshift, PetscBool symmetric, PetscBool blockcompressed, PetscInt *nn, const PetscInt *inia[], const PetscInt *inja[], PetscBool *done)
1418: {
1419: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
1420: PetscInt i, j, n = a->mbs, nz = a->i[n], *tia, *tja, bs = A->rmap->bs, k, l, cnt;
1421: PetscInt **ia = (PetscInt **)inia, **ja = (PetscInt **)inja;
1423: PetscFunctionBegin;
1424: *nn = n;
1425: if (!ia) PetscFunctionReturn(PETSC_SUCCESS);
1426: if (symmetric) {
1427: PetscCall(MatToSymmetricIJ_SeqAIJ(n, a->i, a->j, PETSC_TRUE, 0, 0, &tia, &tja));
1428: nz = tia[n];
1429: } else {
1430: tia = a->i;
1431: tja = a->j;
1432: }
1434: if (!blockcompressed && bs > 1) {
1435: (*nn) *= bs;
1436: /* malloc & create the natural set of indices */
1437: PetscCall(PetscMalloc1((n + 1) * bs, ia));
1438: if (n) {
1439: (*ia)[0] = oshift;
1440: for (j = 1; j < bs; j++) (*ia)[j] = (tia[1] - tia[0]) * bs + (*ia)[j - 1];
1441: }
1443: for (i = 1; i < n; i++) {
1444: (*ia)[i * bs] = (tia[i] - tia[i - 1]) * bs + (*ia)[i * bs - 1];
1445: for (j = 1; j < bs; j++) (*ia)[i * bs + j] = (tia[i + 1] - tia[i]) * bs + (*ia)[i * bs + j - 1];
1446: }
1447: if (n) (*ia)[n * bs] = (tia[n] - tia[n - 1]) * bs + (*ia)[n * bs - 1];
1449: if (inja) {
1450: PetscCall(PetscMalloc1(nz * bs * bs, ja));
1451: cnt = 0;
1452: for (i = 0; i < n; i++) {
1453: for (j = 0; j < bs; j++) {
1454: for (k = tia[i]; k < tia[i + 1]; k++) {
1455: for (l = 0; l < bs; l++) (*ja)[cnt++] = bs * tja[k] + l;
1456: }
1457: }
1458: }
1459: }
1461: if (symmetric) { /* deallocate memory allocated in MatToSymmetricIJ_SeqAIJ() */
1462: PetscCall(PetscFree(tia));
1463: PetscCall(PetscFree(tja));
1464: }
1465: } else if (oshift == 1) {
1466: if (symmetric) {
1467: nz = tia[A->rmap->n / bs];
1468: /* add 1 to i and j indices */
1469: for (i = 0; i < A->rmap->n / bs + 1; i++) tia[i] = tia[i] + 1;
1470: *ia = tia;
1471: if (ja) {
1472: for (i = 0; i < nz; i++) tja[i] = tja[i] + 1;
1473: *ja = tja;
1474: }
1475: } else {
1476: nz = a->i[A->rmap->n / bs];
1477: /* malloc space and add 1 to i and j indices */
1478: PetscCall(PetscMalloc1(A->rmap->n / bs + 1, ia));
1479: for (i = 0; i < A->rmap->n / bs + 1; i++) (*ia)[i] = a->i[i] + 1;
1480: if (ja) {
1481: PetscCall(PetscMalloc1(nz, ja));
1482: for (i = 0; i < nz; i++) (*ja)[i] = a->j[i] + 1;
1483: }
1484: }
1485: } else {
1486: *ia = tia;
1487: if (ja) *ja = tja;
1488: }
1489: PetscFunctionReturn(PETSC_SUCCESS);
1490: }
1492: static PetscErrorCode MatRestoreRowIJ_SeqBAIJ(Mat A, PetscInt oshift, PetscBool symmetric, PetscBool blockcompressed, PetscInt *nn, const PetscInt *ia[], const PetscInt *ja[], PetscBool *done)
1493: {
1494: PetscFunctionBegin;
1495: if (!ia) PetscFunctionReturn(PETSC_SUCCESS);
1496: if ((!blockcompressed && A->rmap->bs > 1) || (symmetric || oshift == 1)) {
1497: PetscCall(PetscFree(*ia));
1498: if (ja) PetscCall(PetscFree(*ja));
1499: }
1500: PetscFunctionReturn(PETSC_SUCCESS);
1501: }
1503: PetscErrorCode MatDestroy_SeqBAIJ(Mat A)
1504: {
1505: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
1507: PetscFunctionBegin;
1508: if (A->hash_active) {
1509: PetscInt bs;
1510: A->ops[0] = a->cops;
1511: PetscCall(PetscHMapIJVDestroy(&a->ht));
1512: PetscCall(MatGetBlockSize(A, &bs));
1513: if (bs > 1) PetscCall(PetscHSetIJDestroy(&a->bht));
1514: PetscCall(PetscFree(a->dnz));
1515: PetscCall(PetscFree(a->bdnz));
1516: A->hash_active = PETSC_FALSE;
1517: }
1518: PetscCall(PetscLogObjectState((PetscObject)A, "Rows=%" PetscInt_FMT ", Cols=%" PetscInt_FMT ", NZ=%" PetscInt_FMT, A->rmap->N, A->cmap->n, a->nz));
1519: PetscCall(MatSeqXAIJFreeAIJ(A, &a->a, &a->j, &a->i));
1520: PetscCall(ISDestroy(&a->row));
1521: PetscCall(ISDestroy(&a->col));
1522: PetscCall(PetscFree(a->diag));
1523: PetscCall(PetscFree(a->idiag));
1524: if (a->free_imax_ilen) PetscCall(PetscFree2(a->imax, a->ilen));
1525: PetscCall(PetscFree(a->solve_work));
1526: PetscCall(PetscFree(a->mult_work));
1527: PetscCall(PetscFree(a->sor_workt));
1528: PetscCall(PetscFree(a->sor_work));
1529: PetscCall(ISDestroy(&a->icol));
1530: PetscCall(PetscFree(a->saved_values));
1531: PetscCall(PetscFree2(a->compressedrow.i, a->compressedrow.rindex));
1533: PetscCall(MatDestroy(&a->sbaijMat));
1534: PetscCall(MatDestroy(&a->parent));
1535: PetscCall(PetscFree(A->data));
1537: PetscCall(PetscObjectChangeTypeName((PetscObject)A, NULL));
1538: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSeqBAIJGetArray_C", NULL));
1539: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSeqBAIJRestoreArray_C", NULL));
1540: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatStoreValues_C", NULL));
1541: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatRetrieveValues_C", NULL));
1542: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSeqBAIJSetColumnIndices_C", NULL));
1543: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqbaij_seqaij_C", NULL));
1544: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqbaij_seqsbaij_C", NULL));
1545: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSeqBAIJSetPreallocation_C", NULL));
1546: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSeqBAIJSetPreallocationCSR_C", NULL));
1547: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqbaij_seqbstrm_C", NULL));
1548: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatIsTranspose_C", NULL));
1549: #if PetscDefined(HAVE_HYPRE)
1550: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqbaij_hypre_C", NULL));
1551: #endif
1552: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqbaij_is_C", NULL));
1553: #if PetscDefined(HAVE_LIBXSMM)
1554: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqbaij_seqbaijlibxsmm_C", NULL));
1555: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_seqbaijlibxsmm_seqdense_C", NULL));
1556: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqbaijlibxsmm_seqbaij_C", NULL));
1557: #endif
1558: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatFactorGetSolverType_C", NULL));
1559: PetscFunctionReturn(PETSC_SUCCESS);
1560: }
1562: static PetscErrorCode MatSetOption_SeqBAIJ(Mat A, MatOption op, PetscBool flg)
1563: {
1564: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
1566: PetscFunctionBegin;
1567: switch (op) {
1568: case MAT_ROW_ORIENTED:
1569: a->roworiented = flg;
1570: break;
1571: case MAT_KEEP_NONZERO_PATTERN:
1572: a->keepnonzeropattern = flg;
1573: break;
1574: case MAT_NEW_NONZERO_LOCATIONS:
1575: a->nonew = (flg ? 0 : 1);
1576: break;
1577: case MAT_NEW_NONZERO_LOCATION_ERR:
1578: a->nonew = (flg ? -1 : 0);
1579: break;
1580: case MAT_NEW_NONZERO_ALLOCATION_ERR:
1581: a->nonew = (flg ? -2 : 0);
1582: break;
1583: case MAT_UNUSED_NONZERO_LOCATION_ERR:
1584: a->nounused = (flg ? -1 : 0);
1585: break;
1586: case MAT_STRUCTURE_ONLY:
1587: if (flg) {
1588: PetscCall(MatXAIJDeallocatea(A, &a->a));
1589: a->a = NULL;
1590: }
1591: break;
1592: default:
1593: break;
1594: }
1595: PetscFunctionReturn(PETSC_SUCCESS);
1596: }
1598: /* used for both SeqBAIJ and SeqSBAIJ matrices */
1599: PetscErrorCode MatGetRow_SeqBAIJ_private(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v, PetscInt *ai, PetscInt *aj, PetscScalar *aa)
1600: {
1601: PetscInt itmp, i, j, k, M, bn, bp, *idx_i, bs, bs2;
1602: MatScalar *aa_i;
1603: PetscScalar *v_i;
1605: PetscFunctionBegin;
1606: bs = A->rmap->bs;
1607: bs2 = bs * bs;
1608: PetscCheck(row >= 0 && row < A->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row %" PetscInt_FMT " out of range", row);
1610: bn = row / bs; /* Block number */
1611: bp = row % bs; /* Block Position */
1612: M = ai[bn + 1] - ai[bn];
1613: *nz = bs * M;
1615: if (v) {
1616: *v = NULL;
1617: if (*nz) {
1618: PetscCall(PetscMalloc1(*nz, v));
1619: for (i = 0; i < M; i++) { /* for each block in the block row */
1620: v_i = *v + i * bs;
1621: aa_i = aa + bs2 * (ai[bn] + i);
1622: for (j = bp, k = 0; j < bs2; j += bs, k++) v_i[k] = aa_i[j];
1623: }
1624: }
1625: }
1627: if (idx) {
1628: *idx = NULL;
1629: if (*nz) {
1630: PetscCall(PetscMalloc1(*nz, idx));
1631: for (i = 0; i < M; i++) { /* for each block in the block row */
1632: idx_i = *idx + i * bs;
1633: itmp = bs * aj[ai[bn] + i];
1634: for (j = 0; j < bs; j++) idx_i[j] = itmp++;
1635: }
1636: }
1637: }
1638: PetscFunctionReturn(PETSC_SUCCESS);
1639: }
1641: PetscErrorCode MatGetRow_SeqBAIJ(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
1642: {
1643: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
1645: PetscFunctionBegin;
1646: PetscCall(MatGetRow_SeqBAIJ_private(A, row, nz, idx, v, a->i, a->j, a->a));
1647: PetscFunctionReturn(PETSC_SUCCESS);
1648: }
1650: PetscErrorCode MatRestoreRow_SeqBAIJ(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
1651: {
1652: PetscFunctionBegin;
1653: if (idx) PetscCall(PetscFree(*idx));
1654: if (v) PetscCall(PetscFree(*v));
1655: PetscFunctionReturn(PETSC_SUCCESS);
1656: }
1658: static PetscErrorCode MatTranspose_SeqBAIJ(Mat A, MatReuse reuse, Mat *B)
1659: {
1660: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *at;
1661: Mat C;
1662: PetscInt i, j, k, *aj = a->j, *ai = a->i, bs = A->rmap->bs, mbs = a->mbs, nbs = a->nbs, *atfill;
1663: PetscInt bs2 = a->bs2, *ati, *atj, anzj, kr;
1664: MatScalar *ata, *aa = a->a;
1666: PetscFunctionBegin;
1667: if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *B));
1668: PetscCall(PetscCalloc1(1 + nbs, &atfill));
1669: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
1670: for (i = 0; i < ai[mbs]; i++) atfill[aj[i]] += 1; /* count num of non-zeros in row aj[i] */
1672: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C));
1673: PetscCall(MatSetSizes(C, A->cmap->n, A->rmap->N, A->cmap->n, A->rmap->N));
1674: PetscCall(MatSetType(C, ((PetscObject)A)->type_name));
1675: PetscCall(MatSeqBAIJSetPreallocation(C, bs, 0, atfill));
1677: at = (Mat_SeqBAIJ *)C->data;
1678: ati = at->i;
1679: for (i = 0; i < nbs; i++) at->ilen[i] = at->imax[i] = ati[i + 1] - ati[i];
1680: } else {
1681: C = *B;
1682: at = (Mat_SeqBAIJ *)C->data;
1683: ati = at->i;
1684: }
1686: atj = at->j;
1687: ata = at->a;
1689: /* Copy ati into atfill so we have locations of the next free space in atj */
1690: PetscCall(PetscArraycpy(atfill, ati, nbs));
1692: /* Walk through A row-wise and mark nonzero entries of A^T. */
1693: for (i = 0; i < mbs; i++) {
1694: anzj = ai[i + 1] - ai[i];
1695: for (j = 0; j < anzj; j++) {
1696: atj[atfill[*aj]] = i;
1697: for (kr = 0; kr < bs; kr++) {
1698: for (k = 0; k < bs; k++) ata[bs2 * atfill[*aj] + k * bs + kr] = *aa++;
1699: }
1700: atfill[*aj++] += 1;
1701: }
1702: }
1703: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
1704: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
1706: /* Clean up temporary space and complete requests. */
1707: PetscCall(PetscFree(atfill));
1709: if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
1710: PetscCall(MatSetBlockSizes(C, A->cmap->bs, A->rmap->bs));
1711: *B = C;
1712: } else {
1713: PetscCall(MatHeaderMerge(A, &C));
1714: }
1715: PetscFunctionReturn(PETSC_SUCCESS);
1716: }
1718: static PetscErrorCode MatCompare_SeqBAIJ_Private(Mat A, Mat B, PetscReal tol, PetscBool *flg)
1719: {
1720: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)B->data;
1722: PetscFunctionBegin;
1723: /* If the matrix/block dimensions are not equal, or no of nonzeros or shift */
1724: if (A->rmap->N != B->rmap->N || A->cmap->n != B->cmap->n || A->rmap->bs != B->rmap->bs || a->nz != b->nz) {
1725: *flg = PETSC_FALSE;
1726: PetscFunctionReturn(PETSC_SUCCESS);
1727: }
1729: /* if the a->i are the same */
1730: PetscCall(PetscArraycmp(a->i, b->i, a->mbs + 1, flg));
1731: if (!*flg) PetscFunctionReturn(PETSC_SUCCESS);
1733: /* if a->j are the same */
1734: PetscCall(PetscArraycmp(a->j, b->j, a->nz, flg));
1735: if (!*flg) PetscFunctionReturn(PETSC_SUCCESS);
1737: if (tol == 0.0) PetscCall(PetscArraycmp(a->a, b->a, a->nz * A->rmap->bs * A->rmap->bs, flg)); /* if a->a are the same */
1738: else {
1739: *flg = PETSC_TRUE;
1740: for (PetscInt i = 0; (i < a->nz * A->rmap->bs * A->rmap->bs) && *flg; ++i)
1741: if (PetscAbsScalar(a->a[i] - b->a[i]) > tol) *flg = PETSC_FALSE;
1742: }
1743: PetscFunctionReturn(PETSC_SUCCESS);
1744: }
1746: static PetscErrorCode MatIsTranspose_SeqBAIJ(Mat A, Mat B, PetscReal tol, PetscBool *f)
1747: {
1748: Mat Btrans;
1750: PetscFunctionBegin;
1751: PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &Btrans));
1752: PetscCall(MatCompare_SeqBAIJ_Private(A, Btrans, tol, f));
1753: PetscCall(MatDestroy(&Btrans));
1754: PetscFunctionReturn(PETSC_SUCCESS);
1755: }
1757: static PetscErrorCode MatEqual_SeqBAIJ(Mat A, Mat B, PetscBool *flg)
1758: {
1759: PetscFunctionBegin;
1760: PetscCall(MatCompare_SeqBAIJ_Private(A, B, 0.0, flg));
1761: PetscFunctionReturn(PETSC_SUCCESS);
1762: }
1764: /* Used for both SeqBAIJ and SeqSBAIJ matrices */
1765: PetscErrorCode MatView_SeqBAIJ_Binary(Mat mat, PetscViewer viewer)
1766: {
1767: Mat_SeqBAIJ *A = (Mat_SeqBAIJ *)mat->data;
1768: PetscInt header[4], M, N, m, bs, nz, cnt, i, j, k, l;
1769: PetscInt *rowlens, *colidxs;
1770: PetscScalar *matvals;
1772: PetscFunctionBegin;
1773: PetscCall(PetscViewerSetUp(viewer));
1775: M = mat->rmap->N;
1776: N = mat->cmap->N;
1777: m = mat->rmap->n;
1778: bs = mat->rmap->bs;
1779: nz = bs * bs * A->nz;
1781: /* write matrix header */
1782: header[0] = MAT_FILE_CLASSID;
1783: header[1] = M;
1784: header[2] = N;
1785: header[3] = nz;
1786: PetscCall(PetscViewerBinaryWrite(viewer, header, 4, PETSC_INT));
1788: /* store row lengths */
1789: PetscCall(PetscMalloc1(m, &rowlens));
1790: for (cnt = 0, i = 0; i < A->mbs; i++)
1791: for (j = 0; j < bs; j++) rowlens[cnt++] = bs * (A->i[i + 1] - A->i[i]);
1792: PetscCall(PetscViewerBinaryWrite(viewer, rowlens, m, PETSC_INT));
1793: PetscCall(PetscFree(rowlens));
1795: /* store column indices */
1796: PetscCall(PetscMalloc1(nz, &colidxs));
1797: for (cnt = 0, i = 0; i < A->mbs; i++)
1798: for (k = 0; k < bs; k++)
1799: for (j = A->i[i]; j < A->i[i + 1]; j++)
1800: for (l = 0; l < bs; l++) colidxs[cnt++] = bs * A->j[j] + l;
1801: PetscCheck(cnt == nz, PETSC_COMM_SELF, PETSC_ERR_LIB, "Internal PETSc error: cnt = %" PetscInt_FMT " nz = %" PetscInt_FMT, cnt, nz);
1802: PetscCall(PetscViewerBinaryWrite(viewer, colidxs, nz, PETSC_INT));
1803: PetscCall(PetscFree(colidxs));
1805: /* store nonzero values */
1806: PetscCall(PetscMalloc1(nz, &matvals));
1807: for (cnt = 0, i = 0; i < A->mbs; i++)
1808: for (k = 0; k < bs; k++)
1809: for (j = A->i[i]; j < A->i[i + 1]; j++)
1810: for (l = 0; l < bs; l++) matvals[cnt++] = A->a[bs * (bs * j + l) + k];
1811: PetscCheck(cnt == nz, PETSC_COMM_SELF, PETSC_ERR_LIB, "Internal PETSc error: cnt = %" PetscInt_FMT " nz = %" PetscInt_FMT, cnt, nz);
1812: PetscCall(PetscViewerBinaryWrite(viewer, matvals, nz, PETSC_SCALAR));
1813: PetscCall(PetscFree(matvals));
1815: /* write block size option to the viewer's .info file */
1816: PetscCall(MatView_Binary_BlockSizes(mat, viewer));
1817: PetscFunctionReturn(PETSC_SUCCESS);
1818: }
1820: static PetscErrorCode MatView_SeqBAIJ_ASCII_structonly(Mat A, PetscViewer viewer)
1821: {
1822: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
1823: PetscInt i, bs = A->rmap->bs, k;
1825: PetscFunctionBegin;
1826: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
1827: for (i = 0; i < a->mbs; i++) {
1828: PetscCall(PetscViewerASCIIPrintf(viewer, "row %" PetscInt_FMT "-%" PetscInt_FMT ":", i * bs, i * bs + bs - 1));
1829: for (k = a->i[i]; k < a->i[i + 1]; k++) PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT "-%" PetscInt_FMT ") ", bs * a->j[k], bs * a->j[k] + bs - 1));
1830: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
1831: }
1832: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
1833: PetscFunctionReturn(PETSC_SUCCESS);
1834: }
1836: static PetscErrorCode MatView_SeqBAIJ_ASCII(Mat A, PetscViewer viewer)
1837: {
1838: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
1839: PetscInt i, j, bs = A->rmap->bs, k, l, bs2 = a->bs2;
1840: PetscViewerFormat format;
1842: PetscFunctionBegin;
1843: if (A->structure_only) {
1844: PetscCall(MatView_SeqBAIJ_ASCII_structonly(A, viewer));
1845: PetscFunctionReturn(PETSC_SUCCESS);
1846: }
1848: PetscCall(PetscViewerGetFormat(viewer, &format));
1849: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
1850: } else if (format == PETSC_VIEWER_ASCII_MATLAB) {
1851: const char *matname;
1852: Mat aij;
1853: PetscCall(MatConvert(A, MATSEQAIJ, MAT_INITIAL_MATRIX, &aij));
1854: PetscCall(PetscObjectGetName((PetscObject)A, &matname));
1855: PetscCall(PetscObjectSetName((PetscObject)aij, matname));
1856: PetscCall(MatView(aij, viewer));
1857: PetscCall(MatDestroy(&aij));
1858: } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) {
1859: PetscFunctionReturn(PETSC_SUCCESS);
1860: } else if (format == PETSC_VIEWER_ASCII_COMMON) {
1861: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
1862: for (i = 0; i < a->mbs; i++) {
1863: for (j = 0; j < bs; j++) {
1864: PetscCall(PetscViewerASCIIPrintf(viewer, "row %" PetscInt_FMT ":", i * bs + j));
1865: for (k = a->i[i]; k < a->i[i + 1]; k++) {
1866: for (l = 0; l < bs; l++) {
1867: if (PetscDefined(USE_COMPLEX) && PetscImaginaryPart(a->a[bs2 * k + l * bs + j]) > 0.0 && PetscRealPart(a->a[bs2 * k + l * bs + j]) != 0.0) {
1868: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g + %gi) ", bs * a->j[k] + l, (double)PetscRealPart(a->a[bs2 * k + l * bs + j]), (double)PetscImaginaryPart(a->a[bs2 * k + l * bs + j])));
1869: } else if (PetscDefined(USE_COMPLEX) && PetscImaginaryPart(a->a[bs2 * k + l * bs + j]) < 0.0 && PetscRealPart(a->a[bs2 * k + l * bs + j]) != 0.0) {
1870: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g - %gi) ", bs * a->j[k] + l, (double)PetscRealPart(a->a[bs2 * k + l * bs + j]), -(double)PetscImaginaryPart(a->a[bs2 * k + l * bs + j])));
1871: } else if (PetscRealPart(a->a[bs2 * k + l * bs + j]) != 0.0) {
1872: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", bs * a->j[k] + l, (double)PetscRealPart(a->a[bs2 * k + l * bs + j])));
1873: }
1874: }
1875: }
1876: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
1877: }
1878: }
1879: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
1880: } else {
1881: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
1882: for (i = 0; i < a->mbs; i++) {
1883: for (j = 0; j < bs; j++) {
1884: PetscCall(PetscViewerASCIIPrintf(viewer, "row %" PetscInt_FMT ":", i * bs + j));
1885: for (k = a->i[i]; k < a->i[i + 1]; k++) {
1886: for (l = 0; l < bs; l++) {
1887: if (PetscDefined(USE_COMPLEX) && PetscImaginaryPart(a->a[bs2 * k + l * bs + j]) > 0.0) {
1888: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g + %g i) ", bs * a->j[k] + l, (double)PetscRealPart(a->a[bs2 * k + l * bs + j]), (double)PetscImaginaryPart(a->a[bs2 * k + l * bs + j])));
1889: } else if (PetscDefined(USE_COMPLEX) && PetscImaginaryPart(a->a[bs2 * k + l * bs + j]) < 0.0) {
1890: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g - %g i) ", bs * a->j[k] + l, (double)PetscRealPart(a->a[bs2 * k + l * bs + j]), -(double)PetscImaginaryPart(a->a[bs2 * k + l * bs + j])));
1891: } else {
1892: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", bs * a->j[k] + l, (double)PetscRealPart(a->a[bs2 * k + l * bs + j])));
1893: }
1894: }
1895: }
1896: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
1897: }
1898: }
1899: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
1900: }
1901: PetscCall(PetscViewerFlush(viewer));
1902: PetscFunctionReturn(PETSC_SUCCESS);
1903: }
1905: #include <petscdraw.h>
1906: #if defined(__GNUC__) && !defined(__clang__)
1907: #pragma GCC diagnostic push
1908: #pragma GCC diagnostic ignored "-Wclobbered"
1909: #endif
1910: static PetscErrorCode MatView_SeqBAIJ_Draw_Zoom(PetscDraw draw, void *Aa)
1911: {
1912: Mat A = (Mat)Aa;
1913: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
1914: PetscInt row, i, j, k, l, mbs = a->mbs, bs = A->rmap->bs, bs2 = a->bs2;
1915: PetscReal xl, yl, xr, yr, x_l, x_r, y_l, y_r;
1916: MatScalar *aa;
1917: PetscViewer viewer;
1918: PetscViewerFormat format;
1919: int color;
1921: PetscFunctionBegin;
1922: PetscCall(PetscObjectQuery((PetscObject)A, "Zoomviewer", (PetscObject *)&viewer));
1923: PetscCall(PetscViewerGetFormat(viewer, &format));
1924: PetscCall(PetscDrawGetCoordinates(draw, &xl, &yl, &xr, &yr));
1926: /* loop over matrix elements drawing boxes */
1928: if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1929: PetscDrawCollectiveBegin(draw);
1930: /* Blue for negative, Cyan for zero and Red for positive */
1931: color = PETSC_DRAW_BLUE;
1932: for (i = 0, row = 0; i < mbs; i++, row += bs) {
1933: for (j = a->i[i]; j < a->i[i + 1]; j++) {
1934: y_l = A->rmap->N - row - 1.0;
1935: y_r = y_l + 1.0;
1936: x_l = a->j[j] * bs;
1937: x_r = x_l + 1.0;
1938: aa = a->a + j * bs2;
1939: for (k = 0; k < bs; k++) {
1940: for (l = 0; l < bs; l++) {
1941: if (PetscRealPart(*aa++) >= 0.) continue;
1942: PetscCall(PetscDrawRectangle(draw, x_l + k, y_l - l, x_r + k, y_r - l, color, color, color, color));
1943: }
1944: }
1945: }
1946: }
1947: color = PETSC_DRAW_CYAN;
1948: for (i = 0, row = 0; i < mbs; i++, row += bs) {
1949: for (j = a->i[i]; j < a->i[i + 1]; j++) {
1950: y_l = A->rmap->N - row - 1.0;
1951: y_r = y_l + 1.0;
1952: x_l = a->j[j] * bs;
1953: x_r = x_l + 1.0;
1954: aa = a->a + j * bs2;
1955: for (k = 0; k < bs; k++) {
1956: for (l = 0; l < bs; l++) {
1957: if (PetscRealPart(*aa++) != 0.) continue;
1958: PetscCall(PetscDrawRectangle(draw, x_l + k, y_l - l, x_r + k, y_r - l, color, color, color, color));
1959: }
1960: }
1961: }
1962: }
1963: color = PETSC_DRAW_RED;
1964: for (i = 0, row = 0; i < mbs; i++, row += bs) {
1965: for (j = a->i[i]; j < a->i[i + 1]; j++) {
1966: y_l = A->rmap->N - row - 1.0;
1967: y_r = y_l + 1.0;
1968: x_l = a->j[j] * bs;
1969: x_r = x_l + 1.0;
1970: aa = a->a + j * bs2;
1971: for (k = 0; k < bs; k++) {
1972: for (l = 0; l < bs; l++) {
1973: if (PetscRealPart(*aa++) <= 0.) continue;
1974: PetscCall(PetscDrawRectangle(draw, x_l + k, y_l - l, x_r + k, y_r - l, color, color, color, color));
1975: }
1976: }
1977: }
1978: }
1979: PetscDrawCollectiveEnd(draw);
1980: } else {
1981: /* use contour shading to indicate magnitude of values */
1982: /* first determine max of all nonzero values */
1983: PetscReal minv = 0.0, maxv = 0.0;
1984: PetscDraw popup;
1986: for (i = 0; i < a->nz * a->bs2; i++) {
1987: if (PetscAbsScalar(a->a[i]) > maxv) maxv = PetscAbsScalar(a->a[i]);
1988: }
1989: if (minv >= maxv) maxv = minv + PETSC_SMALL;
1990: PetscCall(PetscDrawGetPopup(draw, &popup));
1991: PetscCall(PetscDrawScalePopup(popup, 0.0, maxv));
1993: PetscDrawCollectiveBegin(draw);
1994: for (i = 0, row = 0; i < mbs; i++, row += bs) {
1995: for (j = a->i[i]; j < a->i[i + 1]; j++) {
1996: y_l = A->rmap->N - row - 1.0;
1997: y_r = y_l + 1.0;
1998: x_l = a->j[j] * bs;
1999: x_r = x_l + 1.0;
2000: aa = a->a + j * bs2;
2001: for (k = 0; k < bs; k++) {
2002: for (l = 0; l < bs; l++) {
2003: MatScalar v = *aa++;
2004: color = PetscDrawRealToColor(PetscAbsScalar(v), minv, maxv);
2005: PetscCall(PetscDrawRectangle(draw, x_l + k, y_l - l, x_r + k, y_r - l, color, color, color, color));
2006: }
2007: }
2008: }
2009: }
2010: PetscDrawCollectiveEnd(draw);
2011: }
2012: PetscFunctionReturn(PETSC_SUCCESS);
2013: }
2014: #if defined(__GNUC__) && !defined(__clang__)
2015: #pragma GCC diagnostic pop
2016: #endif
2018: static PetscErrorCode MatView_SeqBAIJ_Draw(Mat A, PetscViewer viewer)
2019: {
2020: PetscReal xl, yl, xr, yr, w, h;
2021: PetscDraw draw;
2022: PetscBool isnull;
2024: PetscFunctionBegin;
2025: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
2026: PetscCall(PetscDrawIsNull(draw, &isnull));
2027: if (isnull) PetscFunctionReturn(PETSC_SUCCESS);
2029: xr = A->cmap->n;
2030: yr = A->rmap->N;
2031: h = yr / 10.0;
2032: w = xr / 10.0;
2033: xr += w;
2034: yr += h;
2035: xl = -w;
2036: yl = -h;
2037: PetscCall(PetscDrawSetCoordinates(draw, xl, yl, xr, yr));
2038: PetscCall(PetscObjectCompose((PetscObject)A, "Zoomviewer", (PetscObject)viewer));
2039: PetscCall(PetscDrawZoom(draw, MatView_SeqBAIJ_Draw_Zoom, A));
2040: PetscCall(PetscObjectCompose((PetscObject)A, "Zoomviewer", NULL));
2041: PetscCall(PetscDrawSave(draw));
2042: PetscFunctionReturn(PETSC_SUCCESS);
2043: }
2045: PetscErrorCode MatView_SeqBAIJ(Mat A, PetscViewer viewer)
2046: {
2047: PetscBool isascii, isbinary, isdraw;
2049: PetscFunctionBegin;
2050: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2051: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
2052: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
2053: if (isascii) {
2054: PetscCall(MatView_SeqBAIJ_ASCII(A, viewer));
2055: } else if (isbinary) {
2056: PetscCall(MatView_SeqBAIJ_Binary(A, viewer));
2057: } else if (isdraw) {
2058: PetscCall(MatView_SeqBAIJ_Draw(A, viewer));
2059: } else {
2060: Mat B;
2061: PetscCall(MatConvert(A, MATSEQAIJ, MAT_INITIAL_MATRIX, &B));
2062: PetscCall(MatView(B, viewer));
2063: PetscCall(MatDestroy(&B));
2064: }
2065: PetscFunctionReturn(PETSC_SUCCESS);
2066: }
2068: PetscErrorCode MatGetValues_SeqBAIJ(Mat A, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], PetscScalar v[])
2069: {
2070: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
2071: PetscInt *rp, k, low, high, t, row, nrow, i, col, l, *aj = a->j;
2072: PetscInt *ai = a->i, *ailen = a->ilen;
2073: PetscInt brow, bcol, ridx, cidx, bs = A->rmap->bs, bs2 = a->bs2;
2074: MatScalar *ap, *aa = a->a;
2075: PetscBool roworiented = a->roworiented;
2076: PetscScalar *value;
2078: PetscFunctionBegin;
2079: for (k = 0; k < m; k++) { /* loop over rows */
2080: row = im[k];
2081: if (row < 0) continue; /* negative row */
2082: brow = row / bs;
2083: PetscCheck(row < A->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row %" PetscInt_FMT " too large", row);
2084: rp = PetscSafePointerPlusOffset(aj, ai[brow]);
2085: ap = PetscSafePointerPlusOffset(aa, bs2 * ai[brow]);
2086: nrow = ailen[brow];
2087: for (l = 0; l < n; l++) { /* loop over columns */
2088: if (in[l] < 0) continue; /* negative column */
2089: PetscCheck(in[l] < A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column %" PetscInt_FMT " too large", in[l]);
2090: value = roworiented ? &v[l + k * n] : &v[k + l * m];
2091: col = in[l];
2092: bcol = col / bs;
2093: cidx = col % bs;
2094: ridx = row % bs;
2095: high = nrow;
2096: low = 0; /* assume unsorted */
2097: while (high - low > 5) {
2098: t = (low + high) / 2;
2099: if (rp[t] > bcol) high = t;
2100: else low = t;
2101: }
2102: for (i = low; i < high; i++) {
2103: if (rp[i] > bcol) break;
2104: if (rp[i] == bcol) {
2105: *value = ap[bs2 * i + bs * cidx + ridx];
2106: goto finished;
2107: }
2108: }
2109: *value = 0.0;
2110: finished:;
2111: }
2112: }
2113: PetscFunctionReturn(PETSC_SUCCESS);
2114: }
2116: PetscErrorCode MatSetValuesBlocked_SeqBAIJ(Mat A, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], const PetscScalar v[], InsertMode is)
2117: {
2118: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
2119: PetscInt *rp, k, low, high, t, ii, jj, row, nrow, i, col, l, rmax, N, lastcol = -1;
2120: PetscInt *imax = a->imax, *ai = a->i, *ailen = a->ilen;
2121: PetscInt *aj = a->j, nonew = a->nonew, bs2 = a->bs2, bs = A->rmap->bs, stepval;
2122: PetscBool roworiented = a->roworiented;
2123: const PetscScalar *value = v;
2124: MatScalar *ap = NULL, *aa = a->a, *bap;
2126: PetscFunctionBegin;
2127: if (roworiented) {
2128: stepval = (n - 1) * bs;
2129: } else {
2130: stepval = (m - 1) * bs;
2131: }
2132: for (k = 0; k < m; k++) { /* loop over added rows */
2133: row = im[k];
2134: if (row < 0) continue;
2135: PetscCheck(row < a->mbs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Block row index too large %" PetscInt_FMT " max %" PetscInt_FMT, row, a->mbs - 1);
2136: rp = aj + ai[row];
2137: if (!A->structure_only) ap = aa + bs2 * ai[row];
2138: rmax = imax[row];
2139: nrow = ailen[row];
2140: low = 0;
2141: high = nrow;
2142: for (l = 0; l < n; l++) { /* loop over added columns */
2143: if (in[l] < 0) continue;
2144: PetscCheck(in[l] < a->nbs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Block column index too large %" PetscInt_FMT " max %" PetscInt_FMT, in[l], a->nbs - 1);
2145: col = in[l];
2146: if (!A->structure_only) {
2147: if (roworiented) {
2148: value = v + (k * (stepval + bs) + l) * bs;
2149: } else {
2150: value = v + (l * (stepval + bs) + k) * bs;
2151: }
2152: }
2153: if (col <= lastcol) low = 0;
2154: else high = nrow;
2155: lastcol = col;
2156: while (high - low > 7) {
2157: t = (low + high) / 2;
2158: if (rp[t] > col) high = t;
2159: else low = t;
2160: }
2161: for (i = low; i < high; i++) {
2162: if (rp[i] > col) break;
2163: if (rp[i] == col) {
2164: if (A->structure_only) goto noinsert2;
2165: bap = ap + bs2 * i;
2166: if (roworiented) {
2167: if (is == ADD_VALUES) {
2168: for (ii = 0; ii < bs; ii++, value += stepval) {
2169: for (jj = ii; jj < bs2; jj += bs) bap[jj] += *value++;
2170: }
2171: } else {
2172: for (ii = 0; ii < bs; ii++, value += stepval) {
2173: for (jj = ii; jj < bs2; jj += bs) bap[jj] = *value++;
2174: }
2175: }
2176: } else {
2177: if (is == ADD_VALUES) {
2178: for (ii = 0; ii < bs; ii++, value += bs + stepval) {
2179: for (jj = 0; jj < bs; jj++) bap[jj] += value[jj];
2180: bap += bs;
2181: }
2182: } else {
2183: for (ii = 0; ii < bs; ii++, value += bs + stepval) {
2184: for (jj = 0; jj < bs; jj++) bap[jj] = value[jj];
2185: bap += bs;
2186: }
2187: }
2188: }
2189: goto noinsert2;
2190: }
2191: }
2192: if (nonew == 1) goto noinsert2;
2193: PetscCheck(nonew != -1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Inserting a new blocked index new nonzero block (%" PetscInt_FMT ", %" PetscInt_FMT ") in the matrix", row, col);
2194: if (A->structure_only) {
2195: MatSeqXAIJReallocateAIJ_structure_only(A, a->mbs, bs2, nrow, row, col, rmax, ai, aj, rp, imax, nonew, MatScalar);
2196: } else {
2197: MatSeqXAIJReallocateAIJ(A, a->mbs, bs2, nrow, row, col, rmax, aa, ai, aj, rp, ap, imax, nonew, MatScalar);
2198: }
2199: N = nrow++ - 1;
2200: high++;
2201: /* shift up all the later entries in this row */
2202: PetscCall(PetscArraymove(rp + i + 1, rp + i, N - i + 1));
2203: rp[i] = col;
2204: if (!A->structure_only) {
2205: PetscCall(PetscArraymove(ap + bs2 * (i + 1), ap + bs2 * i, bs2 * (N - i + 1)));
2206: bap = ap + bs2 * i;
2207: if (roworiented) {
2208: for (ii = 0; ii < bs; ii++, value += stepval) {
2209: for (jj = ii; jj < bs2; jj += bs) bap[jj] = *value++;
2210: }
2211: } else {
2212: for (ii = 0; ii < bs; ii++, value += stepval) {
2213: for (jj = 0; jj < bs; jj++) *bap++ = *value++;
2214: }
2215: }
2216: }
2217: noinsert2:;
2218: low = i;
2219: }
2220: ailen[row] = nrow;
2221: }
2222: PetscFunctionReturn(PETSC_SUCCESS);
2223: }
2225: PetscErrorCode MatAssemblyEnd_SeqBAIJ(Mat A, MatAssemblyType mode)
2226: {
2227: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
2228: PetscInt fshift = 0, i, *ai = a->i, *aj = a->j, *imax = a->imax;
2229: PetscInt m = A->rmap->N, *ip, N, *ailen = a->ilen;
2230: PetscInt mbs = a->mbs, bs2 = a->bs2, rmax = 0;
2231: MatScalar *aa = a->a, *ap;
2232: PetscReal ratio = 0.6;
2234: PetscFunctionBegin;
2235: if (mode == MAT_FLUSH_ASSEMBLY || (A->was_assembled && A->ass_nonzerostate == A->nonzerostate)) PetscFunctionReturn(PETSC_SUCCESS);
2237: if (m) rmax = ailen[0];
2238: for (i = 1; i < mbs; i++) {
2239: /* move each row back by the amount of empty slots (fshift) before it*/
2240: fshift += imax[i - 1] - ailen[i - 1];
2241: rmax = PetscMax(rmax, ailen[i]);
2242: if (fshift) {
2243: ip = aj + ai[i];
2244: N = ailen[i];
2245: PetscCall(PetscArraymove(ip - fshift, ip, N));
2246: if (!A->structure_only) {
2247: ap = aa + bs2 * ai[i];
2248: PetscCall(PetscArraymove(ap - bs2 * fshift, ap, bs2 * N));
2249: }
2250: }
2251: ai[i] = ai[i - 1] + ailen[i - 1];
2252: }
2253: if (mbs) {
2254: fshift += imax[mbs - 1] - ailen[mbs - 1];
2255: ai[mbs] = ai[mbs - 1] + ailen[mbs - 1];
2256: }
2258: /* reset ilen and imax for each row */
2259: a->nonzerorowcnt = 0;
2260: for (i = 0; i < mbs; i++) {
2261: ailen[i] = imax[i] = ai[i + 1] - ai[i];
2262: a->nonzerorowcnt += (ailen[i] > 0);
2263: }
2264: a->nz = ai[mbs];
2266: if (fshift) PetscCheck(a->nounused != -1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unused space detected in matrix: %" PetscInt_FMT " X %" PetscInt_FMT " block size %" PetscInt_FMT ", %" PetscInt_FMT " unneeded", m, A->cmap->n, A->rmap->bs, fshift * bs2);
2267: PetscCall(PetscInfo(A, "Matrix size: %" PetscInt_FMT " X %" PetscInt_FMT ", block size %" PetscInt_FMT "; storage space: %" PetscInt_FMT " unneeded, %" PetscInt_FMT " used\n", m, A->cmap->n, A->rmap->bs, fshift * bs2, a->nz * bs2));
2268: PetscCall(PetscInfo(A, "Number of mallocs during MatSetValues is %" PetscInt_FMT "\n", a->reallocs));
2269: PetscCall(PetscInfo(A, "Most nonzeros blocks in any row is %" PetscInt_FMT "\n", rmax));
2271: A->info.mallocs += a->reallocs;
2272: a->reallocs = 0;
2273: A->info.nz_unneeded = (PetscReal)fshift * bs2;
2274: a->rmax = rmax;
2276: if (!A->structure_only) PetscCall(MatCheckCompressedRow(A, a->nonzerorowcnt, &a->compressedrow, a->i, mbs, ratio));
2277: PetscFunctionReturn(PETSC_SUCCESS);
2278: }
2280: /*
2281: This function returns an array of flags which indicate the locations of contiguous
2282: blocks that should be zeroed. for eg: if bs = 3 and is = [0,1,2,3,5,6,7,8,9]
2283: then the resulting sizes = [3,1,1,3,1] corresponding to sets [(0,1,2),(3),(5),(6,7,8),(9)]
2284: Assume: sizes should be long enough to hold all the values.
2285: */
2286: static PetscErrorCode MatZeroRows_SeqBAIJ_Check_Blocks(PetscInt idx[], PetscInt n, PetscInt bs, PetscInt sizes[], PetscInt *bs_max)
2287: {
2288: PetscInt j = 0;
2290: PetscFunctionBegin;
2291: for (PetscInt i = 0; i < n; j++) {
2292: PetscInt row = idx[i];
2293: if (row % bs != 0) { /* Not the beginning of a block */
2294: sizes[j] = 1;
2295: i++;
2296: } else if (i + bs > n) { /* complete block doesn't exist (at idx end) */
2297: sizes[j] = 1; /* Also makes sure at least 'bs' values exist for next else */
2298: i++;
2299: } else { /* Beginning of the block, so check if the complete block exists */
2300: PetscBool flg = PETSC_TRUE;
2301: for (PetscInt k = 1; k < bs; k++) {
2302: if (row + k != idx[i + k]) { /* break in the block */
2303: flg = PETSC_FALSE;
2304: break;
2305: }
2306: }
2307: if (flg) { /* No break in the bs */
2308: sizes[j] = bs;
2309: i += bs;
2310: } else {
2311: sizes[j] = 1;
2312: i++;
2313: }
2314: }
2315: }
2316: *bs_max = j;
2317: PetscFunctionReturn(PETSC_SUCCESS);
2318: }
2320: PetscErrorCode MatZeroRows_SeqBAIJ(Mat A, PetscInt is_n, const PetscInt is_idx[], PetscScalar diag, Vec x, Vec b)
2321: {
2322: Mat_SeqBAIJ *baij = (Mat_SeqBAIJ *)A->data;
2323: PetscInt i, j, k, count, *rows;
2324: PetscInt bs = A->rmap->bs, bs2 = baij->bs2, *sizes, row, bs_max;
2325: PetscScalar zero = 0.0;
2326: MatScalar *aa;
2327: const PetscScalar *xx;
2328: PetscScalar *bb;
2330: PetscFunctionBegin;
2331: /* fix right-hand side if needed */
2332: if (x && b) {
2333: PetscCall(VecGetArrayRead(x, &xx));
2334: PetscCall(VecGetArray(b, &bb));
2335: for (i = 0; i < is_n; i++) bb[is_idx[i]] = diag * xx[is_idx[i]];
2336: PetscCall(VecRestoreArrayRead(x, &xx));
2337: PetscCall(VecRestoreArray(b, &bb));
2338: }
2340: /* Make a copy of the IS and sort it */
2341: /* allocate memory for rows,sizes */
2342: PetscCall(PetscMalloc2(is_n, &rows, 2 * is_n, &sizes));
2344: /* copy IS values to rows, and sort them */
2345: for (i = 0; i < is_n; i++) rows[i] = is_idx[i];
2346: PetscCall(PetscSortInt(is_n, rows));
2348: if (baij->keepnonzeropattern) {
2349: for (i = 0; i < is_n; i++) sizes[i] = 1;
2350: bs_max = is_n;
2351: } else {
2352: PetscCall(MatZeroRows_SeqBAIJ_Check_Blocks(rows, is_n, bs, sizes, &bs_max));
2353: A->nonzerostate++;
2354: }
2356: for (i = 0, j = 0; i < bs_max; j += sizes[i], i++) {
2357: row = rows[j];
2358: PetscCheck(row >= 0 && row <= A->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "row %" PetscInt_FMT " out of range", row);
2359: count = (baij->i[row / bs + 1] - baij->i[row / bs]) * bs;
2360: aa = PetscSafePointerPlusOffset(baij->a, baij->i[row / bs] * bs2 + (row % bs));
2361: if (sizes[i] == bs && !baij->keepnonzeropattern) {
2362: if (diag != (PetscScalar)0.0) {
2363: if (baij->ilen[row / bs] > 0) {
2364: baij->ilen[row / bs] = 1;
2365: baij->j[baij->i[row / bs]] = row / bs;
2367: PetscCall(PetscArrayzero(aa, count * bs));
2368: }
2369: /* Now insert all the diagonal values for this bs */
2370: for (k = 0; k < bs; k++) PetscUseTypeMethod(A, setvalues, 1, rows + j + k, 1, rows + j + k, &diag, INSERT_VALUES);
2371: } else { /* (diag == 0.0) */
2372: baij->ilen[row / bs] = 0;
2373: } /* end (diag == 0.0) */
2374: } else { /* (sizes[i] != bs) */
2375: PetscAssert(sizes[i] == 1, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Internal Error. Value should be 1");
2376: for (k = 0; k < count; k++) {
2377: aa[0] = zero;
2378: aa += bs;
2379: }
2380: if (diag != (PetscScalar)0.0) PetscUseTypeMethod(A, setvalues, 1, rows + j, 1, rows + j, &diag, INSERT_VALUES);
2381: }
2382: }
2384: PetscCall(PetscFree2(rows, sizes));
2385: PetscCall(MatAssemblyEnd_SeqBAIJ(A, MAT_FINAL_ASSEMBLY));
2386: PetscFunctionReturn(PETSC_SUCCESS);
2387: }
2389: static PetscErrorCode MatZeroRowsColumns_SeqBAIJ(Mat A, PetscInt is_n, const PetscInt is_idx[], PetscScalar diag, Vec x, Vec b)
2390: {
2391: Mat_SeqBAIJ *baij = (Mat_SeqBAIJ *)A->data;
2392: PetscInt i, j, k, count;
2393: PetscInt bs = A->rmap->bs, bs2 = baij->bs2, row, col;
2394: PetscScalar zero = 0.0;
2395: MatScalar *aa;
2396: const PetscScalar *xx;
2397: PetscScalar *bb;
2398: PetscBool *zeroed, vecs = PETSC_FALSE;
2400: PetscFunctionBegin;
2401: /* fix right-hand side if needed */
2402: if (x && b) {
2403: PetscCall(VecGetArrayRead(x, &xx));
2404: PetscCall(VecGetArray(b, &bb));
2405: vecs = PETSC_TRUE;
2406: }
2408: /* zero the columns */
2409: PetscCall(PetscCalloc1(A->rmap->n, &zeroed));
2410: for (i = 0; i < is_n; i++) {
2411: PetscCheck(is_idx[i] >= 0 && is_idx[i] < A->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "row %" PetscInt_FMT " out of range", is_idx[i]);
2412: zeroed[is_idx[i]] = PETSC_TRUE;
2413: }
2414: for (i = 0; i < A->rmap->N; i++) {
2415: if (!zeroed[i]) {
2416: row = i / bs;
2417: for (j = baij->i[row]; j < baij->i[row + 1]; j++) {
2418: for (k = 0; k < bs; k++) {
2419: col = bs * baij->j[j] + k;
2420: if (zeroed[col]) {
2421: aa = baij->a + j * bs2 + (i % bs) + bs * k;
2422: if (vecs) bb[i] -= aa[0] * xx[col];
2423: aa[0] = 0.0;
2424: }
2425: }
2426: }
2427: } else if (vecs) bb[i] = diag * xx[i];
2428: }
2429: PetscCall(PetscFree(zeroed));
2430: if (vecs) {
2431: PetscCall(VecRestoreArrayRead(x, &xx));
2432: PetscCall(VecRestoreArray(b, &bb));
2433: }
2435: /* zero the rows */
2436: for (i = 0; i < is_n; i++) {
2437: row = is_idx[i];
2438: count = (baij->i[row / bs + 1] - baij->i[row / bs]) * bs;
2439: aa = PetscSafePointerPlusOffset(baij->a, baij->i[row / bs] * bs2 + (row % bs));
2440: for (k = 0; k < count; k++) {
2441: aa[0] = zero;
2442: aa += bs;
2443: }
2444: if (diag != (PetscScalar)0.0) PetscUseTypeMethod(A, setvalues, 1, &row, 1, &row, &diag, INSERT_VALUES);
2445: }
2446: PetscCall(MatAssemblyEnd_SeqBAIJ(A, MAT_FINAL_ASSEMBLY));
2447: PetscFunctionReturn(PETSC_SUCCESS);
2448: }
2450: PetscErrorCode MatSetValues_SeqBAIJ(Mat A, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], const PetscScalar v[], InsertMode is)
2451: {
2452: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
2453: PetscInt *rp, k, low, high, t, ii, row, nrow, i, col, l, rmax, N, lastcol = -1;
2454: PetscInt *imax = a->imax, *ai = a->i, *ailen = a->ilen;
2455: PetscInt *aj = a->j, nonew = a->nonew, bs = A->rmap->bs, brow, bcol;
2456: PetscInt ridx, cidx, bs2 = a->bs2;
2457: PetscBool roworiented = a->roworiented;
2458: MatScalar *ap = NULL, value = 0.0, *aa = a->a, *bap;
2460: PetscFunctionBegin;
2461: for (k = 0; k < m; k++) { /* loop over added rows */
2462: row = im[k];
2463: brow = row / bs;
2464: if (row < 0) continue;
2465: PetscCheck(row < A->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, row, A->rmap->N - 1);
2466: rp = PetscSafePointerPlusOffset(aj, ai[brow]);
2467: if (!A->structure_only) ap = PetscSafePointerPlusOffset(aa, bs2 * ai[brow]);
2468: rmax = imax[brow];
2469: nrow = ailen[brow];
2470: low = 0;
2471: high = nrow;
2472: for (l = 0; l < n; l++) { /* loop over added columns */
2473: if (in[l] < 0) continue;
2474: PetscCheck(in[l] < A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT, in[l], A->cmap->n - 1);
2475: col = in[l];
2476: bcol = col / bs;
2477: ridx = row % bs;
2478: cidx = col % bs;
2479: if (!A->structure_only) {
2480: if (roworiented) {
2481: value = v[l + k * n];
2482: } else {
2483: value = v[k + l * m];
2484: }
2485: }
2486: if (col <= lastcol) low = 0;
2487: else high = nrow;
2488: lastcol = col;
2489: while (high - low > 7) {
2490: t = (low + high) / 2;
2491: if (rp[t] > bcol) high = t;
2492: else low = t;
2493: }
2494: for (i = low; i < high; i++) {
2495: if (rp[i] > bcol) break;
2496: if (rp[i] == bcol) {
2497: bap = PetscSafePointerPlusOffset(ap, bs2 * i + bs * cidx + ridx);
2498: if (!A->structure_only) {
2499: if (is == ADD_VALUES) *bap += value;
2500: else *bap = value;
2501: }
2502: goto noinsert1;
2503: }
2504: }
2505: if (nonew == 1) goto noinsert1;
2506: PetscCheck(nonew != -1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Inserting a new nonzero (%" PetscInt_FMT ", %" PetscInt_FMT ") in the matrix", row, col);
2507: if (A->structure_only) {
2508: MatSeqXAIJReallocateAIJ_structure_only(A, a->mbs, bs2, nrow, brow, bcol, rmax, ai, aj, rp, imax, nonew, MatScalar);
2509: } else {
2510: MatSeqXAIJReallocateAIJ(A, a->mbs, bs2, nrow, brow, bcol, rmax, aa, ai, aj, rp, ap, imax, nonew, MatScalar);
2511: }
2512: N = nrow++ - 1;
2513: high++;
2514: /* shift up all the later entries in this row */
2515: PetscCall(PetscArraymove(rp + i + 1, rp + i, N - i + 1));
2516: rp[i] = bcol;
2517: if (!A->structure_only) {
2518: PetscCall(PetscArraymove(ap + bs2 * (i + 1), ap + bs2 * i, bs2 * (N - i + 1)));
2519: PetscCall(PetscArrayzero(ap + bs2 * i, bs2));
2520: ap[bs2 * i + bs * cidx + ridx] = value;
2521: }
2522: a->nz++;
2523: noinsert1:;
2524: low = i;
2525: }
2526: ailen[brow] = nrow;
2527: }
2528: PetscFunctionReturn(PETSC_SUCCESS);
2529: }
2531: static PetscErrorCode MatILUFactor_SeqBAIJ(Mat inA, IS row, IS col, const MatFactorInfo *info)
2532: {
2533: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)inA->data;
2534: Mat outA;
2535: PetscBool row_identity, col_identity;
2537: PetscFunctionBegin;
2538: PetscCheck(info->levels == 0, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only levels = 0 supported for in-place ILU");
2539: PetscCall(ISIdentity(row, &row_identity));
2540: PetscCall(ISIdentity(col, &col_identity));
2541: PetscCheck(row_identity && col_identity, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Row and column permutations must be identity for in-place ILU");
2543: outA = inA;
2544: inA->factortype = MAT_FACTOR_LU;
2545: PetscCall(PetscFree(inA->solvertype));
2546: PetscCall(PetscStrallocpy(MATSOLVERPETSC, &inA->solvertype));
2548: PetscCall(PetscObjectReference((PetscObject)row));
2549: PetscCall(ISDestroy(&a->row));
2550: a->row = row;
2551: PetscCall(PetscObjectReference((PetscObject)col));
2552: PetscCall(ISDestroy(&a->col));
2553: a->col = col;
2555: /* Create the invert permutation so that it can be used in MatLUFactorNumeric() */
2556: PetscCall(ISDestroy(&a->icol));
2557: PetscCall(ISInvertPermutation(col, PETSC_DECIDE, &a->icol));
2559: PetscCall(MatSeqBAIJSetNumericFactorization_inplace(inA, (PetscBool)(row_identity && col_identity)));
2560: if (!a->solve_work) PetscCall(PetscMalloc1(inA->rmap->N + inA->rmap->bs, &a->solve_work));
2561: PetscCall(MatLUFactorNumeric(outA, inA, info));
2562: PetscFunctionReturn(PETSC_SUCCESS);
2563: }
2565: static PetscErrorCode MatSeqBAIJSetColumnIndices_SeqBAIJ(Mat mat, const PetscInt *indices)
2566: {
2567: Mat_SeqBAIJ *baij = (Mat_SeqBAIJ *)mat->data;
2569: PetscFunctionBegin;
2570: baij->nz = baij->maxnz;
2571: PetscCall(PetscArraycpy(baij->j, indices, baij->nz));
2572: PetscCall(PetscArraycpy(baij->ilen, baij->imax, baij->mbs));
2573: PetscFunctionReturn(PETSC_SUCCESS);
2574: }
2576: /*@
2577: MatSeqBAIJSetColumnIndices - Set the column indices for all the block rows in the matrix.
2579: Input Parameters:
2580: + mat - the `MATSEQBAIJ` matrix
2581: - indices - the block column indices
2583: Level: advanced
2585: Notes:
2586: This can be called if you have precomputed the nonzero structure of the
2587: matrix and want to provide it to the matrix object to improve the performance
2588: of the `MatSetValues()` operation.
2590: You MUST have set the correct numbers of nonzeros per row in the call to
2591: `MatCreateSeqBAIJ()`, and the columns indices MUST be sorted.
2593: MUST be called before any calls to `MatSetValues()`
2595: .seealso: [](ch_matrices), `Mat`, `MATSEQBAIJ`, `MatSetValues()`
2596: @*/
2597: PetscErrorCode MatSeqBAIJSetColumnIndices(Mat mat, PetscInt *indices)
2598: {
2599: PetscFunctionBegin;
2601: PetscAssertPointer(indices, 2);
2602: PetscUseMethod(mat, "MatSeqBAIJSetColumnIndices_C", (Mat, const PetscInt *), (mat, (const PetscInt *)indices));
2603: PetscFunctionReturn(PETSC_SUCCESS);
2604: }
2606: static PetscErrorCode MatGetRowMaxAbs_SeqBAIJ(Mat A, Vec v, PetscInt idx[])
2607: {
2608: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
2609: PetscInt i, j, n, row, bs, *ai, *aj, mbs;
2610: PetscReal atmp;
2611: PetscScalar *x, zero = 0.0;
2612: MatScalar *aa;
2613: PetscInt ncols, brow, krow, kcol;
2615: PetscFunctionBegin;
2616: PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
2617: bs = A->rmap->bs;
2618: aa = a->a;
2619: ai = a->i;
2620: aj = a->j;
2621: mbs = a->mbs;
2623: PetscCall(VecSet(v, zero));
2624: PetscCall(VecGetArray(v, &x));
2625: PetscCall(VecGetLocalSize(v, &n));
2626: PetscCheck(n == A->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector");
2627: for (i = 0; i < mbs; i++) {
2628: ncols = ai[1] - ai[0];
2629: ai++;
2630: brow = bs * i;
2631: for (j = 0; j < ncols; j++) {
2632: for (kcol = 0; kcol < bs; kcol++) {
2633: for (krow = 0; krow < bs; krow++) {
2634: atmp = PetscAbsScalar(*aa);
2635: aa++;
2636: row = brow + krow; /* row index */
2637: if (PetscAbsScalar(x[row]) < atmp) {
2638: x[row] = atmp;
2639: if (idx) idx[row] = bs * (*aj) + kcol;
2640: }
2641: }
2642: }
2643: aj++;
2644: }
2645: }
2646: PetscCall(VecRestoreArray(v, &x));
2647: PetscFunctionReturn(PETSC_SUCCESS);
2648: }
2650: static PetscErrorCode MatGetRowSumAbs_SeqBAIJ(Mat A, Vec v)
2651: {
2652: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
2653: PetscInt i, j, n, row, bs, *ai, mbs;
2654: PetscReal atmp;
2655: PetscScalar *x, zero = 0.0;
2656: MatScalar *aa;
2657: PetscInt ncols, brow, krow, kcol;
2659: PetscFunctionBegin;
2660: PetscCheck(!A->factortype, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
2661: bs = A->rmap->bs;
2662: aa = a->a;
2663: ai = a->i;
2664: mbs = a->mbs;
2666: PetscCall(VecSet(v, zero));
2667: PetscCall(VecGetArrayWrite(v, &x));
2668: PetscCall(VecGetLocalSize(v, &n));
2669: PetscCheck(n == A->rmap->N, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector");
2670: for (i = 0; i < mbs; i++) {
2671: ncols = ai[1] - ai[0];
2672: ai++;
2673: brow = bs * i;
2674: for (j = 0; j < ncols; j++) {
2675: for (kcol = 0; kcol < bs; kcol++) {
2676: for (krow = 0; krow < bs; krow++) {
2677: atmp = PetscAbsScalar(*aa);
2678: aa++;
2679: row = brow + krow; /* row index */
2680: x[row] += atmp;
2681: }
2682: }
2683: }
2684: }
2685: PetscCall(VecRestoreArrayWrite(v, &x));
2686: PetscFunctionReturn(PETSC_SUCCESS);
2687: }
2689: static PetscErrorCode MatCopy_SeqBAIJ(Mat A, Mat B, MatStructure str)
2690: {
2691: PetscFunctionBegin;
2692: /* If the two matrices have the same copy implementation, use fast copy. */
2693: if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) {
2694: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
2695: Mat_SeqBAIJ *b = (Mat_SeqBAIJ *)B->data;
2696: PetscInt ambs = a->mbs, bmbs = b->mbs, abs = A->rmap->bs, bbs = B->rmap->bs, bs2 = abs * abs;
2698: PetscCheck(a->i[ambs] == b->i[bmbs], PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of nonzero blocks in matrices A %" PetscInt_FMT " and B %" PetscInt_FMT " are different", a->i[ambs], b->i[bmbs]);
2699: PetscCheck(abs == bbs, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Block size A %" PetscInt_FMT " and B %" PetscInt_FMT " are different", abs, bbs);
2700: PetscCall(PetscArraycpy(b->a, a->a, bs2 * a->i[ambs]));
2701: PetscCall(PetscObjectStateIncrease((PetscObject)B));
2702: } else {
2703: PetscCall(MatCopy_Basic(A, B, str));
2704: }
2705: PetscFunctionReturn(PETSC_SUCCESS);
2706: }
2708: static PetscErrorCode MatSeqBAIJGetArray_SeqBAIJ(Mat A, PetscScalar *array[])
2709: {
2710: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
2712: PetscFunctionBegin;
2713: *array = a->a;
2714: PetscFunctionReturn(PETSC_SUCCESS);
2715: }
2717: static PetscErrorCode MatSeqBAIJRestoreArray_SeqBAIJ(Mat A, PetscScalar *array[])
2718: {
2719: PetscFunctionBegin;
2720: *array = NULL;
2721: PetscFunctionReturn(PETSC_SUCCESS);
2722: }
2724: PetscErrorCode MatAXPYGetPreallocation_SeqBAIJ(Mat Y, Mat X, PetscInt *nnz)
2725: {
2726: PetscInt bs = Y->rmap->bs, mbs = Y->rmap->N / bs;
2727: Mat_SeqBAIJ *x = (Mat_SeqBAIJ *)X->data;
2728: Mat_SeqBAIJ *y = (Mat_SeqBAIJ *)Y->data;
2730: PetscFunctionBegin;
2731: /* Set the number of nonzeros in the new matrix */
2732: PetscCall(MatAXPYGetPreallocation_SeqX_private(mbs, x->i, x->j, y->i, y->j, nnz));
2733: PetscFunctionReturn(PETSC_SUCCESS);
2734: }
2736: PetscErrorCode MatAXPY_SeqBAIJ(Mat Y, PetscScalar a, Mat X, MatStructure str)
2737: {
2738: Mat_SeqBAIJ *x = (Mat_SeqBAIJ *)X->data, *y = (Mat_SeqBAIJ *)Y->data;
2739: PetscInt bs = Y->rmap->bs, bs2 = bs * bs;
2740: PetscBLASInt one = 1;
2742: PetscFunctionBegin;
2743: if (str == UNKNOWN_NONZERO_PATTERN || (PetscDefined(USE_DEBUG) && str == SAME_NONZERO_PATTERN)) {
2744: PetscBool e = x->nz == y->nz && x->mbs == y->mbs && bs == X->rmap->bs ? PETSC_TRUE : PETSC_FALSE;
2745: if (e) {
2746: PetscCall(PetscArraycmp(x->i, y->i, x->mbs + 1, &e));
2747: if (e) {
2748: PetscCall(PetscArraycmp(x->j, y->j, x->i[x->mbs], &e));
2749: if (e) str = SAME_NONZERO_PATTERN;
2750: }
2751: }
2752: if (!e) PetscCheck(str != SAME_NONZERO_PATTERN, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "MatStructure is not SAME_NONZERO_PATTERN");
2753: }
2754: if (str == SAME_NONZERO_PATTERN) {
2755: PetscScalar alpha = a;
2756: PetscBLASInt bnz;
2757: PetscCall(PetscBLASIntCast(x->nz * bs2, &bnz));
2758: PetscCallBLAS("BLASaxpy", BLASaxpy_(&bnz, &alpha, x->a, &one, y->a, &one));
2759: PetscCall(PetscObjectStateIncrease((PetscObject)Y));
2760: } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */
2761: PetscCall(MatAXPY_Basic(Y, a, X, str));
2762: } else {
2763: Mat B;
2764: PetscInt *nnz;
2765: PetscCheck(bs == X->rmap->bs, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Matrices must have same block size");
2766: PetscCall(PetscMalloc1(Y->rmap->N, &nnz));
2767: PetscCall(MatCreate(PetscObjectComm((PetscObject)Y), &B));
2768: PetscCall(PetscObjectSetName((PetscObject)B, ((PetscObject)Y)->name));
2769: PetscCall(MatSetSizes(B, Y->rmap->n, Y->cmap->n, Y->rmap->N, Y->cmap->N));
2770: PetscCall(MatSetBlockSizesFromMats(B, Y, Y));
2771: PetscCall(MatSetType(B, (MatType)((PetscObject)Y)->type_name));
2772: PetscCall(MatAXPYGetPreallocation_SeqBAIJ(Y, X, nnz));
2773: PetscCall(MatSeqBAIJSetPreallocation(B, bs, 0, nnz));
2774: PetscCall(MatAXPY_BasicWithPreallocation(B, Y, a, X, str));
2775: PetscCall(MatHeaderMerge(Y, &B));
2776: PetscCall(PetscFree(nnz));
2777: }
2778: PetscFunctionReturn(PETSC_SUCCESS);
2779: }
2781: PETSC_INTERN PetscErrorCode MatConjugate_SeqBAIJ(Mat A)
2782: {
2783: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
2784: PetscInt i, nz = a->bs2 * a->i[a->mbs];
2785: MatScalar *aa = a->a;
2787: PetscFunctionBegin;
2788: for (i = 0; i < nz; i++) aa[i] = PetscConj(aa[i]);
2789: PetscFunctionReturn(PETSC_SUCCESS);
2790: }
2792: static PetscErrorCode MatRealPart_SeqBAIJ(Mat A)
2793: {
2794: #if PetscDefined(USE_COMPLEX)
2795: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
2796: PetscInt i, nz = a->bs2 * a->i[a->mbs];
2797: MatScalar *aa = a->a;
2799: PetscFunctionBegin;
2800: for (i = 0; i < nz; i++) aa[i] = PetscRealPart(aa[i]);
2801: PetscFunctionReturn(PETSC_SUCCESS);
2802: #else
2803: (void)A;
2804: return PETSC_SUCCESS;
2805: #endif
2806: }
2808: static PetscErrorCode MatImaginaryPart_SeqBAIJ(Mat A)
2809: {
2810: #if PetscDefined(USE_COMPLEX)
2811: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
2812: PetscInt i, nz = a->bs2 * a->i[a->mbs];
2813: MatScalar *aa = a->a;
2815: PetscFunctionBegin;
2816: for (i = 0; i < nz; i++) aa[i] = PetscImaginaryPart(aa[i]);
2817: PetscFunctionReturn(PETSC_SUCCESS);
2818: #else
2819: (void)A;
2820: return PETSC_SUCCESS;
2821: #endif
2822: }
2824: /*
2825: Code almost identical to MatGetColumnIJ_SeqAIJ() should share common code
2826: */
2827: static PetscErrorCode MatGetColumnIJ_SeqBAIJ(Mat A, PetscInt oshift, PetscBool symmetric, PetscBool inodecompressed, PetscInt *nn, const PetscInt *ia[], const PetscInt *ja[], PetscBool *done)
2828: {
2829: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
2830: PetscInt bs = A->rmap->bs, i, *collengths, *cia, *cja, n = A->cmap->n / bs, m = A->rmap->n / bs;
2831: PetscInt nz = a->i[m], row, *jj, mr, col;
2833: PetscFunctionBegin;
2834: *nn = n;
2835: if (!ia) PetscFunctionReturn(PETSC_SUCCESS);
2836: PetscCheck(!symmetric, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not for BAIJ matrices");
2837: PetscCall(PetscCalloc1(n, &collengths));
2838: PetscCall(PetscMalloc1(n + 1, &cia));
2839: PetscCall(PetscMalloc1(nz, &cja));
2840: jj = a->j;
2841: for (i = 0; i < nz; i++) collengths[jj[i]]++;
2842: cia[0] = oshift;
2843: for (i = 0; i < n; i++) cia[i + 1] = cia[i] + collengths[i];
2844: PetscCall(PetscArrayzero(collengths, n));
2845: jj = a->j;
2846: for (row = 0; row < m; row++) {
2847: mr = a->i[row + 1] - a->i[row];
2848: for (i = 0; i < mr; i++) {
2849: col = *jj++;
2851: cja[cia[col] + collengths[col]++ - oshift] = row + oshift;
2852: }
2853: }
2854: PetscCall(PetscFree(collengths));
2855: *ia = cia;
2856: *ja = cja;
2857: PetscFunctionReturn(PETSC_SUCCESS);
2858: }
2860: static PetscErrorCode MatRestoreColumnIJ_SeqBAIJ(Mat A, PetscInt oshift, PetscBool symmetric, PetscBool inodecompressed, PetscInt *n, const PetscInt *ia[], const PetscInt *ja[], PetscBool *done)
2861: {
2862: PetscFunctionBegin;
2863: if (!ia) PetscFunctionReturn(PETSC_SUCCESS);
2864: PetscCall(PetscFree(*ia));
2865: PetscCall(PetscFree(*ja));
2866: PetscFunctionReturn(PETSC_SUCCESS);
2867: }
2869: /*
2870: MatGetColumnIJ_SeqBAIJ_Color() and MatRestoreColumnIJ_SeqBAIJ_Color() are customized from
2871: MatGetColumnIJ_SeqBAIJ() and MatRestoreColumnIJ_SeqBAIJ() by adding an output
2872: spidx[], index of a->a, to be used in MatTransposeColoringCreate() and MatFDColoringCreate()
2873: */
2874: PetscErrorCode MatGetColumnIJ_SeqBAIJ_Color(Mat A, PetscInt oshift, PetscBool symmetric, PetscBool inodecompressed, PetscInt *nn, const PetscInt *ia[], const PetscInt *ja[], PetscInt *spidx[], PetscBool *done)
2875: {
2876: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
2877: PetscInt i, *collengths, *cia, *cja, n = a->nbs, m = a->mbs;
2878: PetscInt nz = a->i[m], row, *jj, mr, col;
2879: PetscInt *cspidx;
2881: PetscFunctionBegin;
2882: *nn = n;
2883: if (!ia) PetscFunctionReturn(PETSC_SUCCESS);
2885: PetscCall(PetscCalloc1(n, &collengths));
2886: PetscCall(PetscMalloc1(n + 1, &cia));
2887: PetscCall(PetscMalloc1(nz, &cja));
2888: PetscCall(PetscMalloc1(nz, &cspidx));
2889: jj = a->j;
2890: for (i = 0; i < nz; i++) collengths[jj[i]]++;
2891: cia[0] = oshift;
2892: for (i = 0; i < n; i++) cia[i + 1] = cia[i] + collengths[i];
2893: PetscCall(PetscArrayzero(collengths, n));
2894: jj = a->j;
2895: for (row = 0; row < m; row++) {
2896: mr = a->i[row + 1] - a->i[row];
2897: for (i = 0; i < mr; i++) {
2898: col = *jj++;
2899: cspidx[cia[col] + collengths[col] - oshift] = a->i[row] + i; /* index of a->j */
2900: cja[cia[col] + collengths[col]++ - oshift] = row + oshift;
2901: }
2902: }
2903: PetscCall(PetscFree(collengths));
2904: *ia = cia;
2905: *ja = cja;
2906: *spidx = cspidx;
2907: PetscFunctionReturn(PETSC_SUCCESS);
2908: }
2910: PetscErrorCode MatRestoreColumnIJ_SeqBAIJ_Color(Mat A, PetscInt oshift, PetscBool symmetric, PetscBool inodecompressed, PetscInt *n, const PetscInt *ia[], const PetscInt *ja[], PetscInt *spidx[], PetscBool *done)
2911: {
2912: PetscFunctionBegin;
2913: PetscCall(MatRestoreColumnIJ_SeqBAIJ(A, oshift, symmetric, inodecompressed, n, ia, ja, done));
2914: PetscCall(PetscFree(*spidx));
2915: PetscFunctionReturn(PETSC_SUCCESS);
2916: }
2918: static PetscErrorCode MatShift_SeqBAIJ(Mat Y, PetscScalar a)
2919: {
2920: Mat_SeqBAIJ *aij = (Mat_SeqBAIJ *)Y->data;
2922: PetscFunctionBegin;
2923: if (!Y->preallocated || !aij->nz) PetscCall(MatSeqBAIJSetPreallocation(Y, Y->rmap->bs, 1, NULL));
2924: PetscCall(MatShift_Basic(Y, a));
2925: PetscFunctionReturn(PETSC_SUCCESS);
2926: }
2928: PetscErrorCode MatEliminateZeros_SeqBAIJ(Mat A, PetscBool keep)
2929: {
2930: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
2931: PetscInt fshift = 0, fshift_prev = 0, i, *ai = a->i, *aj = a->j, *imax = a->imax, j, k;
2932: PetscInt m = A->rmap->N, *ailen = a->ilen;
2933: PetscInt mbs = a->mbs, bs2 = a->bs2, rmax = 0;
2934: MatScalar *aa = a->a, *ap;
2935: PetscBool zero;
2937: PetscFunctionBegin;
2938: PetscCheck(A->assembled, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Cannot eliminate zeros for unassembled matrix");
2939: if (m) rmax = ailen[0];
2940: for (i = 1, a->nonzerorowcnt = 0; i <= mbs; i++) {
2941: for (k = ai[i - 1]; k < ai[i]; k++) {
2942: zero = PETSC_TRUE;
2943: ap = aa + bs2 * k;
2944: for (j = 0; j < bs2 && zero; j++) {
2945: if (ap[j] != 0.0) zero = PETSC_FALSE;
2946: }
2947: if (zero && (aj[k] != i - 1 || !keep)) fshift++;
2948: else {
2949: if (zero && aj[k] == i - 1) PetscCall(PetscInfo(A, "Keep the diagonal block at row %" PetscInt_FMT "\n", i - 1));
2950: aj[k - fshift] = aj[k];
2951: PetscCall(PetscArraymove(ap - bs2 * fshift, ap, bs2));
2952: }
2953: }
2954: ai[i - 1] -= fshift_prev;
2955: fshift_prev = fshift;
2956: ailen[i - 1] = imax[i - 1] = ai[i] - fshift - ai[i - 1];
2957: a->nonzerorowcnt += ((ai[i] - fshift - ai[i - 1]) > 0);
2958: rmax = PetscMax(rmax, ailen[i - 1]);
2959: }
2960: if (fshift) {
2961: if (mbs) {
2962: ai[mbs] -= fshift;
2963: a->nz = ai[mbs];
2964: }
2965: PetscCall(PetscInfo(A, "Matrix size: %" PetscInt_FMT " X %" PetscInt_FMT "; zeros eliminated: %" PetscInt_FMT "; nonzeros left: %" PetscInt_FMT "\n", m, A->cmap->n, fshift, a->nz));
2966: A->nonzerostate++;
2967: A->info.nz_unneeded += (PetscReal)fshift;
2968: a->rmax = rmax;
2969: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
2970: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
2971: }
2972: PetscFunctionReturn(PETSC_SUCCESS);
2973: }
2975: static struct _MatOps MatOps_Values = {MatSetValues_SeqBAIJ,
2976: MatGetRow_SeqBAIJ,
2977: MatRestoreRow_SeqBAIJ,
2978: MatMult_SeqBAIJ_N,
2979: /* 4*/ MatMultAdd_SeqBAIJ_N,
2980: MatMultTranspose_SeqBAIJ,
2981: MatMultTransposeAdd_SeqBAIJ,
2982: NULL,
2983: NULL,
2984: NULL,
2985: /* 10*/ NULL,
2986: MatLUFactor_SeqBAIJ,
2987: NULL,
2988: NULL,
2989: MatTranspose_SeqBAIJ,
2990: /* 15*/ MatGetInfo_SeqBAIJ,
2991: MatEqual_SeqBAIJ,
2992: MatGetDiagonal_SeqBAIJ,
2993: MatDiagonalScale_SeqBAIJ,
2994: MatNorm_SeqBAIJ,
2995: /* 20*/ NULL,
2996: MatAssemblyEnd_SeqBAIJ,
2997: MatSetOption_SeqBAIJ,
2998: MatZeroEntries_SeqBAIJ,
2999: /* 24*/ MatZeroRows_SeqBAIJ,
3000: NULL,
3001: NULL,
3002: NULL,
3003: NULL,
3004: /* 29*/ MatSetUp_Seq_Hash,
3005: NULL,
3006: NULL,
3007: NULL,
3008: NULL,
3009: /* 34*/ MatDuplicate_SeqBAIJ,
3010: NULL,
3011: NULL,
3012: MatILUFactor_SeqBAIJ,
3013: NULL,
3014: /* 39*/ MatAXPY_SeqBAIJ,
3015: MatCreateSubMatrices_SeqBAIJ,
3016: MatIncreaseOverlap_SeqBAIJ,
3017: MatGetValues_SeqBAIJ,
3018: MatCopy_SeqBAIJ,
3019: /* 44*/ NULL,
3020: MatScale_SeqBAIJ,
3021: MatShift_SeqBAIJ,
3022: NULL,
3023: MatZeroRowsColumns_SeqBAIJ,
3024: /* 49*/ NULL,
3025: MatGetRowIJ_SeqBAIJ,
3026: MatRestoreRowIJ_SeqBAIJ,
3027: MatGetColumnIJ_SeqBAIJ,
3028: MatRestoreColumnIJ_SeqBAIJ,
3029: /* 54*/ MatFDColoringCreate_SeqXAIJ,
3030: NULL,
3031: NULL,
3032: NULL,
3033: MatSetValuesBlocked_SeqBAIJ,
3034: /* 59*/ MatCreateSubMatrix_SeqBAIJ,
3035: MatDestroy_SeqBAIJ,
3036: MatView_SeqBAIJ,
3037: NULL,
3038: NULL,
3039: /* 64*/ NULL,
3040: NULL,
3041: NULL,
3042: NULL,
3043: MatGetRowMaxAbs_SeqBAIJ,
3044: /* 69*/ NULL,
3045: MatConvert_Basic,
3046: NULL,
3047: MatFDColoringApply_BAIJ,
3048: NULL,
3049: /* 74*/ NULL,
3050: NULL,
3051: NULL,
3052: NULL,
3053: MatLoad_SeqBAIJ,
3054: /* 79*/ NULL,
3055: NULL,
3056: NULL,
3057: NULL,
3058: NULL,
3059: /* 84*/ NULL,
3060: NULL,
3061: NULL,
3062: NULL,
3063: NULL,
3064: /* 89*/ NULL,
3065: NULL,
3066: NULL,
3067: NULL,
3068: MatConjugate_SeqBAIJ,
3069: /* 94*/ NULL,
3070: NULL,
3071: MatRealPart_SeqBAIJ,
3072: MatImaginaryPart_SeqBAIJ,
3073: NULL,
3074: /* 99*/ NULL,
3075: NULL,
3076: NULL,
3077: NULL,
3078: NULL,
3079: /*104*/ NULL,
3080: NULL,
3081: NULL,
3082: NULL,
3083: NULL,
3084: /*109*/ NULL,
3085: NULL,
3086: MatMultHermitianTranspose_SeqBAIJ,
3087: MatMultHermitianTransposeAdd_SeqBAIJ,
3088: NULL,
3089: /*114*/ NULL,
3090: MatGetColumnReductions_SeqBAIJ,
3091: MatInvertBlockDiagonal_SeqBAIJ,
3092: NULL,
3093: NULL,
3094: /*119*/ NULL,
3095: NULL,
3096: NULL,
3097: NULL,
3098: NULL,
3099: /*124*/ NULL,
3100: MatSetBlockSizes_Default,
3101: NULL,
3102: MatFDColoringSetUp_SeqXAIJ,
3103: NULL,
3104: /*129*/ MatCreateMPIMatConcatenateSeqMat_SeqBAIJ,
3105: MatDestroySubMatrices_SeqBAIJ,
3106: NULL,
3107: NULL,
3108: NULL,
3109: /*134*/ NULL,
3110: MatEliminateZeros_SeqBAIJ,
3111: MatGetRowSumAbs_SeqBAIJ,
3112: NULL,
3113: NULL,
3114: /*139*/ NULL,
3115: MatCopyHashToXAIJ_Seq_Hash,
3116: NULL,
3117: NULL,
3118: NULL,
3119: /*144*/ NULL,
3120: NULL,
3121: NULL,
3122: NULL};
3124: static PetscErrorCode MatStoreValues_SeqBAIJ(Mat mat)
3125: {
3126: Mat_SeqBAIJ *aij = (Mat_SeqBAIJ *)mat->data;
3127: PetscInt nz = aij->i[aij->mbs] * aij->bs2;
3129: PetscFunctionBegin;
3130: PetscCheck(aij->nonew == 1, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first");
3132: /* allocate space for values if not already there */
3133: if (!aij->saved_values) PetscCall(PetscMalloc1(nz + 1, &aij->saved_values));
3135: /* copy values over */
3136: PetscCall(PetscArraycpy(aij->saved_values, aij->a, nz));
3137: PetscFunctionReturn(PETSC_SUCCESS);
3138: }
3140: static PetscErrorCode MatRetrieveValues_SeqBAIJ(Mat mat)
3141: {
3142: Mat_SeqBAIJ *aij = (Mat_SeqBAIJ *)mat->data;
3143: PetscInt nz = aij->i[aij->mbs] * aij->bs2;
3145: PetscFunctionBegin;
3146: PetscCheck(aij->nonew == 1, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first");
3147: PetscCheck(aij->saved_values, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call MatStoreValues(A);first");
3149: /* copy values over */
3150: PetscCall(PetscArraycpy(aij->a, aij->saved_values, nz));
3151: PetscFunctionReturn(PETSC_SUCCESS);
3152: }
3154: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqAIJ(Mat, MatType, MatReuse, Mat *);
3155: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqSBAIJ(Mat, MatType, MatReuse, Mat *);
3157: PetscErrorCode MatSeqBAIJSetPreallocation_SeqBAIJ(Mat B, PetscInt bs, PetscInt nz, const PetscInt nnz[])
3158: {
3159: Mat_SeqBAIJ *b = (Mat_SeqBAIJ *)B->data;
3160: PetscInt i, mbs, nbs, bs2;
3161: PetscBool flg = PETSC_FALSE, skipallocation = PETSC_FALSE, realalloc = PETSC_FALSE;
3163: PetscFunctionBegin;
3164: if (B->hash_active) {
3165: PetscInt bs;
3166: B->ops[0] = b->cops;
3167: PetscCall(PetscHMapIJVDestroy(&b->ht));
3168: PetscCall(MatGetBlockSize(B, &bs));
3169: if (bs > 1) PetscCall(PetscHSetIJDestroy(&b->bht));
3170: PetscCall(PetscFree(b->dnz));
3171: PetscCall(PetscFree(b->bdnz));
3172: B->hash_active = PETSC_FALSE;
3173: }
3174: if (nz >= 0 || nnz) realalloc = PETSC_TRUE;
3175: if (nz == MAT_SKIP_ALLOCATION) {
3176: skipallocation = PETSC_TRUE;
3177: nz = 0;
3178: }
3180: PetscCall(MatSetBlockSize(B, bs));
3181: PetscCall(PetscLayoutSetUp(B->rmap));
3182: PetscCall(PetscLayoutSetUp(B->cmap));
3183: PetscCall(PetscLayoutGetBlockSize(B->rmap, &bs));
3185: B->preallocated = PETSC_TRUE;
3187: mbs = B->rmap->n / bs;
3188: nbs = B->cmap->n / bs;
3189: bs2 = bs * bs;
3191: PetscCheck(mbs * bs == B->rmap->n && nbs * bs == B->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number rows %" PetscInt_FMT ", cols %" PetscInt_FMT " must be divisible by blocksize %" PetscInt_FMT, B->rmap->N, B->cmap->n, bs);
3193: if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5;
3194: PetscCheck(nz >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "nz cannot be less than 0: value %" PetscInt_FMT, nz);
3195: if (nnz) {
3196: for (i = 0; i < mbs; i++) {
3197: PetscCheck(nnz[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "nnz cannot be less than 0: local row %" PetscInt_FMT " value %" PetscInt_FMT, i, nnz[i]);
3198: PetscCheck(nnz[i] <= nbs, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "nnz cannot be greater than block row length: local row %" PetscInt_FMT " value %" PetscInt_FMT " rowlength %" PetscInt_FMT, i, nnz[i], nbs);
3199: }
3200: }
3202: PetscOptionsBegin(PetscObjectComm((PetscObject)B), NULL, "Optimize options for SEQBAIJ matrix 2 ", "Mat");
3203: PetscCall(PetscOptionsBool("-mat_no_unroll", "Do not optimize for block size (slow)", NULL, flg, &flg, NULL));
3204: PetscOptionsEnd();
3206: if (!flg) {
3207: switch (bs) {
3208: case 1:
3209: B->ops->mult = MatMult_SeqBAIJ_1;
3210: B->ops->multadd = MatMultAdd_SeqBAIJ_1;
3211: break;
3212: case 2:
3213: B->ops->mult = MatMult_SeqBAIJ_2;
3214: B->ops->multadd = MatMultAdd_SeqBAIJ_2;
3215: break;
3216: case 3:
3217: B->ops->mult = MatMult_SeqBAIJ_3;
3218: B->ops->multadd = MatMultAdd_SeqBAIJ_3;
3219: break;
3220: case 4:
3221: B->ops->mult = MatMult_SeqBAIJ_4;
3222: B->ops->multadd = MatMultAdd_SeqBAIJ_4;
3223: break;
3224: case 5:
3225: B->ops->mult = MatMult_SeqBAIJ_5;
3226: B->ops->multadd = MatMultAdd_SeqBAIJ_5;
3227: break;
3228: case 6:
3229: B->ops->mult = MatMult_SeqBAIJ_6;
3230: B->ops->multadd = MatMultAdd_SeqBAIJ_6;
3231: break;
3232: case 7:
3233: B->ops->mult = MatMult_SeqBAIJ_7;
3234: B->ops->multadd = MatMultAdd_SeqBAIJ_7;
3235: break;
3236: case 9: {
3237: PetscInt version = 1;
3238: PetscCall(PetscOptionsGetInt(NULL, ((PetscObject)B)->prefix, "-mat_baij_mult_version", &version, NULL));
3239: switch (version) {
3240: #if PetscDefined(HAVE_IMMINTRIN_H) && defined(__AVX2__) && defined(__FMA__) && PetscDefined(USE_REAL_DOUBLE) && !PetscDefined(USE_COMPLEX) && !PetscDefined(USE_64BIT_INDICES)
3241: case 1:
3242: B->ops->mult = MatMult_SeqBAIJ_9_AVX2;
3243: B->ops->multadd = MatMultAdd_SeqBAIJ_9_AVX2;
3244: PetscCall(PetscInfo(B, "Using AVX2 for MatMult for BAIJ for blocksize %" PetscInt_FMT "\n", bs));
3245: break;
3246: #endif
3247: default:
3248: B->ops->mult = MatMult_SeqBAIJ_N;
3249: B->ops->multadd = MatMultAdd_SeqBAIJ_N;
3250: PetscCall(PetscInfo(B, "Using BLAS for MatMult for BAIJ for blocksize %" PetscInt_FMT "\n", bs));
3251: break;
3252: }
3253: break;
3254: }
3255: case 11:
3256: B->ops->mult = MatMult_SeqBAIJ_11;
3257: B->ops->multadd = MatMultAdd_SeqBAIJ_11;
3258: break;
3259: case 12: {
3260: PetscInt version = 1;
3261: PetscCall(PetscOptionsGetInt(NULL, ((PetscObject)B)->prefix, "-mat_baij_mult_version", &version, NULL));
3262: switch (version) {
3263: case 1:
3264: B->ops->mult = MatMult_SeqBAIJ_12_ver1;
3265: B->ops->multadd = MatMultAdd_SeqBAIJ_12_ver1;
3266: PetscCall(PetscInfo(B, "Using version %" PetscInt_FMT " of MatMult for BAIJ for blocksize %" PetscInt_FMT "\n", version, bs));
3267: break;
3268: case 2:
3269: B->ops->mult = MatMult_SeqBAIJ_12_ver2;
3270: B->ops->multadd = MatMultAdd_SeqBAIJ_12_ver2;
3271: PetscCall(PetscInfo(B, "Using version %" PetscInt_FMT " of MatMult for BAIJ for blocksize %" PetscInt_FMT "\n", version, bs));
3272: break;
3273: #if PetscDefined(HAVE_IMMINTRIN_H) && defined(__AVX2__) && defined(__FMA__) && PetscDefined(USE_REAL_DOUBLE) && !PetscDefined(USE_COMPLEX) && !PetscDefined(USE_64BIT_INDICES)
3274: case 3:
3275: B->ops->mult = MatMult_SeqBAIJ_12_AVX2;
3276: B->ops->multadd = MatMultAdd_SeqBAIJ_12_ver1;
3277: PetscCall(PetscInfo(B, "Using AVX2 for MatMult for BAIJ for blocksize %" PetscInt_FMT "\n", bs));
3278: break;
3279: #endif
3280: default:
3281: B->ops->mult = MatMult_SeqBAIJ_N;
3282: B->ops->multadd = MatMultAdd_SeqBAIJ_N;
3283: PetscCall(PetscInfo(B, "Using BLAS for MatMult for BAIJ for blocksize %" PetscInt_FMT "\n", bs));
3284: break;
3285: }
3286: break;
3287: }
3288: case 15: {
3289: PetscInt version = 1;
3290: PetscCall(PetscOptionsGetInt(NULL, ((PetscObject)B)->prefix, "-mat_baij_mult_version", &version, NULL));
3291: switch (version) {
3292: case 1:
3293: B->ops->mult = MatMult_SeqBAIJ_15_ver1;
3294: PetscCall(PetscInfo(B, "Using version %" PetscInt_FMT " of MatMult for BAIJ for blocksize %" PetscInt_FMT "\n", version, bs));
3295: break;
3296: case 2:
3297: B->ops->mult = MatMult_SeqBAIJ_15_ver2;
3298: PetscCall(PetscInfo(B, "Using version %" PetscInt_FMT " of MatMult for BAIJ for blocksize %" PetscInt_FMT "\n", version, bs));
3299: break;
3300: case 3:
3301: B->ops->mult = MatMult_SeqBAIJ_15_ver3;
3302: PetscCall(PetscInfo(B, "Using version %" PetscInt_FMT " of MatMult for BAIJ for blocksize %" PetscInt_FMT "\n", version, bs));
3303: break;
3304: case 4:
3305: B->ops->mult = MatMult_SeqBAIJ_15_ver4;
3306: PetscCall(PetscInfo(B, "Using version %" PetscInt_FMT " of MatMult for BAIJ for blocksize %" PetscInt_FMT "\n", version, bs));
3307: break;
3308: default:
3309: B->ops->mult = MatMult_SeqBAIJ_N;
3310: PetscCall(PetscInfo(B, "Using BLAS for MatMult for BAIJ for blocksize %" PetscInt_FMT "\n", bs));
3311: break;
3312: }
3313: B->ops->multadd = MatMultAdd_SeqBAIJ_N;
3314: break;
3315: }
3316: default:
3317: B->ops->mult = MatMult_SeqBAIJ_N;
3318: B->ops->multadd = MatMultAdd_SeqBAIJ_N;
3319: PetscCall(PetscInfo(B, "Using BLAS for MatMult for BAIJ for blocksize %" PetscInt_FMT "\n", bs));
3320: break;
3321: }
3322: }
3323: B->ops->sor = MatSOR_SeqBAIJ;
3324: b->mbs = mbs;
3325: b->nbs = nbs;
3326: if (!skipallocation) {
3327: if (!b->imax) {
3328: PetscCall(PetscMalloc2(mbs, &b->imax, mbs, &b->ilen));
3330: b->free_imax_ilen = PETSC_TRUE;
3331: }
3332: /* b->ilen will count nonzeros in each block row so far. */
3333: for (i = 0; i < mbs; i++) b->ilen[i] = 0;
3334: if (!nnz) {
3335: if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5;
3336: else if (nz < 0) nz = 1;
3337: nz = PetscMin(nz, nbs);
3338: for (i = 0; i < mbs; i++) b->imax[i] = nz;
3339: PetscCall(PetscIntMultError(nz, mbs, &nz));
3340: } else {
3341: PetscInt64 nz64 = 0;
3342: for (i = 0; i < mbs; i++) {
3343: b->imax[i] = nnz[i];
3344: nz64 += nnz[i];
3345: }
3346: PetscCall(PetscIntCast(nz64, &nz));
3347: }
3349: /* allocate the matrix space */
3350: PetscCall(MatSeqXAIJFreeAIJ(B, &b->a, &b->j, &b->i));
3351: PetscCall(PetscShmgetAllocateArray(nz, sizeof(PetscInt), (void **)&b->j));
3352: PetscCall(PetscShmgetAllocateArray(B->rmap->N + 1, sizeof(PetscInt), (void **)&b->i));
3353: if (B->structure_only) {
3354: b->free_a = PETSC_FALSE;
3355: } else {
3356: PetscInt nzbs2 = 0;
3357: PetscCall(PetscIntMultError(nz, bs2, &nzbs2));
3358: PetscCall(PetscShmgetAllocateArray(nzbs2, sizeof(PetscScalar), (void **)&b->a));
3359: b->free_a = PETSC_TRUE;
3360: PetscCall(PetscArrayzero(b->a, nzbs2));
3361: }
3362: b->free_ij = PETSC_TRUE;
3363: PetscCall(PetscArrayzero(b->j, nz));
3365: b->i[0] = 0;
3366: for (i = 1; i < mbs + 1; i++) b->i[i] = b->i[i - 1] + b->imax[i - 1];
3367: } else {
3368: b->free_a = PETSC_FALSE;
3369: b->free_ij = PETSC_FALSE;
3370: }
3372: b->bs2 = bs2;
3373: b->mbs = mbs;
3374: b->nz = 0;
3375: b->maxnz = nz;
3376: B->info.nz_unneeded = (PetscReal)b->maxnz * bs2;
3377: B->was_assembled = PETSC_FALSE;
3378: B->assembled = PETSC_FALSE;
3379: if (realalloc) PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
3380: PetscFunctionReturn(PETSC_SUCCESS);
3381: }
3383: static PetscErrorCode MatSeqBAIJSetPreallocationCSR_SeqBAIJ(Mat B, PetscInt bs, const PetscInt ii[], const PetscInt jj[], const PetscScalar V[])
3384: {
3385: PetscInt i, m, nz, nz_max = 0, *nnz;
3386: PetscScalar *values = NULL;
3387: PetscBool roworiented = ((Mat_SeqBAIJ *)B->data)->roworiented;
3389: PetscFunctionBegin;
3390: PetscCheck(bs >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid block size specified, must be positive but it is %" PetscInt_FMT, bs);
3391: PetscCall(PetscLayoutSetBlockSize(B->rmap, bs));
3392: PetscCall(PetscLayoutSetBlockSize(B->cmap, bs));
3393: PetscCall(PetscLayoutSetUp(B->rmap));
3394: PetscCall(PetscLayoutSetUp(B->cmap));
3395: PetscCall(PetscLayoutGetBlockSize(B->rmap, &bs));
3396: m = B->rmap->n / bs;
3398: PetscCheck(ii[0] == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "ii[0] must be 0 but it is %" PetscInt_FMT, ii[0]);
3399: PetscCall(PetscMalloc1(m + 1, &nnz));
3400: for (i = 0; i < m; i++) {
3401: nz = ii[i + 1] - ii[i];
3402: PetscCheck(nz >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local row %" PetscInt_FMT " has a negative number of columns %" PetscInt_FMT, i, nz);
3403: nz_max = PetscMax(nz_max, nz);
3404: nnz[i] = nz;
3405: }
3406: PetscCall(MatSeqBAIJSetPreallocation(B, bs, 0, nnz));
3407: PetscCall(PetscFree(nnz));
3409: values = (PetscScalar *)V;
3410: if (!values) PetscCall(PetscCalloc1(bs * bs * (nz_max + 1), &values));
3411: for (i = 0; i < m; i++) {
3412: PetscInt ncols = ii[i + 1] - ii[i];
3413: const PetscInt *icols = jj + ii[i];
3414: if (bs == 1 || !roworiented) {
3415: const PetscScalar *svals = values + (V ? (bs * bs * ii[i]) : 0);
3416: PetscCall(MatSetValuesBlocked_SeqBAIJ(B, 1, &i, ncols, icols, svals, INSERT_VALUES));
3417: } else {
3418: for (PetscInt j = 0; j < ncols; j++) {
3419: const PetscScalar *svals = values + (V ? (bs * bs * (ii[i] + j)) : 0);
3420: PetscCall(MatSetValuesBlocked_SeqBAIJ(B, 1, &i, 1, &icols[j], svals, INSERT_VALUES));
3421: }
3422: }
3423: }
3424: if (!V) PetscCall(PetscFree(values));
3425: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
3426: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
3427: PetscCall(MatSetOption(B, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
3428: PetscFunctionReturn(PETSC_SUCCESS);
3429: }
3431: /*@
3432: MatSeqBAIJGetArray - gives read/write access to the array where the data for a `MATSEQBAIJ` matrix is stored
3434: Not Collective
3436: Input Parameter:
3437: . A - a `MATSEQBAIJ` matrix
3439: Output Parameter:
3440: . array - pointer to the data
3442: Level: intermediate
3444: .seealso: [](ch_matrices), `Mat`, `MATSEQBAIJ`, `MatSeqBAIJRestoreArray()`, `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArray()`
3445: @*/
3446: PetscErrorCode MatSeqBAIJGetArray(Mat A, PetscScalar *array[])
3447: {
3448: PetscFunctionBegin;
3449: PetscUseMethod(A, "MatSeqBAIJGetArray_C", (Mat, PetscScalar **), (A, array));
3450: PetscFunctionReturn(PETSC_SUCCESS);
3451: }
3453: /*@
3454: MatSeqBAIJRestoreArray - returns access to the array where the data for a `MATSEQBAIJ` matrix is stored obtained by `MatSeqBAIJGetArray()`
3456: Not Collective
3458: Input Parameters:
3459: + A - a `MATSEQBAIJ` matrix
3460: - array - pointer to the data
3462: Level: intermediate
3464: .seealso: [](ch_matrices), `Mat`, `MatSeqBAIJGetArray()`, `MatSeqAIJGetArray()`, `MatSeqAIJRestoreArray()`
3465: @*/
3466: PetscErrorCode MatSeqBAIJRestoreArray(Mat A, PetscScalar *array[])
3467: {
3468: PetscFunctionBegin;
3469: PetscUseMethod(A, "MatSeqBAIJRestoreArray_C", (Mat, PetscScalar **), (A, array));
3470: PetscCall(PetscObjectStateIncrease((PetscObject)A));
3471: PetscFunctionReturn(PETSC_SUCCESS);
3472: }
3474: /*MC
3475: MATSEQBAIJ - MATSEQBAIJ = "seqbaij" - A matrix type to be used for sequential block sparse matrices, based on
3476: block sparse compressed row format.
3478: Options Database Keys:
3479: + -mat_type seqbaij - sets the matrix type to `MATSEQBAIJ` during a call to `MatSetFromOptions()`
3480: - -mat_baij_mult_version version - indicate the version of the matrix-vector product to use (0 often indicates using BLAS)
3482: Level: beginner
3484: Notes:
3485: Call `MatSetOption(A, MAT_STRUCTURE_ONLY, PETSC_TRUE)` before preallocation or `MatSetUp()` to store only the nonzero pattern.
3486: The assembled matrix has no numerical value array. Row and column indices supplied during insertion are retained, while numerical values are ignored.
3487: Such matrices can be used for structural operations, but not for numerical operations.
3489: Run with `-info` to see what version of the matrix-vector product is being used
3491: .seealso: [](ch_matrices), `Mat`, `MatCreateSeqBAIJ()`
3492: M*/
3494: PETSC_EXTERN PetscErrorCode MatCreate_SeqBAIJ(Mat B)
3495: {
3496: PetscMPIInt size;
3497: Mat_SeqBAIJ *b;
3499: PetscFunctionBegin;
3500: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B), &size));
3501: PetscCheck(size == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Comm must be of size 1");
3503: PetscCall(PetscNew(&b));
3504: B->data = (void *)b;
3505: B->ops[0] = MatOps_Values;
3507: b->row = NULL;
3508: b->col = NULL;
3509: b->icol = NULL;
3510: b->reallocs = 0;
3511: b->saved_values = NULL;
3513: b->roworiented = PETSC_TRUE;
3514: b->nonew = 0;
3515: b->diag = NULL;
3516: B->spptr = NULL;
3517: B->info.nz_unneeded = (PetscReal)b->maxnz * b->bs2;
3518: b->keepnonzeropattern = PETSC_FALSE;
3520: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqBAIJGetArray_C", MatSeqBAIJGetArray_SeqBAIJ));
3521: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqBAIJRestoreArray_C", MatSeqBAIJRestoreArray_SeqBAIJ));
3522: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatStoreValues_C", MatStoreValues_SeqBAIJ));
3523: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatRetrieveValues_C", MatRetrieveValues_SeqBAIJ));
3524: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqBAIJSetColumnIndices_C", MatSeqBAIJSetColumnIndices_SeqBAIJ));
3525: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqbaij_seqaij_C", MatConvert_SeqBAIJ_SeqAIJ));
3526: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqbaij_seqsbaij_C", MatConvert_SeqBAIJ_SeqSBAIJ));
3527: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqBAIJSetPreallocation_C", MatSeqBAIJSetPreallocation_SeqBAIJ));
3528: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqBAIJSetPreallocationCSR_C", MatSeqBAIJSetPreallocationCSR_SeqBAIJ));
3529: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatIsTranspose_C", MatIsTranspose_SeqBAIJ));
3530: #if PetscDefined(HAVE_HYPRE)
3531: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqbaij_hypre_C", MatConvert_AIJ_HYPRE));
3532: #endif
3533: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqbaij_is_C", MatConvert_XAIJ_IS));
3534: #if PetscDefined(HAVE_LIBXSMM)
3535: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqbaij_seqbaijlibxsmm_C", MatConvert_SeqBAIJ_SeqBAIJLIBXSMM));
3536: #endif
3537: PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATSEQBAIJ));
3538: PetscFunctionReturn(PETSC_SUCCESS);
3539: }
3541: PETSC_INTERN PetscErrorCode MatDuplicateNoCreate_SeqBAIJ(Mat C, Mat A, MatDuplicateOption cpvalues, PetscBool mallocmatspace)
3542: {
3543: Mat_SeqBAIJ *c = (Mat_SeqBAIJ *)C->data, *a = (Mat_SeqBAIJ *)A->data;
3544: PetscInt i, mbs = a->mbs, nz = a->nz, bs2 = a->bs2;
3546: PetscFunctionBegin;
3547: PetscCheck(A->assembled, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Cannot duplicate unassembled matrix");
3548: PetscCheck(a->i[mbs] == nz, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Corrupt matrix");
3549: PetscCall(MatSetOption(C, MAT_STRUCTURE_ONLY, A->structure_only));
3551: if (cpvalues == MAT_SHARE_NONZERO_PATTERN) {
3552: c->imax = a->imax;
3553: c->ilen = a->ilen;
3554: c->free_imax_ilen = PETSC_FALSE;
3555: } else {
3556: PetscCall(PetscMalloc2(mbs, &c->imax, mbs, &c->ilen));
3557: for (i = 0; i < mbs; i++) {
3558: c->imax[i] = a->imax[i];
3559: c->ilen[i] = a->ilen[i];
3560: }
3561: c->free_imax_ilen = PETSC_TRUE;
3562: }
3564: /* allocate the matrix space */
3565: if (mallocmatspace) {
3566: if (cpvalues == MAT_SHARE_NONZERO_PATTERN) {
3567: if (!A->structure_only) {
3568: PetscCall(PetscShmgetAllocateArray(bs2 * nz, sizeof(PetscScalar), (void **)&c->a));
3569: PetscCall(PetscArrayzero(c->a, bs2 * nz));
3570: }
3571: c->free_a = PETSC_TRUE;
3572: c->i = a->i;
3573: c->j = a->j;
3574: c->free_ij = PETSC_FALSE;
3575: c->parent = A;
3576: C->preallocated = PETSC_TRUE;
3577: C->assembled = PETSC_TRUE;
3579: PetscCall(PetscObjectReference((PetscObject)A));
3580: PetscCall(MatSetOption(A, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
3581: PetscCall(MatSetOption(C, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
3582: } else {
3583: if (!A->structure_only) PetscCall(PetscShmgetAllocateArray(bs2 * nz, sizeof(PetscScalar), (void **)&c->a));
3584: PetscCall(PetscShmgetAllocateArray(nz, sizeof(PetscInt), (void **)&c->j));
3585: PetscCall(PetscShmgetAllocateArray(mbs + 1, sizeof(PetscInt), (void **)&c->i));
3586: c->free_a = PETSC_TRUE;
3587: c->free_ij = PETSC_TRUE;
3589: PetscCall(PetscArraycpy(c->i, a->i, mbs + 1));
3590: if (mbs > 0) {
3591: PetscCall(PetscArraycpy(c->j, a->j, nz));
3592: if (!A->structure_only) {
3593: if (cpvalues == MAT_COPY_VALUES) PetscCall(PetscArraycpy(c->a, a->a, bs2 * nz));
3594: else PetscCall(PetscArrayzero(c->a, bs2 * nz));
3595: }
3596: }
3597: C->preallocated = PETSC_TRUE;
3598: C->assembled = PETSC_TRUE;
3599: }
3600: }
3602: c->roworiented = a->roworiented;
3603: c->nonew = a->nonew;
3605: PetscCall(PetscLayoutReference(A->rmap, &C->rmap));
3606: PetscCall(PetscLayoutReference(A->cmap, &C->cmap));
3608: c->bs2 = a->bs2;
3609: c->mbs = a->mbs;
3610: c->nbs = a->nbs;
3611: c->nz = a->nz;
3612: c->maxnz = a->nz; /* Since we allocate exactly the right amount */
3613: c->solve_work = NULL;
3614: c->mult_work = NULL;
3615: c->sor_workt = NULL;
3616: c->sor_work = NULL;
3618: c->compressedrow.use = a->compressedrow.use;
3619: c->compressedrow.nrows = a->compressedrow.nrows;
3620: if (a->compressedrow.use) {
3621: i = a->compressedrow.nrows;
3622: PetscCall(PetscMalloc2(i + 1, &c->compressedrow.i, i + 1, &c->compressedrow.rindex));
3623: PetscCall(PetscArraycpy(c->compressedrow.i, a->compressedrow.i, i + 1));
3624: PetscCall(PetscArraycpy(c->compressedrow.rindex, a->compressedrow.rindex, i));
3625: } else {
3626: c->compressedrow.use = PETSC_FALSE;
3627: c->compressedrow.i = NULL;
3628: c->compressedrow.rindex = NULL;
3629: }
3630: c->nonzerorowcnt = a->nonzerorowcnt;
3631: C->nonzerostate = A->nonzerostate;
3633: PetscCall(PetscFunctionListDuplicate(((PetscObject)A)->qlist, &((PetscObject)C)->qlist));
3634: PetscFunctionReturn(PETSC_SUCCESS);
3635: }
3637: PetscErrorCode MatDuplicate_SeqBAIJ(Mat A, MatDuplicateOption cpvalues, Mat *B)
3638: {
3639: PetscFunctionBegin;
3640: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B));
3641: PetscCall(MatSetSizes(*B, A->rmap->N, A->cmap->n, A->rmap->N, A->cmap->n));
3642: PetscCall(MatSetType(*B, MATSEQBAIJ));
3643: PetscCall(MatDuplicateNoCreate_SeqBAIJ(*B, A, cpvalues, PETSC_TRUE));
3644: PetscFunctionReturn(PETSC_SUCCESS);
3645: }
3647: /* Used for both SeqBAIJ and SeqSBAIJ matrices */
3648: PetscErrorCode MatLoad_SeqBAIJ_Binary(Mat mat, PetscViewer viewer)
3649: {
3650: PetscInt header[4], M, N, nz, bs, m, n, mbs, nbs, rows, cols, sum, i, j, k;
3651: PetscInt *rowidxs, *colidxs;
3652: PetscScalar *matvals;
3654: PetscFunctionBegin;
3655: PetscCall(PetscViewerSetUp(viewer));
3657: /* read matrix header */
3658: PetscCall(PetscViewerBinaryRead(viewer, header, 4, NULL, PETSC_INT));
3659: PetscCheck(header[0] == MAT_FILE_CLASSID, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Not a matrix object in file");
3660: M = header[1];
3661: N = header[2];
3662: nz = header[3];
3663: PetscCheck(M >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix row size (%" PetscInt_FMT ") in file is negative", M);
3664: PetscCheck(N >= 0, PetscObjectComm((PetscObject)viewer), PETSC_ERR_FILE_UNEXPECTED, "Matrix column size (%" PetscInt_FMT ") in file is negative", N);
3665: PetscCheck(nz >= 0, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Matrix stored in special format on disk, cannot load as SeqBAIJ");
3667: /* set block sizes from the viewer's .info file */
3668: PetscCall(MatLoad_Binary_BlockSizes(mat, viewer));
3669: /* set local and global sizes if not set already */
3670: if (mat->rmap->n < 0) mat->rmap->n = M;
3671: if (mat->cmap->n < 0) mat->cmap->n = N;
3672: if (mat->rmap->N < 0) mat->rmap->N = M;
3673: if (mat->cmap->N < 0) mat->cmap->N = N;
3674: PetscCall(PetscLayoutSetUp(mat->rmap));
3675: PetscCall(PetscLayoutSetUp(mat->cmap));
3677: /* check if the matrix sizes are correct */
3678: PetscCall(MatGetSize(mat, &rows, &cols));
3679: PetscCheck(M == rows && N == cols, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different sizes (%" PetscInt_FMT ", %" PetscInt_FMT ") than the input matrix (%" PetscInt_FMT ", %" PetscInt_FMT ")", M, N, rows, cols);
3680: PetscCall(MatGetBlockSize(mat, &bs));
3681: PetscCall(MatGetLocalSize(mat, &m, &n));
3682: mbs = m / bs;
3683: nbs = n / bs;
3685: /* read in row lengths, column indices and nonzero values */
3686: PetscCall(PetscMalloc1(m + 1, &rowidxs));
3687: PetscCall(PetscViewerBinaryRead(viewer, rowidxs + 1, m, NULL, PETSC_INT));
3688: rowidxs[0] = 0;
3689: for (i = 0; i < m; i++) rowidxs[i + 1] += rowidxs[i];
3690: sum = rowidxs[m];
3691: PetscCheck(sum == nz, PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Inconsistent matrix data in file: nonzeros = %" PetscInt_FMT ", sum-row-lengths = %" PetscInt_FMT, nz, sum);
3693: /* read in column indices and nonzero values */
3694: PetscCall(PetscMalloc2(rowidxs[m], &colidxs, nz, &matvals));
3695: PetscCall(PetscViewerBinaryRead(viewer, colidxs, rowidxs[m], NULL, PETSC_INT));
3696: PetscCall(PetscViewerBinaryRead(viewer, matvals, rowidxs[m], NULL, PETSC_SCALAR));
3698: { /* preallocate matrix storage */
3699: PetscBT bt; /* helper bit set to count nonzeros */
3700: PetscInt *nnz;
3701: PetscBool sbaij;
3703: PetscCall(PetscBTCreate(nbs, &bt));
3704: PetscCall(PetscCalloc1(mbs, &nnz));
3705: PetscCall(PetscObjectTypeCompare((PetscObject)mat, MATSEQSBAIJ, &sbaij));
3706: for (i = 0; i < mbs; i++) {
3707: PetscCall(PetscBTMemzero(nbs, bt));
3708: for (k = 0; k < bs; k++) {
3709: PetscInt row = bs * i + k;
3710: for (j = rowidxs[row]; j < rowidxs[row + 1]; j++) {
3711: PetscInt col = colidxs[j];
3712: if (!sbaij || col >= row)
3713: if (!PetscBTLookupSet(bt, col / bs)) nnz[i]++;
3714: }
3715: }
3716: }
3717: PetscCall(PetscBTDestroy(&bt));
3718: PetscCall(MatSeqBAIJSetPreallocation(mat, bs, 0, nnz));
3719: PetscCall(MatSeqSBAIJSetPreallocation(mat, bs, 0, nnz));
3720: PetscCall(PetscFree(nnz));
3721: }
3723: /* store matrix values */
3724: for (i = 0; i < m; i++) {
3725: PetscInt row = i, s = rowidxs[i], e = rowidxs[i + 1];
3726: PetscUseTypeMethod(mat, setvalues, 1, &row, e - s, colidxs + s, matvals + s, INSERT_VALUES);
3727: }
3729: PetscCall(PetscFree(rowidxs));
3730: PetscCall(PetscFree2(colidxs, matvals));
3731: PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
3732: PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY));
3733: PetscFunctionReturn(PETSC_SUCCESS);
3734: }
3736: PetscErrorCode MatLoad_SeqBAIJ(Mat mat, PetscViewer viewer)
3737: {
3738: PetscBool isbinary;
3740: PetscFunctionBegin;
3741: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
3742: PetscCheck(isbinary, PetscObjectComm((PetscObject)viewer), PETSC_ERR_SUP, "Viewer type %s not yet supported for reading %s matrices", ((PetscObject)viewer)->type_name, ((PetscObject)mat)->type_name);
3743: PetscCall(MatLoad_SeqBAIJ_Binary(mat, viewer));
3744: PetscFunctionReturn(PETSC_SUCCESS);
3745: }
3747: /*@
3748: MatCreateSeqBAIJ - Creates a sparse matrix in `MATSEQAIJ` (block
3749: compressed row) format. For good matrix assembly performance the
3750: user should preallocate the matrix storage by setting the parameter `nz`
3751: (or the array `nnz`).
3753: Collective
3755: Input Parameters:
3756: + comm - MPI communicator, set to `PETSC_COMM_SELF`
3757: . bs - size of block, the blocks are ALWAYS square. One can use `MatSetBlockSizes()` to set a different row and column blocksize but the row
3758: blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with `MatCreateVecs()`
3759: . m - number of rows
3760: . n - number of columns
3761: . nz - number of nonzero blocks per block row (same for all rows)
3762: - nnz - array containing the number of nonzero blocks in the various block rows
3763: (possibly different for each block row) or `NULL`
3765: Output Parameter:
3766: . A - the matrix
3768: Options Database Keys:
3769: + -mat_no_unroll - uses code that does not unroll the loops in the block calculations (much slower)
3770: - -mat_block_size - size of the blocks to use
3772: Level: intermediate
3774: Notes:
3775: It is recommended that one use `MatCreateFromOptions()` or the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`,
3776: MatXXXXSetPreallocation() paradigm instead of this routine directly.
3777: [MatXXXXSetPreallocation() is, for example, `MatSeqAIJSetPreallocation()`]
3779: The number of rows and columns must be divisible by blocksize.
3781: If the `nnz` parameter is given then the `nz` parameter is ignored
3783: A nonzero block is any block that as 1 or more nonzeros in it
3785: The `MATSEQBAIJ` format is fully compatible with standard Fortran
3786: storage. That is, the stored row and column indices can begin at
3787: either one (as in Fortran) or zero.
3789: Specify the preallocated storage with either `nz` or `nnz` (not both).
3790: Set `nz` = `PETSC_DEFAULT` and `nnz` = `NULL` for PETSc to control dynamic memory
3791: allocation. See [Sparse Matrices](sec_matsparse) for details.
3792: matrices.
3794: .seealso: [](ch_matrices), `Mat`, [Sparse Matrices](sec_matsparse), `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatCreateBAIJ()`
3795: @*/
3796: PetscErrorCode MatCreateSeqBAIJ(MPI_Comm comm, PetscInt bs, PetscInt m, PetscInt n, PetscInt nz, const PetscInt nnz[], Mat *A)
3797: {
3798: PetscFunctionBegin;
3799: PetscCall(MatCreate(comm, A));
3800: PetscCall(MatSetSizes(*A, m, n, m, n));
3801: PetscCall(MatSetType(*A, MATSEQBAIJ));
3802: PetscCall(MatSeqBAIJSetPreallocation(*A, bs, nz, (PetscInt *)nnz));
3803: PetscFunctionReturn(PETSC_SUCCESS);
3804: }
3806: /*@
3807: MatSeqBAIJSetPreallocation - Sets the block size and expected nonzeros
3808: per row in the matrix. For good matrix assembly performance the
3809: user should preallocate the matrix storage by setting the parameter `nz`
3810: (or the array `nnz`).
3812: Collective
3814: Input Parameters:
3815: + B - the matrix
3816: . bs - size of block, the blocks are ALWAYS square. One can use `MatSetBlockSizes()` to set a different row and column blocksize but the row
3817: blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with `MatCreateVecs()`
3818: . nz - number of block nonzeros per block row (same for all rows)
3819: - nnz - array containing the number of block nonzeros in the various block rows
3820: (possibly different for each block row) or `NULL`
3822: Options Database Keys:
3823: + -mat_no_unroll - uses code that does not unroll the loops in the block calculations (much slower)
3824: - -mat_block_size - size of the blocks to use
3826: Level: intermediate
3828: Notes:
3829: If the `nnz` parameter is given then the `nz` parameter is ignored
3831: You can call `MatGetInfo()` to get information on how effective the preallocation was;
3832: for example the fields mallocs,nz_allocated,nz_used,nz_unneeded;
3833: You can also run with the option `-info` and look for messages with the string
3834: malloc in them to see if additional memory allocation was needed.
3836: The `MATSEQBAIJ` format is fully compatible with standard Fortran
3837: storage. That is, the stored row and column indices can begin at
3838: either one (as in Fortran) or zero.
3840: Specify the preallocated storage with either `nz` or `nnz` (not both).
3841: Set `nz` = `PETSC_DEFAULT` and `nnz` = `NULL` for PETSc to control dynamic memory
3842: allocation. See [Sparse Matrices](sec_matsparse) for details.
3844: .seealso: [](ch_matrices), `Mat`, [Sparse Matrices](sec_matsparse), `MatCreate()`, `MatCreateSeqAIJ()`, `MatSetValues()`, `MatCreateBAIJ()`, `MatGetInfo()`
3845: @*/
3846: PetscErrorCode MatSeqBAIJSetPreallocation(Mat B, PetscInt bs, PetscInt nz, const PetscInt nnz[])
3847: {
3848: PetscFunctionBegin;
3852: PetscTryMethod(B, "MatSeqBAIJSetPreallocation_C", (Mat, PetscInt, PetscInt, const PetscInt[]), (B, bs, nz, nnz));
3853: PetscFunctionReturn(PETSC_SUCCESS);
3854: }
3856: /*@
3857: MatSeqBAIJSetPreallocationCSR - Creates a sparse sequential matrix in `MATSEQBAIJ` format using the given nonzero structure and (optional) numerical values
3859: Collective
3861: Input Parameters:
3862: + B - the matrix
3863: . bs - the blocksize
3864: . i - the indices into `j` for the start of each local row (indices start with zero)
3865: . j - the column indices for each local row (indices start with zero) these must be sorted for each row
3866: - v - optional values in the matrix, use `NULL` if not provided
3868: Level: advanced
3870: Notes:
3871: The `i`,`j`,`v` values are COPIED with this routine; to avoid the copy use `MatCreateSeqBAIJWithArrays()`
3873: The order of the entries in values is specified by the `MatOption` `MAT_ROW_ORIENTED`. For example, C programs
3874: may want to use the default `MAT_ROW_ORIENTED` of `PETSC_TRUE` and use an array v[nnz][bs][bs] where the second index is
3875: over rows within a block and the last index is over columns within a block row. Fortran programs will likely set
3876: `MAT_ROW_ORIENTED` of `PETSC_FALSE` and use a Fortran array v(bs,bs,nnz) in which the first index is over rows within a
3877: block column and the second index is over columns within a block.
3879: Though this routine has Preallocation() in the name it also sets the exact nonzero locations of the matrix entries and usually the numerical values as well
3881: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatCreateSeqBAIJ()`, `MatSetValues()`, `MatSeqBAIJSetPreallocation()`, `MATSEQBAIJ`
3882: @*/
3883: PetscErrorCode MatSeqBAIJSetPreallocationCSR(Mat B, PetscInt bs, const PetscInt i[], const PetscInt j[], const PetscScalar v[])
3884: {
3885: PetscFunctionBegin;
3889: PetscTryMethod(B, "MatSeqBAIJSetPreallocationCSR_C", (Mat, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[]), (B, bs, i, j, v));
3890: PetscFunctionReturn(PETSC_SUCCESS);
3891: }
3893: /*@
3894: MatCreateSeqBAIJWithArrays - Creates a `MATSEQBAIJ` matrix using matrix elements provided by the user.
3896: Collective
3898: Input Parameters:
3899: + comm - must be an MPI communicator of size 1
3900: . bs - size of block
3901: . m - number of rows
3902: . n - number of columns
3903: . i - row indices; that is i[0] = 0, i[row] = i[row-1] + number of elements in that row block row of the matrix
3904: . j - column indices
3905: - a - matrix values
3907: Output Parameter:
3908: . mat - the matrix
3910: Level: advanced
3912: Notes:
3913: The `i`, `j`, and `a` arrays are not copied by this routine, the user must free these arrays
3914: once the matrix is destroyed
3916: You cannot set new nonzero locations into this matrix, that will generate an error.
3918: The `i` and `j` indices are 0 based
3920: When block size is greater than 1 the matrix values must be stored using the `MATSEQBAIJ` storage format
3922: The order of the entries in values is the same as the block compressed sparse row storage format; that is, it is
3923: the same as a three dimensional array in Fortran values(bs,bs,nnz) that contains the first column of the first
3924: block, followed by the second column of the first block etc etc. That is, the blocks are contiguous in memory
3925: with column-major ordering within blocks.
3927: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatCreateBAIJ()`, `MatCreateSeqBAIJ()`
3928: @*/
3929: PetscErrorCode MatCreateSeqBAIJWithArrays(MPI_Comm comm, PetscInt bs, PetscInt m, PetscInt n, PetscInt i[], PetscInt j[], PetscScalar a[], Mat *mat)
3930: {
3931: Mat_SeqBAIJ *baij;
3933: PetscFunctionBegin;
3934: PetscCheck(bs == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "block size %" PetscInt_FMT " > 1 is not supported yet", bs);
3935: if (m > 0) PetscCheck(i[0] == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "i (row indices) must start with 0");
3937: PetscCall(MatCreate(comm, mat));
3938: PetscCall(MatSetSizes(*mat, m, n, m, n));
3939: PetscCall(MatSetType(*mat, MATSEQBAIJ));
3940: PetscCall(MatSeqBAIJSetPreallocation(*mat, bs, MAT_SKIP_ALLOCATION, NULL));
3941: baij = (Mat_SeqBAIJ *)(*mat)->data;
3942: PetscCall(PetscMalloc2(m, &baij->imax, m, &baij->ilen));
3944: baij->i = i;
3945: baij->j = j;
3946: baij->a = a;
3948: baij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/
3949: baij->free_a = PETSC_FALSE;
3950: baij->free_ij = PETSC_FALSE;
3951: baij->free_imax_ilen = PETSC_TRUE;
3953: for (PetscInt ii = 0; ii < m; ii++) {
3954: const PetscInt row_len = i[ii + 1] - i[ii];
3956: baij->ilen[ii] = baij->imax[ii] = row_len;
3957: PetscCheck(row_len >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative row length in i (row indices) row = %" PetscInt_FMT " length = %" PetscInt_FMT, ii, row_len);
3958: }
3959: if (PetscDefined(USE_DEBUG)) {
3960: for (PetscInt ii = 0; ii < baij->i[m]; ii++) {
3961: PetscCheck(j[ii] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative column index at location = %" PetscInt_FMT " index = %" PetscInt_FMT, ii, j[ii]);
3962: PetscCheck(j[ii] <= n - 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column index to large at location = %" PetscInt_FMT " index = %" PetscInt_FMT, ii, j[ii]);
3963: }
3964: }
3966: PetscCall(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY));
3967: PetscCall(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY));
3968: PetscFunctionReturn(PETSC_SUCCESS);
3969: }
3971: PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqBAIJ(MPI_Comm comm, Mat inmat, PetscInt n, MatReuse scall, Mat *outmat)
3972: {
3973: PetscFunctionBegin;
3974: PetscCall(MatCreateMPIMatConcatenateSeqMat_MPIBAIJ(comm, inmat, n, scall, outmat));
3975: PetscFunctionReturn(PETSC_SUCCESS);
3976: }