1 #include <../src/mat/impls/sbaij/seq/sbaij.h>
2 #include <petsc/private/kernels/blockinvert.h>
3
4 /* Version for when blocks are 4 by 4 */
MatCholeskyFactorNumeric_SeqSBAIJ_4(Mat C,Mat A,const MatFactorInfo * info)5 PetscErrorCode MatCholeskyFactorNumeric_SeqSBAIJ_4(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, *diag, *rtmp, *rtmp_ptr;
13 PetscBool pivotinblocks = b->pivotinblocks;
14 PetscReal shift = info->shiftamount;
15 PetscBool allowzeropivot, zeropivotdetected = PETSC_FALSE;
16
17 PetscFunctionBegin;
18 /* initialization */
19 allowzeropivot = PetscNot(A->erroriffailure);
20 PetscCall(PetscCalloc1(16 * mbs, &rtmp));
21 PetscCall(PetscMalloc2(mbs, &il, mbs, &jl));
22 il[0] = 0;
23 for (i = 0; i < mbs; i++) jl[i] = mbs;
24
25 PetscCall(PetscMalloc2(16, &dk, 16, &uik));
26 PetscCall(ISGetIndices(perm, &perm_ptr));
27
28 /* check permutation */
29 if (!a->permute) {
30 ai = a->i;
31 aj = a->j;
32 aa = a->a;
33 } else {
34 ai = a->inew;
35 aj = a->jnew;
36 PetscCall(PetscMalloc1(16 * ai[mbs], &aa));
37 PetscCall(PetscArraycpy(aa, a->a, 16 * ai[mbs]));
38 PetscCall(PetscMalloc1(ai[mbs], &a2anew));
39 PetscCall(PetscArraycpy(a2anew, a->a2anew, ai[mbs]));
40
41 for (i = 0; i < mbs; i++) {
42 jmin = ai[i];
43 jmax = ai[i + 1];
44 for (j = jmin; j < jmax; j++) {
45 while (a2anew[j] != j) {
46 k = a2anew[j];
47 a2anew[j] = a2anew[k];
48 a2anew[k] = k;
49 for (k1 = 0; k1 < 16; k1++) {
50 dk[k1] = aa[k * 16 + k1];
51 aa[k * 16 + k1] = aa[j * 16 + k1];
52 aa[j * 16 + k1] = dk[k1];
53 }
54 }
55 /* transform column-oriented blocks that lie in the lower triangle to row-oriented blocks */
56 if (i > aj[j]) {
57 /* printf("change orientation, row: %d, col: %d\n",i,aj[j]); */
58 ap = aa + j * 16; /* ptr to the beginning of j-th block of aa */
59 for (k = 0; k < 16; k++) dk[k] = ap[k]; /* dk <- j-th block of aa */
60 for (k = 0; k < 4; k++) { /* j-th block of aa <- dk^T */
61 for (k1 = 0; k1 < 4; k1++) *ap++ = dk[k + 4 * k1];
62 }
63 }
64 }
65 }
66 PetscCall(PetscFree(a2anew));
67 }
68
69 /* for each row k */
70 for (k = 0; k < mbs; k++) {
71 /*initialize k-th row with elements nonzero in row perm(k) of A */
72 jmin = ai[perm_ptr[k]];
73 jmax = ai[perm_ptr[k] + 1];
74 if (jmin < jmax) {
75 ap = aa + jmin * 16;
76 for (j = jmin; j < jmax; j++) {
77 vj = perm_ptr[aj[j]]; /* block col. index */
78 rtmp_ptr = rtmp + vj * 16;
79 for (i = 0; i < 16; i++) *rtmp_ptr++ = *ap++;
80 }
81 }
82
83 /* modify k-th row by adding in those rows i with U(i,k) != 0 */
84 PetscCall(PetscArraycpy(dk, rtmp + k * 16, 16));
85 i = jl[k]; /* first row to be added to k_th row */
86
87 while (i < mbs) {
88 nexti = jl[i]; /* next row to be added to k_th row */
89
90 /* compute multiplier */
91 ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */
92
93 /* uik = -inv(Di)*U_bar(i,k) */
94 diag = ba + i * 16;
95 u = ba + ili * 16;
96
97 uik[0] = -(diag[0] * u[0] + diag[4] * u[1] + diag[8] * u[2] + diag[12] * u[3]);
98 uik[1] = -(diag[1] * u[0] + diag[5] * u[1] + diag[9] * u[2] + diag[13] * u[3]);
99 uik[2] = -(diag[2] * u[0] + diag[6] * u[1] + diag[10] * u[2] + diag[14] * u[3]);
100 uik[3] = -(diag[3] * u[0] + diag[7] * u[1] + diag[11] * u[2] + diag[15] * u[3]);
101
102 uik[4] = -(diag[0] * u[4] + diag[4] * u[5] + diag[8] * u[6] + diag[12] * u[7]);
103 uik[5] = -(diag[1] * u[4] + diag[5] * u[5] + diag[9] * u[6] + diag[13] * u[7]);
104 uik[6] = -(diag[2] * u[4] + diag[6] * u[5] + diag[10] * u[6] + diag[14] * u[7]);
105 uik[7] = -(diag[3] * u[4] + diag[7] * u[5] + diag[11] * u[6] + diag[15] * u[7]);
106
107 uik[8] = -(diag[0] * u[8] + diag[4] * u[9] + diag[8] * u[10] + diag[12] * u[11]);
108 uik[9] = -(diag[1] * u[8] + diag[5] * u[9] + diag[9] * u[10] + diag[13] * u[11]);
109 uik[10] = -(diag[2] * u[8] + diag[6] * u[9] + diag[10] * u[10] + diag[14] * u[11]);
110 uik[11] = -(diag[3] * u[8] + diag[7] * u[9] + diag[11] * u[10] + diag[15] * u[11]);
111
112 uik[12] = -(diag[0] * u[12] + diag[4] * u[13] + diag[8] * u[14] + diag[12] * u[15]);
113 uik[13] = -(diag[1] * u[12] + diag[5] * u[13] + diag[9] * u[14] + diag[13] * u[15]);
114 uik[14] = -(diag[2] * u[12] + diag[6] * u[13] + diag[10] * u[14] + diag[14] * u[15]);
115 uik[15] = -(diag[3] * u[12] + diag[7] * u[13] + diag[11] * u[14] + diag[15] * u[15]);
116
117 /* update D(k) += -U(i,k)^T * U_bar(i,k) */
118 dk[0] += uik[0] * u[0] + uik[1] * u[1] + uik[2] * u[2] + uik[3] * u[3];
119 dk[1] += uik[4] * u[0] + uik[5] * u[1] + uik[6] * u[2] + uik[7] * u[3];
120 dk[2] += uik[8] * u[0] + uik[9] * u[1] + uik[10] * u[2] + uik[11] * u[3];
121 dk[3] += uik[12] * u[0] + uik[13] * u[1] + uik[14] * u[2] + uik[15] * u[3];
122
123 dk[4] += uik[0] * u[4] + uik[1] * u[5] + uik[2] * u[6] + uik[3] * u[7];
124 dk[5] += uik[4] * u[4] + uik[5] * u[5] + uik[6] * u[6] + uik[7] * u[7];
125 dk[6] += uik[8] * u[4] + uik[9] * u[5] + uik[10] * u[6] + uik[11] * u[7];
126 dk[7] += uik[12] * u[4] + uik[13] * u[5] + uik[14] * u[6] + uik[15] * u[7];
127
128 dk[8] += uik[0] * u[8] + uik[1] * u[9] + uik[2] * u[10] + uik[3] * u[11];
129 dk[9] += uik[4] * u[8] + uik[5] * u[9] + uik[6] * u[10] + uik[7] * u[11];
130 dk[10] += uik[8] * u[8] + uik[9] * u[9] + uik[10] * u[10] + uik[11] * u[11];
131 dk[11] += uik[12] * u[8] + uik[13] * u[9] + uik[14] * u[10] + uik[15] * u[11];
132
133 dk[12] += uik[0] * u[12] + uik[1] * u[13] + uik[2] * u[14] + uik[3] * u[15];
134 dk[13] += uik[4] * u[12] + uik[5] * u[13] + uik[6] * u[14] + uik[7] * u[15];
135 dk[14] += uik[8] * u[12] + uik[9] * u[13] + uik[10] * u[14] + uik[11] * u[15];
136 dk[15] += uik[12] * u[12] + uik[13] * u[13] + uik[14] * u[14] + uik[15] * u[15];
137
138 PetscCall(PetscLogFlops(64.0 * 4.0));
139
140 /* update -U(i,k) */
141 PetscCall(PetscArraycpy(ba + ili * 16, uik, 16));
142
143 /* add multiple of row i to k-th row ... */
144 jmin = ili + 1;
145 jmax = bi[i + 1];
146 if (jmin < jmax) {
147 for (j = jmin; j < jmax; j++) {
148 /* rtmp += -U(i,k)^T * U_bar(i,j) */
149 rtmp_ptr = rtmp + bj[j] * 16;
150 u = ba + j * 16;
151 rtmp_ptr[0] += uik[0] * u[0] + uik[1] * u[1] + uik[2] * u[2] + uik[3] * u[3];
152 rtmp_ptr[1] += uik[4] * u[0] + uik[5] * u[1] + uik[6] * u[2] + uik[7] * u[3];
153 rtmp_ptr[2] += uik[8] * u[0] + uik[9] * u[1] + uik[10] * u[2] + uik[11] * u[3];
154 rtmp_ptr[3] += uik[12] * u[0] + uik[13] * u[1] + uik[14] * u[2] + uik[15] * u[3];
155
156 rtmp_ptr[4] += uik[0] * u[4] + uik[1] * u[5] + uik[2] * u[6] + uik[3] * u[7];
157 rtmp_ptr[5] += uik[4] * u[4] + uik[5] * u[5] + uik[6] * u[6] + uik[7] * u[7];
158 rtmp_ptr[6] += uik[8] * u[4] + uik[9] * u[5] + uik[10] * u[6] + uik[11] * u[7];
159 rtmp_ptr[7] += uik[12] * u[4] + uik[13] * u[5] + uik[14] * u[6] + uik[15] * u[7];
160
161 rtmp_ptr[8] += uik[0] * u[8] + uik[1] * u[9] + uik[2] * u[10] + uik[3] * u[11];
162 rtmp_ptr[9] += uik[4] * u[8] + uik[5] * u[9] + uik[6] * u[10] + uik[7] * u[11];
163 rtmp_ptr[10] += uik[8] * u[8] + uik[9] * u[9] + uik[10] * u[10] + uik[11] * u[11];
164 rtmp_ptr[11] += uik[12] * u[8] + uik[13] * u[9] + uik[14] * u[10] + uik[15] * u[11];
165
166 rtmp_ptr[12] += uik[0] * u[12] + uik[1] * u[13] + uik[2] * u[14] + uik[3] * u[15];
167 rtmp_ptr[13] += uik[4] * u[12] + uik[5] * u[13] + uik[6] * u[14] + uik[7] * u[15];
168 rtmp_ptr[14] += uik[8] * u[12] + uik[9] * u[13] + uik[10] * u[14] + uik[11] * u[15];
169 rtmp_ptr[15] += uik[12] * u[12] + uik[13] * u[13] + uik[14] * u[14] + uik[15] * u[15];
170 }
171 PetscCall(PetscLogFlops(2.0 * 64.0 * (jmax - jmin)));
172
173 /* ... add i to row list for next nonzero entry */
174 il[i] = jmin; /* update il(i) in column k+1, ... mbs-1 */
175 j = bj[jmin];
176 jl[i] = jl[j];
177 jl[j] = i; /* update jl */
178 }
179 i = nexti;
180 }
181
182 /* save nonzero entries in k-th row of U ... */
183
184 /* invert diagonal block */
185 diag = ba + k * 16;
186 PetscCall(PetscArraycpy(diag, dk, 16));
187
188 if (pivotinblocks) {
189 PetscCall(PetscKernel_A_gets_inverse_A_4(diag, shift, allowzeropivot, &zeropivotdetected));
190 if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
191 } else {
192 PetscCall(PetscKernel_A_gets_inverse_A_4_nopivot(diag));
193 }
194
195 jmin = bi[k];
196 jmax = bi[k + 1];
197 if (jmin < jmax) {
198 for (j = jmin; j < jmax; j++) {
199 vj = bj[j]; /* block col. index of U */
200 u = ba + j * 16;
201 rtmp_ptr = rtmp + vj * 16;
202 for (k1 = 0; k1 < 16; k1++) {
203 *u++ = *rtmp_ptr;
204 *rtmp_ptr++ = 0.0;
205 }
206 }
207
208 /* ... add k to row list for first nonzero entry in k-th row */
209 il[k] = jmin;
210 i = bj[jmin];
211 jl[k] = jl[i];
212 jl[i] = k;
213 }
214 }
215
216 PetscCall(PetscFree(rtmp));
217 PetscCall(PetscFree2(il, jl));
218 PetscCall(PetscFree2(dk, uik));
219 if (a->permute) PetscCall(PetscFree(aa));
220
221 PetscCall(ISRestoreIndices(perm, &perm_ptr));
222
223 C->ops->solve = MatSolve_SeqSBAIJ_4_inplace;
224 C->ops->solvetranspose = MatSolve_SeqSBAIJ_4_inplace;
225 C->assembled = PETSC_TRUE;
226 C->preallocated = PETSC_TRUE;
227
228 PetscCall(PetscLogFlops(1.3333 * 64 * b->mbs)); /* from inverting diagonal blocks */
229 PetscFunctionReturn(PETSC_SUCCESS);
230 }
231