Actual source code: sbaijfact9.c
1: #include <../src/mat/impls/sbaij/seq/sbaij.h>
2: #include <petsc/private/kernels/blockinvert.h>
4: /* Version for when blocks are 6 by 6 */
5: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_6(Mat C, Mat A, const MatFactorInfo *info)
6: {
7: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data, *b = (Mat_SeqSBAIJ *)C->data;
8: IS perm = b->row;
9: const PetscInt *ai, *aj, *perm_ptr, mbs = a->mbs, *bi = b->i, *bj = b->j;
10: PetscInt i, j, *a2anew, k, k1, jmin, jmax, *jl, *il, vj, nexti, ili;
11: MatScalar *ba = b->a, *aa, *ap, *dk, *uik;
12: MatScalar *u, *d, *w, *wp, u0, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10, u11, u12;
13: MatScalar u13, u14, u15, u16, u17, u18, u19, u20, u21, u22, u23, u24, u25, u26, u27;
14: MatScalar u28, u29, u30, u31, u32, u33, u34, u35;
15: PetscReal shift = info->shiftamount;
16: PetscBool allowzeropivot, zeropivotdetected;
18: PetscFunctionBegin;
19: /* initialization */
20: allowzeropivot = PetscNot(A->erroriffailure);
21: PetscCall(PetscCalloc1(36 * mbs, &w));
22: PetscCall(PetscMalloc2(mbs, &il, mbs, &jl));
23: il[0] = 0;
24: for (i = 0; i < mbs; i++) jl[i] = mbs;
26: PetscCall(PetscMalloc2(36, &dk, 36, &uik));
27: PetscCall(ISGetIndices(perm, &perm_ptr));
29: /* check permutation */
30: if (!a->permute) {
31: ai = a->i;
32: aj = a->j;
33: aa = a->a;
34: } else {
35: ai = a->inew;
36: aj = a->jnew;
37: PetscCall(PetscMalloc1(36 * ai[mbs], &aa));
38: PetscCall(PetscArraycpy(aa, a->a, 36 * ai[mbs]));
39: PetscCall(PetscMalloc1(ai[mbs], &a2anew));
40: PetscCall(PetscArraycpy(a2anew, a->a2anew, ai[mbs]));
42: for (i = 0; i < mbs; i++) {
43: jmin = ai[i];
44: jmax = ai[i + 1];
45: for (j = jmin; j < jmax; j++) {
46: while (a2anew[j] != j) {
47: k = a2anew[j];
48: a2anew[j] = a2anew[k];
49: a2anew[k] = k;
50: for (k1 = 0; k1 < 36; k1++) {
51: dk[k1] = aa[k * 36 + k1];
52: aa[k * 36 + k1] = aa[j * 36 + k1];
53: aa[j * 36 + k1] = dk[k1];
54: }
55: }
56: /* transform column-oriented blocks that lie in the lower triangle to row-oriented blocks */
57: if (i > aj[j]) {
58: /* printf("change orientation, row: %d, col: %d\n",i,aj[j]); */
59: ap = aa + j * 36; /* ptr to the beginning of j-th block of aa */
60: for (k = 0; k < 36; k++) dk[k] = ap[k]; /* dk <- j-th block of aa */
61: for (k = 0; k < 6; k++) { /* j-th block of aa <- dk^T */
62: for (k1 = 0; k1 < 6; k1++) *ap++ = dk[k + 6 * k1];
63: }
64: }
65: }
66: }
67: PetscCall(PetscFree(a2anew));
68: }
70: /* for each row k */
71: for (k = 0; k < mbs; k++) {
72: /*initialize k-th row with elements nonzero in row perm(k) of A */
73: jmin = ai[perm_ptr[k]];
74: jmax = ai[perm_ptr[k] + 1];
75: if (jmin < jmax) {
76: ap = aa + jmin * 36;
77: for (j = jmin; j < jmax; j++) {
78: vj = perm_ptr[aj[j]]; /* block col. index */
79: wp = w + vj * 36;
80: for (i = 0; i < 36; i++) *wp++ = *ap++;
81: }
82: }
84: /* modify k-th row by adding in those rows i with U(i,k) != 0 */
85: PetscCall(PetscArraycpy(dk, w + k * 36, 36));
86: i = jl[k]; /* first row to be added to k_th row */
88: while (i < mbs) {
89: nexti = jl[i]; /* next row to be added to k_th row */
91: /* compute multiplier */
92: ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
94: /* uik = -inv(Di)*U_bar(i,k) */
95: d = ba + i * 36;
96: u = ba + ili * 36;
98: u0 = u[0];
99: u1 = u[1];
100: u2 = u[2];
101: u3 = u[3];
102: u4 = u[4];
103: u5 = u[5];
104: u6 = u[6];
105: u7 = u[7];
106: u8 = u[8];
107: u9 = u[9];
108: u10 = u[10];
109: u11 = u[11];
110: u12 = u[12];
111: u13 = u[13];
112: u14 = u[14];
113: u15 = u[15];
114: u16 = u[16];
115: u17 = u[17];
116: u18 = u[18];
117: u19 = u[19];
118: u20 = u[20];
119: u21 = u[21];
120: u22 = u[22];
121: u23 = u[23];
122: u24 = u[24];
123: u25 = u[25];
124: u26 = u[26];
125: u27 = u[27];
126: u28 = u[28];
127: u29 = u[29];
128: u30 = u[30];
129: u31 = u[31];
130: u32 = u[32];
131: u33 = u[33];
132: u34 = u[34];
133: u35 = u[35];
135: uik[0] = -(d[0] * u0 + d[6] * u1 + d[12] * u2 + d[18] * u3 + d[24] * u4 + d[30] * u5);
136: uik[1] = -(d[1] * u0 + d[7] * u1 + d[13] * u2 + d[19] * u3 + d[25] * u4 + d[31] * u5);
137: uik[2] = -(d[2] * u0 + d[8] * u1 + d[14] * u2 + d[20] * u3 + d[26] * u4 + d[32] * u5);
138: uik[3] = -(d[3] * u0 + d[9] * u1 + d[15] * u2 + d[21] * u3 + d[27] * u4 + d[33] * u5);
139: uik[4] = -(d[4] * u0 + d[10] * u1 + d[16] * u2 + d[22] * u3 + d[28] * u4 + d[34] * u5);
140: uik[5] = -(d[5] * u0 + d[11] * u1 + d[17] * u2 + d[23] * u3 + d[29] * u4 + d[35] * u5);
142: uik[6] = -(d[0] * u6 + d[6] * u7 + d[12] * u8 + d[18] * u9 + d[24] * u10 + d[30] * u11);
143: uik[7] = -(d[1] * u6 + d[7] * u7 + d[13] * u8 + d[19] * u9 + d[25] * u10 + d[31] * u11);
144: uik[8] = -(d[2] * u6 + d[8] * u7 + d[14] * u8 + d[20] * u9 + d[26] * u10 + d[32] * u11);
145: uik[9] = -(d[3] * u6 + d[9] * u7 + d[15] * u8 + d[21] * u9 + d[27] * u10 + d[33] * u11);
146: uik[10] = -(d[4] * u6 + d[10] * u7 + d[16] * u8 + d[22] * u9 + d[28] * u10 + d[34] * u11);
147: uik[11] = -(d[5] * u6 + d[11] * u7 + d[17] * u8 + d[23] * u9 + d[29] * u10 + d[35] * u11);
149: uik[12] = -(d[0] * u12 + d[6] * u13 + d[12] * u14 + d[18] * u15 + d[24] * u16 + d[30] * u17);
150: uik[13] = -(d[1] * u12 + d[7] * u13 + d[13] * u14 + d[19] * u15 + d[25] * u16 + d[31] * u17);
151: uik[14] = -(d[2] * u12 + d[8] * u13 + d[14] * u14 + d[20] * u15 + d[26] * u16 + d[32] * u17);
152: uik[15] = -(d[3] * u12 + d[9] * u13 + d[15] * u14 + d[21] * u15 + d[27] * u16 + d[33] * u17);
153: uik[16] = -(d[4] * u12 + d[10] * u13 + d[16] * u14 + d[22] * u15 + d[28] * u16 + d[34] * u17);
154: uik[17] = -(d[5] * u12 + d[11] * u13 + d[17] * u14 + d[23] * u15 + d[29] * u16 + d[35] * u17);
156: uik[18] = -(d[0] * u18 + d[6] * u19 + d[12] * u20 + d[18] * u21 + d[24] * u22 + d[30] * u23);
157: uik[19] = -(d[1] * u18 + d[7] * u19 + d[13] * u20 + d[19] * u21 + d[25] * u22 + d[31] * u23);
158: uik[20] = -(d[2] * u18 + d[8] * u19 + d[14] * u20 + d[20] * u21 + d[26] * u22 + d[32] * u23);
159: uik[21] = -(d[3] * u18 + d[9] * u19 + d[15] * u20 + d[21] * u21 + d[27] * u22 + d[33] * u23);
160: uik[22] = -(d[4] * u18 + d[10] * u19 + d[16] * u20 + d[22] * u21 + d[28] * u22 + d[34] * u23);
161: uik[23] = -(d[5] * u18 + d[11] * u19 + d[17] * u20 + d[23] * u21 + d[29] * u22 + d[35] * u23);
163: uik[24] = -(d[0] * u24 + d[6] * u25 + d[12] * u26 + d[18] * u27 + d[24] * u28 + d[30] * u29);
164: uik[25] = -(d[1] * u24 + d[7] * u25 + d[13] * u26 + d[19] * u27 + d[25] * u28 + d[31] * u29);
165: uik[26] = -(d[2] * u24 + d[8] * u25 + d[14] * u26 + d[20] * u27 + d[26] * u28 + d[32] * u29);
166: uik[27] = -(d[3] * u24 + d[9] * u25 + d[15] * u26 + d[21] * u27 + d[27] * u28 + d[33] * u29);
167: uik[28] = -(d[4] * u24 + d[10] * u25 + d[16] * u26 + d[22] * u27 + d[28] * u28 + d[34] * u29);
168: uik[29] = -(d[5] * u24 + d[11] * u25 + d[17] * u26 + d[23] * u27 + d[29] * u28 + d[35] * u29);
170: uik[30] = -(d[0] * u30 + d[6] * u31 + d[12] * u32 + d[18] * u33 + d[24] * u34 + d[30] * u35);
171: uik[31] = -(d[1] * u30 + d[7] * u31 + d[13] * u32 + d[19] * u33 + d[25] * u34 + d[31] * u35);
172: uik[32] = -(d[2] * u30 + d[8] * u31 + d[14] * u32 + d[20] * u33 + d[26] * u34 + d[32] * u35);
173: uik[33] = -(d[3] * u30 + d[9] * u31 + d[15] * u32 + d[21] * u33 + d[27] * u34 + d[33] * u35);
174: uik[34] = -(d[4] * u30 + d[10] * u31 + d[16] * u32 + d[22] * u33 + d[28] * u34 + d[34] * u35);
175: uik[35] = -(d[5] * u30 + d[11] * u31 + d[17] * u32 + d[23] * u33 + d[29] * u34 + d[35] * u35);
177: /* update D(k) += -U(i,k)^T * U_bar(i,k) */
178: dk[0] += uik[0] * u0 + uik[1] * u1 + uik[2] * u2 + uik[3] * u3 + uik[4] * u4 + uik[5] * u5;
179: dk[1] += uik[6] * u0 + uik[7] * u1 + uik[8] * u2 + uik[9] * u3 + uik[10] * u4 + uik[11] * u5;
180: dk[2] += uik[12] * u0 + uik[13] * u1 + uik[14] * u2 + uik[15] * u3 + uik[16] * u4 + uik[17] * u5;
181: dk[3] += uik[18] * u0 + uik[19] * u1 + uik[20] * u2 + uik[21] * u3 + uik[22] * u4 + uik[23] * u5;
182: dk[4] += uik[24] * u0 + uik[25] * u1 + uik[26] * u2 + uik[27] * u3 + uik[28] * u4 + uik[29] * u5;
183: dk[5] += uik[30] * u0 + uik[31] * u1 + uik[32] * u2 + uik[33] * u3 + uik[34] * u4 + uik[35] * u5;
185: dk[6] += uik[0] * u6 + uik[1] * u7 + uik[2] * u8 + uik[3] * u9 + uik[4] * u10 + uik[5] * u11;
186: dk[7] += uik[6] * u6 + uik[7] * u7 + uik[8] * u8 + uik[9] * u9 + uik[10] * u10 + uik[11] * u11;
187: dk[8] += uik[12] * u6 + uik[13] * u7 + uik[14] * u8 + uik[15] * u9 + uik[16] * u10 + uik[17] * u11;
188: dk[9] += uik[18] * u6 + uik[19] * u7 + uik[20] * u8 + uik[21] * u9 + uik[22] * u10 + uik[23] * u11;
189: dk[10] += uik[24] * u6 + uik[25] * u7 + uik[26] * u8 + uik[27] * u9 + uik[28] * u10 + uik[29] * u11;
190: dk[11] += uik[30] * u6 + uik[31] * u7 + uik[32] * u8 + uik[33] * u9 + uik[34] * u10 + uik[35] * u11;
192: dk[12] += uik[0] * u12 + uik[1] * u13 + uik[2] * u14 + uik[3] * u15 + uik[4] * u16 + uik[5] * u17;
193: dk[13] += uik[6] * u12 + uik[7] * u13 + uik[8] * u14 + uik[9] * u15 + uik[10] * u16 + uik[11] * u17;
194: dk[14] += uik[12] * u12 + uik[13] * u13 + uik[14] * u14 + uik[15] * u15 + uik[16] * u16 + uik[17] * u17;
195: dk[15] += uik[18] * u12 + uik[19] * u13 + uik[20] * u14 + uik[21] * u15 + uik[22] * u16 + uik[23] * u17;
196: dk[16] += uik[24] * u12 + uik[25] * u13 + uik[26] * u14 + uik[27] * u15 + uik[28] * u16 + uik[29] * u17;
197: dk[17] += uik[30] * u12 + uik[31] * u13 + uik[32] * u14 + uik[33] * u15 + uik[34] * u16 + uik[35] * u17;
199: dk[18] += uik[0] * u18 + uik[1] * u19 + uik[2] * u20 + uik[3] * u21 + uik[4] * u22 + uik[5] * u23;
200: dk[19] += uik[6] * u18 + uik[7] * u19 + uik[8] * u20 + uik[9] * u21 + uik[10] * u22 + uik[11] * u23;
201: dk[20] += uik[12] * u18 + uik[13] * u19 + uik[14] * u20 + uik[15] * u21 + uik[16] * u22 + uik[17] * u23;
202: dk[21] += uik[18] * u18 + uik[19] * u19 + uik[20] * u20 + uik[21] * u21 + uik[22] * u22 + uik[23] * u23;
203: dk[22] += uik[24] * u18 + uik[25] * u19 + uik[26] * u20 + uik[27] * u21 + uik[28] * u22 + uik[29] * u23;
204: dk[23] += uik[30] * u18 + uik[31] * u19 + uik[32] * u20 + uik[33] * u21 + uik[34] * u22 + uik[35] * u23;
206: dk[24] += uik[0] * u24 + uik[1] * u25 + uik[2] * u26 + uik[3] * u27 + uik[4] * u28 + uik[5] * u29;
207: dk[25] += uik[6] * u24 + uik[7] * u25 + uik[8] * u26 + uik[9] * u27 + uik[10] * u28 + uik[11] * u29;
208: dk[26] += uik[12] * u24 + uik[13] * u25 + uik[14] * u26 + uik[15] * u27 + uik[16] * u28 + uik[17] * u29;
209: dk[27] += uik[18] * u24 + uik[19] * u25 + uik[20] * u26 + uik[21] * u27 + uik[22] * u28 + uik[23] * u29;
210: dk[28] += uik[24] * u24 + uik[25] * u25 + uik[26] * u26 + uik[27] * u27 + uik[28] * u28 + uik[29] * u29;
211: dk[29] += uik[30] * u24 + uik[31] * u25 + uik[32] * u26 + uik[33] * u27 + uik[34] * u28 + uik[35] * u29;
213: dk[30] += uik[0] * u30 + uik[1] * u31 + uik[2] * u32 + uik[3] * u33 + uik[4] * u34 + uik[5] * u35;
214: dk[31] += uik[6] * u30 + uik[7] * u31 + uik[8] * u32 + uik[9] * u33 + uik[10] * u34 + uik[11] * u35;
215: dk[32] += uik[12] * u30 + uik[13] * u31 + uik[14] * u32 + uik[15] * u33 + uik[16] * u34 + uik[17] * u35;
216: dk[33] += uik[18] * u30 + uik[19] * u31 + uik[20] * u32 + uik[21] * u33 + uik[22] * u34 + uik[23] * u35;
217: dk[34] += uik[24] * u30 + uik[25] * u31 + uik[26] * u32 + uik[27] * u33 + uik[28] * u34 + uik[29] * u35;
218: dk[35] += uik[30] * u30 + uik[31] * u31 + uik[32] * u32 + uik[33] * u33 + uik[34] * u34 + uik[35] * u35;
220: PetscCall(PetscLogFlops(216.0 * 4.0));
222: /* update -U(i,k) */
223: PetscCall(PetscArraycpy(ba + ili * 36, uik, 36));
225: /* add multiple of row i to k-th row ... */
226: jmin = ili + 1;
227: jmax = bi[i + 1];
228: if (jmin < jmax) {
229: for (j = jmin; j < jmax; j++) {
230: /* w += -U(i,k)^T * U_bar(i,j) */
231: wp = w + bj[j] * 36;
232: u = ba + j * 36;
234: u0 = u[0];
235: u1 = u[1];
236: u2 = u[2];
237: u3 = u[3];
238: u4 = u[4];
239: u5 = u[5];
240: u6 = u[6];
241: u7 = u[7];
242: u8 = u[8];
243: u9 = u[9];
244: u10 = u[10];
245: u11 = u[11];
246: u12 = u[12];
247: u13 = u[13];
248: u14 = u[14];
249: u15 = u[15];
250: u16 = u[16];
251: u17 = u[17];
252: u18 = u[18];
253: u19 = u[19];
254: u20 = u[20];
255: u21 = u[21];
256: u22 = u[22];
257: u23 = u[23];
258: u24 = u[24];
259: u25 = u[25];
260: u26 = u[26];
261: u27 = u[27];
262: u28 = u[28];
263: u29 = u[29];
264: u30 = u[30];
265: u31 = u[31];
266: u32 = u[32];
267: u33 = u[33];
268: u34 = u[34];
269: u35 = u[35];
271: wp[0] += uik[0] * u0 + uik[1] * u1 + uik[2] * u2 + uik[3] * u3 + uik[4] * u4 + uik[5] * u5;
272: wp[1] += uik[6] * u0 + uik[7] * u1 + uik[8] * u2 + uik[9] * u3 + uik[10] * u4 + uik[11] * u5;
273: wp[2] += uik[12] * u0 + uik[13] * u1 + uik[14] * u2 + uik[15] * u3 + uik[16] * u4 + uik[17] * u5;
274: wp[3] += uik[18] * u0 + uik[19] * u1 + uik[20] * u2 + uik[21] * u3 + uik[22] * u4 + uik[23] * u5;
275: wp[4] += uik[24] * u0 + uik[25] * u1 + uik[26] * u2 + uik[27] * u3 + uik[28] * u4 + uik[29] * u5;
276: wp[5] += uik[30] * u0 + uik[31] * u1 + uik[32] * u2 + uik[33] * u3 + uik[34] * u4 + uik[35] * u5;
278: wp[6] += uik[0] * u6 + uik[1] * u7 + uik[2] * u8 + uik[3] * u9 + uik[4] * u10 + uik[5] * u11;
279: wp[7] += uik[6] * u6 + uik[7] * u7 + uik[8] * u8 + uik[9] * u9 + uik[10] * u10 + uik[11] * u11;
280: wp[8] += uik[12] * u6 + uik[13] * u7 + uik[14] * u8 + uik[15] * u9 + uik[16] * u10 + uik[17] * u11;
281: wp[9] += uik[18] * u6 + uik[19] * u7 + uik[20] * u8 + uik[21] * u9 + uik[22] * u10 + uik[23] * u11;
282: wp[10] += uik[24] * u6 + uik[25] * u7 + uik[26] * u8 + uik[27] * u9 + uik[28] * u10 + uik[29] * u11;
283: wp[11] += uik[30] * u6 + uik[31] * u7 + uik[32] * u8 + uik[33] * u9 + uik[34] * u10 + uik[35] * u11;
285: wp[12] += uik[0] * u12 + uik[1] * u13 + uik[2] * u14 + uik[3] * u15 + uik[4] * u16 + uik[5] * u17;
286: wp[13] += uik[6] * u12 + uik[7] * u13 + uik[8] * u14 + uik[9] * u15 + uik[10] * u16 + uik[11] * u17;
287: wp[14] += uik[12] * u12 + uik[13] * u13 + uik[14] * u14 + uik[15] * u15 + uik[16] * u16 + uik[17] * u17;
288: wp[15] += uik[18] * u12 + uik[19] * u13 + uik[20] * u14 + uik[21] * u15 + uik[22] * u16 + uik[23] * u17;
289: wp[16] += uik[24] * u12 + uik[25] * u13 + uik[26] * u14 + uik[27] * u15 + uik[28] * u16 + uik[29] * u17;
290: wp[17] += uik[30] * u12 + uik[31] * u13 + uik[32] * u14 + uik[33] * u15 + uik[34] * u16 + uik[35] * u17;
292: wp[18] += uik[0] * u18 + uik[1] * u19 + uik[2] * u20 + uik[3] * u21 + uik[4] * u22 + uik[5] * u23;
293: wp[19] += uik[6] * u18 + uik[7] * u19 + uik[8] * u20 + uik[9] * u21 + uik[10] * u22 + uik[11] * u23;
294: wp[20] += uik[12] * u18 + uik[13] * u19 + uik[14] * u20 + uik[15] * u21 + uik[16] * u22 + uik[17] * u23;
295: wp[21] += uik[18] * u18 + uik[19] * u19 + uik[20] * u20 + uik[21] * u21 + uik[22] * u22 + uik[23] * u23;
296: wp[22] += uik[24] * u18 + uik[25] * u19 + uik[26] * u20 + uik[27] * u21 + uik[28] * u22 + uik[29] * u23;
297: wp[23] += uik[30] * u18 + uik[31] * u19 + uik[32] * u20 + uik[33] * u21 + uik[34] * u22 + uik[35] * u23;
299: wp[24] += uik[0] * u24 + uik[1] * u25 + uik[2] * u26 + uik[3] * u27 + uik[4] * u28 + uik[5] * u29;
300: wp[25] += uik[6] * u24 + uik[7] * u25 + uik[8] * u26 + uik[9] * u27 + uik[10] * u28 + uik[11] * u29;
301: wp[26] += uik[12] * u24 + uik[13] * u25 + uik[14] * u26 + uik[15] * u27 + uik[16] * u28 + uik[17] * u29;
302: wp[27] += uik[18] * u24 + uik[19] * u25 + uik[20] * u26 + uik[21] * u27 + uik[22] * u28 + uik[23] * u29;
303: wp[28] += uik[24] * u24 + uik[25] * u25 + uik[26] * u26 + uik[27] * u27 + uik[28] * u28 + uik[29] * u29;
304: wp[29] += uik[30] * u24 + uik[31] * u25 + uik[32] * u26 + uik[33] * u27 + uik[34] * u28 + uik[35] * u29;
306: wp[30] += uik[0] * u30 + uik[1] * u31 + uik[2] * u32 + uik[3] * u33 + uik[4] * u34 + uik[5] * u35;
307: wp[31] += uik[6] * u30 + uik[7] * u31 + uik[8] * u32 + uik[9] * u33 + uik[10] * u34 + uik[11] * u35;
308: wp[32] += uik[12] * u30 + uik[13] * u31 + uik[14] * u32 + uik[15] * u33 + uik[16] * u34 + uik[17] * u35;
309: wp[33] += uik[18] * u30 + uik[19] * u31 + uik[20] * u32 + uik[21] * u33 + uik[22] * u34 + uik[23] * u35;
310: wp[34] += uik[24] * u30 + uik[25] * u31 + uik[26] * u32 + uik[27] * u33 + uik[28] * u34 + uik[29] * u35;
311: wp[35] += uik[30] * u30 + uik[31] * u31 + uik[32] * u32 + uik[33] * u33 + uik[34] * u34 + uik[35] * u35;
312: }
313: PetscCall(PetscLogFlops(2.0 * 216.0 * (jmax - jmin)));
315: /* ... add i to row list for next nonzero entry */
316: il[i] = jmin; /* update il(i) in column k+1, ... mbs-1 */
317: j = bj[jmin];
318: jl[i] = jl[j];
319: jl[j] = i; /* update jl */
320: }
321: i = nexti;
322: }
324: /* save nonzero entries in k-th row of U ... */
326: /* invert diagonal block */
327: d = ba + k * 36;
328: PetscCall(PetscArraycpy(d, dk, 36));
329: PetscCall(PetscKernel_A_gets_inverse_A_6(d, shift, allowzeropivot, &zeropivotdetected));
330: if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
332: jmin = bi[k];
333: jmax = bi[k + 1];
334: if (jmin < jmax) {
335: for (j = jmin; j < jmax; j++) {
336: vj = bj[j]; /* block col. index of U */
337: u = ba + j * 36;
338: wp = w + vj * 36;
339: for (k1 = 0; k1 < 36; k1++) {
340: *u++ = *wp;
341: *wp++ = 0.0;
342: }
343: }
345: /* ... add k to row list for first nonzero entry in k-th row */
346: il[k] = jmin;
347: i = bj[jmin];
348: jl[k] = jl[i];
349: jl[i] = k;
350: }
351: }
353: PetscCall(PetscFree(w));
354: PetscCall(PetscFree2(il, jl));
355: PetscCall(PetscFree2(dk, uik));
356: if (a->permute) PetscCall(PetscFree(aa));
358: PetscCall(ISRestoreIndices(perm, &perm_ptr));
360: C->ops->solve = MatSolve_SeqSBAIJ_6_inplace;
361: C->ops->solvetranspose = MatSolve_SeqSBAIJ_6_inplace;
362: C->assembled = PETSC_TRUE;
363: C->preallocated = PETSC_TRUE;
365: PetscCall(PetscLogFlops(1.3333 * 216 * b->mbs)); /* from inverting diagonal blocks */
366: PetscFunctionReturn(PETSC_SUCCESS);
367: }