Actual source code: sbaijfact9.c


  2: #include <../src/mat/impls/sbaij/seq/sbaij.h>
  3: #include <petsc/private/kernels/blockinvert.h>

  5: /* Version for when blocks are 6 by 6 */
  6: PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_6(Mat C, Mat A, const MatFactorInfo *info)
  7: {
  8:   Mat_SeqSBAIJ   *a = (Mat_SeqSBAIJ *)A->data, *b = (Mat_SeqSBAIJ *)C->data;
  9:   IS              perm = b->row;
 10:   const PetscInt *ai, *aj, *perm_ptr, mbs = a->mbs, *bi = b->i, *bj = b->j;
 11:   PetscInt        i, j, *a2anew, 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;
 16:   PetscReal       shift = info->shiftamount;
 17:   PetscBool       allowzeropivot, zeropivotdetected;

 19:   PetscFunctionBegin;
 20:   /* initialization */
 21:   allowzeropivot = PetscNot(A->erroriffailure);
 22:   PetscCall(PetscCalloc1(36 * mbs, &w));
 23:   PetscCall(PetscMalloc2(mbs, &il, mbs, &jl));
 24:   il[0] = 0;
 25:   for (i = 0; i < mbs; i++) jl[i] = mbs;

 27:   PetscCall(PetscMalloc2(36, &dk, 36, &uik));
 28:   PetscCall(ISGetIndices(perm, &perm_ptr));

 30:   /* check permutation */
 31:   if (!a->permute) {
 32:     ai = a->i;
 33:     aj = a->j;
 34:     aa = a->a;
 35:   } else {
 36:     ai = a->inew;
 37:     aj = a->jnew;
 38:     PetscCall(PetscMalloc1(36 * ai[mbs], &aa));
 39:     PetscCall(PetscArraycpy(aa, a->a, 36 * ai[mbs]));
 40:     PetscCall(PetscMalloc1(ai[mbs], &a2anew));
 41:     PetscCall(PetscArraycpy(a2anew, a->a2anew, ai[mbs]));

 43:     for (i = 0; i < mbs; i++) {
 44:       jmin = ai[i];
 45:       jmax = ai[i + 1];
 46:       for (j = jmin; j < jmax; j++) {
 47:         while (a2anew[j] != j) {
 48:           k         = a2anew[j];
 49:           a2anew[j] = a2anew[k];
 50:           a2anew[k] = k;
 51:           for (k1 = 0; k1 < 36; k1++) {
 52:             dk[k1]          = aa[k * 36 + k1];
 53:             aa[k * 36 + k1] = aa[j * 36 + k1];
 54:             aa[j * 36 + k1] = dk[k1];
 55:           }
 56:         }
 57:         /* transform columnoriented blocks that lie in the lower triangle to roworiented blocks */
 58:         if (i > aj[j]) {
 59:           /* printf("change orientation, row: %d, col: %d\n",i,aj[j]); */
 60:           ap = aa + j * 36;                       /* ptr to the beginning of j-th block of aa */
 61:           for (k = 0; k < 36; k++) dk[k] = ap[k]; /* dk <- j-th block of aa */
 62:           for (k = 0; k < 6; k++) {               /* j-th block of aa <- dk^T */
 63:             for (k1 = 0; k1 < 6; k1++) *ap++ = dk[k + 6 * k1];
 64:           }
 65:         }
 66:       }
 67:     }
 68:     PetscCall(PetscFree(a2anew));
 69:   }

 71:   /* for each row k */
 72:   for (k = 0; k < mbs; k++) {
 73:     /*initialize k-th row with elements nonzero in row perm(k) of A */
 74:     jmin = ai[perm_ptr[k]];
 75:     jmax = ai[perm_ptr[k] + 1];
 76:     if (jmin < jmax) {
 77:       ap = aa + jmin * 36;
 78:       for (j = jmin; j < jmax; j++) {
 79:         vj = perm_ptr[aj[j]]; /* block col. index */
 80:         wp = w + vj * 36;
 81:         for (i = 0; i < 36; i++) *wp++ = *ap++;
 82:       }
 83:     }

 85:     /* modify k-th row by adding in those rows i with U(i,k) != 0 */
 86:     PetscCall(PetscArraycpy(dk, w + k * 36, 36));
 87:     i = jl[k]; /* first row to be added to k_th row  */

 89:     while (i < mbs) {
 90:       nexti = jl[i]; /* next row to be added to k_th row */

 92:       /* compute multiplier */
 93:       ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */

 95:       /* uik = -inv(Di)*U_bar(i,k) */
 96:       d = ba + i * 36;
 97:       u = ba + ili * 36;

 99:       u0  = u[0];
100:       u1  = u[1];
101:       u2  = u[2];
102:       u3  = u[3];
103:       u4  = u[4];
104:       u5  = u[5];
105:       u6  = u[6];
106:       u7  = u[7];
107:       u8  = u[8];
108:       u9  = u[9];
109:       u10 = u[10];
110:       u11 = u[11];
111:       u12 = u[12];
112:       u13 = u[13];
113:       u14 = u[14];
114:       u15 = u[15];
115:       u16 = u[16];
116:       u17 = u[17];
117:       u18 = u[18];
118:       u19 = u[19];
119:       u20 = u[20];
120:       u21 = u[21];
121:       u22 = u[22];
122:       u23 = u[23];
123:       u24 = u[24];
124:       u25 = u[25];
125:       u26 = u[26];
126:       u27 = u[27];
127:       u28 = u[28];
128:       u29 = u[29];
129:       u30 = u[30];
130:       u31 = u[31];
131:       u32 = u[32];
132:       u33 = u[33];
133:       u34 = u[34];
134:       u35 = u[35];

136:       uik[0] = -(d[0] * u0 + d[6] * u1 + d[12] * u2 + d[18] * u3 + d[24] * u4 + d[30] * u5);
137:       uik[1] = -(d[1] * u0 + d[7] * u1 + d[13] * u2 + d[19] * u3 + d[25] * u4 + d[31] * u5);
138:       uik[2] = -(d[2] * u0 + d[8] * u1 + d[14] * u2 + d[20] * u3 + d[26] * u4 + d[32] * u5);
139:       uik[3] = -(d[3] * u0 + d[9] * u1 + d[15] * u2 + d[21] * u3 + d[27] * u4 + d[33] * u5);
140:       uik[4] = -(d[4] * u0 + d[10] * u1 + d[16] * u2 + d[22] * u3 + d[28] * u4 + d[34] * u5);
141:       uik[5] = -(d[5] * u0 + d[11] * u1 + d[17] * u2 + d[23] * u3 + d[29] * u4 + d[35] * u5);

143:       uik[6]  = -(d[0] * u6 + d[6] * u7 + d[12] * u8 + d[18] * u9 + d[24] * u10 + d[30] * u11);
144:       uik[7]  = -(d[1] * u6 + d[7] * u7 + d[13] * u8 + d[19] * u9 + d[25] * u10 + d[31] * u11);
145:       uik[8]  = -(d[2] * u6 + d[8] * u7 + d[14] * u8 + d[20] * u9 + d[26] * u10 + d[32] * u11);
146:       uik[9]  = -(d[3] * u6 + d[9] * u7 + d[15] * u8 + d[21] * u9 + d[27] * u10 + d[33] * u11);
147:       uik[10] = -(d[4] * u6 + d[10] * u7 + d[16] * u8 + d[22] * u9 + d[28] * u10 + d[34] * u11);
148:       uik[11] = -(d[5] * u6 + d[11] * u7 + d[17] * u8 + d[23] * u9 + d[29] * u10 + d[35] * u11);

150:       uik[12] = -(d[0] * u12 + d[6] * u13 + d[12] * u14 + d[18] * u15 + d[24] * u16 + d[30] * u17);
151:       uik[13] = -(d[1] * u12 + d[7] * u13 + d[13] * u14 + d[19] * u15 + d[25] * u16 + d[31] * u17);
152:       uik[14] = -(d[2] * u12 + d[8] * u13 + d[14] * u14 + d[20] * u15 + d[26] * u16 + d[32] * u17);
153:       uik[15] = -(d[3] * u12 + d[9] * u13 + d[15] * u14 + d[21] * u15 + d[27] * u16 + d[33] * u17);
154:       uik[16] = -(d[4] * u12 + d[10] * u13 + d[16] * u14 + d[22] * u15 + d[28] * u16 + d[34] * u17);
155:       uik[17] = -(d[5] * u12 + d[11] * u13 + d[17] * u14 + d[23] * u15 + d[29] * u16 + d[35] * u17);

157:       uik[18] = -(d[0] * u18 + d[6] * u19 + d[12] * u20 + d[18] * u21 + d[24] * u22 + d[30] * u23);
158:       uik[19] = -(d[1] * u18 + d[7] * u19 + d[13] * u20 + d[19] * u21 + d[25] * u22 + d[31] * u23);
159:       uik[20] = -(d[2] * u18 + d[8] * u19 + d[14] * u20 + d[20] * u21 + d[26] * u22 + d[32] * u23);
160:       uik[21] = -(d[3] * u18 + d[9] * u19 + d[15] * u20 + d[21] * u21 + d[27] * u22 + d[33] * u23);
161:       uik[22] = -(d[4] * u18 + d[10] * u19 + d[16] * u20 + d[22] * u21 + d[28] * u22 + d[34] * u23);
162:       uik[23] = -(d[5] * u18 + d[11] * u19 + d[17] * u20 + d[23] * u21 + d[29] * u22 + d[35] * u23);

164:       uik[24] = -(d[0] * u24 + d[6] * u25 + d[12] * u26 + d[18] * u27 + d[24] * u28 + d[30] * u29);
165:       uik[25] = -(d[1] * u24 + d[7] * u25 + d[13] * u26 + d[19] * u27 + d[25] * u28 + d[31] * u29);
166:       uik[26] = -(d[2] * u24 + d[8] * u25 + d[14] * u26 + d[20] * u27 + d[26] * u28 + d[32] * u29);
167:       uik[27] = -(d[3] * u24 + d[9] * u25 + d[15] * u26 + d[21] * u27 + d[27] * u28 + d[33] * u29);
168:       uik[28] = -(d[4] * u24 + d[10] * u25 + d[16] * u26 + d[22] * u27 + d[28] * u28 + d[34] * u29);
169:       uik[29] = -(d[5] * u24 + d[11] * u25 + d[17] * u26 + d[23] * u27 + d[29] * u28 + d[35] * u29);

171:       uik[30] = -(d[0] * u30 + d[6] * u31 + d[12] * u32 + d[18] * u33 + d[24] * u34 + d[30] * u35);
172:       uik[31] = -(d[1] * u30 + d[7] * u31 + d[13] * u32 + d[19] * u33 + d[25] * u34 + d[31] * u35);
173:       uik[32] = -(d[2] * u30 + d[8] * u31 + d[14] * u32 + d[20] * u33 + d[26] * u34 + d[32] * u35);
174:       uik[33] = -(d[3] * u30 + d[9] * u31 + d[15] * u32 + d[21] * u33 + d[27] * u34 + d[33] * u35);
175:       uik[34] = -(d[4] * u30 + d[10] * u31 + d[16] * u32 + d[22] * u33 + d[28] * u34 + d[34] * u35);
176:       uik[35] = -(d[5] * u30 + d[11] * u31 + d[17] * u32 + d[23] * u33 + d[29] * u34 + d[35] * u35);

178:       /* update D(k) += -U(i,k)^T * U_bar(i,k) */
179:       dk[0] += uik[0] * u0 + uik[1] * u1 + uik[2] * u2 + uik[3] * u3 + uik[4] * u4 + uik[5] * u5;
180:       dk[1] += uik[6] * u0 + uik[7] * u1 + uik[8] * u2 + uik[9] * u3 + uik[10] * u4 + uik[11] * u5;
181:       dk[2] += uik[12] * u0 + uik[13] * u1 + uik[14] * u2 + uik[15] * u3 + uik[16] * u4 + uik[17] * u5;
182:       dk[3] += uik[18] * u0 + uik[19] * u1 + uik[20] * u2 + uik[21] * u3 + uik[22] * u4 + uik[23] * u5;
183:       dk[4] += uik[24] * u0 + uik[25] * u1 + uik[26] * u2 + uik[27] * u3 + uik[28] * u4 + uik[29] * u5;
184:       dk[5] += uik[30] * u0 + uik[31] * u1 + uik[32] * u2 + uik[33] * u3 + uik[34] * u4 + uik[35] * u5;

186:       dk[6] += uik[0] * u6 + uik[1] * u7 + uik[2] * u8 + uik[3] * u9 + uik[4] * u10 + uik[5] * u11;
187:       dk[7] += uik[6] * u6 + uik[7] * u7 + uik[8] * u8 + uik[9] * u9 + uik[10] * u10 + uik[11] * u11;
188:       dk[8] += uik[12] * u6 + uik[13] * u7 + uik[14] * u8 + uik[15] * u9 + uik[16] * u10 + uik[17] * u11;
189:       dk[9] += uik[18] * u6 + uik[19] * u7 + uik[20] * u8 + uik[21] * u9 + uik[22] * u10 + uik[23] * u11;
190:       dk[10] += uik[24] * u6 + uik[25] * u7 + uik[26] * u8 + uik[27] * u9 + uik[28] * u10 + uik[29] * u11;
191:       dk[11] += uik[30] * u6 + uik[31] * u7 + uik[32] * u8 + uik[33] * u9 + uik[34] * u10 + uik[35] * u11;

193:       dk[12] += uik[0] * u12 + uik[1] * u13 + uik[2] * u14 + uik[3] * u15 + uik[4] * u16 + uik[5] * u17;
194:       dk[13] += uik[6] * u12 + uik[7] * u13 + uik[8] * u14 + uik[9] * u15 + uik[10] * u16 + uik[11] * u17;
195:       dk[14] += uik[12] * u12 + uik[13] * u13 + uik[14] * u14 + uik[15] * u15 + uik[16] * u16 + uik[17] * u17;
196:       dk[15] += uik[18] * u12 + uik[19] * u13 + uik[20] * u14 + uik[21] * u15 + uik[22] * u16 + uik[23] * u17;
197:       dk[16] += uik[24] * u12 + uik[25] * u13 + uik[26] * u14 + uik[27] * u15 + uik[28] * u16 + uik[29] * u17;
198:       dk[17] += uik[30] * u12 + uik[31] * u13 + uik[32] * u14 + uik[33] * u15 + uik[34] * u16 + uik[35] * u17;

200:       dk[18] += uik[0] * u18 + uik[1] * u19 + uik[2] * u20 + uik[3] * u21 + uik[4] * u22 + uik[5] * u23;
201:       dk[19] += uik[6] * u18 + uik[7] * u19 + uik[8] * u20 + uik[9] * u21 + uik[10] * u22 + uik[11] * u23;
202:       dk[20] += uik[12] * u18 + uik[13] * u19 + uik[14] * u20 + uik[15] * u21 + uik[16] * u22 + uik[17] * u23;
203:       dk[21] += uik[18] * u18 + uik[19] * u19 + uik[20] * u20 + uik[21] * u21 + uik[22] * u22 + uik[23] * u23;
204:       dk[22] += uik[24] * u18 + uik[25] * u19 + uik[26] * u20 + uik[27] * u21 + uik[28] * u22 + uik[29] * u23;
205:       dk[23] += uik[30] * u18 + uik[31] * u19 + uik[32] * u20 + uik[33] * u21 + uik[34] * u22 + uik[35] * u23;

207:       dk[24] += uik[0] * u24 + uik[1] * u25 + uik[2] * u26 + uik[3] * u27 + uik[4] * u28 + uik[5] * u29;
208:       dk[25] += uik[6] * u24 + uik[7] * u25 + uik[8] * u26 + uik[9] * u27 + uik[10] * u28 + uik[11] * u29;
209:       dk[26] += uik[12] * u24 + uik[13] * u25 + uik[14] * u26 + uik[15] * u27 + uik[16] * u28 + uik[17] * u29;
210:       dk[27] += uik[18] * u24 + uik[19] * u25 + uik[20] * u26 + uik[21] * u27 + uik[22] * u28 + uik[23] * u29;
211:       dk[28] += uik[24] * u24 + uik[25] * u25 + uik[26] * u26 + uik[27] * u27 + uik[28] * u28 + uik[29] * u29;
212:       dk[29] += uik[30] * u24 + uik[31] * u25 + uik[32] * u26 + uik[33] * u27 + uik[34] * u28 + uik[35] * u29;

214:       dk[30] += uik[0] * u30 + uik[1] * u31 + uik[2] * u32 + uik[3] * u33 + uik[4] * u34 + uik[5] * u35;
215:       dk[31] += uik[6] * u30 + uik[7] * u31 + uik[8] * u32 + uik[9] * u33 + uik[10] * u34 + uik[11] * u35;
216:       dk[32] += uik[12] * u30 + uik[13] * u31 + uik[14] * u32 + uik[15] * u33 + uik[16] * u34 + uik[17] * u35;
217:       dk[33] += uik[18] * u30 + uik[19] * u31 + uik[20] * u32 + uik[21] * u33 + uik[22] * u34 + uik[23] * u35;
218:       dk[34] += uik[24] * u30 + uik[25] * u31 + uik[26] * u32 + uik[27] * u33 + uik[28] * u34 + uik[29] * u35;
219:       dk[35] += uik[30] * u30 + uik[31] * u31 + uik[32] * u32 + uik[33] * u33 + uik[34] * u34 + uik[35] * u35;

221:       PetscCall(PetscLogFlops(216.0 * 4.0));

223:       /* update -U(i,k) */
224:       PetscCall(PetscArraycpy(ba + ili * 36, uik, 36));

226:       /* add multiple of row i to k-th row ... */
227:       jmin = ili + 1;
228:       jmax = bi[i + 1];
229:       if (jmin < jmax) {
230:         for (j = jmin; j < jmax; j++) {
231:           /* w += -U(i,k)^T * U_bar(i,j) */
232:           wp = w + bj[j] * 36;
233:           u  = ba + j * 36;

235:           u0  = u[0];
236:           u1  = u[1];
237:           u2  = u[2];
238:           u3  = u[3];
239:           u4  = u[4];
240:           u5  = u[5];
241:           u6  = u[6];
242:           u7  = u[7];
243:           u8  = u[8];
244:           u9  = u[9];
245:           u10 = u[10];
246:           u11 = u[11];
247:           u12 = u[12];
248:           u13 = u[13];
249:           u14 = u[14];
250:           u15 = u[15];
251:           u16 = u[16];
252:           u17 = u[17];
253:           u18 = u[18];
254:           u19 = u[19];
255:           u20 = u[20];
256:           u21 = u[21];
257:           u22 = u[22];
258:           u23 = u[23];
259:           u24 = u[24];
260:           u25 = u[25];
261:           u26 = u[26];
262:           u27 = u[27];
263:           u28 = u[28];
264:           u29 = u[29];
265:           u30 = u[30];
266:           u31 = u[31];
267:           u32 = u[32];
268:           u33 = u[33];
269:           u34 = u[34];
270:           u35 = u[35];

272:           wp[0] += uik[0] * u0 + uik[1] * u1 + uik[2] * u2 + uik[3] * u3 + uik[4] * u4 + uik[5] * u5;
273:           wp[1] += uik[6] * u0 + uik[7] * u1 + uik[8] * u2 + uik[9] * u3 + uik[10] * u4 + uik[11] * u5;
274:           wp[2] += uik[12] * u0 + uik[13] * u1 + uik[14] * u2 + uik[15] * u3 + uik[16] * u4 + uik[17] * u5;
275:           wp[3] += uik[18] * u0 + uik[19] * u1 + uik[20] * u2 + uik[21] * u3 + uik[22] * u4 + uik[23] * u5;
276:           wp[4] += uik[24] * u0 + uik[25] * u1 + uik[26] * u2 + uik[27] * u3 + uik[28] * u4 + uik[29] * u5;
277:           wp[5] += uik[30] * u0 + uik[31] * u1 + uik[32] * u2 + uik[33] * u3 + uik[34] * u4 + uik[35] * u5;

279:           wp[6] += uik[0] * u6 + uik[1] * u7 + uik[2] * u8 + uik[3] * u9 + uik[4] * u10 + uik[5] * u11;
280:           wp[7] += uik[6] * u6 + uik[7] * u7 + uik[8] * u8 + uik[9] * u9 + uik[10] * u10 + uik[11] * u11;
281:           wp[8] += uik[12] * u6 + uik[13] * u7 + uik[14] * u8 + uik[15] * u9 + uik[16] * u10 + uik[17] * u11;
282:           wp[9] += uik[18] * u6 + uik[19] * u7 + uik[20] * u8 + uik[21] * u9 + uik[22] * u10 + uik[23] * u11;
283:           wp[10] += uik[24] * u6 + uik[25] * u7 + uik[26] * u8 + uik[27] * u9 + uik[28] * u10 + uik[29] * u11;
284:           wp[11] += uik[30] * u6 + uik[31] * u7 + uik[32] * u8 + uik[33] * u9 + uik[34] * u10 + uik[35] * u11;

286:           wp[12] += uik[0] * u12 + uik[1] * u13 + uik[2] * u14 + uik[3] * u15 + uik[4] * u16 + uik[5] * u17;
287:           wp[13] += uik[6] * u12 + uik[7] * u13 + uik[8] * u14 + uik[9] * u15 + uik[10] * u16 + uik[11] * u17;
288:           wp[14] += uik[12] * u12 + uik[13] * u13 + uik[14] * u14 + uik[15] * u15 + uik[16] * u16 + uik[17] * u17;
289:           wp[15] += uik[18] * u12 + uik[19] * u13 + uik[20] * u14 + uik[21] * u15 + uik[22] * u16 + uik[23] * u17;
290:           wp[16] += uik[24] * u12 + uik[25] * u13 + uik[26] * u14 + uik[27] * u15 + uik[28] * u16 + uik[29] * u17;
291:           wp[17] += uik[30] * u12 + uik[31] * u13 + uik[32] * u14 + uik[33] * u15 + uik[34] * u16 + uik[35] * u17;

293:           wp[18] += uik[0] * u18 + uik[1] * u19 + uik[2] * u20 + uik[3] * u21 + uik[4] * u22 + uik[5] * u23;
294:           wp[19] += uik[6] * u18 + uik[7] * u19 + uik[8] * u20 + uik[9] * u21 + uik[10] * u22 + uik[11] * u23;
295:           wp[20] += uik[12] * u18 + uik[13] * u19 + uik[14] * u20 + uik[15] * u21 + uik[16] * u22 + uik[17] * u23;
296:           wp[21] += uik[18] * u18 + uik[19] * u19 + uik[20] * u20 + uik[21] * u21 + uik[22] * u22 + uik[23] * u23;
297:           wp[22] += uik[24] * u18 + uik[25] * u19 + uik[26] * u20 + uik[27] * u21 + uik[28] * u22 + uik[29] * u23;
298:           wp[23] += uik[30] * u18 + uik[31] * u19 + uik[32] * u20 + uik[33] * u21 + uik[34] * u22 + uik[35] * u23;

300:           wp[24] += uik[0] * u24 + uik[1] * u25 + uik[2] * u26 + uik[3] * u27 + uik[4] * u28 + uik[5] * u29;
301:           wp[25] += uik[6] * u24 + uik[7] * u25 + uik[8] * u26 + uik[9] * u27 + uik[10] * u28 + uik[11] * u29;
302:           wp[26] += uik[12] * u24 + uik[13] * u25 + uik[14] * u26 + uik[15] * u27 + uik[16] * u28 + uik[17] * u29;
303:           wp[27] += uik[18] * u24 + uik[19] * u25 + uik[20] * u26 + uik[21] * u27 + uik[22] * u28 + uik[23] * u29;
304:           wp[28] += uik[24] * u24 + uik[25] * u25 + uik[26] * u26 + uik[27] * u27 + uik[28] * u28 + uik[29] * u29;
305:           wp[29] += uik[30] * u24 + uik[31] * u25 + uik[32] * u26 + uik[33] * u27 + uik[34] * u28 + uik[35] * u29;

307:           wp[30] += uik[0] * u30 + uik[1] * u31 + uik[2] * u32 + uik[3] * u33 + uik[4] * u34 + uik[5] * u35;
308:           wp[31] += uik[6] * u30 + uik[7] * u31 + uik[8] * u32 + uik[9] * u33 + uik[10] * u34 + uik[11] * u35;
309:           wp[32] += uik[12] * u30 + uik[13] * u31 + uik[14] * u32 + uik[15] * u33 + uik[16] * u34 + uik[17] * u35;
310:           wp[33] += uik[18] * u30 + uik[19] * u31 + uik[20] * u32 + uik[21] * u33 + uik[22] * u34 + uik[23] * u35;
311:           wp[34] += uik[24] * u30 + uik[25] * u31 + uik[26] * u32 + uik[27] * u33 + uik[28] * u34 + uik[29] * u35;
312:           wp[35] += uik[30] * u30 + uik[31] * u31 + uik[32] * u32 + uik[33] * u33 + uik[34] * u34 + uik[35] * u35;
313:         }
314:         PetscCall(PetscLogFlops(2.0 * 216.0 * (jmax - jmin)));

316:         /* ... add i to row list for next nonzero entry */
317:         il[i] = jmin; /* update il(i) in column k+1, ... mbs-1 */
318:         j     = bj[jmin];
319:         jl[i] = jl[j];
320:         jl[j] = i; /* update jl */
321:       }
322:       i = nexti;
323:     }

325:     /* save nonzero entries in k-th row of U ... */

327:     /* invert diagonal block */
328:     d = ba + k * 36;
329:     PetscCall(PetscArraycpy(d, dk, 36));
330:     PetscCall(PetscKernel_A_gets_inverse_A_6(d, shift, allowzeropivot, &zeropivotdetected));
331:     if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;

333:     jmin = bi[k];
334:     jmax = bi[k + 1];
335:     if (jmin < jmax) {
336:       for (j = jmin; j < jmax; j++) {
337:         vj = bj[j]; /* block col. index of U */
338:         u  = ba + j * 36;
339:         wp = w + vj * 36;
340:         for (k1 = 0; k1 < 36; k1++) {
341:           *u++  = *wp;
342:           *wp++ = 0.0;
343:         }
344:       }

346:       /* ... add k to row list for first nonzero entry in k-th row */
347:       il[k] = jmin;
348:       i     = bj[jmin];
349:       jl[k] = jl[i];
350:       jl[i] = k;
351:     }
352:   }

354:   PetscCall(PetscFree(w));
355:   PetscCall(PetscFree2(il, jl));
356:   PetscCall(PetscFree2(dk, uik));
357:   if (a->permute) PetscCall(PetscFree(aa));

359:   PetscCall(ISRestoreIndices(perm, &perm_ptr));

361:   C->ops->solve          = MatSolve_SeqSBAIJ_6_inplace;
362:   C->ops->solvetranspose = MatSolve_SeqSBAIJ_6_inplace;
363:   C->assembled           = PETSC_TRUE;
364:   C->preallocated        = PETSC_TRUE;

366:   PetscCall(PetscLogFlops(1.3333 * 216 * b->mbs)); /* from inverting diagonal blocks */
367:   PetscFunctionReturn(PETSC_SUCCESS);
368: }