1 #include <../src/mat/impls/sbaij/seq/sbaij.h>
2 #include <petsc/private/kernels/blockinvert.h>
3
4 /*
5 Version for when blocks are 7 by 7 Using natural ordering
6 */
MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering(Mat C,Mat A,const MatFactorInfo * info)7 PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering(Mat C, Mat A, const MatFactorInfo *info)
8 {
9 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data, *b = (Mat_SeqSBAIJ *)C->data;
10 PetscInt i, j, mbs = a->mbs, *bi = b->i, *bj = b->j;
11 PetscInt *ai, *aj, k, k1, jmin, jmax, *jl, *il, vj, nexti, ili;
12 MatScalar *ba = b->a, *aa, *ap, *dk, *uik;
13 MatScalar *u, *d, *w, *wp, u0, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10, u11, u12;
14 MatScalar u13, u14, u15, u16, u17, u18, u19, u20, u21, u22, u23, u24, u25, u26, u27;
15 MatScalar u28, u29, u30, u31, u32, u33, u34, u35, u36, u37, u38, u39, u40, u41;
16 MatScalar u42, u43, u44, u45, u46, u47, u48;
17 PetscReal shift = info->shiftamount;
18 PetscBool allowzeropivot, zeropivotdetected;
19
20 PetscFunctionBegin;
21 /* initialization */
22 allowzeropivot = PetscNot(A->erroriffailure);
23 PetscCall(PetscCalloc1(49 * mbs, &w));
24 PetscCall(PetscMalloc2(mbs, &il, mbs, &jl));
25 il[0] = 0;
26 for (i = 0; i < mbs; i++) jl[i] = mbs;
27
28 PetscCall(PetscMalloc2(49, &dk, 49, &uik));
29 ai = a->i;
30 aj = a->j;
31 aa = a->a;
32
33 /* for each row k */
34 for (k = 0; k < mbs; k++) {
35 /*initialize k-th row with elements nonzero in row k of A */
36 jmin = ai[k];
37 jmax = ai[k + 1];
38 if (jmin < jmax) {
39 ap = aa + jmin * 49;
40 for (j = jmin; j < jmax; j++) {
41 vj = aj[j]; /* block col. index */
42 wp = w + vj * 49;
43 PetscCall(PetscArraycpy(wp, ap, 49));
44 ap += 49;
45 }
46 }
47
48 /* modify k-th row by adding in those rows i with U(i,k) != 0 */
49 PetscCall(PetscArraycpy(dk, w + k * 49, 49));
50 i = jl[k]; /* first row to be added to k_th row */
51
52 while (i < mbs) {
53 nexti = jl[i]; /* next row to be added to k_th row */
54
55 /* compute multiplier */
56 ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
57
58 /* uik = -inv(Di)*U_bar(i,k) */
59 d = ba + i * 49;
60 u = ba + ili * 49;
61
62 u0 = u[0];
63 u1 = u[1];
64 u2 = u[2];
65 u3 = u[3];
66 u4 = u[4];
67 u5 = u[5];
68 u6 = u[6];
69 u7 = u[7];
70 u8 = u[8];
71 u9 = u[9];
72 u10 = u[10];
73 u11 = u[11];
74 u12 = u[12];
75 u13 = u[13];
76 u14 = u[14];
77 u15 = u[15];
78 u16 = u[16];
79 u17 = u[17];
80 u18 = u[18];
81 u19 = u[19];
82 u20 = u[20];
83 u21 = u[21];
84 u22 = u[22];
85 u23 = u[23];
86 u24 = u[24];
87 u25 = u[25];
88 u26 = u[26];
89 u27 = u[27];
90 u28 = u[28];
91 u29 = u[29];
92 u30 = u[30];
93 u31 = u[31];
94 u32 = u[32];
95 u33 = u[33];
96 u34 = u[34];
97 u35 = u[35];
98 u36 = u[36];
99 u37 = u[37];
100 u38 = u[38];
101 u39 = u[39];
102 u40 = u[40];
103 u41 = u[41];
104 u42 = u[42];
105 u43 = u[43];
106 u44 = u[44];
107 u45 = u[45];
108 u46 = u[46];
109 u47 = u[47];
110 u48 = u[48];
111
112 uik[0] = -(d[0] * u0 + d[7] * u1 + d[14] * u2 + d[21] * u3 + d[28] * u4 + d[35] * u5 + d[42] * u6);
113 uik[1] = -(d[1] * u0 + d[8] * u1 + d[15] * u2 + d[22] * u3 + d[29] * u4 + d[36] * u5 + d[43] * u6);
114 uik[2] = -(d[2] * u0 + d[9] * u1 + d[16] * u2 + d[23] * u3 + d[30] * u4 + d[37] * u5 + d[44] * u6);
115 uik[3] = -(d[3] * u0 + d[10] * u1 + d[17] * u2 + d[24] * u3 + d[31] * u4 + d[38] * u5 + d[45] * u6);
116 uik[4] = -(d[4] * u0 + d[11] * u1 + d[18] * u2 + d[25] * u3 + d[32] * u4 + d[39] * u5 + d[46] * u6);
117 uik[5] = -(d[5] * u0 + d[12] * u1 + d[19] * u2 + d[26] * u3 + d[33] * u4 + d[40] * u5 + d[47] * u6);
118 uik[6] = -(d[6] * u0 + d[13] * u1 + d[20] * u2 + d[27] * u3 + d[34] * u4 + d[41] * u5 + d[48] * u6);
119
120 uik[7] = -(d[0] * u7 + d[7] * u8 + d[14] * u9 + d[21] * u10 + d[28] * u11 + d[35] * u12 + d[42] * u13);
121 uik[8] = -(d[1] * u7 + d[8] * u8 + d[15] * u9 + d[22] * u10 + d[29] * u11 + d[36] * u12 + d[43] * u13);
122 uik[9] = -(d[2] * u7 + d[9] * u8 + d[16] * u9 + d[23] * u10 + d[30] * u11 + d[37] * u12 + d[44] * u13);
123 uik[10] = -(d[3] * u7 + d[10] * u8 + d[17] * u9 + d[24] * u10 + d[31] * u11 + d[38] * u12 + d[45] * u13);
124 uik[11] = -(d[4] * u7 + d[11] * u8 + d[18] * u9 + d[25] * u10 + d[32] * u11 + d[39] * u12 + d[46] * u13);
125 uik[12] = -(d[5] * u7 + d[12] * u8 + d[19] * u9 + d[26] * u10 + d[33] * u11 + d[40] * u12 + d[47] * u13);
126 uik[13] = -(d[6] * u7 + d[13] * u8 + d[20] * u9 + d[27] * u10 + d[34] * u11 + d[41] * u12 + d[48] * u13);
127
128 uik[14] = -(d[0] * u14 + d[7] * u15 + d[14] * u16 + d[21] * u17 + d[28] * u18 + d[35] * u19 + d[42] * u20);
129 uik[15] = -(d[1] * u14 + d[8] * u15 + d[15] * u16 + d[22] * u17 + d[29] * u18 + d[36] * u19 + d[43] * u20);
130 uik[16] = -(d[2] * u14 + d[9] * u15 + d[16] * u16 + d[23] * u17 + d[30] * u18 + d[37] * u19 + d[44] * u20);
131 uik[17] = -(d[3] * u14 + d[10] * u15 + d[17] * u16 + d[24] * u17 + d[31] * u18 + d[38] * u19 + d[45] * u20);
132 uik[18] = -(d[4] * u14 + d[11] * u15 + d[18] * u16 + d[25] * u17 + d[32] * u18 + d[39] * u19 + d[46] * u20);
133 uik[19] = -(d[5] * u14 + d[12] * u15 + d[19] * u16 + d[26] * u17 + d[33] * u18 + d[40] * u19 + d[47] * u20);
134 uik[20] = -(d[6] * u14 + d[13] * u15 + d[20] * u16 + d[27] * u17 + d[34] * u18 + d[41] * u19 + d[48] * u20);
135
136 uik[21] = -(d[0] * u21 + d[7] * u22 + d[14] * u23 + d[21] * u24 + d[28] * u25 + d[35] * u26 + d[42] * u27);
137 uik[22] = -(d[1] * u21 + d[8] * u22 + d[15] * u23 + d[22] * u24 + d[29] * u25 + d[36] * u26 + d[43] * u27);
138 uik[23] = -(d[2] * u21 + d[9] * u22 + d[16] * u23 + d[23] * u24 + d[30] * u25 + d[37] * u26 + d[44] * u27);
139 uik[24] = -(d[3] * u21 + d[10] * u22 + d[17] * u23 + d[24] * u24 + d[31] * u25 + d[38] * u26 + d[45] * u27);
140 uik[25] = -(d[4] * u21 + d[11] * u22 + d[18] * u23 + d[25] * u24 + d[32] * u25 + d[39] * u26 + d[46] * u27);
141 uik[26] = -(d[5] * u21 + d[12] * u22 + d[19] * u23 + d[26] * u24 + d[33] * u25 + d[40] * u26 + d[47] * u27);
142 uik[27] = -(d[6] * u21 + d[13] * u22 + d[20] * u23 + d[27] * u24 + d[34] * u25 + d[41] * u26 + d[48] * u27);
143
144 uik[28] = -(d[0] * u28 + d[7] * u29 + d[14] * u30 + d[21] * u31 + d[28] * u32 + d[35] * u33 + d[42] * u34);
145 uik[29] = -(d[1] * u28 + d[8] * u29 + d[15] * u30 + d[22] * u31 + d[29] * u32 + d[36] * u33 + d[43] * u34);
146 uik[30] = -(d[2] * u28 + d[9] * u29 + d[16] * u30 + d[23] * u31 + d[30] * u32 + d[37] * u33 + d[44] * u34);
147 uik[31] = -(d[3] * u28 + d[10] * u29 + d[17] * u30 + d[24] * u31 + d[31] * u32 + d[38] * u33 + d[45] * u34);
148 uik[32] = -(d[4] * u28 + d[11] * u29 + d[18] * u30 + d[25] * u31 + d[32] * u32 + d[39] * u33 + d[46] * u34);
149 uik[33] = -(d[5] * u28 + d[12] * u29 + d[19] * u30 + d[26] * u31 + d[33] * u32 + d[40] * u33 + d[47] * u34);
150 uik[34] = -(d[6] * u28 + d[13] * u29 + d[20] * u30 + d[27] * u31 + d[34] * u32 + d[41] * u33 + d[48] * u34);
151
152 uik[35] = -(d[0] * u35 + d[7] * u36 + d[14] * u37 + d[21] * u38 + d[28] * u39 + d[35] * u40 + d[42] * u41);
153 uik[36] = -(d[1] * u35 + d[8] * u36 + d[15] * u37 + d[22] * u38 + d[29] * u39 + d[36] * u40 + d[43] * u41);
154 uik[37] = -(d[2] * u35 + d[9] * u36 + d[16] * u37 + d[23] * u38 + d[30] * u39 + d[37] * u40 + d[44] * u41);
155 uik[38] = -(d[3] * u35 + d[10] * u36 + d[17] * u37 + d[24] * u38 + d[31] * u39 + d[38] * u40 + d[45] * u41);
156 uik[39] = -(d[4] * u35 + d[11] * u36 + d[18] * u37 + d[25] * u38 + d[32] * u39 + d[39] * u40 + d[46] * u41);
157 uik[40] = -(d[5] * u35 + d[12] * u36 + d[19] * u37 + d[26] * u38 + d[33] * u39 + d[40] * u40 + d[47] * u41);
158 uik[41] = -(d[6] * u35 + d[13] * u36 + d[20] * u37 + d[27] * u38 + d[34] * u39 + d[41] * u40 + d[48] * u41);
159
160 uik[42] = -(d[0] * u42 + d[7] * u43 + d[14] * u44 + d[21] * u45 + d[28] * u46 + d[35] * u47 + d[42] * u48);
161 uik[43] = -(d[1] * u42 + d[8] * u43 + d[15] * u44 + d[22] * u45 + d[29] * u46 + d[36] * u47 + d[43] * u48);
162 uik[44] = -(d[2] * u42 + d[9] * u43 + d[16] * u44 + d[23] * u45 + d[30] * u46 + d[37] * u47 + d[44] * u48);
163 uik[45] = -(d[3] * u42 + d[10] * u43 + d[17] * u44 + d[24] * u45 + d[31] * u46 + d[38] * u47 + d[45] * u48);
164 uik[46] = -(d[4] * u42 + d[11] * u43 + d[18] * u44 + d[25] * u45 + d[32] * u46 + d[39] * u47 + d[46] * u48);
165 uik[47] = -(d[5] * u42 + d[12] * u43 + d[19] * u44 + d[26] * u45 + d[33] * u46 + d[40] * u47 + d[47] * u48);
166 uik[48] = -(d[6] * u42 + d[13] * u43 + d[20] * u44 + d[27] * u45 + d[34] * u46 + d[41] * u47 + d[48] * u48);
167
168 /* update D(k) += -U(i,k)^T * U_bar(i,k) */
169 dk[0] += uik[0] * u0 + uik[1] * u1 + uik[2] * u2 + uik[3] * u3 + uik[4] * u4 + uik[5] * u5 + uik[6] * u6;
170 dk[1] += uik[7] * u0 + uik[8] * u1 + uik[9] * u2 + uik[10] * u3 + uik[11] * u4 + uik[12] * u5 + uik[13] * u6;
171 dk[2] += uik[14] * u0 + uik[15] * u1 + uik[16] * u2 + uik[17] * u3 + uik[18] * u4 + uik[19] * u5 + uik[20] * u6;
172 dk[3] += uik[21] * u0 + uik[22] * u1 + uik[23] * u2 + uik[24] * u3 + uik[25] * u4 + uik[26] * u5 + uik[27] * u6;
173 dk[4] += uik[28] * u0 + uik[29] * u1 + uik[30] * u2 + uik[31] * u3 + uik[32] * u4 + uik[33] * u5 + uik[34] * u6;
174 dk[5] += uik[35] * u0 + uik[36] * u1 + uik[37] * u2 + uik[38] * u3 + uik[39] * u4 + uik[40] * u5 + uik[41] * u6;
175 dk[6] += uik[42] * u0 + uik[43] * u1 + uik[44] * u2 + uik[45] * u3 + uik[46] * u4 + uik[47] * u5 + uik[48] * u6;
176
177 dk[7] += uik[0] * u7 + uik[1] * u8 + uik[2] * u9 + uik[3] * u10 + uik[4] * u11 + uik[5] * u12 + uik[6] * u13;
178 dk[8] += uik[7] * u7 + uik[8] * u8 + uik[9] * u9 + uik[10] * u10 + uik[11] * u11 + uik[12] * u12 + uik[13] * u13;
179 dk[9] += uik[14] * u7 + uik[15] * u8 + uik[16] * u9 + uik[17] * u10 + uik[18] * u11 + uik[19] * u12 + uik[20] * u13;
180 dk[10] += uik[21] * u7 + uik[22] * u8 + uik[23] * u9 + uik[24] * u10 + uik[25] * u11 + uik[26] * u12 + uik[27] * u13;
181 dk[11] += uik[28] * u7 + uik[29] * u8 + uik[30] * u9 + uik[31] * u10 + uik[32] * u11 + uik[33] * u12 + uik[34] * u13;
182 dk[12] += uik[35] * u7 + uik[36] * u8 + uik[37] * u9 + uik[38] * u10 + uik[39] * u11 + uik[40] * u12 + uik[41] * u13;
183 dk[13] += uik[42] * u7 + uik[43] * u8 + uik[44] * u9 + uik[45] * u10 + uik[46] * u11 + uik[47] * u12 + uik[48] * u13;
184
185 dk[14] += uik[0] * u14 + uik[1] * u15 + uik[2] * u16 + uik[3] * u17 + uik[4] * u18 + uik[5] * u19 + uik[6] * u20;
186 dk[15] += uik[7] * u14 + uik[8] * u15 + uik[9] * u16 + uik[10] * u17 + uik[11] * u18 + uik[12] * u19 + uik[13] * u20;
187 dk[16] += uik[14] * u14 + uik[15] * u15 + uik[16] * u16 + uik[17] * u17 + uik[18] * u18 + uik[19] * u19 + uik[20] * u20;
188 dk[17] += uik[21] * u14 + uik[22] * u15 + uik[23] * u16 + uik[24] * u17 + uik[25] * u18 + uik[26] * u19 + uik[27] * u20;
189 dk[18] += uik[28] * u14 + uik[29] * u15 + uik[30] * u16 + uik[31] * u17 + uik[32] * u18 + uik[33] * u19 + uik[34] * u20;
190 dk[19] += uik[35] * u14 + uik[36] * u15 + uik[37] * u16 + uik[38] * u17 + uik[39] * u18 + uik[40] * u19 + uik[41] * u20;
191 dk[20] += uik[42] * u14 + uik[43] * u15 + uik[44] * u16 + uik[45] * u17 + uik[46] * u18 + uik[47] * u19 + uik[48] * u20;
192
193 dk[21] += uik[0] * u21 + uik[1] * u22 + uik[2] * u23 + uik[3] * u24 + uik[4] * u25 + uik[5] * u26 + uik[6] * u27;
194 dk[22] += uik[7] * u21 + uik[8] * u22 + uik[9] * u23 + uik[10] * u24 + uik[11] * u25 + uik[12] * u26 + uik[13] * u27;
195 dk[23] += uik[14] * u21 + uik[15] * u22 + uik[16] * u23 + uik[17] * u24 + uik[18] * u25 + uik[19] * u26 + uik[20] * u27;
196 dk[24] += uik[21] * u21 + uik[22] * u22 + uik[23] * u23 + uik[24] * u24 + uik[25] * u25 + uik[26] * u26 + uik[27] * u27;
197 dk[25] += uik[28] * u21 + uik[29] * u22 + uik[30] * u23 + uik[31] * u24 + uik[32] * u25 + uik[33] * u26 + uik[34] * u27;
198 dk[26] += uik[35] * u21 + uik[36] * u22 + uik[37] * u23 + uik[38] * u24 + uik[39] * u25 + uik[40] * u26 + uik[41] * u27;
199 dk[27] += uik[42] * u21 + uik[43] * u22 + uik[44] * u23 + uik[45] * u24 + uik[46] * u25 + uik[47] * u26 + uik[48] * u27;
200
201 dk[28] += uik[0] * u28 + uik[1] * u29 + uik[2] * u30 + uik[3] * u31 + uik[4] * u32 + uik[5] * u33 + uik[6] * u34;
202 dk[29] += uik[7] * u28 + uik[8] * u29 + uik[9] * u30 + uik[10] * u31 + uik[11] * u32 + uik[12] * u33 + uik[13] * u34;
203 dk[30] += uik[14] * u28 + uik[15] * u29 + uik[16] * u30 + uik[17] * u31 + uik[18] * u32 + uik[19] * u33 + uik[20] * u34;
204 dk[31] += uik[21] * u28 + uik[22] * u29 + uik[23] * u30 + uik[24] * u31 + uik[25] * u32 + uik[26] * u33 + uik[27] * u34;
205 dk[32] += uik[28] * u28 + uik[29] * u29 + uik[30] * u30 + uik[31] * u31 + uik[32] * u32 + uik[33] * u33 + uik[34] * u34;
206 dk[33] += uik[35] * u28 + uik[36] * u29 + uik[37] * u30 + uik[38] * u31 + uik[39] * u32 + uik[40] * u33 + uik[41] * u34;
207 dk[34] += uik[42] * u28 + uik[43] * u29 + uik[44] * u30 + uik[45] * u31 + uik[46] * u32 + uik[47] * u33 + uik[48] * u34;
208
209 dk[35] += uik[0] * u35 + uik[1] * u36 + uik[2] * u37 + uik[3] * u38 + uik[4] * u39 + uik[5] * u40 + uik[6] * u41;
210 dk[36] += uik[7] * u35 + uik[8] * u36 + uik[9] * u37 + uik[10] * u38 + uik[11] * u39 + uik[12] * u40 + uik[13] * u41;
211 dk[37] += uik[14] * u35 + uik[15] * u36 + uik[16] * u37 + uik[17] * u38 + uik[18] * u39 + uik[19] * u40 + uik[20] * u41;
212 dk[38] += uik[21] * u35 + uik[22] * u36 + uik[23] * u37 + uik[24] * u38 + uik[25] * u39 + uik[26] * u40 + uik[27] * u41;
213 dk[39] += uik[28] * u35 + uik[29] * u36 + uik[30] * u37 + uik[31] * u38 + uik[32] * u39 + uik[33] * u40 + uik[34] * u41;
214 dk[40] += uik[35] * u35 + uik[36] * u36 + uik[37] * u37 + uik[38] * u38 + uik[39] * u39 + uik[40] * u40 + uik[41] * u41;
215 dk[41] += uik[42] * u35 + uik[43] * u36 + uik[44] * u37 + uik[45] * u38 + uik[46] * u39 + uik[47] * u40 + uik[48] * u41;
216
217 dk[42] += uik[0] * u42 + uik[1] * u43 + uik[2] * u44 + uik[3] * u45 + uik[4] * u46 + uik[5] * u47 + uik[6] * u48;
218 dk[43] += uik[7] * u42 + uik[8] * u43 + uik[9] * u44 + uik[10] * u45 + uik[11] * u46 + uik[12] * u47 + uik[13] * u48;
219 dk[44] += uik[14] * u42 + uik[15] * u43 + uik[16] * u44 + uik[17] * u45 + uik[18] * u46 + uik[19] * u47 + uik[20] * u48;
220 dk[45] += uik[21] * u42 + uik[22] * u43 + uik[23] * u44 + uik[24] * u45 + uik[25] * u46 + uik[26] * u47 + uik[27] * u48;
221 dk[46] += uik[28] * u42 + uik[29] * u43 + uik[30] * u44 + uik[31] * u45 + uik[32] * u46 + uik[33] * u47 + uik[34] * u48;
222 dk[47] += uik[35] * u42 + uik[36] * u43 + uik[37] * u44 + uik[38] * u45 + uik[39] * u46 + uik[40] * u47 + uik[41] * u48;
223 dk[48] += uik[42] * u42 + uik[43] * u43 + uik[44] * u44 + uik[45] * u45 + uik[46] * u46 + uik[47] * u47 + uik[48] * u48;
224
225 PetscCall(PetscLogFlops(343.0 * 4.0));
226
227 /* update -U(i,k) */
228 PetscCall(PetscArraycpy(ba + ili * 49, uik, 49));
229
230 /* add multiple of row i to k-th row ... */
231 jmin = ili + 1;
232 jmax = bi[i + 1];
233 if (jmin < jmax) {
234 for (j = jmin; j < jmax; j++) {
235 /* w += -U(i,k)^T * U_bar(i,j) */
236 wp = w + bj[j] * 49;
237 u = ba + j * 49;
238
239 u0 = u[0];
240 u1 = u[1];
241 u2 = u[2];
242 u3 = u[3];
243 u4 = u[4];
244 u5 = u[5];
245 u6 = u[6];
246 u7 = u[7];
247 u8 = u[8];
248 u9 = u[9];
249 u10 = u[10];
250 u11 = u[11];
251 u12 = u[12];
252 u13 = u[13];
253 u14 = u[14];
254 u15 = u[15];
255 u16 = u[16];
256 u17 = u[17];
257 u18 = u[18];
258 u19 = u[19];
259 u20 = u[20];
260 u21 = u[21];
261 u22 = u[22];
262 u23 = u[23];
263 u24 = u[24];
264 u25 = u[25];
265 u26 = u[26];
266 u27 = u[27];
267 u28 = u[28];
268 u29 = u[29];
269 u30 = u[30];
270 u31 = u[31];
271 u32 = u[32];
272 u33 = u[33];
273 u34 = u[34];
274 u35 = u[35];
275 u36 = u[36];
276 u37 = u[37];
277 u38 = u[38];
278 u39 = u[39];
279 u40 = u[40];
280 u41 = u[41];
281 u42 = u[42];
282 u43 = u[43];
283 u44 = u[44];
284 u45 = u[45];
285 u46 = u[46];
286 u47 = u[47];
287 u48 = u[48];
288
289 wp[0] += uik[0] * u0 + uik[1] * u1 + uik[2] * u2 + uik[3] * u3 + uik[4] * u4 + uik[5] * u5 + uik[6] * u6;
290 wp[1] += uik[7] * u0 + uik[8] * u1 + uik[9] * u2 + uik[10] * u3 + uik[11] * u4 + uik[12] * u5 + uik[13] * u6;
291 wp[2] += uik[14] * u0 + uik[15] * u1 + uik[16] * u2 + uik[17] * u3 + uik[18] * u4 + uik[19] * u5 + uik[20] * u6;
292 wp[3] += uik[21] * u0 + uik[22] * u1 + uik[23] * u2 + uik[24] * u3 + uik[25] * u4 + uik[26] * u5 + uik[27] * u6;
293 wp[4] += uik[28] * u0 + uik[29] * u1 + uik[30] * u2 + uik[31] * u3 + uik[32] * u4 + uik[33] * u5 + uik[34] * u6;
294 wp[5] += uik[35] * u0 + uik[36] * u1 + uik[37] * u2 + uik[38] * u3 + uik[39] * u4 + uik[40] * u5 + uik[41] * u6;
295 wp[6] += uik[42] * u0 + uik[43] * u1 + uik[44] * u2 + uik[45] * u3 + uik[46] * u4 + uik[47] * u5 + uik[48] * u6;
296
297 wp[7] += uik[0] * u7 + uik[1] * u8 + uik[2] * u9 + uik[3] * u10 + uik[4] * u11 + uik[5] * u12 + uik[6] * u13;
298 wp[8] += uik[7] * u7 + uik[8] * u8 + uik[9] * u9 + uik[10] * u10 + uik[11] * u11 + uik[12] * u12 + uik[13] * u13;
299 wp[9] += uik[14] * u7 + uik[15] * u8 + uik[16] * u9 + uik[17] * u10 + uik[18] * u11 + uik[19] * u12 + uik[20] * u13;
300 wp[10] += uik[21] * u7 + uik[22] * u8 + uik[23] * u9 + uik[24] * u10 + uik[25] * u11 + uik[26] * u12 + uik[27] * u13;
301 wp[11] += uik[28] * u7 + uik[29] * u8 + uik[30] * u9 + uik[31] * u10 + uik[32] * u11 + uik[33] * u12 + uik[34] * u13;
302 wp[12] += uik[35] * u7 + uik[36] * u8 + uik[37] * u9 + uik[38] * u10 + uik[39] * u11 + uik[40] * u12 + uik[41] * u13;
303 wp[13] += uik[42] * u7 + uik[43] * u8 + uik[44] * u9 + uik[45] * u10 + uik[46] * u11 + uik[47] * u12 + uik[48] * u13;
304
305 wp[14] += uik[0] * u14 + uik[1] * u15 + uik[2] * u16 + uik[3] * u17 + uik[4] * u18 + uik[5] * u19 + uik[6] * u20;
306 wp[15] += uik[7] * u14 + uik[8] * u15 + uik[9] * u16 + uik[10] * u17 + uik[11] * u18 + uik[12] * u19 + uik[13] * u20;
307 wp[16] += uik[14] * u14 + uik[15] * u15 + uik[16] * u16 + uik[17] * u17 + uik[18] * u18 + uik[19] * u19 + uik[20] * u20;
308 wp[17] += uik[21] * u14 + uik[22] * u15 + uik[23] * u16 + uik[24] * u17 + uik[25] * u18 + uik[26] * u19 + uik[27] * u20;
309 wp[18] += uik[28] * u14 + uik[29] * u15 + uik[30] * u16 + uik[31] * u17 + uik[32] * u18 + uik[33] * u19 + uik[34] * u20;
310 wp[19] += uik[35] * u14 + uik[36] * u15 + uik[37] * u16 + uik[38] * u17 + uik[39] * u18 + uik[40] * u19 + uik[41] * u20;
311 wp[20] += uik[42] * u14 + uik[43] * u15 + uik[44] * u16 + uik[45] * u17 + uik[46] * u18 + uik[47] * u19 + uik[48] * u20;
312
313 wp[21] += uik[0] * u21 + uik[1] * u22 + uik[2] * u23 + uik[3] * u24 + uik[4] * u25 + uik[5] * u26 + uik[6] * u27;
314 wp[22] += uik[7] * u21 + uik[8] * u22 + uik[9] * u23 + uik[10] * u24 + uik[11] * u25 + uik[12] * u26 + uik[13] * u27;
315 wp[23] += uik[14] * u21 + uik[15] * u22 + uik[16] * u23 + uik[17] * u24 + uik[18] * u25 + uik[19] * u26 + uik[20] * u27;
316 wp[24] += uik[21] * u21 + uik[22] * u22 + uik[23] * u23 + uik[24] * u24 + uik[25] * u25 + uik[26] * u26 + uik[27] * u27;
317 wp[25] += uik[28] * u21 + uik[29] * u22 + uik[30] * u23 + uik[31] * u24 + uik[32] * u25 + uik[33] * u26 + uik[34] * u27;
318 wp[26] += uik[35] * u21 + uik[36] * u22 + uik[37] * u23 + uik[38] * u24 + uik[39] * u25 + uik[40] * u26 + uik[41] * u27;
319 wp[27] += uik[42] * u21 + uik[43] * u22 + uik[44] * u23 + uik[45] * u24 + uik[46] * u25 + uik[47] * u26 + uik[48] * u27;
320
321 wp[28] += uik[0] * u28 + uik[1] * u29 + uik[2] * u30 + uik[3] * u31 + uik[4] * u32 + uik[5] * u33 + uik[6] * u34;
322 wp[29] += uik[7] * u28 + uik[8] * u29 + uik[9] * u30 + uik[10] * u31 + uik[11] * u32 + uik[12] * u33 + uik[13] * u34;
323 wp[30] += uik[14] * u28 + uik[15] * u29 + uik[16] * u30 + uik[17] * u31 + uik[18] * u32 + uik[19] * u33 + uik[20] * u34;
324 wp[31] += uik[21] * u28 + uik[22] * u29 + uik[23] * u30 + uik[24] * u31 + uik[25] * u32 + uik[26] * u33 + uik[27] * u34;
325 wp[32] += uik[28] * u28 + uik[29] * u29 + uik[30] * u30 + uik[31] * u31 + uik[32] * u32 + uik[33] * u33 + uik[34] * u34;
326 wp[33] += uik[35] * u28 + uik[36] * u29 + uik[37] * u30 + uik[38] * u31 + uik[39] * u32 + uik[40] * u33 + uik[41] * u34;
327 wp[34] += uik[42] * u28 + uik[43] * u29 + uik[44] * u30 + uik[45] * u31 + uik[46] * u32 + uik[47] * u33 + uik[48] * u34;
328
329 wp[35] += uik[0] * u35 + uik[1] * u36 + uik[2] * u37 + uik[3] * u38 + uik[4] * u39 + uik[5] * u40 + uik[6] * u41;
330 wp[36] += uik[7] * u35 + uik[8] * u36 + uik[9] * u37 + uik[10] * u38 + uik[11] * u39 + uik[12] * u40 + uik[13] * u41;
331 wp[37] += uik[14] * u35 + uik[15] * u36 + uik[16] * u37 + uik[17] * u38 + uik[18] * u39 + uik[19] * u40 + uik[20] * u41;
332 wp[38] += uik[21] * u35 + uik[22] * u36 + uik[23] * u37 + uik[24] * u38 + uik[25] * u39 + uik[26] * u40 + uik[27] * u41;
333 wp[39] += uik[28] * u35 + uik[29] * u36 + uik[30] * u37 + uik[31] * u38 + uik[32] * u39 + uik[33] * u40 + uik[34] * u41;
334 wp[40] += uik[35] * u35 + uik[36] * u36 + uik[37] * u37 + uik[38] * u38 + uik[39] * u39 + uik[40] * u40 + uik[41] * u41;
335 wp[41] += uik[42] * u35 + uik[43] * u36 + uik[44] * u37 + uik[45] * u38 + uik[46] * u39 + uik[47] * u40 + uik[48] * u41;
336
337 wp[42] += uik[0] * u42 + uik[1] * u43 + uik[2] * u44 + uik[3] * u45 + uik[4] * u46 + uik[5] * u47 + uik[6] * u48;
338 wp[43] += uik[7] * u42 + uik[8] * u43 + uik[9] * u44 + uik[10] * u45 + uik[11] * u46 + uik[12] * u47 + uik[13] * u48;
339 wp[44] += uik[14] * u42 + uik[15] * u43 + uik[16] * u44 + uik[17] * u45 + uik[18] * u46 + uik[19] * u47 + uik[20] * u48;
340 wp[45] += uik[21] * u42 + uik[22] * u43 + uik[23] * u44 + uik[24] * u45 + uik[25] * u46 + uik[26] * u47 + uik[27] * u48;
341 wp[46] += uik[28] * u42 + uik[29] * u43 + uik[30] * u44 + uik[31] * u45 + uik[32] * u46 + uik[33] * u47 + uik[34] * u48;
342 wp[47] += uik[35] * u42 + uik[36] * u43 + uik[37] * u44 + uik[38] * u45 + uik[39] * u46 + uik[40] * u47 + uik[41] * u48;
343 wp[48] += uik[42] * u42 + uik[43] * u43 + uik[44] * u44 + uik[45] * u45 + uik[46] * u46 + uik[47] * u47 + uik[48] * u48;
344 }
345 PetscCall(PetscLogFlops(2.0 * 343.0 * (jmax - jmin)));
346
347 /* ... add i to row list for next nonzero entry */
348 il[i] = jmin; /* update il(i) in column k+1, ... mbs-1 */
349 j = bj[jmin];
350 jl[i] = jl[j];
351 jl[j] = i; /* update jl */
352 }
353 i = nexti;
354 }
355
356 /* save nonzero entries in k-th row of U ... */
357
358 /* invert diagonal block */
359 d = ba + k * 49;
360 PetscCall(PetscArraycpy(d, dk, 49));
361 PetscCall(PetscKernel_A_gets_inverse_A_7(d, shift, allowzeropivot, &zeropivotdetected));
362 if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
363
364 jmin = bi[k];
365 jmax = bi[k + 1];
366 if (jmin < jmax) {
367 for (j = jmin; j < jmax; j++) {
368 vj = bj[j]; /* block col. index of U */
369 u = ba + j * 49;
370 wp = w + vj * 49;
371 for (k1 = 0; k1 < 49; k1++) {
372 *u++ = *wp;
373 *wp++ = 0.0;
374 }
375 }
376
377 /* ... add k to row list for first nonzero entry in k-th row */
378 il[k] = jmin;
379 i = bj[jmin];
380 jl[k] = jl[i];
381 jl[i] = k;
382 }
383 }
384
385 PetscCall(PetscFree(w));
386 PetscCall(PetscFree2(il, jl));
387 PetscCall(PetscFree2(dk, uik));
388
389 C->ops->solve = MatSolve_SeqSBAIJ_7_NaturalOrdering_inplace;
390 C->ops->solvetranspose = MatSolve_SeqSBAIJ_7_NaturalOrdering_inplace;
391 C->ops->forwardsolve = MatForwardSolve_SeqSBAIJ_7_NaturalOrdering_inplace;
392 C->ops->backwardsolve = MatBackwardSolve_SeqSBAIJ_7_NaturalOrdering_inplace;
393 C->assembled = PETSC_TRUE;
394 C->preallocated = PETSC_TRUE;
395
396 PetscCall(PetscLogFlops(1.3333 * 343 * b->mbs)); /* from inverting diagonal blocks */
397 PetscFunctionReturn(PETSC_SUCCESS);
398 }
399