Actual source code: baijfact9.c
1: /*
2: Factorization code for BAIJ format.
3: */
4: #include <../src/mat/impls/baij/seq/baij.h>
5: #include <petsc/private/kernels/blockinvert.h>
7: /*
8: Version for when blocks are 5 by 5
9: */
10: PetscErrorCode MatILUFactorNumeric_SeqBAIJ_5_inplace(Mat C, Mat A, const MatFactorInfo *info)
11: {
12: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)C->data;
13: IS isrow = b->row, isicol = b->icol;
14: const PetscInt *r, *ic;
15: PetscInt *bi = b->i, *bj = b->j, *ajtmpold, *ajtmp;
16: PetscInt i, j, n = a->mbs, nz, row, idx, ipvt[5];
17: const PetscInt *diag_offset;
18: PetscInt *ai = a->i, *aj = a->j, *pj;
19: MatScalar *w, *pv, *rtmp, *x, *pc;
20: const MatScalar *v, *aa = a->a;
21: MatScalar p1, p2, p3, p4, m1, m2, m3, m4, m5, m6, m7, m8, m9, x1, x2, x3, x4;
22: MatScalar p5, p6, p7, p8, p9, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15, x16;
23: MatScalar x17, x18, x19, x20, x21, x22, x23, x24, x25, p10, p11, p12, p13, p14;
24: MatScalar p15, p16, p17, p18, p19, p20, p21, p22, p23, p24, p25, m10, m11, m12;
25: MatScalar m13, m14, m15, m16, m17, m18, m19, m20, m21, m22, m23, m24, m25;
26: MatScalar *ba = b->a, work[25];
27: PetscReal shift = info->shiftamount;
28: PetscBool allowzeropivot, zeropivotdetected;
30: PetscFunctionBegin;
31: /* Since A is C and C is labeled as a factored matrix we need to lie to MatGetDiagonalMarkers_SeqBAIJ() to get it to compute the diagonals */
32: A->factortype = MAT_FACTOR_NONE;
33: PetscCall(MatGetDiagonalMarkers_SeqBAIJ(A, &diag_offset, NULL));
34: A->factortype = MAT_FACTOR_ILU;
35: allowzeropivot = PetscNot(A->erroriffailure);
36: PetscCall(ISGetIndices(isrow, &r));
37: PetscCall(ISGetIndices(isicol, &ic));
38: PetscCall(PetscMalloc1(25 * (n + 1), &rtmp));
40: #define PETSC_USE_MEMZERO 1
41: #define PETSC_USE_MEMCPY 1
43: for (i = 0; i < n; i++) {
44: nz = bi[i + 1] - bi[i];
45: ajtmp = bj + bi[i];
46: for (j = 0; j < nz; j++) {
47: if (PetscDefined(USE_MEMZERO)) PetscCall(PetscArrayzero(rtmp + 25 * ajtmp[j], 25));
48: else {
49: x = rtmp + 25 * ajtmp[j];
50: x[0] = x[1] = x[2] = x[3] = x[4] = x[5] = x[6] = x[7] = x[8] = x[9] = 0.0;
51: x[10] = x[11] = x[12] = x[13] = x[14] = x[15] = x[16] = x[17] = 0.0;
52: x[18] = x[19] = x[20] = x[21] = x[22] = x[23] = x[24] = 0.0;
53: }
54: }
55: /* load in initial (unfactored row) */
56: idx = r[i];
57: nz = ai[idx + 1] - ai[idx];
58: ajtmpold = aj + ai[idx];
59: v = aa + 25 * ai[idx];
60: for (j = 0; j < nz; j++) {
61: if (PetscDefined(USE_MEMCPY)) PetscCall(PetscArraycpy(rtmp + 25 * ic[ajtmpold[j]], v, 25));
62: else {
63: x = rtmp + 25 * ic[ajtmpold[j]];
64: x[0] = v[0];
65: x[1] = v[1];
66: x[2] = v[2];
67: x[3] = v[3];
68: x[4] = v[4];
69: x[5] = v[5];
70: x[6] = v[6];
71: x[7] = v[7];
72: x[8] = v[8];
73: x[9] = v[9];
74: x[10] = v[10];
75: x[11] = v[11];
76: x[12] = v[12];
77: x[13] = v[13];
78: x[14] = v[14];
79: x[15] = v[15];
80: x[16] = v[16];
81: x[17] = v[17];
82: x[18] = v[18];
83: x[19] = v[19];
84: x[20] = v[20];
85: x[21] = v[21];
86: x[22] = v[22];
87: x[23] = v[23];
88: x[24] = v[24];
89: }
90: v += 25;
91: }
92: row = *ajtmp++;
93: while (row < i) {
94: pc = rtmp + 25 * row;
95: p1 = pc[0];
96: p2 = pc[1];
97: p3 = pc[2];
98: p4 = pc[3];
99: p5 = pc[4];
100: p6 = pc[5];
101: p7 = pc[6];
102: p8 = pc[7];
103: p9 = pc[8];
104: p10 = pc[9];
105: p11 = pc[10];
106: p12 = pc[11];
107: p13 = pc[12];
108: p14 = pc[13];
109: p15 = pc[14];
110: p16 = pc[15];
111: p17 = pc[16];
112: p18 = pc[17];
113: p19 = pc[18];
114: p20 = pc[19];
115: p21 = pc[20];
116: p22 = pc[21];
117: p23 = pc[22];
118: p24 = pc[23];
119: p25 = pc[24];
120: if (p1 != 0.0 || p2 != 0.0 || p3 != 0.0 || p4 != 0.0 || p5 != 0.0 || p6 != 0.0 || p7 != 0.0 || p8 != 0.0 || p9 != 0.0 || p10 != 0.0 || p11 != 0.0 || p12 != 0.0 || p13 != 0.0 || p14 != 0.0 || p15 != 0.0 || p16 != 0.0 || p17 != 0.0 || p18 != 0.0 || p19 != 0.0 || p20 != 0.0 || p21 != 0.0 || p22 != 0.0 || p23 != 0.0 || p24 != 0.0 || p25 != 0.0) {
121: pv = ba + 25 * diag_offset[row];
122: pj = bj + diag_offset[row] + 1;
123: x1 = pv[0];
124: x2 = pv[1];
125: x3 = pv[2];
126: x4 = pv[3];
127: x5 = pv[4];
128: x6 = pv[5];
129: x7 = pv[6];
130: x8 = pv[7];
131: x9 = pv[8];
132: x10 = pv[9];
133: x11 = pv[10];
134: x12 = pv[11];
135: x13 = pv[12];
136: x14 = pv[13];
137: x15 = pv[14];
138: x16 = pv[15];
139: x17 = pv[16];
140: x18 = pv[17];
141: x19 = pv[18];
142: x20 = pv[19];
143: x21 = pv[20];
144: x22 = pv[21];
145: x23 = pv[22];
146: x24 = pv[23];
147: x25 = pv[24];
148: pc[0] = m1 = p1 * x1 + p6 * x2 + p11 * x3 + p16 * x4 + p21 * x5;
149: pc[1] = m2 = p2 * x1 + p7 * x2 + p12 * x3 + p17 * x4 + p22 * x5;
150: pc[2] = m3 = p3 * x1 + p8 * x2 + p13 * x3 + p18 * x4 + p23 * x5;
151: pc[3] = m4 = p4 * x1 + p9 * x2 + p14 * x3 + p19 * x4 + p24 * x5;
152: pc[4] = m5 = p5 * x1 + p10 * x2 + p15 * x3 + p20 * x4 + p25 * x5;
154: pc[5] = m6 = p1 * x6 + p6 * x7 + p11 * x8 + p16 * x9 + p21 * x10;
155: pc[6] = m7 = p2 * x6 + p7 * x7 + p12 * x8 + p17 * x9 + p22 * x10;
156: pc[7] = m8 = p3 * x6 + p8 * x7 + p13 * x8 + p18 * x9 + p23 * x10;
157: pc[8] = m9 = p4 * x6 + p9 * x7 + p14 * x8 + p19 * x9 + p24 * x10;
158: pc[9] = m10 = p5 * x6 + p10 * x7 + p15 * x8 + p20 * x9 + p25 * x10;
160: pc[10] = m11 = p1 * x11 + p6 * x12 + p11 * x13 + p16 * x14 + p21 * x15;
161: pc[11] = m12 = p2 * x11 + p7 * x12 + p12 * x13 + p17 * x14 + p22 * x15;
162: pc[12] = m13 = p3 * x11 + p8 * x12 + p13 * x13 + p18 * x14 + p23 * x15;
163: pc[13] = m14 = p4 * x11 + p9 * x12 + p14 * x13 + p19 * x14 + p24 * x15;
164: pc[14] = m15 = p5 * x11 + p10 * x12 + p15 * x13 + p20 * x14 + p25 * x15;
166: pc[15] = m16 = p1 * x16 + p6 * x17 + p11 * x18 + p16 * x19 + p21 * x20;
167: pc[16] = m17 = p2 * x16 + p7 * x17 + p12 * x18 + p17 * x19 + p22 * x20;
168: pc[17] = m18 = p3 * x16 + p8 * x17 + p13 * x18 + p18 * x19 + p23 * x20;
169: pc[18] = m19 = p4 * x16 + p9 * x17 + p14 * x18 + p19 * x19 + p24 * x20;
170: pc[19] = m20 = p5 * x16 + p10 * x17 + p15 * x18 + p20 * x19 + p25 * x20;
172: pc[20] = m21 = p1 * x21 + p6 * x22 + p11 * x23 + p16 * x24 + p21 * x25;
173: pc[21] = m22 = p2 * x21 + p7 * x22 + p12 * x23 + p17 * x24 + p22 * x25;
174: pc[22] = m23 = p3 * x21 + p8 * x22 + p13 * x23 + p18 * x24 + p23 * x25;
175: pc[23] = m24 = p4 * x21 + p9 * x22 + p14 * x23 + p19 * x24 + p24 * x25;
176: pc[24] = m25 = p5 * x21 + p10 * x22 + p15 * x23 + p20 * x24 + p25 * x25;
178: nz = bi[row + 1] - diag_offset[row] - 1;
179: pv += 25;
180: for (j = 0; j < nz; j++) {
181: x1 = pv[0];
182: x2 = pv[1];
183: x3 = pv[2];
184: x4 = pv[3];
185: x5 = pv[4];
186: x6 = pv[5];
187: x7 = pv[6];
188: x8 = pv[7];
189: x9 = pv[8];
190: x10 = pv[9];
191: x11 = pv[10];
192: x12 = pv[11];
193: x13 = pv[12];
194: x14 = pv[13];
195: x15 = pv[14];
196: x16 = pv[15];
197: x17 = pv[16];
198: x18 = pv[17];
199: x19 = pv[18];
200: x20 = pv[19];
201: x21 = pv[20];
202: x22 = pv[21];
203: x23 = pv[22];
204: x24 = pv[23];
205: x25 = pv[24];
206: x = rtmp + 25 * pj[j];
207: x[0] -= m1 * x1 + m6 * x2 + m11 * x3 + m16 * x4 + m21 * x5;
208: x[1] -= m2 * x1 + m7 * x2 + m12 * x3 + m17 * x4 + m22 * x5;
209: x[2] -= m3 * x1 + m8 * x2 + m13 * x3 + m18 * x4 + m23 * x5;
210: x[3] -= m4 * x1 + m9 * x2 + m14 * x3 + m19 * x4 + m24 * x5;
211: x[4] -= m5 * x1 + m10 * x2 + m15 * x3 + m20 * x4 + m25 * x5;
213: x[5] -= m1 * x6 + m6 * x7 + m11 * x8 + m16 * x9 + m21 * x10;
214: x[6] -= m2 * x6 + m7 * x7 + m12 * x8 + m17 * x9 + m22 * x10;
215: x[7] -= m3 * x6 + m8 * x7 + m13 * x8 + m18 * x9 + m23 * x10;
216: x[8] -= m4 * x6 + m9 * x7 + m14 * x8 + m19 * x9 + m24 * x10;
217: x[9] -= m5 * x6 + m10 * x7 + m15 * x8 + m20 * x9 + m25 * x10;
219: x[10] -= m1 * x11 + m6 * x12 + m11 * x13 + m16 * x14 + m21 * x15;
220: x[11] -= m2 * x11 + m7 * x12 + m12 * x13 + m17 * x14 + m22 * x15;
221: x[12] -= m3 * x11 + m8 * x12 + m13 * x13 + m18 * x14 + m23 * x15;
222: x[13] -= m4 * x11 + m9 * x12 + m14 * x13 + m19 * x14 + m24 * x15;
223: x[14] -= m5 * x11 + m10 * x12 + m15 * x13 + m20 * x14 + m25 * x15;
225: x[15] -= m1 * x16 + m6 * x17 + m11 * x18 + m16 * x19 + m21 * x20;
226: x[16] -= m2 * x16 + m7 * x17 + m12 * x18 + m17 * x19 + m22 * x20;
227: x[17] -= m3 * x16 + m8 * x17 + m13 * x18 + m18 * x19 + m23 * x20;
228: x[18] -= m4 * x16 + m9 * x17 + m14 * x18 + m19 * x19 + m24 * x20;
229: x[19] -= m5 * x16 + m10 * x17 + m15 * x18 + m20 * x19 + m25 * x20;
231: x[20] -= m1 * x21 + m6 * x22 + m11 * x23 + m16 * x24 + m21 * x25;
232: x[21] -= m2 * x21 + m7 * x22 + m12 * x23 + m17 * x24 + m22 * x25;
233: x[22] -= m3 * x21 + m8 * x22 + m13 * x23 + m18 * x24 + m23 * x25;
234: x[23] -= m4 * x21 + m9 * x22 + m14 * x23 + m19 * x24 + m24 * x25;
235: x[24] -= m5 * x21 + m10 * x22 + m15 * x23 + m20 * x24 + m25 * x25;
237: pv += 25;
238: }
239: PetscCall(PetscLogFlops(250.0 * nz + 225.0));
240: }
241: row = *ajtmp++;
242: }
243: /* finished row so stick it into b->a */
244: pv = ba + 25 * bi[i];
245: pj = bj + bi[i];
246: nz = bi[i + 1] - bi[i];
247: for (j = 0; j < nz; j++) {
248: if (PetscDefined(USE_MEMCPY)) PetscCall(PetscArraycpy(pv, rtmp + 25 * pj[j], 25));
249: else {
250: x = rtmp + 25 * pj[j];
251: pv[0] = x[0];
252: pv[1] = x[1];
253: pv[2] = x[2];
254: pv[3] = x[3];
255: pv[4] = x[4];
256: pv[5] = x[5];
257: pv[6] = x[6];
258: pv[7] = x[7];
259: pv[8] = x[8];
260: pv[9] = x[9];
261: pv[10] = x[10];
262: pv[11] = x[11];
263: pv[12] = x[12];
264: pv[13] = x[13];
265: pv[14] = x[14];
266: pv[15] = x[15];
267: pv[16] = x[16];
268: pv[17] = x[17];
269: pv[18] = x[18];
270: pv[19] = x[19];
271: pv[20] = x[20];
272: pv[21] = x[21];
273: pv[22] = x[22];
274: pv[23] = x[23];
275: pv[24] = x[24];
276: }
277: pv += 25;
278: }
279: /* invert diagonal block */
280: w = ba + 25 * diag_offset[i];
281: PetscCall(PetscKernel_A_gets_inverse_A_5(w, ipvt, work, shift, allowzeropivot, &zeropivotdetected));
282: if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
283: }
285: PetscCall(PetscFree(rtmp));
286: PetscCall(ISRestoreIndices(isicol, &ic));
287: PetscCall(ISRestoreIndices(isrow, &r));
289: C->ops->solve = MatSolve_SeqBAIJ_5_inplace;
290: C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_5_inplace;
291: C->assembled = PETSC_TRUE;
293: PetscCall(PetscLogFlops(1.333333333333 * 5 * 5 * 5 * b->mbs)); /* from inverting diagonal blocks */
294: PetscFunctionReturn(PETSC_SUCCESS);
295: }
297: /* MatLUFactorNumeric_SeqBAIJ_5 -
298: copied from MatLUFactorNumeric_SeqBAIJ_N_inplace() and manually re-implemented
299: PetscKernel_A_gets_A_times_B()
300: PetscKernel_A_gets_A_minus_B_times_C()
301: PetscKernel_A_gets_inverse_A()
302: */
304: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_5(Mat B, Mat A, const MatFactorInfo *info)
305: {
306: Mat C = B;
307: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)C->data;
308: IS isrow = b->row, isicol = b->icol;
309: const PetscInt *r, *ic;
310: PetscInt i, j, k, nz, nzL, row;
311: const PetscInt n = a->mbs, *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j;
312: const PetscInt *ajtmp, *bjtmp, *bdiag = b->diag, *pj, bs2 = a->bs2;
313: MatScalar *rtmp, *pc, *mwork, *v, *pv, *aa = a->a, work[25];
314: PetscInt flg, ipvt[5];
315: PetscReal shift = info->shiftamount;
316: PetscBool allowzeropivot, zeropivotdetected;
318: PetscFunctionBegin;
319: allowzeropivot = PetscNot(A->erroriffailure);
320: PetscCall(ISGetIndices(isrow, &r));
321: PetscCall(ISGetIndices(isicol, &ic));
323: /* generate work space needed by the factorization */
324: PetscCall(PetscMalloc2(bs2 * n, &rtmp, bs2, &mwork));
325: PetscCall(PetscArrayzero(rtmp, bs2 * n));
327: for (i = 0; i < n; i++) {
328: /* zero rtmp */
329: /* L part */
330: nz = bi[i + 1] - bi[i];
331: bjtmp = bj + bi[i];
332: for (j = 0; j < nz; j++) PetscCall(PetscArrayzero(rtmp + bs2 * bjtmp[j], bs2));
334: /* U part */
335: nz = bdiag[i] - bdiag[i + 1];
336: bjtmp = bj + bdiag[i + 1] + 1;
337: for (j = 0; j < nz; j++) PetscCall(PetscArrayzero(rtmp + bs2 * bjtmp[j], bs2));
339: /* load in initial (unfactored row) */
340: nz = ai[r[i] + 1] - ai[r[i]];
341: ajtmp = aj + ai[r[i]];
342: v = aa + bs2 * ai[r[i]];
343: for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(rtmp + bs2 * ic[ajtmp[j]], v + bs2 * j, bs2));
345: /* elimination */
346: bjtmp = bj + bi[i];
347: nzL = bi[i + 1] - bi[i];
348: for (k = 0; k < nzL; k++) {
349: row = bjtmp[k];
350: pc = rtmp + bs2 * row;
351: for (flg = 0, j = 0; j < bs2; j++) {
352: if (pc[j] != 0.0) {
353: flg = 1;
354: break;
355: }
356: }
357: if (flg) {
358: pv = b->a + bs2 * bdiag[row];
359: /* PetscKernel_A_gets_A_times_B(bs,pc,pv,mwork); *pc = *pc * (*pv); */
360: PetscCall(PetscKernel_A_gets_A_times_B_5(pc, pv, mwork));
362: pj = b->j + bdiag[row + 1] + 1; /* beginning of U(row,:) */
363: pv = b->a + bs2 * (bdiag[row + 1] + 1);
364: nz = bdiag[row] - bdiag[row + 1] - 1; /* num of entries inU(row,:), excluding diag */
365: for (j = 0; j < nz; j++) {
366: /* PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j); */
367: /* rtmp+bs2*pj[j] = rtmp+bs2*pj[j] - (*pc)*(pv+bs2*j) */
368: v = rtmp + bs2 * pj[j];
369: PetscCall(PetscKernel_A_gets_A_minus_B_times_C_5(v, pc, pv));
370: pv += bs2;
371: }
372: PetscCall(PetscLogFlops(250.0 * nz + 225)); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */
373: }
374: }
376: /* finished row so stick it into b->a */
377: /* L part */
378: pv = b->a + bs2 * bi[i];
379: pj = b->j + bi[i];
380: nz = bi[i + 1] - bi[i];
381: for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2));
383: /* Mark diagonal and invert diagonal for simpler triangular solves */
384: pv = b->a + bs2 * bdiag[i];
385: pj = b->j + bdiag[i];
386: PetscCall(PetscArraycpy(pv, rtmp + bs2 * pj[0], bs2));
387: PetscCall(PetscKernel_A_gets_inverse_A_5(pv, ipvt, work, shift, allowzeropivot, &zeropivotdetected));
388: if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
390: /* U part */
391: pv = b->a + bs2 * (bdiag[i + 1] + 1);
392: pj = b->j + bdiag[i + 1] + 1;
393: nz = bdiag[i] - bdiag[i + 1] - 1;
394: for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2));
395: }
397: PetscCall(PetscFree2(rtmp, mwork));
398: PetscCall(ISRestoreIndices(isicol, &ic));
399: PetscCall(ISRestoreIndices(isrow, &r));
401: C->ops->solve = MatSolve_SeqBAIJ_5;
402: C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_5;
403: C->assembled = PETSC_TRUE;
405: PetscCall(PetscLogFlops(1.333333333333 * 5 * 5 * 5 * n)); /* from inverting diagonal blocks */
406: PetscFunctionReturn(PETSC_SUCCESS);
407: }
409: /*
410: Version for when blocks are 5 by 5 Using natural ordering
411: */
412: PetscErrorCode MatILUFactorNumeric_SeqBAIJ_5_NaturalOrdering_inplace(Mat C, Mat A, const MatFactorInfo *info)
413: {
414: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)C->data;
415: PetscInt i, j, n = a->mbs, *bi = b->i, *bj = b->j, ipvt[5];
416: PetscInt *ajtmpold, *ajtmp, nz, row;
417: PetscInt *diag_offset = b->diag, *ai = a->i, *aj = a->j, *pj;
418: MatScalar *pv, *v, *rtmp, *pc, *w, *x;
419: MatScalar x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15;
420: MatScalar x16, x17, x18, x19, x20, x21, x22, x23, x24, x25;
421: MatScalar p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12, p13, p14, p15;
422: MatScalar p16, p17, p18, p19, p20, p21, p22, p23, p24, p25;
423: MatScalar m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, m11, m12, m13, m14, m15;
424: MatScalar m16, m17, m18, m19, m20, m21, m22, m23, m24, m25;
425: MatScalar *ba = b->a, *aa = a->a, work[25];
426: PetscReal shift = info->shiftamount;
427: PetscBool allowzeropivot, zeropivotdetected;
429: PetscFunctionBegin;
430: allowzeropivot = PetscNot(A->erroriffailure);
431: PetscCall(PetscMalloc1(25 * (n + 1), &rtmp));
432: for (i = 0; i < n; i++) {
433: nz = bi[i + 1] - bi[i];
434: ajtmp = bj + bi[i];
435: for (j = 0; j < nz; j++) {
436: x = rtmp + 25 * ajtmp[j];
437: x[0] = x[1] = x[2] = x[3] = x[4] = x[5] = x[6] = x[7] = x[8] = x[9] = 0.0;
438: x[10] = x[11] = x[12] = x[13] = x[14] = x[15] = 0.0;
439: x[16] = x[17] = x[18] = x[19] = x[20] = x[21] = x[22] = x[23] = x[24] = 0.0;
440: }
441: /* load in initial (unfactored row) */
442: nz = ai[i + 1] - ai[i];
443: ajtmpold = aj + ai[i];
444: v = aa + 25 * ai[i];
445: for (j = 0; j < nz; j++) {
446: x = rtmp + 25 * ajtmpold[j];
447: x[0] = v[0];
448: x[1] = v[1];
449: x[2] = v[2];
450: x[3] = v[3];
451: x[4] = v[4];
452: x[5] = v[5];
453: x[6] = v[6];
454: x[7] = v[7];
455: x[8] = v[8];
456: x[9] = v[9];
457: x[10] = v[10];
458: x[11] = v[11];
459: x[12] = v[12];
460: x[13] = v[13];
461: x[14] = v[14];
462: x[15] = v[15];
463: x[16] = v[16];
464: x[17] = v[17];
465: x[18] = v[18];
466: x[19] = v[19];
467: x[20] = v[20];
468: x[21] = v[21];
469: x[22] = v[22];
470: x[23] = v[23];
471: x[24] = v[24];
472: v += 25;
473: }
474: row = *ajtmp++;
475: while (row < i) {
476: pc = rtmp + 25 * row;
477: p1 = pc[0];
478: p2 = pc[1];
479: p3 = pc[2];
480: p4 = pc[3];
481: p5 = pc[4];
482: p6 = pc[5];
483: p7 = pc[6];
484: p8 = pc[7];
485: p9 = pc[8];
486: p10 = pc[9];
487: p11 = pc[10];
488: p12 = pc[11];
489: p13 = pc[12];
490: p14 = pc[13];
491: p15 = pc[14];
492: p16 = pc[15];
493: p17 = pc[16];
494: p18 = pc[17];
495: p19 = pc[18];
496: p20 = pc[19];
497: p21 = pc[20];
498: p22 = pc[21];
499: p23 = pc[22];
500: p24 = pc[23];
501: p25 = pc[24];
502: if (p1 != 0.0 || p2 != 0.0 || p3 != 0.0 || p4 != 0.0 || p5 != 0.0 || p6 != 0.0 || p7 != 0.0 || p8 != 0.0 || p9 != 0.0 || p10 != 0.0 || p11 != 0.0 || p12 != 0.0 || p13 != 0.0 || p14 != 0.0 || p15 != 0.0 || p16 != 0.0 || p17 != 0.0 || p18 != 0.0 || p19 != 0.0 || p20 != 0.0 || p21 != 0.0 || p22 != 0.0 || p23 != 0.0 || p24 != 0.0 || p25 != 0.0) {
503: pv = ba + 25 * diag_offset[row];
504: pj = bj + diag_offset[row] + 1;
505: x1 = pv[0];
506: x2 = pv[1];
507: x3 = pv[2];
508: x4 = pv[3];
509: x5 = pv[4];
510: x6 = pv[5];
511: x7 = pv[6];
512: x8 = pv[7];
513: x9 = pv[8];
514: x10 = pv[9];
515: x11 = pv[10];
516: x12 = pv[11];
517: x13 = pv[12];
518: x14 = pv[13];
519: x15 = pv[14];
520: x16 = pv[15];
521: x17 = pv[16];
522: x18 = pv[17];
523: x19 = pv[18];
524: x20 = pv[19];
525: x21 = pv[20];
526: x22 = pv[21];
527: x23 = pv[22];
528: x24 = pv[23];
529: x25 = pv[24];
530: pc[0] = m1 = p1 * x1 + p6 * x2 + p11 * x3 + p16 * x4 + p21 * x5;
531: pc[1] = m2 = p2 * x1 + p7 * x2 + p12 * x3 + p17 * x4 + p22 * x5;
532: pc[2] = m3 = p3 * x1 + p8 * x2 + p13 * x3 + p18 * x4 + p23 * x5;
533: pc[3] = m4 = p4 * x1 + p9 * x2 + p14 * x3 + p19 * x4 + p24 * x5;
534: pc[4] = m5 = p5 * x1 + p10 * x2 + p15 * x3 + p20 * x4 + p25 * x5;
536: pc[5] = m6 = p1 * x6 + p6 * x7 + p11 * x8 + p16 * x9 + p21 * x10;
537: pc[6] = m7 = p2 * x6 + p7 * x7 + p12 * x8 + p17 * x9 + p22 * x10;
538: pc[7] = m8 = p3 * x6 + p8 * x7 + p13 * x8 + p18 * x9 + p23 * x10;
539: pc[8] = m9 = p4 * x6 + p9 * x7 + p14 * x8 + p19 * x9 + p24 * x10;
540: pc[9] = m10 = p5 * x6 + p10 * x7 + p15 * x8 + p20 * x9 + p25 * x10;
542: pc[10] = m11 = p1 * x11 + p6 * x12 + p11 * x13 + p16 * x14 + p21 * x15;
543: pc[11] = m12 = p2 * x11 + p7 * x12 + p12 * x13 + p17 * x14 + p22 * x15;
544: pc[12] = m13 = p3 * x11 + p8 * x12 + p13 * x13 + p18 * x14 + p23 * x15;
545: pc[13] = m14 = p4 * x11 + p9 * x12 + p14 * x13 + p19 * x14 + p24 * x15;
546: pc[14] = m15 = p5 * x11 + p10 * x12 + p15 * x13 + p20 * x14 + p25 * x15;
548: pc[15] = m16 = p1 * x16 + p6 * x17 + p11 * x18 + p16 * x19 + p21 * x20;
549: pc[16] = m17 = p2 * x16 + p7 * x17 + p12 * x18 + p17 * x19 + p22 * x20;
550: pc[17] = m18 = p3 * x16 + p8 * x17 + p13 * x18 + p18 * x19 + p23 * x20;
551: pc[18] = m19 = p4 * x16 + p9 * x17 + p14 * x18 + p19 * x19 + p24 * x20;
552: pc[19] = m20 = p5 * x16 + p10 * x17 + p15 * x18 + p20 * x19 + p25 * x20;
554: pc[20] = m21 = p1 * x21 + p6 * x22 + p11 * x23 + p16 * x24 + p21 * x25;
555: pc[21] = m22 = p2 * x21 + p7 * x22 + p12 * x23 + p17 * x24 + p22 * x25;
556: pc[22] = m23 = p3 * x21 + p8 * x22 + p13 * x23 + p18 * x24 + p23 * x25;
557: pc[23] = m24 = p4 * x21 + p9 * x22 + p14 * x23 + p19 * x24 + p24 * x25;
558: pc[24] = m25 = p5 * x21 + p10 * x22 + p15 * x23 + p20 * x24 + p25 * x25;
560: nz = bi[row + 1] - diag_offset[row] - 1;
561: pv += 25;
562: for (j = 0; j < nz; j++) {
563: x1 = pv[0];
564: x2 = pv[1];
565: x3 = pv[2];
566: x4 = pv[3];
567: x5 = pv[4];
568: x6 = pv[5];
569: x7 = pv[6];
570: x8 = pv[7];
571: x9 = pv[8];
572: x10 = pv[9];
573: x11 = pv[10];
574: x12 = pv[11];
575: x13 = pv[12];
576: x14 = pv[13];
577: x15 = pv[14];
578: x16 = pv[15];
579: x17 = pv[16];
580: x18 = pv[17];
581: x19 = pv[18];
582: x20 = pv[19];
583: x21 = pv[20];
584: x22 = pv[21];
585: x23 = pv[22];
586: x24 = pv[23];
587: x25 = pv[24];
588: x = rtmp + 25 * pj[j];
589: x[0] -= m1 * x1 + m6 * x2 + m11 * x3 + m16 * x4 + m21 * x5;
590: x[1] -= m2 * x1 + m7 * x2 + m12 * x3 + m17 * x4 + m22 * x5;
591: x[2] -= m3 * x1 + m8 * x2 + m13 * x3 + m18 * x4 + m23 * x5;
592: x[3] -= m4 * x1 + m9 * x2 + m14 * x3 + m19 * x4 + m24 * x5;
593: x[4] -= m5 * x1 + m10 * x2 + m15 * x3 + m20 * x4 + m25 * x5;
595: x[5] -= m1 * x6 + m6 * x7 + m11 * x8 + m16 * x9 + m21 * x10;
596: x[6] -= m2 * x6 + m7 * x7 + m12 * x8 + m17 * x9 + m22 * x10;
597: x[7] -= m3 * x6 + m8 * x7 + m13 * x8 + m18 * x9 + m23 * x10;
598: x[8] -= m4 * x6 + m9 * x7 + m14 * x8 + m19 * x9 + m24 * x10;
599: x[9] -= m5 * x6 + m10 * x7 + m15 * x8 + m20 * x9 + m25 * x10;
601: x[10] -= m1 * x11 + m6 * x12 + m11 * x13 + m16 * x14 + m21 * x15;
602: x[11] -= m2 * x11 + m7 * x12 + m12 * x13 + m17 * x14 + m22 * x15;
603: x[12] -= m3 * x11 + m8 * x12 + m13 * x13 + m18 * x14 + m23 * x15;
604: x[13] -= m4 * x11 + m9 * x12 + m14 * x13 + m19 * x14 + m24 * x15;
605: x[14] -= m5 * x11 + m10 * x12 + m15 * x13 + m20 * x14 + m25 * x15;
607: x[15] -= m1 * x16 + m6 * x17 + m11 * x18 + m16 * x19 + m21 * x20;
608: x[16] -= m2 * x16 + m7 * x17 + m12 * x18 + m17 * x19 + m22 * x20;
609: x[17] -= m3 * x16 + m8 * x17 + m13 * x18 + m18 * x19 + m23 * x20;
610: x[18] -= m4 * x16 + m9 * x17 + m14 * x18 + m19 * x19 + m24 * x20;
611: x[19] -= m5 * x16 + m10 * x17 + m15 * x18 + m20 * x19 + m25 * x20;
613: x[20] -= m1 * x21 + m6 * x22 + m11 * x23 + m16 * x24 + m21 * x25;
614: x[21] -= m2 * x21 + m7 * x22 + m12 * x23 + m17 * x24 + m22 * x25;
615: x[22] -= m3 * x21 + m8 * x22 + m13 * x23 + m18 * x24 + m23 * x25;
616: x[23] -= m4 * x21 + m9 * x22 + m14 * x23 + m19 * x24 + m24 * x25;
617: x[24] -= m5 * x21 + m10 * x22 + m15 * x23 + m20 * x24 + m25 * x25;
618: pv += 25;
619: }
620: PetscCall(PetscLogFlops(250.0 * nz + 225.0));
621: }
622: row = *ajtmp++;
623: }
624: /* finished row so stick it into b->a */
625: pv = ba + 25 * bi[i];
626: pj = bj + bi[i];
627: nz = bi[i + 1] - bi[i];
628: for (j = 0; j < nz; j++) {
629: x = rtmp + 25 * pj[j];
630: pv[0] = x[0];
631: pv[1] = x[1];
632: pv[2] = x[2];
633: pv[3] = x[3];
634: pv[4] = x[4];
635: pv[5] = x[5];
636: pv[6] = x[6];
637: pv[7] = x[7];
638: pv[8] = x[8];
639: pv[9] = x[9];
640: pv[10] = x[10];
641: pv[11] = x[11];
642: pv[12] = x[12];
643: pv[13] = x[13];
644: pv[14] = x[14];
645: pv[15] = x[15];
646: pv[16] = x[16];
647: pv[17] = x[17];
648: pv[18] = x[18];
649: pv[19] = x[19];
650: pv[20] = x[20];
651: pv[21] = x[21];
652: pv[22] = x[22];
653: pv[23] = x[23];
654: pv[24] = x[24];
655: pv += 25;
656: }
657: /* invert diagonal block */
658: w = ba + 25 * diag_offset[i];
659: PetscCall(PetscKernel_A_gets_inverse_A_5(w, ipvt, work, shift, allowzeropivot, &zeropivotdetected));
660: if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
661: }
663: PetscCall(PetscFree(rtmp));
665: C->ops->solve = MatSolve_SeqBAIJ_5_NaturalOrdering_inplace;
666: C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_5_NaturalOrdering_inplace;
667: C->assembled = PETSC_TRUE;
669: PetscCall(PetscLogFlops(1.333333333333 * 5 * 5 * 5 * b->mbs)); /* from inverting diagonal blocks */
670: PetscFunctionReturn(PETSC_SUCCESS);
671: }
673: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_5_NaturalOrdering(Mat B, Mat A, const MatFactorInfo *info)
674: {
675: Mat C = B;
676: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)C->data;
677: PetscInt i, j, k, nz, nzL, row;
678: const PetscInt n = a->mbs, *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j;
679: const PetscInt *ajtmp, *bjtmp, *bdiag = b->diag, *pj, bs2 = a->bs2;
680: MatScalar *rtmp, *pc, *mwork, *v, *vv, *pv, *aa = a->a, work[25];
681: PetscInt flg, ipvt[5];
682: PetscReal shift = info->shiftamount;
683: PetscBool allowzeropivot, zeropivotdetected;
685: PetscFunctionBegin;
686: allowzeropivot = PetscNot(A->erroriffailure);
688: /* generate work space needed by the factorization */
689: PetscCall(PetscMalloc2(bs2 * n, &rtmp, bs2, &mwork));
690: PetscCall(PetscArrayzero(rtmp, bs2 * n));
692: for (i = 0; i < n; i++) {
693: /* zero rtmp */
694: /* L part */
695: nz = bi[i + 1] - bi[i];
696: bjtmp = bj + bi[i];
697: for (j = 0; j < nz; j++) PetscCall(PetscArrayzero(rtmp + bs2 * bjtmp[j], bs2));
699: /* U part */
700: nz = bdiag[i] - bdiag[i + 1];
701: bjtmp = bj + bdiag[i + 1] + 1;
702: for (j = 0; j < nz; j++) PetscCall(PetscArrayzero(rtmp + bs2 * bjtmp[j], bs2));
704: /* load in initial (unfactored row) */
705: nz = ai[i + 1] - ai[i];
706: ajtmp = aj + ai[i];
707: v = aa + bs2 * ai[i];
708: for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(rtmp + bs2 * ajtmp[j], v + bs2 * j, bs2));
710: /* elimination */
711: bjtmp = bj + bi[i];
712: nzL = bi[i + 1] - bi[i];
713: for (k = 0; k < nzL; k++) {
714: row = bjtmp[k];
715: pc = rtmp + bs2 * row;
716: for (flg = 0, j = 0; j < bs2; j++) {
717: if (pc[j] != 0.0) {
718: flg = 1;
719: break;
720: }
721: }
722: if (flg) {
723: pv = b->a + bs2 * bdiag[row];
724: /* PetscKernel_A_gets_A_times_B(bs,pc,pv,mwork); *pc = *pc * (*pv); */
725: PetscCall(PetscKernel_A_gets_A_times_B_5(pc, pv, mwork));
727: pj = b->j + bdiag[row + 1] + 1; /* beginning of U(row,:) */
728: pv = b->a + bs2 * (bdiag[row + 1] + 1);
729: nz = bdiag[row] - bdiag[row + 1] - 1; /* num of entries inU(row,:), excluding diag */
730: for (j = 0; j < nz; j++) {
731: /* PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j); */
732: /* rtmp+bs2*pj[j] = rtmp+bs2*pj[j] - (*pc)*(pv+bs2*j) */
733: vv = rtmp + bs2 * pj[j];
734: PetscCall(PetscKernel_A_gets_A_minus_B_times_C_5(vv, pc, pv));
735: pv += bs2;
736: }
737: PetscCall(PetscLogFlops(250.0 * nz + 225)); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */
738: }
739: }
741: /* finished row so stick it into b->a */
742: /* L part */
743: pv = b->a + bs2 * bi[i];
744: pj = b->j + bi[i];
745: nz = bi[i + 1] - bi[i];
746: for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2));
748: /* Mark diagonal and invert diagonal for simpler triangular solves */
749: pv = b->a + bs2 * bdiag[i];
750: pj = b->j + bdiag[i];
751: PetscCall(PetscArraycpy(pv, rtmp + bs2 * pj[0], bs2));
752: PetscCall(PetscKernel_A_gets_inverse_A_5(pv, ipvt, work, shift, allowzeropivot, &zeropivotdetected));
753: if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
755: /* U part */
756: pv = b->a + bs2 * (bdiag[i + 1] + 1);
757: pj = b->j + bdiag[i + 1] + 1;
758: nz = bdiag[i] - bdiag[i + 1] - 1;
759: for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2));
760: }
761: PetscCall(PetscFree2(rtmp, mwork));
763: C->ops->solve = MatSolve_SeqBAIJ_5_NaturalOrdering;
764: C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_5_NaturalOrdering;
765: C->assembled = PETSC_TRUE;
767: PetscCall(PetscLogFlops(1.333333333333 * 5 * 5 * 5 * n)); /* from inverting diagonal blocks */
768: PetscFunctionReturn(PETSC_SUCCESS);
769: }
771: /*
772: Version for when blocks are 9 by 9
773: */
774: #if PetscDefined(HAVE_IMMINTRIN_H) && defined(__AVX2__) && defined(__FMA__) && PetscDefined(USE_REAL_DOUBLE) && !PetscDefined(USE_COMPLEX) && !PetscDefined(USE_64BIT_INDICES)
775: #include <immintrin.h>
776: PetscErrorCode MatLUFactorNumeric_SeqBAIJ_9_NaturalOrdering(Mat B, Mat A, const MatFactorInfo *info)
777: {
778: Mat C = B;
779: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)C->data;
780: PetscInt i, j, k, nz, nzL, row;
781: const PetscInt n = a->mbs, *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j;
782: const PetscInt *ajtmp, *bjtmp, *bdiag = b->diag, *pj, bs2 = a->bs2;
783: MatScalar *rtmp, *pc, *mwork, *v, *pv, *aa = a->a;
784: PetscInt flg;
785: PetscReal shift = info->shiftamount;
786: PetscBool allowzeropivot, zeropivotdetected;
788: PetscFunctionBegin;
789: allowzeropivot = PetscNot(A->erroriffailure);
791: /* generate work space needed by the factorization */
792: PetscCall(PetscMalloc2(bs2 * n, &rtmp, bs2, &mwork));
793: PetscCall(PetscArrayzero(rtmp, bs2 * n));
795: for (i = 0; i < n; i++) {
796: /* zero rtmp */
797: /* L part */
798: nz = bi[i + 1] - bi[i];
799: bjtmp = bj + bi[i];
800: for (j = 0; j < nz; j++) PetscCall(PetscArrayzero(rtmp + bs2 * bjtmp[j], bs2));
802: /* U part */
803: nz = bdiag[i] - bdiag[i + 1];
804: bjtmp = bj + bdiag[i + 1] + 1;
805: for (j = 0; j < nz; j++) PetscCall(PetscArrayzero(rtmp + bs2 * bjtmp[j], bs2));
807: /* load in initial (unfactored row) */
808: nz = ai[i + 1] - ai[i];
809: ajtmp = aj + ai[i];
810: v = aa + bs2 * ai[i];
811: for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(rtmp + bs2 * ajtmp[j], v + bs2 * j, bs2));
813: /* elimination */
814: bjtmp = bj + bi[i];
815: nzL = bi[i + 1] - bi[i];
816: for (k = 0; k < nzL; k++) {
817: row = bjtmp[k];
818: pc = rtmp + bs2 * row;
819: for (flg = 0, j = 0; j < bs2; j++) {
820: if (pc[j] != 0.0) {
821: flg = 1;
822: break;
823: }
824: }
825: if (flg) {
826: pv = b->a + bs2 * bdiag[row];
827: /* PetscKernel_A_gets_A_times_B(bs,pc,pv,mwork); *pc = *pc * (*pv); */
828: PetscCall(PetscKernel_A_gets_A_times_B_9(pc, pv, mwork));
830: pj = b->j + bdiag[row + 1] + 1; /* beginning of U(row,:) */
831: pv = b->a + bs2 * (bdiag[row + 1] + 1);
832: nz = bdiag[row] - bdiag[row + 1] - 1; /* num of entries inU(row,:), excluding diag */
833: for (j = 0; j < nz; j++) {
834: /* PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j); */
835: /* rtmp+bs2*pj[j] = rtmp+bs2*pj[j] - (*pc)*(pv+bs2*j) */
836: v = rtmp + bs2 * pj[j];
837: PetscCall(PetscKernel_A_gets_A_minus_B_times_C_9(v, pc, pv + 81 * j));
838: /* pv incremented in PetscKernel_A_gets_A_minus_B_times_C_9 */
839: }
840: PetscCall(PetscLogFlops(1458 * nz + 1377)); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */
841: }
842: }
844: /* finished row so stick it into b->a */
845: /* L part */
846: pv = b->a + bs2 * bi[i];
847: pj = b->j + bi[i];
848: nz = bi[i + 1] - bi[i];
849: for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2));
851: /* Mark diagonal and invert diagonal for simpler triangular solves */
852: pv = b->a + bs2 * bdiag[i];
853: pj = b->j + bdiag[i];
854: PetscCall(PetscArraycpy(pv, rtmp + bs2 * pj[0], bs2));
855: PetscCall(PetscKernel_A_gets_inverse_A_9(pv, shift, allowzeropivot, &zeropivotdetected));
856: if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
858: /* U part */
859: pv = b->a + bs2 * (bdiag[i + 1] + 1);
860: pj = b->j + bdiag[i + 1] + 1;
861: nz = bdiag[i] - bdiag[i + 1] - 1;
862: for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2));
863: }
864: PetscCall(PetscFree2(rtmp, mwork));
866: C->ops->solve = MatSolve_SeqBAIJ_9_NaturalOrdering;
867: C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_N;
868: C->assembled = PETSC_TRUE;
870: PetscCall(PetscLogFlops(1.333333333333 * 9 * 9 * 9 * n)); /* from inverting diagonal blocks */
871: PetscFunctionReturn(PETSC_SUCCESS);
872: }
874: PetscErrorCode MatSolve_SeqBAIJ_9_NaturalOrdering(Mat A, Vec bb, Vec xx)
875: {
876: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
877: const PetscInt *ai = a->i, *aj = a->j, *adiag = a->diag, *vi;
878: PetscInt i, k, n = a->mbs;
879: PetscInt nz, bs = A->rmap->bs, bs2 = a->bs2;
880: const MatScalar *aa = a->a, *v;
881: PetscScalar *x, *s, *t, *ls;
882: const PetscScalar *b;
883: __m256d a0, a1, a2, a3, a4, a5, w0, w1, w2, w3, s0, s1, s2, v0, v1, v2, v3;
885: PetscFunctionBegin;
886: PetscCall(VecGetArrayRead(bb, &b));
887: PetscCall(VecGetArray(xx, &x));
888: t = a->solve_work;
890: /* forward solve the lower triangular */
891: PetscCall(PetscArraycpy(t, b, bs)); /* copy 1st block of b to t */
893: for (i = 1; i < n; i++) {
894: v = aa + bs2 * ai[i];
895: vi = aj + ai[i];
896: nz = ai[i + 1] - ai[i];
897: s = t + bs * i;
898: PetscCall(PetscArraycpy(s, b + bs * i, bs)); /* copy i_th block of b to t */
900: __m256d s0, s1, s2;
901: s0 = _mm256_loadu_pd(s + 0);
902: s1 = _mm256_loadu_pd(s + 4);
903: s2 = _mm256_maskload_pd(s + 8, _mm256_set_epi64x(0LL, 0LL, 0LL, 1LL << 63));
905: for (k = 0; k < nz; k++) {
906: w0 = _mm256_set1_pd((t + bs * vi[k])[0]);
907: a0 = _mm256_loadu_pd(&v[0]);
908: s0 = _mm256_fnmadd_pd(a0, w0, s0);
909: a1 = _mm256_loadu_pd(&v[4]);
910: s1 = _mm256_fnmadd_pd(a1, w0, s1);
911: a2 = _mm256_loadu_pd(&v[8]);
912: s2 = _mm256_fnmadd_pd(a2, w0, s2);
914: w1 = _mm256_set1_pd((t + bs * vi[k])[1]);
915: a3 = _mm256_loadu_pd(&v[9]);
916: s0 = _mm256_fnmadd_pd(a3, w1, s0);
917: a4 = _mm256_loadu_pd(&v[13]);
918: s1 = _mm256_fnmadd_pd(a4, w1, s1);
919: a5 = _mm256_loadu_pd(&v[17]);
920: s2 = _mm256_fnmadd_pd(a5, w1, s2);
922: w2 = _mm256_set1_pd((t + bs * vi[k])[2]);
923: a0 = _mm256_loadu_pd(&v[18]);
924: s0 = _mm256_fnmadd_pd(a0, w2, s0);
925: a1 = _mm256_loadu_pd(&v[22]);
926: s1 = _mm256_fnmadd_pd(a1, w2, s1);
927: a2 = _mm256_loadu_pd(&v[26]);
928: s2 = _mm256_fnmadd_pd(a2, w2, s2);
930: w3 = _mm256_set1_pd((t + bs * vi[k])[3]);
931: a3 = _mm256_loadu_pd(&v[27]);
932: s0 = _mm256_fnmadd_pd(a3, w3, s0);
933: a4 = _mm256_loadu_pd(&v[31]);
934: s1 = _mm256_fnmadd_pd(a4, w3, s1);
935: a5 = _mm256_loadu_pd(&v[35]);
936: s2 = _mm256_fnmadd_pd(a5, w3, s2);
938: w0 = _mm256_set1_pd((t + bs * vi[k])[4]);
939: a0 = _mm256_loadu_pd(&v[36]);
940: s0 = _mm256_fnmadd_pd(a0, w0, s0);
941: a1 = _mm256_loadu_pd(&v[40]);
942: s1 = _mm256_fnmadd_pd(a1, w0, s1);
943: a2 = _mm256_loadu_pd(&v[44]);
944: s2 = _mm256_fnmadd_pd(a2, w0, s2);
946: w1 = _mm256_set1_pd((t + bs * vi[k])[5]);
947: a3 = _mm256_loadu_pd(&v[45]);
948: s0 = _mm256_fnmadd_pd(a3, w1, s0);
949: a4 = _mm256_loadu_pd(&v[49]);
950: s1 = _mm256_fnmadd_pd(a4, w1, s1);
951: a5 = _mm256_loadu_pd(&v[53]);
952: s2 = _mm256_fnmadd_pd(a5, w1, s2);
954: w2 = _mm256_set1_pd((t + bs * vi[k])[6]);
955: a0 = _mm256_loadu_pd(&v[54]);
956: s0 = _mm256_fnmadd_pd(a0, w2, s0);
957: a1 = _mm256_loadu_pd(&v[58]);
958: s1 = _mm256_fnmadd_pd(a1, w2, s1);
959: a2 = _mm256_loadu_pd(&v[62]);
960: s2 = _mm256_fnmadd_pd(a2, w2, s2);
962: w3 = _mm256_set1_pd((t + bs * vi[k])[7]);
963: a3 = _mm256_loadu_pd(&v[63]);
964: s0 = _mm256_fnmadd_pd(a3, w3, s0);
965: a4 = _mm256_loadu_pd(&v[67]);
966: s1 = _mm256_fnmadd_pd(a4, w3, s1);
967: a5 = _mm256_loadu_pd(&v[71]);
968: s2 = _mm256_fnmadd_pd(a5, w3, s2);
970: w0 = _mm256_set1_pd((t + bs * vi[k])[8]);
971: a0 = _mm256_loadu_pd(&v[72]);
972: s0 = _mm256_fnmadd_pd(a0, w0, s0);
973: a1 = _mm256_loadu_pd(&v[76]);
974: s1 = _mm256_fnmadd_pd(a1, w0, s1);
975: a2 = _mm256_maskload_pd(v + 80, _mm256_set_epi64x(0LL, 0LL, 0LL, 1LL << 63));
976: s2 = _mm256_fnmadd_pd(a2, w0, s2);
977: v += bs2;
978: }
979: _mm256_storeu_pd(&s[0], s0);
980: _mm256_storeu_pd(&s[4], s1);
981: _mm256_maskstore_pd(&s[8], _mm256_set_epi64x(0LL, 0LL, 0LL, 1LL << 63), s2);
982: }
984: /* backward solve the upper triangular */
985: ls = a->solve_work + A->cmap->n;
986: for (i = n - 1; i >= 0; i--) {
987: v = aa + bs2 * (adiag[i + 1] + 1);
988: vi = aj + adiag[i + 1] + 1;
989: nz = adiag[i] - adiag[i + 1] - 1;
990: PetscCall(PetscArraycpy(ls, t + i * bs, bs));
992: s0 = _mm256_loadu_pd(ls + 0);
993: s1 = _mm256_loadu_pd(ls + 4);
994: s2 = _mm256_maskload_pd(ls + 8, _mm256_set_epi64x(0LL, 0LL, 0LL, 1LL << 63));
996: for (k = 0; k < nz; k++) {
997: w0 = _mm256_set1_pd((t + bs * vi[k])[0]);
998: a0 = _mm256_loadu_pd(&v[0]);
999: s0 = _mm256_fnmadd_pd(a0, w0, s0);
1000: a1 = _mm256_loadu_pd(&v[4]);
1001: s1 = _mm256_fnmadd_pd(a1, w0, s1);
1002: a2 = _mm256_loadu_pd(&v[8]);
1003: s2 = _mm256_fnmadd_pd(a2, w0, s2);
1005: /* v += 9; */
1006: w1 = _mm256_set1_pd((t + bs * vi[k])[1]);
1007: a3 = _mm256_loadu_pd(&v[9]);
1008: s0 = _mm256_fnmadd_pd(a3, w1, s0);
1009: a4 = _mm256_loadu_pd(&v[13]);
1010: s1 = _mm256_fnmadd_pd(a4, w1, s1);
1011: a5 = _mm256_loadu_pd(&v[17]);
1012: s2 = _mm256_fnmadd_pd(a5, w1, s2);
1014: /* v += 9; */
1015: w2 = _mm256_set1_pd((t + bs * vi[k])[2]);
1016: a0 = _mm256_loadu_pd(&v[18]);
1017: s0 = _mm256_fnmadd_pd(a0, w2, s0);
1018: a1 = _mm256_loadu_pd(&v[22]);
1019: s1 = _mm256_fnmadd_pd(a1, w2, s1);
1020: a2 = _mm256_loadu_pd(&v[26]);
1021: s2 = _mm256_fnmadd_pd(a2, w2, s2);
1023: /* v += 9; */
1024: w3 = _mm256_set1_pd((t + bs * vi[k])[3]);
1025: a3 = _mm256_loadu_pd(&v[27]);
1026: s0 = _mm256_fnmadd_pd(a3, w3, s0);
1027: a4 = _mm256_loadu_pd(&v[31]);
1028: s1 = _mm256_fnmadd_pd(a4, w3, s1);
1029: a5 = _mm256_loadu_pd(&v[35]);
1030: s2 = _mm256_fnmadd_pd(a5, w3, s2);
1032: /* v += 9; */
1033: w0 = _mm256_set1_pd((t + bs * vi[k])[4]);
1034: a0 = _mm256_loadu_pd(&v[36]);
1035: s0 = _mm256_fnmadd_pd(a0, w0, s0);
1036: a1 = _mm256_loadu_pd(&v[40]);
1037: s1 = _mm256_fnmadd_pd(a1, w0, s1);
1038: a2 = _mm256_loadu_pd(&v[44]);
1039: s2 = _mm256_fnmadd_pd(a2, w0, s2);
1041: /* v += 9; */
1042: w1 = _mm256_set1_pd((t + bs * vi[k])[5]);
1043: a3 = _mm256_loadu_pd(&v[45]);
1044: s0 = _mm256_fnmadd_pd(a3, w1, s0);
1045: a4 = _mm256_loadu_pd(&v[49]);
1046: s1 = _mm256_fnmadd_pd(a4, w1, s1);
1047: a5 = _mm256_loadu_pd(&v[53]);
1048: s2 = _mm256_fnmadd_pd(a5, w1, s2);
1050: /* v += 9; */
1051: w2 = _mm256_set1_pd((t + bs * vi[k])[6]);
1052: a0 = _mm256_loadu_pd(&v[54]);
1053: s0 = _mm256_fnmadd_pd(a0, w2, s0);
1054: a1 = _mm256_loadu_pd(&v[58]);
1055: s1 = _mm256_fnmadd_pd(a1, w2, s1);
1056: a2 = _mm256_loadu_pd(&v[62]);
1057: s2 = _mm256_fnmadd_pd(a2, w2, s2);
1059: /* v += 9; */
1060: w3 = _mm256_set1_pd((t + bs * vi[k])[7]);
1061: a3 = _mm256_loadu_pd(&v[63]);
1062: s0 = _mm256_fnmadd_pd(a3, w3, s0);
1063: a4 = _mm256_loadu_pd(&v[67]);
1064: s1 = _mm256_fnmadd_pd(a4, w3, s1);
1065: a5 = _mm256_loadu_pd(&v[71]);
1066: s2 = _mm256_fnmadd_pd(a5, w3, s2);
1068: /* v += 9; */
1069: w0 = _mm256_set1_pd((t + bs * vi[k])[8]);
1070: a0 = _mm256_loadu_pd(&v[72]);
1071: s0 = _mm256_fnmadd_pd(a0, w0, s0);
1072: a1 = _mm256_loadu_pd(&v[76]);
1073: s1 = _mm256_fnmadd_pd(a1, w0, s1);
1074: a2 = _mm256_maskload_pd(v + 80, _mm256_set_epi64x(0LL, 0LL, 0LL, 1LL << 63));
1075: s2 = _mm256_fnmadd_pd(a2, w0, s2);
1076: v += bs2;
1077: }
1079: _mm256_storeu_pd(&ls[0], s0);
1080: _mm256_storeu_pd(&ls[4], s1);
1081: _mm256_maskstore_pd(&ls[8], _mm256_set_epi64x(0LL, 0LL, 0LL, 1LL << 63), s2);
1083: w0 = _mm256_setzero_pd();
1084: w1 = _mm256_setzero_pd();
1085: w2 = _mm256_setzero_pd();
1087: /* first row */
1088: v0 = _mm256_set1_pd(ls[0]);
1089: a0 = _mm256_loadu_pd(&(aa + bs2 * adiag[i])[0]);
1090: w0 = _mm256_fmadd_pd(a0, v0, w0);
1091: a1 = _mm256_loadu_pd(&(aa + bs2 * adiag[i])[4]);
1092: w1 = _mm256_fmadd_pd(a1, v0, w1);
1093: a2 = _mm256_loadu_pd(&(aa + bs2 * adiag[i])[8]);
1094: w2 = _mm256_fmadd_pd(a2, v0, w2);
1096: /* second row */
1097: v1 = _mm256_set1_pd(ls[1]);
1098: a3 = _mm256_loadu_pd(&(aa + bs2 * adiag[i])[9]);
1099: w0 = _mm256_fmadd_pd(a3, v1, w0);
1100: a4 = _mm256_loadu_pd(&(aa + bs2 * adiag[i])[13]);
1101: w1 = _mm256_fmadd_pd(a4, v1, w1);
1102: a5 = _mm256_loadu_pd(&(aa + bs2 * adiag[i])[17]);
1103: w2 = _mm256_fmadd_pd(a5, v1, w2);
1105: /* third row */
1106: v2 = _mm256_set1_pd(ls[2]);
1107: a0 = _mm256_loadu_pd(&(aa + bs2 * adiag[i])[18]);
1108: w0 = _mm256_fmadd_pd(a0, v2, w0);
1109: a1 = _mm256_loadu_pd(&(aa + bs2 * adiag[i])[22]);
1110: w1 = _mm256_fmadd_pd(a1, v2, w1);
1111: a2 = _mm256_loadu_pd(&(aa + bs2 * adiag[i])[26]);
1112: w2 = _mm256_fmadd_pd(a2, v2, w2);
1114: /* fourth row */
1115: v3 = _mm256_set1_pd(ls[3]);
1116: a3 = _mm256_loadu_pd(&(aa + bs2 * adiag[i])[27]);
1117: w0 = _mm256_fmadd_pd(a3, v3, w0);
1118: a4 = _mm256_loadu_pd(&(aa + bs2 * adiag[i])[31]);
1119: w1 = _mm256_fmadd_pd(a4, v3, w1);
1120: a5 = _mm256_loadu_pd(&(aa + bs2 * adiag[i])[35]);
1121: w2 = _mm256_fmadd_pd(a5, v3, w2);
1123: /* fifth row */
1124: v0 = _mm256_set1_pd(ls[4]);
1125: a0 = _mm256_loadu_pd(&(aa + bs2 * adiag[i])[36]);
1126: w0 = _mm256_fmadd_pd(a0, v0, w0);
1127: a1 = _mm256_loadu_pd(&(aa + bs2 * adiag[i])[40]);
1128: w1 = _mm256_fmadd_pd(a1, v0, w1);
1129: a2 = _mm256_loadu_pd(&(aa + bs2 * adiag[i])[44]);
1130: w2 = _mm256_fmadd_pd(a2, v0, w2);
1132: /* sixth row */
1133: v1 = _mm256_set1_pd(ls[5]);
1134: a3 = _mm256_loadu_pd(&(aa + bs2 * adiag[i])[45]);
1135: w0 = _mm256_fmadd_pd(a3, v1, w0);
1136: a4 = _mm256_loadu_pd(&(aa + bs2 * adiag[i])[49]);
1137: w1 = _mm256_fmadd_pd(a4, v1, w1);
1138: a5 = _mm256_loadu_pd(&(aa + bs2 * adiag[i])[53]);
1139: w2 = _mm256_fmadd_pd(a5, v1, w2);
1141: /* seventh row */
1142: v2 = _mm256_set1_pd(ls[6]);
1143: a0 = _mm256_loadu_pd(&(aa + bs2 * adiag[i])[54]);
1144: w0 = _mm256_fmadd_pd(a0, v2, w0);
1145: a1 = _mm256_loadu_pd(&(aa + bs2 * adiag[i])[58]);
1146: w1 = _mm256_fmadd_pd(a1, v2, w1);
1147: a2 = _mm256_loadu_pd(&(aa + bs2 * adiag[i])[62]);
1148: w2 = _mm256_fmadd_pd(a2, v2, w2);
1150: /* eighth row */
1151: v3 = _mm256_set1_pd(ls[7]);
1152: a3 = _mm256_loadu_pd(&(aa + bs2 * adiag[i])[63]);
1153: w0 = _mm256_fmadd_pd(a3, v3, w0);
1154: a4 = _mm256_loadu_pd(&(aa + bs2 * adiag[i])[67]);
1155: w1 = _mm256_fmadd_pd(a4, v3, w1);
1156: a5 = _mm256_loadu_pd(&(aa + bs2 * adiag[i])[71]);
1157: w2 = _mm256_fmadd_pd(a5, v3, w2);
1159: /* ninth row */
1160: v0 = _mm256_set1_pd(ls[8]);
1161: a3 = _mm256_loadu_pd(&(aa + bs2 * adiag[i])[72]);
1162: w0 = _mm256_fmadd_pd(a3, v0, w0);
1163: a4 = _mm256_loadu_pd(&(aa + bs2 * adiag[i])[76]);
1164: w1 = _mm256_fmadd_pd(a4, v0, w1);
1165: a2 = _mm256_maskload_pd(&(aa + bs2 * adiag[i])[80], _mm256_set_epi64x(0LL, 0LL, 0LL, 1LL << 63));
1166: w2 = _mm256_fmadd_pd(a2, v0, w2);
1168: _mm256_storeu_pd(&(t + i * bs)[0], w0);
1169: _mm256_storeu_pd(&(t + i * bs)[4], w1);
1170: _mm256_maskstore_pd(&(t + i * bs)[8], _mm256_set_epi64x(0LL, 0LL, 0LL, 1LL << 63), w2);
1172: PetscCall(PetscArraycpy(x + i * bs, t + i * bs, bs));
1173: }
1175: PetscCall(VecRestoreArrayRead(bb, &b));
1176: PetscCall(VecRestoreArray(xx, &x));
1177: PetscCall(PetscLogFlops(2.0 * (a->bs2) * (a->nz) - A->rmap->bs * A->cmap->n));
1178: PetscFunctionReturn(PETSC_SUCCESS);
1179: }
1180: #endif