sbaij/seq/sbaijfact8.c

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;
12   MatScalar    *ba = b->a, *aa, *ap, *dk, *uik;
13   MatScalar    *u, *d, *rtmp, *rtmp_ptr, work[25];
14   PetscReal     shift = info->shiftamount;
19   allowzeropivot = PetscNot(A->erroriffailure);
26   ai = a->i;
27   aj = a->j;
28   aa = a->a;
32     /*initialize k-th row with elements nonzero in row k of A */
44     /* modify k-th row by adding in those rows i with U(i,k) != 0 */
52       ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
54       /* uik = -inv(Di)*U_bar(i,k) */
55       d = ba + i * 25;
58       uik[0] = -(d[0] * u[0] + d[5] * u[1] + d[10] * u[2] + d[15] * u[3] + d[20] * u[4]);
59       uik[1] = -(d[1] * u[0] + d[6] * u[1] + d[11] * u[2] + d[16] * u[3] + d[21] * u[4]);
60       uik[2] = -(d[2] * u[0] + d[7] * u[1] + d[12] * u[2] + d[17] * u[3] + d[22] * u[4]);
61       uik[3] = -(d[3] * u[0] + d[8] * u[1] + d[13] * u[2] + d[18] * u[3] + d[23] * u[4]);
62       uik[4] = -(d[4] * u[0] + d[9] * u[1] + d[14] * u[2] + d[19] * u[3] + d[24] * u[4]);
64       uik[5] = -(d[0] * u[5] + d[5] * u[6] + d[10] * u[7] + d[15] * u[8] + d[20] * u[9]);
65       uik[6] = -(d[1] * u[5] + d[6] * u[6] + d[11] * u[7] + d[16] * u[8] + d[21] * u[9]);
66       uik[7] = -(d[2] * u[5] + d[7] * u[6] + d[12] * u[7] + d[17] * u[8] + d[22] * u[9]);
67       uik[8] = -(d[3] * u[5] + d[8] * u[6] + d[13] * u[7] + d[18] * u[8] + d[23] * u[9]);
68       uik[9] = -(d[4] * u[5] + d[9] * u[6] + d[14] * u[7] + d[19] * u[8] + d[24] * u[9]);
70       uik[10] = -(d[0] * u[10] + d[5] * u[11] + d[10] * u[12] + d[15] * u[13] + d[20] * u[14]);
71       uik[11] = -(d[1] * u[10] + d[6] * u[11] + d[11] * u[12] + d[16] * u[13] + d[21] * u[14]);
72       uik[12] = -(d[2] * u[10] + d[7] * u[11] + d[12] * u[12] + d[17] * u[13] + d[22] * u[14]);
73       uik[13] = -(d[3] * u[10] + d[8] * u[11] + d[13] * u[12] + d[18] * u[13] + d[23] * u[14]);
74       uik[14] = -(d[4] * u[10] + d[9] * u[11] + d[14] * u[12] + d[19] * u[13] + d[24] * u[14]);
76       uik[15] = -(d[0] * u[15] + d[5] * u[16] + d[10] * u[17] + d[15] * u[18] + d[20] * u[19]);
77       uik[16] = -(d[1] * u[15] + d[6] * u[16] + d[11] * u[17] + d[16] * u[18] + d[21] * u[19]);
78       uik[17] = -(d[2] * u[15] + d[7] * u[16] + d[12] * u[17] + d[17] * u[18] + d[22] * u[19]);
79       uik[18] = -(d[3] * u[15] + d[8] * u[16] + d[13] * u[17] + d[18] * u[18] + d[23] * u[19]);
80       uik[19] = -(d[4] * u[15] + d[9] * u[16] + d[14] * u[17] + d[19] * u[18] + d[24] * u[19]);
82       uik[20] = -(d[0] * u[20] + d[5] * u[21] + d[10] * u[22] + d[15] * u[23] + d[20] * u[24]);
83       uik[21] = -(d[1] * u[20] + d[6] * u[21] + d[11] * u[22] + d[16] * u[23] + d[21] * u[24]);
84       uik[22] = -(d[2] * u[20] + d[7] * u[21] + d[12] * u[22] + d[17] * u[23] + d[22] * u[24]);
85       uik[23] = -(d[3] * u[20] + d[8] * u[21] + d[13] * u[22] + d[18] * u[23] + d[23] * u[24]);
86       uik[24] = -(d[4] * u[20] + d[9] * u[21] + d[14] * u[22] + d[19] * u[23] + d[24] * u[24]);
88       /* update D(k) += -U(i,k)^T * U_bar(i,k) */
121       /* update -U(i,k) */
124       /* add multiple of row i to k-th row ... */
129           /* rtmp += -U(i,k)^T * U_bar(i,j) */
162         PetscCall(PetscLogFlops(2.0 * 125.0 * (jmax - jmin)));
165         il[i] = jmin; /* update il(i) in column k+1, ... mbs-1 */
173     /* save nonzero entries in k-th row of U ... */
176     d = ba + k * 25;
177     PetscCall(PetscArraycpy(d, dk, 25));
178     PetscCall(PetscKernel_A_gets_inverse_A_5(d, ipvt, work, shift, allowzeropivot, &zeropivotdetected));
179     if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
194       /* ... add k to row list for first nonzero entry in k-th row */
206   C->ops->solve          = MatSolve_SeqSBAIJ_5_NaturalOrdering_inplace;
207   C->ops->solvetranspose = MatSolve_SeqSBAIJ_5_NaturalOrdering_inplace;
208   C->ops->forwardsolve   = MatForwardSolve_SeqSBAIJ_5_NaturalOrdering_inplace;
209   C->ops->backwardsolve  = MatBackwardSolve_SeqSBAIJ_5_NaturalOrdering_inplace;
210   C->assembled           = PETSC_TRUE;
211   C->preallocated        = PETSC_TRUE;
213   PetscCall(PetscLogFlops(1.3333 * 125 * b->mbs)); /* from inverting diagonal blocks */