1 /*
2 Factorization code for BAIJ format.
3 */
4 #include <../src/mat/impls/baij/seq/baij.h>
5 #include <petsc/private/kernels/blockinvert.h>
6 /*
7 Version for when blocks are 7 by 7
8 */
MatILUFactorNumeric_SeqBAIJ_7_inplace(Mat C,Mat A,const MatFactorInfo * info)9 PetscErrorCode MatILUFactorNumeric_SeqBAIJ_7_inplace(Mat C, Mat A, const MatFactorInfo *info)
10 {
11 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)C->data;
12 IS isrow = b->row, isicol = b->icol;
13 const PetscInt *r, *ic, *bi = b->i, *bj = b->j, *ajtmp, *ai = a->i, *aj = a->j, *pj, *ajtmpold;
14 const PetscInt *diag_offset;
15 PetscInt i, j, n = a->mbs, nz, row, idx;
16 MatScalar *pv, *v, *rtmp, *pc, *w, *x;
17 MatScalar p1, p2, p3, p4, m1, m2, m3, m4, m5, m6, m7, m8, m9, x1, x2, x3, x4;
18 MatScalar p5, p6, p7, p8, p9, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15, x16;
19 MatScalar x17, x18, x19, x20, x21, x22, x23, x24, x25, p10, p11, p12, p13, p14;
20 MatScalar p15, p16, p17, p18, p19, p20, p21, p22, p23, p24, p25, m10, m11, m12;
21 MatScalar m13, m14, m15, m16, m17, m18, m19, m20, m21, m22, m23, m24, m25;
22 MatScalar p26, p27, p28, p29, p30, p31, p32, p33, p34, p35, p36;
23 MatScalar p37, p38, p39, p40, p41, p42, p43, p44, p45, p46, p47, p48, p49;
24 MatScalar x26, x27, x28, x29, x30, x31, x32, x33, x34, x35, x36;
25 MatScalar x37, x38, x39, x40, x41, x42, x43, x44, x45, x46, x47, x48, x49;
26 MatScalar m26, m27, m28, m29, m30, m31, m32, m33, m34, m35, m36;
27 MatScalar m37, m38, m39, m40, m41, m42, m43, m44, m45, m46, m47, m48, m49;
28 MatScalar *ba = b->a, *aa = a->a;
29 PetscReal shift = info->shiftamount;
30 PetscBool allowzeropivot, zeropivotdetected;
31
32 PetscFunctionBegin;
33 /* 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 */
34 A->factortype = MAT_FACTOR_NONE;
35 PetscCall(MatGetDiagonalMarkers_SeqBAIJ(A, &diag_offset, NULL));
36 A->factortype = MAT_FACTOR_ILU;
37 allowzeropivot = PetscNot(A->erroriffailure);
38 PetscCall(ISGetIndices(isrow, &r));
39 PetscCall(ISGetIndices(isicol, &ic));
40 PetscCall(PetscMalloc1(49 * (n + 1), &rtmp));
41
42 for (i = 0; i < n; i++) {
43 nz = bi[i + 1] - bi[i];
44 ajtmp = bj + bi[i];
45 for (j = 0; j < nz; j++) {
46 x = rtmp + 49 * ajtmp[j];
47 x[0] = x[1] = x[2] = x[3] = x[4] = x[5] = x[6] = x[7] = x[8] = x[9] = 0.0;
48 x[10] = x[11] = x[12] = x[13] = x[14] = x[15] = x[16] = x[17] = 0.0;
49 x[18] = x[19] = x[20] = x[21] = x[22] = x[23] = x[24] = x[25] = 0.0;
50 x[26] = x[27] = x[28] = x[29] = x[30] = x[31] = x[32] = x[33] = 0.0;
51 x[34] = x[35] = x[36] = x[37] = x[38] = x[39] = x[40] = x[41] = 0.0;
52 x[42] = x[43] = x[44] = x[45] = x[46] = x[47] = x[48] = 0.0;
53 }
54 /* load in initial (unfactored row) */
55 idx = r[i];
56 nz = ai[idx + 1] - ai[idx];
57 ajtmpold = aj + ai[idx];
58 v = aa + 49 * ai[idx];
59 for (j = 0; j < nz; j++) {
60 x = rtmp + 49 * ic[ajtmpold[j]];
61 x[0] = v[0];
62 x[1] = v[1];
63 x[2] = v[2];
64 x[3] = v[3];
65 x[4] = v[4];
66 x[5] = v[5];
67 x[6] = v[6];
68 x[7] = v[7];
69 x[8] = v[8];
70 x[9] = v[9];
71 x[10] = v[10];
72 x[11] = v[11];
73 x[12] = v[12];
74 x[13] = v[13];
75 x[14] = v[14];
76 x[15] = v[15];
77 x[16] = v[16];
78 x[17] = v[17];
79 x[18] = v[18];
80 x[19] = v[19];
81 x[20] = v[20];
82 x[21] = v[21];
83 x[22] = v[22];
84 x[23] = v[23];
85 x[24] = v[24];
86 x[25] = v[25];
87 x[26] = v[26];
88 x[27] = v[27];
89 x[28] = v[28];
90 x[29] = v[29];
91 x[30] = v[30];
92 x[31] = v[31];
93 x[32] = v[32];
94 x[33] = v[33];
95 x[34] = v[34];
96 x[35] = v[35];
97 x[36] = v[36];
98 x[37] = v[37];
99 x[38] = v[38];
100 x[39] = v[39];
101 x[40] = v[40];
102 x[41] = v[41];
103 x[42] = v[42];
104 x[43] = v[43];
105 x[44] = v[44];
106 x[45] = v[45];
107 x[46] = v[46];
108 x[47] = v[47];
109 x[48] = v[48];
110 v += 49;
111 }
112 row = *ajtmp++;
113 while (row < i) {
114 pc = rtmp + 49 * row;
115 p1 = pc[0];
116 p2 = pc[1];
117 p3 = pc[2];
118 p4 = pc[3];
119 p5 = pc[4];
120 p6 = pc[5];
121 p7 = pc[6];
122 p8 = pc[7];
123 p9 = pc[8];
124 p10 = pc[9];
125 p11 = pc[10];
126 p12 = pc[11];
127 p13 = pc[12];
128 p14 = pc[13];
129 p15 = pc[14];
130 p16 = pc[15];
131 p17 = pc[16];
132 p18 = pc[17];
133 p19 = pc[18];
134 p20 = pc[19];
135 p21 = pc[20];
136 p22 = pc[21];
137 p23 = pc[22];
138 p24 = pc[23];
139 p25 = pc[24];
140 p26 = pc[25];
141 p27 = pc[26];
142 p28 = pc[27];
143 p29 = pc[28];
144 p30 = pc[29];
145 p31 = pc[30];
146 p32 = pc[31];
147 p33 = pc[32];
148 p34 = pc[33];
149 p35 = pc[34];
150 p36 = pc[35];
151 p37 = pc[36];
152 p38 = pc[37];
153 p39 = pc[38];
154 p40 = pc[39];
155 p41 = pc[40];
156 p42 = pc[41];
157 p43 = pc[42];
158 p44 = pc[43];
159 p45 = pc[44];
160 p46 = pc[45];
161 p47 = pc[46];
162 p48 = pc[47];
163 p49 = pc[48];
164 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 || p26 != 0.0 || p27 != 0.0 || p28 != 0.0 || p29 != 0.0 || p30 != 0.0 || p31 != 0.0 || p32 != 0.0 || p33 != 0.0 || p34 != 0.0 || p35 != 0.0 || p36 != 0.0 || p37 != 0.0 || p38 != 0.0 || p39 != 0.0 || p40 != 0.0 || p41 != 0.0 || p42 != 0.0 || p43 != 0.0 || p44 != 0.0 || p45 != 0.0 || p46 != 0.0 || p47 != 0.0 || p48 != 0.0 || p49 != 0.0) {
165 pv = ba + 49 * diag_offset[row];
166 pj = bj + diag_offset[row] + 1;
167 x1 = pv[0];
168 x2 = pv[1];
169 x3 = pv[2];
170 x4 = pv[3];
171 x5 = pv[4];
172 x6 = pv[5];
173 x7 = pv[6];
174 x8 = pv[7];
175 x9 = pv[8];
176 x10 = pv[9];
177 x11 = pv[10];
178 x12 = pv[11];
179 x13 = pv[12];
180 x14 = pv[13];
181 x15 = pv[14];
182 x16 = pv[15];
183 x17 = pv[16];
184 x18 = pv[17];
185 x19 = pv[18];
186 x20 = pv[19];
187 x21 = pv[20];
188 x22 = pv[21];
189 x23 = pv[22];
190 x24 = pv[23];
191 x25 = pv[24];
192 x26 = pv[25];
193 x27 = pv[26];
194 x28 = pv[27];
195 x29 = pv[28];
196 x30 = pv[29];
197 x31 = pv[30];
198 x32 = pv[31];
199 x33 = pv[32];
200 x34 = pv[33];
201 x35 = pv[34];
202 x36 = pv[35];
203 x37 = pv[36];
204 x38 = pv[37];
205 x39 = pv[38];
206 x40 = pv[39];
207 x41 = pv[40];
208 x42 = pv[41];
209 x43 = pv[42];
210 x44 = pv[43];
211 x45 = pv[44];
212 x46 = pv[45];
213 x47 = pv[46];
214 x48 = pv[47];
215 x49 = pv[48];
216 pc[0] = m1 = p1 * x1 + p8 * x2 + p15 * x3 + p22 * x4 + p29 * x5 + p36 * x6 + p43 * x7;
217 pc[1] = m2 = p2 * x1 + p9 * x2 + p16 * x3 + p23 * x4 + p30 * x5 + p37 * x6 + p44 * x7;
218 pc[2] = m3 = p3 * x1 + p10 * x2 + p17 * x3 + p24 * x4 + p31 * x5 + p38 * x6 + p45 * x7;
219 pc[3] = m4 = p4 * x1 + p11 * x2 + p18 * x3 + p25 * x4 + p32 * x5 + p39 * x6 + p46 * x7;
220 pc[4] = m5 = p5 * x1 + p12 * x2 + p19 * x3 + p26 * x4 + p33 * x5 + p40 * x6 + p47 * x7;
221 pc[5] = m6 = p6 * x1 + p13 * x2 + p20 * x3 + p27 * x4 + p34 * x5 + p41 * x6 + p48 * x7;
222 pc[6] = m7 = p7 * x1 + p14 * x2 + p21 * x3 + p28 * x4 + p35 * x5 + p42 * x6 + p49 * x7;
223
224 pc[7] = m8 = p1 * x8 + p8 * x9 + p15 * x10 + p22 * x11 + p29 * x12 + p36 * x13 + p43 * x14;
225 pc[8] = m9 = p2 * x8 + p9 * x9 + p16 * x10 + p23 * x11 + p30 * x12 + p37 * x13 + p44 * x14;
226 pc[9] = m10 = p3 * x8 + p10 * x9 + p17 * x10 + p24 * x11 + p31 * x12 + p38 * x13 + p45 * x14;
227 pc[10] = m11 = p4 * x8 + p11 * x9 + p18 * x10 + p25 * x11 + p32 * x12 + p39 * x13 + p46 * x14;
228 pc[11] = m12 = p5 * x8 + p12 * x9 + p19 * x10 + p26 * x11 + p33 * x12 + p40 * x13 + p47 * x14;
229 pc[12] = m13 = p6 * x8 + p13 * x9 + p20 * x10 + p27 * x11 + p34 * x12 + p41 * x13 + p48 * x14;
230 pc[13] = m14 = p7 * x8 + p14 * x9 + p21 * x10 + p28 * x11 + p35 * x12 + p42 * x13 + p49 * x14;
231
232 pc[14] = m15 = p1 * x15 + p8 * x16 + p15 * x17 + p22 * x18 + p29 * x19 + p36 * x20 + p43 * x21;
233 pc[15] = m16 = p2 * x15 + p9 * x16 + p16 * x17 + p23 * x18 + p30 * x19 + p37 * x20 + p44 * x21;
234 pc[16] = m17 = p3 * x15 + p10 * x16 + p17 * x17 + p24 * x18 + p31 * x19 + p38 * x20 + p45 * x21;
235 pc[17] = m18 = p4 * x15 + p11 * x16 + p18 * x17 + p25 * x18 + p32 * x19 + p39 * x20 + p46 * x21;
236 pc[18] = m19 = p5 * x15 + p12 * x16 + p19 * x17 + p26 * x18 + p33 * x19 + p40 * x20 + p47 * x21;
237 pc[19] = m20 = p6 * x15 + p13 * x16 + p20 * x17 + p27 * x18 + p34 * x19 + p41 * x20 + p48 * x21;
238 pc[20] = m21 = p7 * x15 + p14 * x16 + p21 * x17 + p28 * x18 + p35 * x19 + p42 * x20 + p49 * x21;
239
240 pc[21] = m22 = p1 * x22 + p8 * x23 + p15 * x24 + p22 * x25 + p29 * x26 + p36 * x27 + p43 * x28;
241 pc[22] = m23 = p2 * x22 + p9 * x23 + p16 * x24 + p23 * x25 + p30 * x26 + p37 * x27 + p44 * x28;
242 pc[23] = m24 = p3 * x22 + p10 * x23 + p17 * x24 + p24 * x25 + p31 * x26 + p38 * x27 + p45 * x28;
243 pc[24] = m25 = p4 * x22 + p11 * x23 + p18 * x24 + p25 * x25 + p32 * x26 + p39 * x27 + p46 * x28;
244 pc[25] = m26 = p5 * x22 + p12 * x23 + p19 * x24 + p26 * x25 + p33 * x26 + p40 * x27 + p47 * x28;
245 pc[26] = m27 = p6 * x22 + p13 * x23 + p20 * x24 + p27 * x25 + p34 * x26 + p41 * x27 + p48 * x28;
246 pc[27] = m28 = p7 * x22 + p14 * x23 + p21 * x24 + p28 * x25 + p35 * x26 + p42 * x27 + p49 * x28;
247
248 pc[28] = m29 = p1 * x29 + p8 * x30 + p15 * x31 + p22 * x32 + p29 * x33 + p36 * x34 + p43 * x35;
249 pc[29] = m30 = p2 * x29 + p9 * x30 + p16 * x31 + p23 * x32 + p30 * x33 + p37 * x34 + p44 * x35;
250 pc[30] = m31 = p3 * x29 + p10 * x30 + p17 * x31 + p24 * x32 + p31 * x33 + p38 * x34 + p45 * x35;
251 pc[31] = m32 = p4 * x29 + p11 * x30 + p18 * x31 + p25 * x32 + p32 * x33 + p39 * x34 + p46 * x35;
252 pc[32] = m33 = p5 * x29 + p12 * x30 + p19 * x31 + p26 * x32 + p33 * x33 + p40 * x34 + p47 * x35;
253 pc[33] = m34 = p6 * x29 + p13 * x30 + p20 * x31 + p27 * x32 + p34 * x33 + p41 * x34 + p48 * x35;
254 pc[34] = m35 = p7 * x29 + p14 * x30 + p21 * x31 + p28 * x32 + p35 * x33 + p42 * x34 + p49 * x35;
255
256 pc[35] = m36 = p1 * x36 + p8 * x37 + p15 * x38 + p22 * x39 + p29 * x40 + p36 * x41 + p43 * x42;
257 pc[36] = m37 = p2 * x36 + p9 * x37 + p16 * x38 + p23 * x39 + p30 * x40 + p37 * x41 + p44 * x42;
258 pc[37] = m38 = p3 * x36 + p10 * x37 + p17 * x38 + p24 * x39 + p31 * x40 + p38 * x41 + p45 * x42;
259 pc[38] = m39 = p4 * x36 + p11 * x37 + p18 * x38 + p25 * x39 + p32 * x40 + p39 * x41 + p46 * x42;
260 pc[39] = m40 = p5 * x36 + p12 * x37 + p19 * x38 + p26 * x39 + p33 * x40 + p40 * x41 + p47 * x42;
261 pc[40] = m41 = p6 * x36 + p13 * x37 + p20 * x38 + p27 * x39 + p34 * x40 + p41 * x41 + p48 * x42;
262 pc[41] = m42 = p7 * x36 + p14 * x37 + p21 * x38 + p28 * x39 + p35 * x40 + p42 * x41 + p49 * x42;
263
264 pc[42] = m43 = p1 * x43 + p8 * x44 + p15 * x45 + p22 * x46 + p29 * x47 + p36 * x48 + p43 * x49;
265 pc[43] = m44 = p2 * x43 + p9 * x44 + p16 * x45 + p23 * x46 + p30 * x47 + p37 * x48 + p44 * x49;
266 pc[44] = m45 = p3 * x43 + p10 * x44 + p17 * x45 + p24 * x46 + p31 * x47 + p38 * x48 + p45 * x49;
267 pc[45] = m46 = p4 * x43 + p11 * x44 + p18 * x45 + p25 * x46 + p32 * x47 + p39 * x48 + p46 * x49;
268 pc[46] = m47 = p5 * x43 + p12 * x44 + p19 * x45 + p26 * x46 + p33 * x47 + p40 * x48 + p47 * x49;
269 pc[47] = m48 = p6 * x43 + p13 * x44 + p20 * x45 + p27 * x46 + p34 * x47 + p41 * x48 + p48 * x49;
270 pc[48] = m49 = p7 * x43 + p14 * x44 + p21 * x45 + p28 * x46 + p35 * x47 + p42 * x48 + p49 * x49;
271
272 nz = bi[row + 1] - diag_offset[row] - 1;
273 pv += 49;
274 for (j = 0; j < nz; j++) {
275 x1 = pv[0];
276 x2 = pv[1];
277 x3 = pv[2];
278 x4 = pv[3];
279 x5 = pv[4];
280 x6 = pv[5];
281 x7 = pv[6];
282 x8 = pv[7];
283 x9 = pv[8];
284 x10 = pv[9];
285 x11 = pv[10];
286 x12 = pv[11];
287 x13 = pv[12];
288 x14 = pv[13];
289 x15 = pv[14];
290 x16 = pv[15];
291 x17 = pv[16];
292 x18 = pv[17];
293 x19 = pv[18];
294 x20 = pv[19];
295 x21 = pv[20];
296 x22 = pv[21];
297 x23 = pv[22];
298 x24 = pv[23];
299 x25 = pv[24];
300 x26 = pv[25];
301 x27 = pv[26];
302 x28 = pv[27];
303 x29 = pv[28];
304 x30 = pv[29];
305 x31 = pv[30];
306 x32 = pv[31];
307 x33 = pv[32];
308 x34 = pv[33];
309 x35 = pv[34];
310 x36 = pv[35];
311 x37 = pv[36];
312 x38 = pv[37];
313 x39 = pv[38];
314 x40 = pv[39];
315 x41 = pv[40];
316 x42 = pv[41];
317 x43 = pv[42];
318 x44 = pv[43];
319 x45 = pv[44];
320 x46 = pv[45];
321 x47 = pv[46];
322 x48 = pv[47];
323 x49 = pv[48];
324 x = rtmp + 49 * pj[j];
325 x[0] -= m1 * x1 + m8 * x2 + m15 * x3 + m22 * x4 + m29 * x5 + m36 * x6 + m43 * x7;
326 x[1] -= m2 * x1 + m9 * x2 + m16 * x3 + m23 * x4 + m30 * x5 + m37 * x6 + m44 * x7;
327 x[2] -= m3 * x1 + m10 * x2 + m17 * x3 + m24 * x4 + m31 * x5 + m38 * x6 + m45 * x7;
328 x[3] -= m4 * x1 + m11 * x2 + m18 * x3 + m25 * x4 + m32 * x5 + m39 * x6 + m46 * x7;
329 x[4] -= m5 * x1 + m12 * x2 + m19 * x3 + m26 * x4 + m33 * x5 + m40 * x6 + m47 * x7;
330 x[5] -= m6 * x1 + m13 * x2 + m20 * x3 + m27 * x4 + m34 * x5 + m41 * x6 + m48 * x7;
331 x[6] -= m7 * x1 + m14 * x2 + m21 * x3 + m28 * x4 + m35 * x5 + m42 * x6 + m49 * x7;
332
333 x[7] -= m1 * x8 + m8 * x9 + m15 * x10 + m22 * x11 + m29 * x12 + m36 * x13 + m43 * x14;
334 x[8] -= m2 * x8 + m9 * x9 + m16 * x10 + m23 * x11 + m30 * x12 + m37 * x13 + m44 * x14;
335 x[9] -= m3 * x8 + m10 * x9 + m17 * x10 + m24 * x11 + m31 * x12 + m38 * x13 + m45 * x14;
336 x[10] -= m4 * x8 + m11 * x9 + m18 * x10 + m25 * x11 + m32 * x12 + m39 * x13 + m46 * x14;
337 x[11] -= m5 * x8 + m12 * x9 + m19 * x10 + m26 * x11 + m33 * x12 + m40 * x13 + m47 * x14;
338 x[12] -= m6 * x8 + m13 * x9 + m20 * x10 + m27 * x11 + m34 * x12 + m41 * x13 + m48 * x14;
339 x[13] -= m7 * x8 + m14 * x9 + m21 * x10 + m28 * x11 + m35 * x12 + m42 * x13 + m49 * x14;
340
341 x[14] -= m1 * x15 + m8 * x16 + m15 * x17 + m22 * x18 + m29 * x19 + m36 * x20 + m43 * x21;
342 x[15] -= m2 * x15 + m9 * x16 + m16 * x17 + m23 * x18 + m30 * x19 + m37 * x20 + m44 * x21;
343 x[16] -= m3 * x15 + m10 * x16 + m17 * x17 + m24 * x18 + m31 * x19 + m38 * x20 + m45 * x21;
344 x[17] -= m4 * x15 + m11 * x16 + m18 * x17 + m25 * x18 + m32 * x19 + m39 * x20 + m46 * x21;
345 x[18] -= m5 * x15 + m12 * x16 + m19 * x17 + m26 * x18 + m33 * x19 + m40 * x20 + m47 * x21;
346 x[19] -= m6 * x15 + m13 * x16 + m20 * x17 + m27 * x18 + m34 * x19 + m41 * x20 + m48 * x21;
347 x[20] -= m7 * x15 + m14 * x16 + m21 * x17 + m28 * x18 + m35 * x19 + m42 * x20 + m49 * x21;
348
349 x[21] -= m1 * x22 + m8 * x23 + m15 * x24 + m22 * x25 + m29 * x26 + m36 * x27 + m43 * x28;
350 x[22] -= m2 * x22 + m9 * x23 + m16 * x24 + m23 * x25 + m30 * x26 + m37 * x27 + m44 * x28;
351 x[23] -= m3 * x22 + m10 * x23 + m17 * x24 + m24 * x25 + m31 * x26 + m38 * x27 + m45 * x28;
352 x[24] -= m4 * x22 + m11 * x23 + m18 * x24 + m25 * x25 + m32 * x26 + m39 * x27 + m46 * x28;
353 x[25] -= m5 * x22 + m12 * x23 + m19 * x24 + m26 * x25 + m33 * x26 + m40 * x27 + m47 * x28;
354 x[26] -= m6 * x22 + m13 * x23 + m20 * x24 + m27 * x25 + m34 * x26 + m41 * x27 + m48 * x28;
355 x[27] -= m7 * x22 + m14 * x23 + m21 * x24 + m28 * x25 + m35 * x26 + m42 * x27 + m49 * x28;
356
357 x[28] -= m1 * x29 + m8 * x30 + m15 * x31 + m22 * x32 + m29 * x33 + m36 * x34 + m43 * x35;
358 x[29] -= m2 * x29 + m9 * x30 + m16 * x31 + m23 * x32 + m30 * x33 + m37 * x34 + m44 * x35;
359 x[30] -= m3 * x29 + m10 * x30 + m17 * x31 + m24 * x32 + m31 * x33 + m38 * x34 + m45 * x35;
360 x[31] -= m4 * x29 + m11 * x30 + m18 * x31 + m25 * x32 + m32 * x33 + m39 * x34 + m46 * x35;
361 x[32] -= m5 * x29 + m12 * x30 + m19 * x31 + m26 * x32 + m33 * x33 + m40 * x34 + m47 * x35;
362 x[33] -= m6 * x29 + m13 * x30 + m20 * x31 + m27 * x32 + m34 * x33 + m41 * x34 + m48 * x35;
363 x[34] -= m7 * x29 + m14 * x30 + m21 * x31 + m28 * x32 + m35 * x33 + m42 * x34 + m49 * x35;
364
365 x[35] -= m1 * x36 + m8 * x37 + m15 * x38 + m22 * x39 + m29 * x40 + m36 * x41 + m43 * x42;
366 x[36] -= m2 * x36 + m9 * x37 + m16 * x38 + m23 * x39 + m30 * x40 + m37 * x41 + m44 * x42;
367 x[37] -= m3 * x36 + m10 * x37 + m17 * x38 + m24 * x39 + m31 * x40 + m38 * x41 + m45 * x42;
368 x[38] -= m4 * x36 + m11 * x37 + m18 * x38 + m25 * x39 + m32 * x40 + m39 * x41 + m46 * x42;
369 x[39] -= m5 * x36 + m12 * x37 + m19 * x38 + m26 * x39 + m33 * x40 + m40 * x41 + m47 * x42;
370 x[40] -= m6 * x36 + m13 * x37 + m20 * x38 + m27 * x39 + m34 * x40 + m41 * x41 + m48 * x42;
371 x[41] -= m7 * x36 + m14 * x37 + m21 * x38 + m28 * x39 + m35 * x40 + m42 * x41 + m49 * x42;
372
373 x[42] -= m1 * x43 + m8 * x44 + m15 * x45 + m22 * x46 + m29 * x47 + m36 * x48 + m43 * x49;
374 x[43] -= m2 * x43 + m9 * x44 + m16 * x45 + m23 * x46 + m30 * x47 + m37 * x48 + m44 * x49;
375 x[44] -= m3 * x43 + m10 * x44 + m17 * x45 + m24 * x46 + m31 * x47 + m38 * x48 + m45 * x49;
376 x[45] -= m4 * x43 + m11 * x44 + m18 * x45 + m25 * x46 + m32 * x47 + m39 * x48 + m46 * x49;
377 x[46] -= m5 * x43 + m12 * x44 + m19 * x45 + m26 * x46 + m33 * x47 + m40 * x48 + m47 * x49;
378 x[47] -= m6 * x43 + m13 * x44 + m20 * x45 + m27 * x46 + m34 * x47 + m41 * x48 + m48 * x49;
379 x[48] -= m7 * x43 + m14 * x44 + m21 * x45 + m28 * x46 + m35 * x47 + m42 * x48 + m49 * x49;
380 pv += 49;
381 }
382 PetscCall(PetscLogFlops(686.0 * nz + 637.0));
383 }
384 row = *ajtmp++;
385 }
386 /* finished row so stick it into b->a */
387 pv = ba + 49 * bi[i];
388 pj = bj + bi[i];
389 nz = bi[i + 1] - bi[i];
390 for (j = 0; j < nz; j++) {
391 x = rtmp + 49 * pj[j];
392 pv[0] = x[0];
393 pv[1] = x[1];
394 pv[2] = x[2];
395 pv[3] = x[3];
396 pv[4] = x[4];
397 pv[5] = x[5];
398 pv[6] = x[6];
399 pv[7] = x[7];
400 pv[8] = x[8];
401 pv[9] = x[9];
402 pv[10] = x[10];
403 pv[11] = x[11];
404 pv[12] = x[12];
405 pv[13] = x[13];
406 pv[14] = x[14];
407 pv[15] = x[15];
408 pv[16] = x[16];
409 pv[17] = x[17];
410 pv[18] = x[18];
411 pv[19] = x[19];
412 pv[20] = x[20];
413 pv[21] = x[21];
414 pv[22] = x[22];
415 pv[23] = x[23];
416 pv[24] = x[24];
417 pv[25] = x[25];
418 pv[26] = x[26];
419 pv[27] = x[27];
420 pv[28] = x[28];
421 pv[29] = x[29];
422 pv[30] = x[30];
423 pv[31] = x[31];
424 pv[32] = x[32];
425 pv[33] = x[33];
426 pv[34] = x[34];
427 pv[35] = x[35];
428 pv[36] = x[36];
429 pv[37] = x[37];
430 pv[38] = x[38];
431 pv[39] = x[39];
432 pv[40] = x[40];
433 pv[41] = x[41];
434 pv[42] = x[42];
435 pv[43] = x[43];
436 pv[44] = x[44];
437 pv[45] = x[45];
438 pv[46] = x[46];
439 pv[47] = x[47];
440 pv[48] = x[48];
441 pv += 49;
442 }
443 /* invert diagonal block */
444 w = ba + 49 * diag_offset[i];
445 PetscCall(PetscKernel_A_gets_inverse_A_7(w, shift, allowzeropivot, &zeropivotdetected));
446 if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
447 }
448
449 PetscCall(PetscFree(rtmp));
450 PetscCall(ISRestoreIndices(isicol, &ic));
451 PetscCall(ISRestoreIndices(isrow, &r));
452
453 C->ops->solve = MatSolve_SeqBAIJ_7_inplace;
454 C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_7_inplace;
455 C->assembled = PETSC_TRUE;
456
457 PetscCall(PetscLogFlops(1.333333333333 * 7 * 7 * 7 * b->mbs)); /* from inverting diagonal blocks */
458 PetscFunctionReturn(PETSC_SUCCESS);
459 }
460
MatLUFactorNumeric_SeqBAIJ_7(Mat B,Mat A,const MatFactorInfo * info)461 PetscErrorCode MatLUFactorNumeric_SeqBAIJ_7(Mat B, Mat A, const MatFactorInfo *info)
462 {
463 Mat C = B;
464 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)C->data;
465 IS isrow = b->row, isicol = b->icol;
466 const PetscInt *r, *ic;
467 PetscInt i, j, k, nz, nzL, row;
468 const PetscInt n = a->mbs, *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j;
469 const PetscInt *ajtmp, *bjtmp, *bdiag = b->diag, *pj, bs2 = a->bs2;
470 MatScalar *rtmp, *pc, *mwork, *v, *pv, *aa = a->a;
471 PetscInt flg;
472 PetscReal shift = info->shiftamount;
473 PetscBool allowzeropivot, zeropivotdetected;
474
475 PetscFunctionBegin;
476 allowzeropivot = PetscNot(A->erroriffailure);
477 PetscCall(ISGetIndices(isrow, &r));
478 PetscCall(ISGetIndices(isicol, &ic));
479
480 /* generate work space needed by the factorization */
481 PetscCall(PetscMalloc2(bs2 * n, &rtmp, bs2, &mwork));
482 PetscCall(PetscArrayzero(rtmp, bs2 * n));
483
484 for (i = 0; i < n; i++) {
485 /* zero rtmp */
486 /* L part */
487 nz = bi[i + 1] - bi[i];
488 bjtmp = bj + bi[i];
489 for (j = 0; j < nz; j++) PetscCall(PetscArrayzero(rtmp + bs2 * bjtmp[j], bs2));
490
491 /* U part */
492 nz = bdiag[i] - bdiag[i + 1];
493 bjtmp = bj + bdiag[i + 1] + 1;
494 for (j = 0; j < nz; j++) PetscCall(PetscArrayzero(rtmp + bs2 * bjtmp[j], bs2));
495
496 /* load in initial (unfactored row) */
497 nz = ai[r[i] + 1] - ai[r[i]];
498 ajtmp = aj + ai[r[i]];
499 v = aa + bs2 * ai[r[i]];
500 for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(rtmp + bs2 * ic[ajtmp[j]], v + bs2 * j, bs2));
501
502 /* elimination */
503 bjtmp = bj + bi[i];
504 nzL = bi[i + 1] - bi[i];
505 for (k = 0; k < nzL; k++) {
506 row = bjtmp[k];
507 pc = rtmp + bs2 * row;
508 for (flg = 0, j = 0; j < bs2; j++) {
509 if (pc[j] != 0.0) {
510 flg = 1;
511 break;
512 }
513 }
514 if (flg) {
515 pv = b->a + bs2 * bdiag[row];
516 /* PetscKernel_A_gets_A_times_B(bs,pc,pv,mwork); *pc = *pc * (*pv); */
517 PetscCall(PetscKernel_A_gets_A_times_B_7(pc, pv, mwork));
518
519 pj = b->j + bdiag[row + 1] + 1; /* beginning of U(row,:) */
520 pv = b->a + bs2 * (bdiag[row + 1] + 1);
521 nz = bdiag[row] - bdiag[row + 1] - 1; /* num of entries inU(row,:), excluding diag */
522 for (j = 0; j < nz; j++) {
523 /* PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j); */
524 /* rtmp+bs2*pj[j] = rtmp+bs2*pj[j] - (*pc)*(pv+bs2*j) */
525 v = rtmp + bs2 * pj[j];
526 PetscCall(PetscKernel_A_gets_A_minus_B_times_C_7(v, pc, pv));
527 pv += bs2;
528 }
529 PetscCall(PetscLogFlops(686.0 * nz + 637)); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */
530 }
531 }
532
533 /* finished row so stick it into b->a */
534 /* L part */
535 pv = b->a + bs2 * bi[i];
536 pj = b->j + bi[i];
537 nz = bi[i + 1] - bi[i];
538 for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2));
539
540 /* Mark diagonal and invert diagonal for simpler triangular solves */
541 pv = b->a + bs2 * bdiag[i];
542 pj = b->j + bdiag[i];
543 PetscCall(PetscArraycpy(pv, rtmp + bs2 * pj[0], bs2));
544 PetscCall(PetscKernel_A_gets_inverse_A_7(pv, shift, allowzeropivot, &zeropivotdetected));
545 if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
546
547 /* U part */
548 pv = b->a + bs2 * (bdiag[i + 1] + 1);
549 pj = b->j + bdiag[i + 1] + 1;
550 nz = bdiag[i] - bdiag[i + 1] - 1;
551 for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2));
552 }
553
554 PetscCall(PetscFree2(rtmp, mwork));
555 PetscCall(ISRestoreIndices(isicol, &ic));
556 PetscCall(ISRestoreIndices(isrow, &r));
557
558 C->ops->solve = MatSolve_SeqBAIJ_7;
559 C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_7;
560 C->assembled = PETSC_TRUE;
561
562 PetscCall(PetscLogFlops(1.333333333333 * 7 * 7 * 7 * n)); /* from inverting diagonal blocks */
563 PetscFunctionReturn(PETSC_SUCCESS);
564 }
565
MatILUFactorNumeric_SeqBAIJ_7_NaturalOrdering_inplace(Mat C,Mat A,const MatFactorInfo * info)566 PetscErrorCode MatILUFactorNumeric_SeqBAIJ_7_NaturalOrdering_inplace(Mat C, Mat A, const MatFactorInfo *info)
567 {
568 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)C->data;
569 PetscInt i, j, n = a->mbs, *bi = b->i, *bj = b->j;
570 PetscInt *ajtmpold, *ajtmp, nz, row;
571 PetscInt *ai = a->i, *aj = a->j, *pj;
572 const PetscInt *diag_offset;
573 MatScalar *pv, *v, *rtmp, *pc, *w, *x;
574 MatScalar x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15;
575 MatScalar x16, x17, x18, x19, x20, x21, x22, x23, x24, x25;
576 MatScalar p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12, p13, p14, p15;
577 MatScalar p16, p17, p18, p19, p20, p21, p22, p23, p24, p25;
578 MatScalar m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, m11, m12, m13, m14, m15;
579 MatScalar m16, m17, m18, m19, m20, m21, m22, m23, m24, m25;
580 MatScalar p26, p27, p28, p29, p30, p31, p32, p33, p34, p35, p36;
581 MatScalar p37, p38, p39, p40, p41, p42, p43, p44, p45, p46, p47, p48, p49;
582 MatScalar x26, x27, x28, x29, x30, x31, x32, x33, x34, x35, x36;
583 MatScalar x37, x38, x39, x40, x41, x42, x43, x44, x45, x46, x47, x48, x49;
584 MatScalar m26, m27, m28, m29, m30, m31, m32, m33, m34, m35, m36;
585 MatScalar m37, m38, m39, m40, m41, m42, m43, m44, m45, m46, m47, m48, m49;
586 MatScalar *ba = b->a, *aa = a->a;
587 PetscReal shift = info->shiftamount;
588 PetscBool allowzeropivot, zeropivotdetected;
589
590 PetscFunctionBegin;
591 /* 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 */
592 A->factortype = MAT_FACTOR_NONE;
593 PetscCall(MatGetDiagonalMarkers_SeqBAIJ(A, &diag_offset, NULL));
594 A->factortype = MAT_FACTOR_ILU;
595 allowzeropivot = PetscNot(A->erroriffailure);
596 PetscCall(PetscMalloc1(49 * (n + 1), &rtmp));
597 for (i = 0; i < n; i++) {
598 nz = bi[i + 1] - bi[i];
599 ajtmp = bj + bi[i];
600 for (j = 0; j < nz; j++) {
601 x = rtmp + 49 * ajtmp[j];
602 x[0] = x[1] = x[2] = x[3] = x[4] = x[5] = x[6] = x[7] = x[8] = x[9] = 0.0;
603 x[10] = x[11] = x[12] = x[13] = x[14] = x[15] = x[16] = x[17] = 0.0;
604 x[18] = x[19] = x[20] = x[21] = x[22] = x[23] = x[24] = x[25] = 0.0;
605 x[26] = x[27] = x[28] = x[29] = x[30] = x[31] = x[32] = x[33] = 0.0;
606 x[34] = x[35] = x[36] = x[37] = x[38] = x[39] = x[40] = x[41] = 0.0;
607 x[42] = x[43] = x[44] = x[45] = x[46] = x[47] = x[48] = 0.0;
608 }
609 /* load in initial (unfactored row) */
610 nz = ai[i + 1] - ai[i];
611 ajtmpold = aj + ai[i];
612 v = aa + 49 * ai[i];
613 for (j = 0; j < nz; j++) {
614 x = rtmp + 49 * ajtmpold[j];
615 x[0] = v[0];
616 x[1] = v[1];
617 x[2] = v[2];
618 x[3] = v[3];
619 x[4] = v[4];
620 x[5] = v[5];
621 x[6] = v[6];
622 x[7] = v[7];
623 x[8] = v[8];
624 x[9] = v[9];
625 x[10] = v[10];
626 x[11] = v[11];
627 x[12] = v[12];
628 x[13] = v[13];
629 x[14] = v[14];
630 x[15] = v[15];
631 x[16] = v[16];
632 x[17] = v[17];
633 x[18] = v[18];
634 x[19] = v[19];
635 x[20] = v[20];
636 x[21] = v[21];
637 x[22] = v[22];
638 x[23] = v[23];
639 x[24] = v[24];
640 x[25] = v[25];
641 x[26] = v[26];
642 x[27] = v[27];
643 x[28] = v[28];
644 x[29] = v[29];
645 x[30] = v[30];
646 x[31] = v[31];
647 x[32] = v[32];
648 x[33] = v[33];
649 x[34] = v[34];
650 x[35] = v[35];
651 x[36] = v[36];
652 x[37] = v[37];
653 x[38] = v[38];
654 x[39] = v[39];
655 x[40] = v[40];
656 x[41] = v[41];
657 x[42] = v[42];
658 x[43] = v[43];
659 x[44] = v[44];
660 x[45] = v[45];
661 x[46] = v[46];
662 x[47] = v[47];
663 x[48] = v[48];
664 v += 49;
665 }
666 row = *ajtmp++;
667 while (row < i) {
668 pc = rtmp + 49 * row;
669 p1 = pc[0];
670 p2 = pc[1];
671 p3 = pc[2];
672 p4 = pc[3];
673 p5 = pc[4];
674 p6 = pc[5];
675 p7 = pc[6];
676 p8 = pc[7];
677 p9 = pc[8];
678 p10 = pc[9];
679 p11 = pc[10];
680 p12 = pc[11];
681 p13 = pc[12];
682 p14 = pc[13];
683 p15 = pc[14];
684 p16 = pc[15];
685 p17 = pc[16];
686 p18 = pc[17];
687 p19 = pc[18];
688 p20 = pc[19];
689 p21 = pc[20];
690 p22 = pc[21];
691 p23 = pc[22];
692 p24 = pc[23];
693 p25 = pc[24];
694 p26 = pc[25];
695 p27 = pc[26];
696 p28 = pc[27];
697 p29 = pc[28];
698 p30 = pc[29];
699 p31 = pc[30];
700 p32 = pc[31];
701 p33 = pc[32];
702 p34 = pc[33];
703 p35 = pc[34];
704 p36 = pc[35];
705 p37 = pc[36];
706 p38 = pc[37];
707 p39 = pc[38];
708 p40 = pc[39];
709 p41 = pc[40];
710 p42 = pc[41];
711 p43 = pc[42];
712 p44 = pc[43];
713 p45 = pc[44];
714 p46 = pc[45];
715 p47 = pc[46];
716 p48 = pc[47];
717 p49 = pc[48];
718 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 || p26 != 0.0 || p27 != 0.0 || p28 != 0.0 || p29 != 0.0 || p30 != 0.0 || p31 != 0.0 || p32 != 0.0 || p33 != 0.0 || p34 != 0.0 || p35 != 0.0 || p36 != 0.0 || p37 != 0.0 || p38 != 0.0 || p39 != 0.0 || p40 != 0.0 || p41 != 0.0 || p42 != 0.0 || p43 != 0.0 || p44 != 0.0 || p45 != 0.0 || p46 != 0.0 || p47 != 0.0 || p48 != 0.0 || p49 != 0.0) {
719 pv = ba + 49 * diag_offset[row];
720 pj = bj + diag_offset[row] + 1;
721 x1 = pv[0];
722 x2 = pv[1];
723 x3 = pv[2];
724 x4 = pv[3];
725 x5 = pv[4];
726 x6 = pv[5];
727 x7 = pv[6];
728 x8 = pv[7];
729 x9 = pv[8];
730 x10 = pv[9];
731 x11 = pv[10];
732 x12 = pv[11];
733 x13 = pv[12];
734 x14 = pv[13];
735 x15 = pv[14];
736 x16 = pv[15];
737 x17 = pv[16];
738 x18 = pv[17];
739 x19 = pv[18];
740 x20 = pv[19];
741 x21 = pv[20];
742 x22 = pv[21];
743 x23 = pv[22];
744 x24 = pv[23];
745 x25 = pv[24];
746 x26 = pv[25];
747 x27 = pv[26];
748 x28 = pv[27];
749 x29 = pv[28];
750 x30 = pv[29];
751 x31 = pv[30];
752 x32 = pv[31];
753 x33 = pv[32];
754 x34 = pv[33];
755 x35 = pv[34];
756 x36 = pv[35];
757 x37 = pv[36];
758 x38 = pv[37];
759 x39 = pv[38];
760 x40 = pv[39];
761 x41 = pv[40];
762 x42 = pv[41];
763 x43 = pv[42];
764 x44 = pv[43];
765 x45 = pv[44];
766 x46 = pv[45];
767 x47 = pv[46];
768 x48 = pv[47];
769 x49 = pv[48];
770 pc[0] = m1 = p1 * x1 + p8 * x2 + p15 * x3 + p22 * x4 + p29 * x5 + p36 * x6 + p43 * x7;
771 pc[1] = m2 = p2 * x1 + p9 * x2 + p16 * x3 + p23 * x4 + p30 * x5 + p37 * x6 + p44 * x7;
772 pc[2] = m3 = p3 * x1 + p10 * x2 + p17 * x3 + p24 * x4 + p31 * x5 + p38 * x6 + p45 * x7;
773 pc[3] = m4 = p4 * x1 + p11 * x2 + p18 * x3 + p25 * x4 + p32 * x5 + p39 * x6 + p46 * x7;
774 pc[4] = m5 = p5 * x1 + p12 * x2 + p19 * x3 + p26 * x4 + p33 * x5 + p40 * x6 + p47 * x7;
775 pc[5] = m6 = p6 * x1 + p13 * x2 + p20 * x3 + p27 * x4 + p34 * x5 + p41 * x6 + p48 * x7;
776 pc[6] = m7 = p7 * x1 + p14 * x2 + p21 * x3 + p28 * x4 + p35 * x5 + p42 * x6 + p49 * x7;
777
778 pc[7] = m8 = p1 * x8 + p8 * x9 + p15 * x10 + p22 * x11 + p29 * x12 + p36 * x13 + p43 * x14;
779 pc[8] = m9 = p2 * x8 + p9 * x9 + p16 * x10 + p23 * x11 + p30 * x12 + p37 * x13 + p44 * x14;
780 pc[9] = m10 = p3 * x8 + p10 * x9 + p17 * x10 + p24 * x11 + p31 * x12 + p38 * x13 + p45 * x14;
781 pc[10] = m11 = p4 * x8 + p11 * x9 + p18 * x10 + p25 * x11 + p32 * x12 + p39 * x13 + p46 * x14;
782 pc[11] = m12 = p5 * x8 + p12 * x9 + p19 * x10 + p26 * x11 + p33 * x12 + p40 * x13 + p47 * x14;
783 pc[12] = m13 = p6 * x8 + p13 * x9 + p20 * x10 + p27 * x11 + p34 * x12 + p41 * x13 + p48 * x14;
784 pc[13] = m14 = p7 * x8 + p14 * x9 + p21 * x10 + p28 * x11 + p35 * x12 + p42 * x13 + p49 * x14;
785
786 pc[14] = m15 = p1 * x15 + p8 * x16 + p15 * x17 + p22 * x18 + p29 * x19 + p36 * x20 + p43 * x21;
787 pc[15] = m16 = p2 * x15 + p9 * x16 + p16 * x17 + p23 * x18 + p30 * x19 + p37 * x20 + p44 * x21;
788 pc[16] = m17 = p3 * x15 + p10 * x16 + p17 * x17 + p24 * x18 + p31 * x19 + p38 * x20 + p45 * x21;
789 pc[17] = m18 = p4 * x15 + p11 * x16 + p18 * x17 + p25 * x18 + p32 * x19 + p39 * x20 + p46 * x21;
790 pc[18] = m19 = p5 * x15 + p12 * x16 + p19 * x17 + p26 * x18 + p33 * x19 + p40 * x20 + p47 * x21;
791 pc[19] = m20 = p6 * x15 + p13 * x16 + p20 * x17 + p27 * x18 + p34 * x19 + p41 * x20 + p48 * x21;
792 pc[20] = m21 = p7 * x15 + p14 * x16 + p21 * x17 + p28 * x18 + p35 * x19 + p42 * x20 + p49 * x21;
793
794 pc[21] = m22 = p1 * x22 + p8 * x23 + p15 * x24 + p22 * x25 + p29 * x26 + p36 * x27 + p43 * x28;
795 pc[22] = m23 = p2 * x22 + p9 * x23 + p16 * x24 + p23 * x25 + p30 * x26 + p37 * x27 + p44 * x28;
796 pc[23] = m24 = p3 * x22 + p10 * x23 + p17 * x24 + p24 * x25 + p31 * x26 + p38 * x27 + p45 * x28;
797 pc[24] = m25 = p4 * x22 + p11 * x23 + p18 * x24 + p25 * x25 + p32 * x26 + p39 * x27 + p46 * x28;
798 pc[25] = m26 = p5 * x22 + p12 * x23 + p19 * x24 + p26 * x25 + p33 * x26 + p40 * x27 + p47 * x28;
799 pc[26] = m27 = p6 * x22 + p13 * x23 + p20 * x24 + p27 * x25 + p34 * x26 + p41 * x27 + p48 * x28;
800 pc[27] = m28 = p7 * x22 + p14 * x23 + p21 * x24 + p28 * x25 + p35 * x26 + p42 * x27 + p49 * x28;
801
802 pc[28] = m29 = p1 * x29 + p8 * x30 + p15 * x31 + p22 * x32 + p29 * x33 + p36 * x34 + p43 * x35;
803 pc[29] = m30 = p2 * x29 + p9 * x30 + p16 * x31 + p23 * x32 + p30 * x33 + p37 * x34 + p44 * x35;
804 pc[30] = m31 = p3 * x29 + p10 * x30 + p17 * x31 + p24 * x32 + p31 * x33 + p38 * x34 + p45 * x35;
805 pc[31] = m32 = p4 * x29 + p11 * x30 + p18 * x31 + p25 * x32 + p32 * x33 + p39 * x34 + p46 * x35;
806 pc[32] = m33 = p5 * x29 + p12 * x30 + p19 * x31 + p26 * x32 + p33 * x33 + p40 * x34 + p47 * x35;
807 pc[33] = m34 = p6 * x29 + p13 * x30 + p20 * x31 + p27 * x32 + p34 * x33 + p41 * x34 + p48 * x35;
808 pc[34] = m35 = p7 * x29 + p14 * x30 + p21 * x31 + p28 * x32 + p35 * x33 + p42 * x34 + p49 * x35;
809
810 pc[35] = m36 = p1 * x36 + p8 * x37 + p15 * x38 + p22 * x39 + p29 * x40 + p36 * x41 + p43 * x42;
811 pc[36] = m37 = p2 * x36 + p9 * x37 + p16 * x38 + p23 * x39 + p30 * x40 + p37 * x41 + p44 * x42;
812 pc[37] = m38 = p3 * x36 + p10 * x37 + p17 * x38 + p24 * x39 + p31 * x40 + p38 * x41 + p45 * x42;
813 pc[38] = m39 = p4 * x36 + p11 * x37 + p18 * x38 + p25 * x39 + p32 * x40 + p39 * x41 + p46 * x42;
814 pc[39] = m40 = p5 * x36 + p12 * x37 + p19 * x38 + p26 * x39 + p33 * x40 + p40 * x41 + p47 * x42;
815 pc[40] = m41 = p6 * x36 + p13 * x37 + p20 * x38 + p27 * x39 + p34 * x40 + p41 * x41 + p48 * x42;
816 pc[41] = m42 = p7 * x36 + p14 * x37 + p21 * x38 + p28 * x39 + p35 * x40 + p42 * x41 + p49 * x42;
817
818 pc[42] = m43 = p1 * x43 + p8 * x44 + p15 * x45 + p22 * x46 + p29 * x47 + p36 * x48 + p43 * x49;
819 pc[43] = m44 = p2 * x43 + p9 * x44 + p16 * x45 + p23 * x46 + p30 * x47 + p37 * x48 + p44 * x49;
820 pc[44] = m45 = p3 * x43 + p10 * x44 + p17 * x45 + p24 * x46 + p31 * x47 + p38 * x48 + p45 * x49;
821 pc[45] = m46 = p4 * x43 + p11 * x44 + p18 * x45 + p25 * x46 + p32 * x47 + p39 * x48 + p46 * x49;
822 pc[46] = m47 = p5 * x43 + p12 * x44 + p19 * x45 + p26 * x46 + p33 * x47 + p40 * x48 + p47 * x49;
823 pc[47] = m48 = p6 * x43 + p13 * x44 + p20 * x45 + p27 * x46 + p34 * x47 + p41 * x48 + p48 * x49;
824 pc[48] = m49 = p7 * x43 + p14 * x44 + p21 * x45 + p28 * x46 + p35 * x47 + p42 * x48 + p49 * x49;
825
826 nz = bi[row + 1] - diag_offset[row] - 1;
827 pv += 49;
828 for (j = 0; j < nz; j++) {
829 x1 = pv[0];
830 x2 = pv[1];
831 x3 = pv[2];
832 x4 = pv[3];
833 x5 = pv[4];
834 x6 = pv[5];
835 x7 = pv[6];
836 x8 = pv[7];
837 x9 = pv[8];
838 x10 = pv[9];
839 x11 = pv[10];
840 x12 = pv[11];
841 x13 = pv[12];
842 x14 = pv[13];
843 x15 = pv[14];
844 x16 = pv[15];
845 x17 = pv[16];
846 x18 = pv[17];
847 x19 = pv[18];
848 x20 = pv[19];
849 x21 = pv[20];
850 x22 = pv[21];
851 x23 = pv[22];
852 x24 = pv[23];
853 x25 = pv[24];
854 x26 = pv[25];
855 x27 = pv[26];
856 x28 = pv[27];
857 x29 = pv[28];
858 x30 = pv[29];
859 x31 = pv[30];
860 x32 = pv[31];
861 x33 = pv[32];
862 x34 = pv[33];
863 x35 = pv[34];
864 x36 = pv[35];
865 x37 = pv[36];
866 x38 = pv[37];
867 x39 = pv[38];
868 x40 = pv[39];
869 x41 = pv[40];
870 x42 = pv[41];
871 x43 = pv[42];
872 x44 = pv[43];
873 x45 = pv[44];
874 x46 = pv[45];
875 x47 = pv[46];
876 x48 = pv[47];
877 x49 = pv[48];
878 x = rtmp + 49 * pj[j];
879 x[0] -= m1 * x1 + m8 * x2 + m15 * x3 + m22 * x4 + m29 * x5 + m36 * x6 + m43 * x7;
880 x[1] -= m2 * x1 + m9 * x2 + m16 * x3 + m23 * x4 + m30 * x5 + m37 * x6 + m44 * x7;
881 x[2] -= m3 * x1 + m10 * x2 + m17 * x3 + m24 * x4 + m31 * x5 + m38 * x6 + m45 * x7;
882 x[3] -= m4 * x1 + m11 * x2 + m18 * x3 + m25 * x4 + m32 * x5 + m39 * x6 + m46 * x7;
883 x[4] -= m5 * x1 + m12 * x2 + m19 * x3 + m26 * x4 + m33 * x5 + m40 * x6 + m47 * x7;
884 x[5] -= m6 * x1 + m13 * x2 + m20 * x3 + m27 * x4 + m34 * x5 + m41 * x6 + m48 * x7;
885 x[6] -= m7 * x1 + m14 * x2 + m21 * x3 + m28 * x4 + m35 * x5 + m42 * x6 + m49 * x7;
886
887 x[7] -= m1 * x8 + m8 * x9 + m15 * x10 + m22 * x11 + m29 * x12 + m36 * x13 + m43 * x14;
888 x[8] -= m2 * x8 + m9 * x9 + m16 * x10 + m23 * x11 + m30 * x12 + m37 * x13 + m44 * x14;
889 x[9] -= m3 * x8 + m10 * x9 + m17 * x10 + m24 * x11 + m31 * x12 + m38 * x13 + m45 * x14;
890 x[10] -= m4 * x8 + m11 * x9 + m18 * x10 + m25 * x11 + m32 * x12 + m39 * x13 + m46 * x14;
891 x[11] -= m5 * x8 + m12 * x9 + m19 * x10 + m26 * x11 + m33 * x12 + m40 * x13 + m47 * x14;
892 x[12] -= m6 * x8 + m13 * x9 + m20 * x10 + m27 * x11 + m34 * x12 + m41 * x13 + m48 * x14;
893 x[13] -= m7 * x8 + m14 * x9 + m21 * x10 + m28 * x11 + m35 * x12 + m42 * x13 + m49 * x14;
894
895 x[14] -= m1 * x15 + m8 * x16 + m15 * x17 + m22 * x18 + m29 * x19 + m36 * x20 + m43 * x21;
896 x[15] -= m2 * x15 + m9 * x16 + m16 * x17 + m23 * x18 + m30 * x19 + m37 * x20 + m44 * x21;
897 x[16] -= m3 * x15 + m10 * x16 + m17 * x17 + m24 * x18 + m31 * x19 + m38 * x20 + m45 * x21;
898 x[17] -= m4 * x15 + m11 * x16 + m18 * x17 + m25 * x18 + m32 * x19 + m39 * x20 + m46 * x21;
899 x[18] -= m5 * x15 + m12 * x16 + m19 * x17 + m26 * x18 + m33 * x19 + m40 * x20 + m47 * x21;
900 x[19] -= m6 * x15 + m13 * x16 + m20 * x17 + m27 * x18 + m34 * x19 + m41 * x20 + m48 * x21;
901 x[20] -= m7 * x15 + m14 * x16 + m21 * x17 + m28 * x18 + m35 * x19 + m42 * x20 + m49 * x21;
902
903 x[21] -= m1 * x22 + m8 * x23 + m15 * x24 + m22 * x25 + m29 * x26 + m36 * x27 + m43 * x28;
904 x[22] -= m2 * x22 + m9 * x23 + m16 * x24 + m23 * x25 + m30 * x26 + m37 * x27 + m44 * x28;
905 x[23] -= m3 * x22 + m10 * x23 + m17 * x24 + m24 * x25 + m31 * x26 + m38 * x27 + m45 * x28;
906 x[24] -= m4 * x22 + m11 * x23 + m18 * x24 + m25 * x25 + m32 * x26 + m39 * x27 + m46 * x28;
907 x[25] -= m5 * x22 + m12 * x23 + m19 * x24 + m26 * x25 + m33 * x26 + m40 * x27 + m47 * x28;
908 x[26] -= m6 * x22 + m13 * x23 + m20 * x24 + m27 * x25 + m34 * x26 + m41 * x27 + m48 * x28;
909 x[27] -= m7 * x22 + m14 * x23 + m21 * x24 + m28 * x25 + m35 * x26 + m42 * x27 + m49 * x28;
910
911 x[28] -= m1 * x29 + m8 * x30 + m15 * x31 + m22 * x32 + m29 * x33 + m36 * x34 + m43 * x35;
912 x[29] -= m2 * x29 + m9 * x30 + m16 * x31 + m23 * x32 + m30 * x33 + m37 * x34 + m44 * x35;
913 x[30] -= m3 * x29 + m10 * x30 + m17 * x31 + m24 * x32 + m31 * x33 + m38 * x34 + m45 * x35;
914 x[31] -= m4 * x29 + m11 * x30 + m18 * x31 + m25 * x32 + m32 * x33 + m39 * x34 + m46 * x35;
915 x[32] -= m5 * x29 + m12 * x30 + m19 * x31 + m26 * x32 + m33 * x33 + m40 * x34 + m47 * x35;
916 x[33] -= m6 * x29 + m13 * x30 + m20 * x31 + m27 * x32 + m34 * x33 + m41 * x34 + m48 * x35;
917 x[34] -= m7 * x29 + m14 * x30 + m21 * x31 + m28 * x32 + m35 * x33 + m42 * x34 + m49 * x35;
918
919 x[35] -= m1 * x36 + m8 * x37 + m15 * x38 + m22 * x39 + m29 * x40 + m36 * x41 + m43 * x42;
920 x[36] -= m2 * x36 + m9 * x37 + m16 * x38 + m23 * x39 + m30 * x40 + m37 * x41 + m44 * x42;
921 x[37] -= m3 * x36 + m10 * x37 + m17 * x38 + m24 * x39 + m31 * x40 + m38 * x41 + m45 * x42;
922 x[38] -= m4 * x36 + m11 * x37 + m18 * x38 + m25 * x39 + m32 * x40 + m39 * x41 + m46 * x42;
923 x[39] -= m5 * x36 + m12 * x37 + m19 * x38 + m26 * x39 + m33 * x40 + m40 * x41 + m47 * x42;
924 x[40] -= m6 * x36 + m13 * x37 + m20 * x38 + m27 * x39 + m34 * x40 + m41 * x41 + m48 * x42;
925 x[41] -= m7 * x36 + m14 * x37 + m21 * x38 + m28 * x39 + m35 * x40 + m42 * x41 + m49 * x42;
926
927 x[42] -= m1 * x43 + m8 * x44 + m15 * x45 + m22 * x46 + m29 * x47 + m36 * x48 + m43 * x49;
928 x[43] -= m2 * x43 + m9 * x44 + m16 * x45 + m23 * x46 + m30 * x47 + m37 * x48 + m44 * x49;
929 x[44] -= m3 * x43 + m10 * x44 + m17 * x45 + m24 * x46 + m31 * x47 + m38 * x48 + m45 * x49;
930 x[45] -= m4 * x43 + m11 * x44 + m18 * x45 + m25 * x46 + m32 * x47 + m39 * x48 + m46 * x49;
931 x[46] -= m5 * x43 + m12 * x44 + m19 * x45 + m26 * x46 + m33 * x47 + m40 * x48 + m47 * x49;
932 x[47] -= m6 * x43 + m13 * x44 + m20 * x45 + m27 * x46 + m34 * x47 + m41 * x48 + m48 * x49;
933 x[48] -= m7 * x43 + m14 * x44 + m21 * x45 + m28 * x46 + m35 * x47 + m42 * x48 + m49 * x49;
934 pv += 49;
935 }
936 PetscCall(PetscLogFlops(686.0 * nz + 637.0));
937 }
938 row = *ajtmp++;
939 }
940 /* finished row so stick it into b->a */
941 pv = ba + 49 * bi[i];
942 pj = bj + bi[i];
943 nz = bi[i + 1] - bi[i];
944 for (j = 0; j < nz; j++) {
945 x = rtmp + 49 * pj[j];
946 pv[0] = x[0];
947 pv[1] = x[1];
948 pv[2] = x[2];
949 pv[3] = x[3];
950 pv[4] = x[4];
951 pv[5] = x[5];
952 pv[6] = x[6];
953 pv[7] = x[7];
954 pv[8] = x[8];
955 pv[9] = x[9];
956 pv[10] = x[10];
957 pv[11] = x[11];
958 pv[12] = x[12];
959 pv[13] = x[13];
960 pv[14] = x[14];
961 pv[15] = x[15];
962 pv[16] = x[16];
963 pv[17] = x[17];
964 pv[18] = x[18];
965 pv[19] = x[19];
966 pv[20] = x[20];
967 pv[21] = x[21];
968 pv[22] = x[22];
969 pv[23] = x[23];
970 pv[24] = x[24];
971 pv[25] = x[25];
972 pv[26] = x[26];
973 pv[27] = x[27];
974 pv[28] = x[28];
975 pv[29] = x[29];
976 pv[30] = x[30];
977 pv[31] = x[31];
978 pv[32] = x[32];
979 pv[33] = x[33];
980 pv[34] = x[34];
981 pv[35] = x[35];
982 pv[36] = x[36];
983 pv[37] = x[37];
984 pv[38] = x[38];
985 pv[39] = x[39];
986 pv[40] = x[40];
987 pv[41] = x[41];
988 pv[42] = x[42];
989 pv[43] = x[43];
990 pv[44] = x[44];
991 pv[45] = x[45];
992 pv[46] = x[46];
993 pv[47] = x[47];
994 pv[48] = x[48];
995 pv += 49;
996 }
997 /* invert diagonal block */
998 w = ba + 49 * diag_offset[i];
999 PetscCall(PetscKernel_A_gets_inverse_A_7(w, shift, allowzeropivot, &zeropivotdetected));
1000 if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
1001 }
1002
1003 PetscCall(PetscFree(rtmp));
1004
1005 C->ops->solve = MatSolve_SeqBAIJ_7_NaturalOrdering_inplace;
1006 C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_7_NaturalOrdering_inplace;
1007 C->assembled = PETSC_TRUE;
1008
1009 PetscCall(PetscLogFlops(1.333333333333 * 7 * 7 * 7 * b->mbs)); /* from inverting diagonal blocks */
1010 PetscFunctionReturn(PETSC_SUCCESS);
1011 }
1012
MatLUFactorNumeric_SeqBAIJ_7_NaturalOrdering(Mat B,Mat A,const MatFactorInfo * info)1013 PetscErrorCode MatLUFactorNumeric_SeqBAIJ_7_NaturalOrdering(Mat B, Mat A, const MatFactorInfo *info)
1014 {
1015 Mat C = B;
1016 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)C->data;
1017 PetscInt i, j, k, nz, nzL, row;
1018 const PetscInt n = a->mbs, *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j;
1019 const PetscInt *ajtmp, *bjtmp, *bdiag = b->diag, *pj, bs2 = a->bs2;
1020 MatScalar *rtmp, *pc, *mwork, *v, *pv, *aa = a->a;
1021 PetscInt flg;
1022 PetscReal shift = info->shiftamount;
1023 PetscBool allowzeropivot, zeropivotdetected;
1024
1025 PetscFunctionBegin;
1026 allowzeropivot = PetscNot(A->erroriffailure);
1027
1028 /* generate work space needed by the factorization */
1029 PetscCall(PetscMalloc2(bs2 * n, &rtmp, bs2, &mwork));
1030 PetscCall(PetscArrayzero(rtmp, bs2 * n));
1031
1032 for (i = 0; i < n; i++) {
1033 /* zero rtmp */
1034 /* L part */
1035 nz = bi[i + 1] - bi[i];
1036 bjtmp = bj + bi[i];
1037 for (j = 0; j < nz; j++) PetscCall(PetscArrayzero(rtmp + bs2 * bjtmp[j], bs2));
1038
1039 /* U part */
1040 nz = bdiag[i] - bdiag[i + 1];
1041 bjtmp = bj + bdiag[i + 1] + 1;
1042 for (j = 0; j < nz; j++) PetscCall(PetscArrayzero(rtmp + bs2 * bjtmp[j], bs2));
1043
1044 /* load in initial (unfactored row) */
1045 nz = ai[i + 1] - ai[i];
1046 ajtmp = aj + ai[i];
1047 v = aa + bs2 * ai[i];
1048 for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(rtmp + bs2 * ajtmp[j], v + bs2 * j, bs2));
1049
1050 /* elimination */
1051 bjtmp = bj + bi[i];
1052 nzL = bi[i + 1] - bi[i];
1053 for (k = 0; k < nzL; k++) {
1054 row = bjtmp[k];
1055 pc = rtmp + bs2 * row;
1056 for (flg = 0, j = 0; j < bs2; j++) {
1057 if (pc[j] != 0.0) {
1058 flg = 1;
1059 break;
1060 }
1061 }
1062 if (flg) {
1063 pv = b->a + bs2 * bdiag[row];
1064 /* PetscKernel_A_gets_A_times_B(bs,pc,pv,mwork); *pc = *pc * (*pv); */
1065 PetscCall(PetscKernel_A_gets_A_times_B_7(pc, pv, mwork));
1066
1067 pj = b->j + bdiag[row + 1] + 1; /* beginning of U(row,:) */
1068 pv = b->a + bs2 * (bdiag[row + 1] + 1);
1069 nz = bdiag[row] - bdiag[row + 1] - 1; /* num of entries inU(row,:), excluding diag */
1070 for (j = 0; j < nz; j++) {
1071 /* PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j); */
1072 /* rtmp+bs2*pj[j] = rtmp+bs2*pj[j] - (*pc)*(pv+bs2*j) */
1073 v = rtmp + bs2 * pj[j];
1074 PetscCall(PetscKernel_A_gets_A_minus_B_times_C_7(v, pc, pv));
1075 pv += bs2;
1076 }
1077 PetscCall(PetscLogFlops(686.0 * nz + 637)); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */
1078 }
1079 }
1080
1081 /* finished row so stick it into b->a */
1082 /* L part */
1083 pv = b->a + bs2 * bi[i];
1084 pj = b->j + bi[i];
1085 nz = bi[i + 1] - bi[i];
1086 for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2));
1087
1088 /* Mark diagonal and invert diagonal for simpler triangular solves */
1089 pv = b->a + bs2 * bdiag[i];
1090 pj = b->j + bdiag[i];
1091 PetscCall(PetscArraycpy(pv, rtmp + bs2 * pj[0], bs2));
1092 PetscCall(PetscKernel_A_gets_inverse_A_7(pv, shift, allowzeropivot, &zeropivotdetected));
1093 if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
1094
1095 /* U part */
1096 pv = b->a + bs2 * (bdiag[i + 1] + 1);
1097 pj = b->j + bdiag[i + 1] + 1;
1098 nz = bdiag[i] - bdiag[i + 1] - 1;
1099 for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2));
1100 }
1101 PetscCall(PetscFree2(rtmp, mwork));
1102
1103 C->ops->solve = MatSolve_SeqBAIJ_7_NaturalOrdering;
1104 C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_7_NaturalOrdering;
1105 C->assembled = PETSC_TRUE;
1106
1107 PetscCall(PetscLogFlops(1.333333333333 * 7 * 7 * 7 * n)); /* from inverting diagonal blocks */
1108 PetscFunctionReturn(PETSC_SUCCESS);
1109 }
1110