sbaij/seq/sbaijfact9.c

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;
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;
15   PetscReal       shift = info->shiftamount;
20   allowzeropivot = PetscNot(A->erroriffailure);
30   if (!a->permute) {
31     ai = a->i;
32     aj = a->j;
33     aa = a->a;
35     ai = a->inew;
36     aj = a->jnew;
38     PetscCall(PetscArraycpy(aa, a->a, 36 * ai[mbs]));
40     PetscCall(PetscArraycpy(a2anew, a->a2anew, ai[mbs]));
56         /* transform column-oriented blocks that lie in the lower triangle to row-oriented blocks */
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 */
72     /*initialize k-th row with elements nonzero in row perm(k) of A */
84     /* modify k-th row by adding in those rows i with U(i,k) != 0 */
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;
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) */
222       /* update -U(i,k) */
225       /* add multiple of row i to k-th row ... */
230           /* w += -U(i,k)^T * U_bar(i,j) */
313         PetscCall(PetscLogFlops(2.0 * 216.0 * (jmax - jmin)));
316         il[i] = jmin; /* update il(i) in column k+1, ... mbs-1 */
324     /* save nonzero entries in k-th row of U ... */
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;
345       /* ... add k to row list for first nonzero entry in k-th row */
356   if (a->permute) PetscCall(PetscFree(aa));
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 */