1 2 /* 3 Factorization code for BAIJ format. 4 */ 5 #include <../src/mat/impls/baij/seq/baij.h> 6 #include <petsc/private/kernels/blockinvert.h> 7 8 /* ------------------------------------------------------------*/ 9 /* 10 Version for when blocks are 6 by 6 11 */ 12 PetscErrorCode MatLUFactorNumeric_SeqBAIJ_6_inplace(Mat C, Mat A, const MatFactorInfo *info) { 13 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)C->data; 14 IS isrow = b->row, isicol = b->icol; 15 const PetscInt *ajtmpold, *ajtmp, *diag_offset = b->diag, *r, *ic, *bi = b->i, *bj = b->j, *ai = a->i, *aj = a->j, *pj; 16 PetscInt nz, row, i, j, n = a->mbs, idx; 17 MatScalar *pv, *v, *rtmp, *pc, *w, *x; 18 MatScalar p1, p2, p3, p4, m1, m2, m3, m4, m5, m6, m7, m8, m9, x1, x2, x3, x4; 19 MatScalar p5, p6, p7, p8, p9, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15, x16; 20 MatScalar x17, x18, x19, x20, x21, x22, x23, x24, x25, p10, p11, p12, p13, p14; 21 MatScalar p15, p16, p17, p18, p19, p20, p21, p22, p23, p24, p25, m10, m11, m12; 22 MatScalar m13, m14, m15, m16, m17, m18, m19, m20, m21, m22, m23, m24, m25; 23 MatScalar p26, p27, p28, p29, p30, p31, p32, p33, p34, p35, p36; 24 MatScalar x26, x27, x28, x29, x30, x31, x32, x33, x34, x35, x36; 25 MatScalar m26, m27, m28, m29, m30, m31, m32, m33, m34, m35, m36; 26 MatScalar *ba = b->a, *aa = a->a; 27 PetscReal shift = info->shiftamount; 28 PetscBool allowzeropivot, zeropivotdetected; 29 30 PetscFunctionBegin; 31 allowzeropivot = PetscNot(A->erroriffailure); 32 PetscCall(ISGetIndices(isrow, &r)); 33 PetscCall(ISGetIndices(isicol, &ic)); 34 PetscCall(PetscMalloc1(36 * (n + 1), &rtmp)); 35 36 for (i = 0; i < n; i++) { 37 nz = bi[i + 1] - bi[i]; 38 ajtmp = bj + bi[i]; 39 for (j = 0; j < nz; j++) { 40 x = rtmp + 36 * ajtmp[j]; 41 x[0] = x[1] = x[2] = x[3] = x[4] = x[5] = x[6] = x[7] = x[8] = x[9] = 0.0; 42 x[10] = x[11] = x[12] = x[13] = x[14] = x[15] = x[16] = x[17] = 0.0; 43 x[18] = x[19] = x[20] = x[21] = x[22] = x[23] = x[24] = x[25] = 0.0; 44 x[26] = x[27] = x[28] = x[29] = x[30] = x[31] = x[32] = x[33] = 0.0; 45 x[34] = x[35] = 0.0; 46 } 47 /* load in initial (unfactored row) */ 48 idx = r[i]; 49 nz = ai[idx + 1] - ai[idx]; 50 ajtmpold = aj + ai[idx]; 51 v = aa + 36 * ai[idx]; 52 for (j = 0; j < nz; j++) { 53 x = rtmp + 36 * ic[ajtmpold[j]]; 54 x[0] = v[0]; 55 x[1] = v[1]; 56 x[2] = v[2]; 57 x[3] = v[3]; 58 x[4] = v[4]; 59 x[5] = v[5]; 60 x[6] = v[6]; 61 x[7] = v[7]; 62 x[8] = v[8]; 63 x[9] = v[9]; 64 x[10] = v[10]; 65 x[11] = v[11]; 66 x[12] = v[12]; 67 x[13] = v[13]; 68 x[14] = v[14]; 69 x[15] = v[15]; 70 x[16] = v[16]; 71 x[17] = v[17]; 72 x[18] = v[18]; 73 x[19] = v[19]; 74 x[20] = v[20]; 75 x[21] = v[21]; 76 x[22] = v[22]; 77 x[23] = v[23]; 78 x[24] = v[24]; 79 x[25] = v[25]; 80 x[26] = v[26]; 81 x[27] = v[27]; 82 x[28] = v[28]; 83 x[29] = v[29]; 84 x[30] = v[30]; 85 x[31] = v[31]; 86 x[32] = v[32]; 87 x[33] = v[33]; 88 x[34] = v[34]; 89 x[35] = v[35]; 90 v += 36; 91 } 92 row = *ajtmp++; 93 while (row < i) { 94 pc = rtmp + 36 * row; 95 p1 = pc[0]; 96 p2 = pc[1]; 97 p3 = pc[2]; 98 p4 = pc[3]; 99 p5 = pc[4]; 100 p6 = pc[5]; 101 p7 = pc[6]; 102 p8 = pc[7]; 103 p9 = pc[8]; 104 p10 = pc[9]; 105 p11 = pc[10]; 106 p12 = pc[11]; 107 p13 = pc[12]; 108 p14 = pc[13]; 109 p15 = pc[14]; 110 p16 = pc[15]; 111 p17 = pc[16]; 112 p18 = pc[17]; 113 p19 = pc[18]; 114 p20 = pc[19]; 115 p21 = pc[20]; 116 p22 = pc[21]; 117 p23 = pc[22]; 118 p24 = pc[23]; 119 p25 = pc[24]; 120 p26 = pc[25]; 121 p27 = pc[26]; 122 p28 = pc[27]; 123 p29 = pc[28]; 124 p30 = pc[29]; 125 p31 = pc[30]; 126 p32 = pc[31]; 127 p33 = pc[32]; 128 p34 = pc[33]; 129 p35 = pc[34]; 130 p36 = pc[35]; 131 if (p1 != 0.0 || p2 != 0.0 || p3 != 0.0 || p4 != 0.0 || p5 != 0.0 || p6 != 0.0 || p7 != 0.0 || p8 != 0.0 || p9 != 0.0 || p10 != 0.0 || p11 != 0.0 || p12 != 0.0 || p13 != 0.0 || p14 != 0.0 || p15 != 0.0 || p16 != 0.0 || p17 != 0.0 || p18 != 0.0 || p19 != 0.0 || p20 != 0.0 || p21 != 0.0 || p22 != 0.0 || p23 != 0.0 || p24 != 0.0 || p25 != 0.0 || p26 != 0.0 || p27 != 0.0 || p28 != 0.0 || p29 != 0.0 || p30 != 0.0 || p31 != 0.0 || p32 != 0.0 || p33 != 0.0 || p34 != 0.0 || p35 != 0.0 || p36 != 0.0) { 132 pv = ba + 36 * diag_offset[row]; 133 pj = bj + diag_offset[row] + 1; 134 x1 = pv[0]; 135 x2 = pv[1]; 136 x3 = pv[2]; 137 x4 = pv[3]; 138 x5 = pv[4]; 139 x6 = pv[5]; 140 x7 = pv[6]; 141 x8 = pv[7]; 142 x9 = pv[8]; 143 x10 = pv[9]; 144 x11 = pv[10]; 145 x12 = pv[11]; 146 x13 = pv[12]; 147 x14 = pv[13]; 148 x15 = pv[14]; 149 x16 = pv[15]; 150 x17 = pv[16]; 151 x18 = pv[17]; 152 x19 = pv[18]; 153 x20 = pv[19]; 154 x21 = pv[20]; 155 x22 = pv[21]; 156 x23 = pv[22]; 157 x24 = pv[23]; 158 x25 = pv[24]; 159 x26 = pv[25]; 160 x27 = pv[26]; 161 x28 = pv[27]; 162 x29 = pv[28]; 163 x30 = pv[29]; 164 x31 = pv[30]; 165 x32 = pv[31]; 166 x33 = pv[32]; 167 x34 = pv[33]; 168 x35 = pv[34]; 169 x36 = pv[35]; 170 pc[0] = m1 = p1 * x1 + p7 * x2 + p13 * x3 + p19 * x4 + p25 * x5 + p31 * x6; 171 pc[1] = m2 = p2 * x1 + p8 * x2 + p14 * x3 + p20 * x4 + p26 * x5 + p32 * x6; 172 pc[2] = m3 = p3 * x1 + p9 * x2 + p15 * x3 + p21 * x4 + p27 * x5 + p33 * x6; 173 pc[3] = m4 = p4 * x1 + p10 * x2 + p16 * x3 + p22 * x4 + p28 * x5 + p34 * x6; 174 pc[4] = m5 = p5 * x1 + p11 * x2 + p17 * x3 + p23 * x4 + p29 * x5 + p35 * x6; 175 pc[5] = m6 = p6 * x1 + p12 * x2 + p18 * x3 + p24 * x4 + p30 * x5 + p36 * x6; 176 177 pc[6] = m7 = p1 * x7 + p7 * x8 + p13 * x9 + p19 * x10 + p25 * x11 + p31 * x12; 178 pc[7] = m8 = p2 * x7 + p8 * x8 + p14 * x9 + p20 * x10 + p26 * x11 + p32 * x12; 179 pc[8] = m9 = p3 * x7 + p9 * x8 + p15 * x9 + p21 * x10 + p27 * x11 + p33 * x12; 180 pc[9] = m10 = p4 * x7 + p10 * x8 + p16 * x9 + p22 * x10 + p28 * x11 + p34 * x12; 181 pc[10] = m11 = p5 * x7 + p11 * x8 + p17 * x9 + p23 * x10 + p29 * x11 + p35 * x12; 182 pc[11] = m12 = p6 * x7 + p12 * x8 + p18 * x9 + p24 * x10 + p30 * x11 + p36 * x12; 183 184 pc[12] = m13 = p1 * x13 + p7 * x14 + p13 * x15 + p19 * x16 + p25 * x17 + p31 * x18; 185 pc[13] = m14 = p2 * x13 + p8 * x14 + p14 * x15 + p20 * x16 + p26 * x17 + p32 * x18; 186 pc[14] = m15 = p3 * x13 + p9 * x14 + p15 * x15 + p21 * x16 + p27 * x17 + p33 * x18; 187 pc[15] = m16 = p4 * x13 + p10 * x14 + p16 * x15 + p22 * x16 + p28 * x17 + p34 * x18; 188 pc[16] = m17 = p5 * x13 + p11 * x14 + p17 * x15 + p23 * x16 + p29 * x17 + p35 * x18; 189 pc[17] = m18 = p6 * x13 + p12 * x14 + p18 * x15 + p24 * x16 + p30 * x17 + p36 * x18; 190 191 pc[18] = m19 = p1 * x19 + p7 * x20 + p13 * x21 + p19 * x22 + p25 * x23 + p31 * x24; 192 pc[19] = m20 = p2 * x19 + p8 * x20 + p14 * x21 + p20 * x22 + p26 * x23 + p32 * x24; 193 pc[20] = m21 = p3 * x19 + p9 * x20 + p15 * x21 + p21 * x22 + p27 * x23 + p33 * x24; 194 pc[21] = m22 = p4 * x19 + p10 * x20 + p16 * x21 + p22 * x22 + p28 * x23 + p34 * x24; 195 pc[22] = m23 = p5 * x19 + p11 * x20 + p17 * x21 + p23 * x22 + p29 * x23 + p35 * x24; 196 pc[23] = m24 = p6 * x19 + p12 * x20 + p18 * x21 + p24 * x22 + p30 * x23 + p36 * x24; 197 198 pc[24] = m25 = p1 * x25 + p7 * x26 + p13 * x27 + p19 * x28 + p25 * x29 + p31 * x30; 199 pc[25] = m26 = p2 * x25 + p8 * x26 + p14 * x27 + p20 * x28 + p26 * x29 + p32 * x30; 200 pc[26] = m27 = p3 * x25 + p9 * x26 + p15 * x27 + p21 * x28 + p27 * x29 + p33 * x30; 201 pc[27] = m28 = p4 * x25 + p10 * x26 + p16 * x27 + p22 * x28 + p28 * x29 + p34 * x30; 202 pc[28] = m29 = p5 * x25 + p11 * x26 + p17 * x27 + p23 * x28 + p29 * x29 + p35 * x30; 203 pc[29] = m30 = p6 * x25 + p12 * x26 + p18 * x27 + p24 * x28 + p30 * x29 + p36 * x30; 204 205 pc[30] = m31 = p1 * x31 + p7 * x32 + p13 * x33 + p19 * x34 + p25 * x35 + p31 * x36; 206 pc[31] = m32 = p2 * x31 + p8 * x32 + p14 * x33 + p20 * x34 + p26 * x35 + p32 * x36; 207 pc[32] = m33 = p3 * x31 + p9 * x32 + p15 * x33 + p21 * x34 + p27 * x35 + p33 * x36; 208 pc[33] = m34 = p4 * x31 + p10 * x32 + p16 * x33 + p22 * x34 + p28 * x35 + p34 * x36; 209 pc[34] = m35 = p5 * x31 + p11 * x32 + p17 * x33 + p23 * x34 + p29 * x35 + p35 * x36; 210 pc[35] = m36 = p6 * x31 + p12 * x32 + p18 * x33 + p24 * x34 + p30 * x35 + p36 * x36; 211 212 nz = bi[row + 1] - diag_offset[row] - 1; 213 pv += 36; 214 for (j = 0; j < nz; j++) { 215 x1 = pv[0]; 216 x2 = pv[1]; 217 x3 = pv[2]; 218 x4 = pv[3]; 219 x5 = pv[4]; 220 x6 = pv[5]; 221 x7 = pv[6]; 222 x8 = pv[7]; 223 x9 = pv[8]; 224 x10 = pv[9]; 225 x11 = pv[10]; 226 x12 = pv[11]; 227 x13 = pv[12]; 228 x14 = pv[13]; 229 x15 = pv[14]; 230 x16 = pv[15]; 231 x17 = pv[16]; 232 x18 = pv[17]; 233 x19 = pv[18]; 234 x20 = pv[19]; 235 x21 = pv[20]; 236 x22 = pv[21]; 237 x23 = pv[22]; 238 x24 = pv[23]; 239 x25 = pv[24]; 240 x26 = pv[25]; 241 x27 = pv[26]; 242 x28 = pv[27]; 243 x29 = pv[28]; 244 x30 = pv[29]; 245 x31 = pv[30]; 246 x32 = pv[31]; 247 x33 = pv[32]; 248 x34 = pv[33]; 249 x35 = pv[34]; 250 x36 = pv[35]; 251 x = rtmp + 36 * pj[j]; 252 x[0] -= m1 * x1 + m7 * x2 + m13 * x3 + m19 * x4 + m25 * x5 + m31 * x6; 253 x[1] -= m2 * x1 + m8 * x2 + m14 * x3 + m20 * x4 + m26 * x5 + m32 * x6; 254 x[2] -= m3 * x1 + m9 * x2 + m15 * x3 + m21 * x4 + m27 * x5 + m33 * x6; 255 x[3] -= m4 * x1 + m10 * x2 + m16 * x3 + m22 * x4 + m28 * x5 + m34 * x6; 256 x[4] -= m5 * x1 + m11 * x2 + m17 * x3 + m23 * x4 + m29 * x5 + m35 * x6; 257 x[5] -= m6 * x1 + m12 * x2 + m18 * x3 + m24 * x4 + m30 * x5 + m36 * x6; 258 259 x[6] -= m1 * x7 + m7 * x8 + m13 * x9 + m19 * x10 + m25 * x11 + m31 * x12; 260 x[7] -= m2 * x7 + m8 * x8 + m14 * x9 + m20 * x10 + m26 * x11 + m32 * x12; 261 x[8] -= m3 * x7 + m9 * x8 + m15 * x9 + m21 * x10 + m27 * x11 + m33 * x12; 262 x[9] -= m4 * x7 + m10 * x8 + m16 * x9 + m22 * x10 + m28 * x11 + m34 * x12; 263 x[10] -= m5 * x7 + m11 * x8 + m17 * x9 + m23 * x10 + m29 * x11 + m35 * x12; 264 x[11] -= m6 * x7 + m12 * x8 + m18 * x9 + m24 * x10 + m30 * x11 + m36 * x12; 265 266 x[12] -= m1 * x13 + m7 * x14 + m13 * x15 + m19 * x16 + m25 * x17 + m31 * x18; 267 x[13] -= m2 * x13 + m8 * x14 + m14 * x15 + m20 * x16 + m26 * x17 + m32 * x18; 268 x[14] -= m3 * x13 + m9 * x14 + m15 * x15 + m21 * x16 + m27 * x17 + m33 * x18; 269 x[15] -= m4 * x13 + m10 * x14 + m16 * x15 + m22 * x16 + m28 * x17 + m34 * x18; 270 x[16] -= m5 * x13 + m11 * x14 + m17 * x15 + m23 * x16 + m29 * x17 + m35 * x18; 271 x[17] -= m6 * x13 + m12 * x14 + m18 * x15 + m24 * x16 + m30 * x17 + m36 * x18; 272 273 x[18] -= m1 * x19 + m7 * x20 + m13 * x21 + m19 * x22 + m25 * x23 + m31 * x24; 274 x[19] -= m2 * x19 + m8 * x20 + m14 * x21 + m20 * x22 + m26 * x23 + m32 * x24; 275 x[20] -= m3 * x19 + m9 * x20 + m15 * x21 + m21 * x22 + m27 * x23 + m33 * x24; 276 x[21] -= m4 * x19 + m10 * x20 + m16 * x21 + m22 * x22 + m28 * x23 + m34 * x24; 277 x[22] -= m5 * x19 + m11 * x20 + m17 * x21 + m23 * x22 + m29 * x23 + m35 * x24; 278 x[23] -= m6 * x19 + m12 * x20 + m18 * x21 + m24 * x22 + m30 * x23 + m36 * x24; 279 280 x[24] -= m1 * x25 + m7 * x26 + m13 * x27 + m19 * x28 + m25 * x29 + m31 * x30; 281 x[25] -= m2 * x25 + m8 * x26 + m14 * x27 + m20 * x28 + m26 * x29 + m32 * x30; 282 x[26] -= m3 * x25 + m9 * x26 + m15 * x27 + m21 * x28 + m27 * x29 + m33 * x30; 283 x[27] -= m4 * x25 + m10 * x26 + m16 * x27 + m22 * x28 + m28 * x29 + m34 * x30; 284 x[28] -= m5 * x25 + m11 * x26 + m17 * x27 + m23 * x28 + m29 * x29 + m35 * x30; 285 x[29] -= m6 * x25 + m12 * x26 + m18 * x27 + m24 * x28 + m30 * x29 + m36 * x30; 286 287 x[30] -= m1 * x31 + m7 * x32 + m13 * x33 + m19 * x34 + m25 * x35 + m31 * x36; 288 x[31] -= m2 * x31 + m8 * x32 + m14 * x33 + m20 * x34 + m26 * x35 + m32 * x36; 289 x[32] -= m3 * x31 + m9 * x32 + m15 * x33 + m21 * x34 + m27 * x35 + m33 * x36; 290 x[33] -= m4 * x31 + m10 * x32 + m16 * x33 + m22 * x34 + m28 * x35 + m34 * x36; 291 x[34] -= m5 * x31 + m11 * x32 + m17 * x33 + m23 * x34 + m29 * x35 + m35 * x36; 292 x[35] -= m6 * x31 + m12 * x32 + m18 * x33 + m24 * x34 + m30 * x35 + m36 * x36; 293 294 pv += 36; 295 } 296 PetscCall(PetscLogFlops(432.0 * nz + 396.0)); 297 } 298 row = *ajtmp++; 299 } 300 /* finished row so stick it into b->a */ 301 pv = ba + 36 * bi[i]; 302 pj = bj + bi[i]; 303 nz = bi[i + 1] - bi[i]; 304 for (j = 0; j < nz; j++) { 305 x = rtmp + 36 * pj[j]; 306 pv[0] = x[0]; 307 pv[1] = x[1]; 308 pv[2] = x[2]; 309 pv[3] = x[3]; 310 pv[4] = x[4]; 311 pv[5] = x[5]; 312 pv[6] = x[6]; 313 pv[7] = x[7]; 314 pv[8] = x[8]; 315 pv[9] = x[9]; 316 pv[10] = x[10]; 317 pv[11] = x[11]; 318 pv[12] = x[12]; 319 pv[13] = x[13]; 320 pv[14] = x[14]; 321 pv[15] = x[15]; 322 pv[16] = x[16]; 323 pv[17] = x[17]; 324 pv[18] = x[18]; 325 pv[19] = x[19]; 326 pv[20] = x[20]; 327 pv[21] = x[21]; 328 pv[22] = x[22]; 329 pv[23] = x[23]; 330 pv[24] = x[24]; 331 pv[25] = x[25]; 332 pv[26] = x[26]; 333 pv[27] = x[27]; 334 pv[28] = x[28]; 335 pv[29] = x[29]; 336 pv[30] = x[30]; 337 pv[31] = x[31]; 338 pv[32] = x[32]; 339 pv[33] = x[33]; 340 pv[34] = x[34]; 341 pv[35] = x[35]; 342 pv += 36; 343 } 344 /* invert diagonal block */ 345 w = ba + 36 * diag_offset[i]; 346 PetscCall(PetscKernel_A_gets_inverse_A_6(w, shift, allowzeropivot, &zeropivotdetected)); 347 if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 348 } 349 350 PetscCall(PetscFree(rtmp)); 351 PetscCall(ISRestoreIndices(isicol, &ic)); 352 PetscCall(ISRestoreIndices(isrow, &r)); 353 354 C->ops->solve = MatSolve_SeqBAIJ_6_inplace; 355 C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_6_inplace; 356 C->assembled = PETSC_TRUE; 357 358 PetscCall(PetscLogFlops(1.333333333333 * 6 * 6 * 6 * b->mbs)); /* from inverting diagonal blocks */ 359 PetscFunctionReturn(0); 360 } 361 362 PetscErrorCode MatLUFactorNumeric_SeqBAIJ_6(Mat B, Mat A, const MatFactorInfo *info) { 363 Mat C = B; 364 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)C->data; 365 IS isrow = b->row, isicol = b->icol; 366 const PetscInt *r, *ic; 367 PetscInt i, j, k, nz, nzL, row; 368 const PetscInt n = a->mbs, *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j; 369 const PetscInt *ajtmp, *bjtmp, *bdiag = b->diag, *pj, bs2 = a->bs2; 370 MatScalar *rtmp, *pc, *mwork, *v, *pv, *aa = a->a; 371 PetscInt flg; 372 PetscReal shift = info->shiftamount; 373 PetscBool allowzeropivot, zeropivotdetected; 374 375 PetscFunctionBegin; 376 allowzeropivot = PetscNot(A->erroriffailure); 377 PetscCall(ISGetIndices(isrow, &r)); 378 PetscCall(ISGetIndices(isicol, &ic)); 379 380 /* generate work space needed by the factorization */ 381 PetscCall(PetscMalloc2(bs2 * n, &rtmp, bs2, &mwork)); 382 PetscCall(PetscArrayzero(rtmp, bs2 * n)); 383 384 for (i = 0; i < n; i++) { 385 /* zero rtmp */ 386 /* L part */ 387 nz = bi[i + 1] - bi[i]; 388 bjtmp = bj + bi[i]; 389 for (j = 0; j < nz; j++) { PetscCall(PetscArrayzero(rtmp + bs2 * bjtmp[j], bs2)); } 390 391 /* U part */ 392 nz = bdiag[i] - bdiag[i + 1]; 393 bjtmp = bj + bdiag[i + 1] + 1; 394 for (j = 0; j < nz; j++) { PetscCall(PetscArrayzero(rtmp + bs2 * bjtmp[j], bs2)); } 395 396 /* load in initial (unfactored row) */ 397 nz = ai[r[i] + 1] - ai[r[i]]; 398 ajtmp = aj + ai[r[i]]; 399 v = aa + bs2 * ai[r[i]]; 400 for (j = 0; j < nz; j++) { PetscCall(PetscArraycpy(rtmp + bs2 * ic[ajtmp[j]], v + bs2 * j, bs2)); } 401 402 /* elimination */ 403 bjtmp = bj + bi[i]; 404 nzL = bi[i + 1] - bi[i]; 405 for (k = 0; k < nzL; k++) { 406 row = bjtmp[k]; 407 pc = rtmp + bs2 * row; 408 for (flg = 0, j = 0; j < bs2; j++) { 409 if (pc[j] != 0.0) { 410 flg = 1; 411 break; 412 } 413 } 414 if (flg) { 415 pv = b->a + bs2 * bdiag[row]; 416 /* PetscKernel_A_gets_A_times_B(bs,pc,pv,mwork); *pc = *pc * (*pv); */ 417 PetscCall(PetscKernel_A_gets_A_times_B_6(pc, pv, mwork)); 418 419 pj = b->j + bdiag[row + 1] + 1; /* beginning of U(row,:) */ 420 pv = b->a + bs2 * (bdiag[row + 1] + 1); 421 nz = bdiag[row] - bdiag[row + 1] - 1; /* num of entries inU(row,:), excluding diag */ 422 for (j = 0; j < nz; j++) { 423 /* PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j); */ 424 /* rtmp+bs2*pj[j] = rtmp+bs2*pj[j] - (*pc)*(pv+bs2*j) */ 425 v = rtmp + bs2 * pj[j]; 426 PetscCall(PetscKernel_A_gets_A_minus_B_times_C_6(v, pc, pv)); 427 pv += bs2; 428 } 429 PetscCall(PetscLogFlops(432.0 * nz + 396)); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */ 430 } 431 } 432 433 /* finished row so stick it into b->a */ 434 /* L part */ 435 pv = b->a + bs2 * bi[i]; 436 pj = b->j + bi[i]; 437 nz = bi[i + 1] - bi[i]; 438 for (j = 0; j < nz; j++) { PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2)); } 439 440 /* Mark diagonal and invert diagonal for simpler triangular solves */ 441 pv = b->a + bs2 * bdiag[i]; 442 pj = b->j + bdiag[i]; 443 PetscCall(PetscArraycpy(pv, rtmp + bs2 * pj[0], bs2)); 444 PetscCall(PetscKernel_A_gets_inverse_A_6(pv, shift, allowzeropivot, &zeropivotdetected)); 445 if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 446 447 /* U part */ 448 pv = b->a + bs2 * (bdiag[i + 1] + 1); 449 pj = b->j + bdiag[i + 1] + 1; 450 nz = bdiag[i] - bdiag[i + 1] - 1; 451 for (j = 0; j < nz; j++) { PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2)); } 452 } 453 454 PetscCall(PetscFree2(rtmp, mwork)); 455 PetscCall(ISRestoreIndices(isicol, &ic)); 456 PetscCall(ISRestoreIndices(isrow, &r)); 457 458 C->ops->solve = MatSolve_SeqBAIJ_6; 459 C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_6; 460 C->assembled = PETSC_TRUE; 461 462 PetscCall(PetscLogFlops(1.333333333333 * 6 * 6 * 6 * n)); /* from inverting diagonal blocks */ 463 PetscFunctionReturn(0); 464 } 465 466 PetscErrorCode MatLUFactorNumeric_SeqBAIJ_6_NaturalOrdering_inplace(Mat C, Mat A, const MatFactorInfo *info) { 467 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)C->data; 468 PetscInt i, j, n = a->mbs, *bi = b->i, *bj = b->j; 469 PetscInt *ajtmpold, *ajtmp, nz, row; 470 PetscInt *diag_offset = b->diag, *ai = a->i, *aj = a->j, *pj; 471 MatScalar *pv, *v, *rtmp, *pc, *w, *x; 472 MatScalar x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15; 473 MatScalar x16, x17, x18, x19, x20, x21, x22, x23, x24, x25; 474 MatScalar p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12, p13, p14, p15; 475 MatScalar p16, p17, p18, p19, p20, p21, p22, p23, p24, p25; 476 MatScalar m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, m11, m12, m13, m14, m15; 477 MatScalar m16, m17, m18, m19, m20, m21, m22, m23, m24, m25; 478 MatScalar p26, p27, p28, p29, p30, p31, p32, p33, p34, p35, p36; 479 MatScalar x26, x27, x28, x29, x30, x31, x32, x33, x34, x35, x36; 480 MatScalar m26, m27, m28, m29, m30, m31, m32, m33, m34, m35, m36; 481 MatScalar *ba = b->a, *aa = a->a; 482 PetscReal shift = info->shiftamount; 483 PetscBool allowzeropivot, zeropivotdetected; 484 485 PetscFunctionBegin; 486 allowzeropivot = PetscNot(A->erroriffailure); 487 PetscCall(PetscMalloc1(36 * (n + 1), &rtmp)); 488 for (i = 0; i < n; i++) { 489 nz = bi[i + 1] - bi[i]; 490 ajtmp = bj + bi[i]; 491 for (j = 0; j < nz; j++) { 492 x = rtmp + 36 * ajtmp[j]; 493 x[0] = x[1] = x[2] = x[3] = x[4] = x[5] = x[6] = x[7] = x[8] = x[9] = 0.0; 494 x[10] = x[11] = x[12] = x[13] = x[14] = x[15] = x[16] = x[17] = 0.0; 495 x[18] = x[19] = x[20] = x[21] = x[22] = x[23] = x[24] = x[25] = 0.0; 496 x[26] = x[27] = x[28] = x[29] = x[30] = x[31] = x[32] = x[33] = 0.0; 497 x[34] = x[35] = 0.0; 498 } 499 /* load in initial (unfactored row) */ 500 nz = ai[i + 1] - ai[i]; 501 ajtmpold = aj + ai[i]; 502 v = aa + 36 * ai[i]; 503 for (j = 0; j < nz; j++) { 504 x = rtmp + 36 * ajtmpold[j]; 505 x[0] = v[0]; 506 x[1] = v[1]; 507 x[2] = v[2]; 508 x[3] = v[3]; 509 x[4] = v[4]; 510 x[5] = v[5]; 511 x[6] = v[6]; 512 x[7] = v[7]; 513 x[8] = v[8]; 514 x[9] = v[9]; 515 x[10] = v[10]; 516 x[11] = v[11]; 517 x[12] = v[12]; 518 x[13] = v[13]; 519 x[14] = v[14]; 520 x[15] = v[15]; 521 x[16] = v[16]; 522 x[17] = v[17]; 523 x[18] = v[18]; 524 x[19] = v[19]; 525 x[20] = v[20]; 526 x[21] = v[21]; 527 x[22] = v[22]; 528 x[23] = v[23]; 529 x[24] = v[24]; 530 x[25] = v[25]; 531 x[26] = v[26]; 532 x[27] = v[27]; 533 x[28] = v[28]; 534 x[29] = v[29]; 535 x[30] = v[30]; 536 x[31] = v[31]; 537 x[32] = v[32]; 538 x[33] = v[33]; 539 x[34] = v[34]; 540 x[35] = v[35]; 541 v += 36; 542 } 543 row = *ajtmp++; 544 while (row < i) { 545 pc = rtmp + 36 * row; 546 p1 = pc[0]; 547 p2 = pc[1]; 548 p3 = pc[2]; 549 p4 = pc[3]; 550 p5 = pc[4]; 551 p6 = pc[5]; 552 p7 = pc[6]; 553 p8 = pc[7]; 554 p9 = pc[8]; 555 p10 = pc[9]; 556 p11 = pc[10]; 557 p12 = pc[11]; 558 p13 = pc[12]; 559 p14 = pc[13]; 560 p15 = pc[14]; 561 p16 = pc[15]; 562 p17 = pc[16]; 563 p18 = pc[17]; 564 p19 = pc[18]; 565 p20 = pc[19]; 566 p21 = pc[20]; 567 p22 = pc[21]; 568 p23 = pc[22]; 569 p24 = pc[23]; 570 p25 = pc[24]; 571 p26 = pc[25]; 572 p27 = pc[26]; 573 p28 = pc[27]; 574 p29 = pc[28]; 575 p30 = pc[29]; 576 p31 = pc[30]; 577 p32 = pc[31]; 578 p33 = pc[32]; 579 p34 = pc[33]; 580 p35 = pc[34]; 581 p36 = pc[35]; 582 if (p1 != 0.0 || p2 != 0.0 || p3 != 0.0 || p4 != 0.0 || p5 != 0.0 || p6 != 0.0 || p7 != 0.0 || p8 != 0.0 || p9 != 0.0 || p10 != 0.0 || p11 != 0.0 || p12 != 0.0 || p13 != 0.0 || p14 != 0.0 || p15 != 0.0 || p16 != 0.0 || p17 != 0.0 || p18 != 0.0 || p19 != 0.0 || p20 != 0.0 || p21 != 0.0 || p22 != 0.0 || p23 != 0.0 || p24 != 0.0 || p25 != 0.0 || p26 != 0.0 || p27 != 0.0 || p28 != 0.0 || p29 != 0.0 || p30 != 0.0 || p31 != 0.0 || p32 != 0.0 || p33 != 0.0 || p34 != 0.0 || p35 != 0.0 || p36 != 0.0) { 583 pv = ba + 36 * diag_offset[row]; 584 pj = bj + diag_offset[row] + 1; 585 x1 = pv[0]; 586 x2 = pv[1]; 587 x3 = pv[2]; 588 x4 = pv[3]; 589 x5 = pv[4]; 590 x6 = pv[5]; 591 x7 = pv[6]; 592 x8 = pv[7]; 593 x9 = pv[8]; 594 x10 = pv[9]; 595 x11 = pv[10]; 596 x12 = pv[11]; 597 x13 = pv[12]; 598 x14 = pv[13]; 599 x15 = pv[14]; 600 x16 = pv[15]; 601 x17 = pv[16]; 602 x18 = pv[17]; 603 x19 = pv[18]; 604 x20 = pv[19]; 605 x21 = pv[20]; 606 x22 = pv[21]; 607 x23 = pv[22]; 608 x24 = pv[23]; 609 x25 = pv[24]; 610 x26 = pv[25]; 611 x27 = pv[26]; 612 x28 = pv[27]; 613 x29 = pv[28]; 614 x30 = pv[29]; 615 x31 = pv[30]; 616 x32 = pv[31]; 617 x33 = pv[32]; 618 x34 = pv[33]; 619 x35 = pv[34]; 620 x36 = pv[35]; 621 pc[0] = m1 = p1 * x1 + p7 * x2 + p13 * x3 + p19 * x4 + p25 * x5 + p31 * x6; 622 pc[1] = m2 = p2 * x1 + p8 * x2 + p14 * x3 + p20 * x4 + p26 * x5 + p32 * x6; 623 pc[2] = m3 = p3 * x1 + p9 * x2 + p15 * x3 + p21 * x4 + p27 * x5 + p33 * x6; 624 pc[3] = m4 = p4 * x1 + p10 * x2 + p16 * x3 + p22 * x4 + p28 * x5 + p34 * x6; 625 pc[4] = m5 = p5 * x1 + p11 * x2 + p17 * x3 + p23 * x4 + p29 * x5 + p35 * x6; 626 pc[5] = m6 = p6 * x1 + p12 * x2 + p18 * x3 + p24 * x4 + p30 * x5 + p36 * x6; 627 628 pc[6] = m7 = p1 * x7 + p7 * x8 + p13 * x9 + p19 * x10 + p25 * x11 + p31 * x12; 629 pc[7] = m8 = p2 * x7 + p8 * x8 + p14 * x9 + p20 * x10 + p26 * x11 + p32 * x12; 630 pc[8] = m9 = p3 * x7 + p9 * x8 + p15 * x9 + p21 * x10 + p27 * x11 + p33 * x12; 631 pc[9] = m10 = p4 * x7 + p10 * x8 + p16 * x9 + p22 * x10 + p28 * x11 + p34 * x12; 632 pc[10] = m11 = p5 * x7 + p11 * x8 + p17 * x9 + p23 * x10 + p29 * x11 + p35 * x12; 633 pc[11] = m12 = p6 * x7 + p12 * x8 + p18 * x9 + p24 * x10 + p30 * x11 + p36 * x12; 634 635 pc[12] = m13 = p1 * x13 + p7 * x14 + p13 * x15 + p19 * x16 + p25 * x17 + p31 * x18; 636 pc[13] = m14 = p2 * x13 + p8 * x14 + p14 * x15 + p20 * x16 + p26 * x17 + p32 * x18; 637 pc[14] = m15 = p3 * x13 + p9 * x14 + p15 * x15 + p21 * x16 + p27 * x17 + p33 * x18; 638 pc[15] = m16 = p4 * x13 + p10 * x14 + p16 * x15 + p22 * x16 + p28 * x17 + p34 * x18; 639 pc[16] = m17 = p5 * x13 + p11 * x14 + p17 * x15 + p23 * x16 + p29 * x17 + p35 * x18; 640 pc[17] = m18 = p6 * x13 + p12 * x14 + p18 * x15 + p24 * x16 + p30 * x17 + p36 * x18; 641 642 pc[18] = m19 = p1 * x19 + p7 * x20 + p13 * x21 + p19 * x22 + p25 * x23 + p31 * x24; 643 pc[19] = m20 = p2 * x19 + p8 * x20 + p14 * x21 + p20 * x22 + p26 * x23 + p32 * x24; 644 pc[20] = m21 = p3 * x19 + p9 * x20 + p15 * x21 + p21 * x22 + p27 * x23 + p33 * x24; 645 pc[21] = m22 = p4 * x19 + p10 * x20 + p16 * x21 + p22 * x22 + p28 * x23 + p34 * x24; 646 pc[22] = m23 = p5 * x19 + p11 * x20 + p17 * x21 + p23 * x22 + p29 * x23 + p35 * x24; 647 pc[23] = m24 = p6 * x19 + p12 * x20 + p18 * x21 + p24 * x22 + p30 * x23 + p36 * x24; 648 649 pc[24] = m25 = p1 * x25 + p7 * x26 + p13 * x27 + p19 * x28 + p25 * x29 + p31 * x30; 650 pc[25] = m26 = p2 * x25 + p8 * x26 + p14 * x27 + p20 * x28 + p26 * x29 + p32 * x30; 651 pc[26] = m27 = p3 * x25 + p9 * x26 + p15 * x27 + p21 * x28 + p27 * x29 + p33 * x30; 652 pc[27] = m28 = p4 * x25 + p10 * x26 + p16 * x27 + p22 * x28 + p28 * x29 + p34 * x30; 653 pc[28] = m29 = p5 * x25 + p11 * x26 + p17 * x27 + p23 * x28 + p29 * x29 + p35 * x30; 654 pc[29] = m30 = p6 * x25 + p12 * x26 + p18 * x27 + p24 * x28 + p30 * x29 + p36 * x30; 655 656 pc[30] = m31 = p1 * x31 + p7 * x32 + p13 * x33 + p19 * x34 + p25 * x35 + p31 * x36; 657 pc[31] = m32 = p2 * x31 + p8 * x32 + p14 * x33 + p20 * x34 + p26 * x35 + p32 * x36; 658 pc[32] = m33 = p3 * x31 + p9 * x32 + p15 * x33 + p21 * x34 + p27 * x35 + p33 * x36; 659 pc[33] = m34 = p4 * x31 + p10 * x32 + p16 * x33 + p22 * x34 + p28 * x35 + p34 * x36; 660 pc[34] = m35 = p5 * x31 + p11 * x32 + p17 * x33 + p23 * x34 + p29 * x35 + p35 * x36; 661 pc[35] = m36 = p6 * x31 + p12 * x32 + p18 * x33 + p24 * x34 + p30 * x35 + p36 * x36; 662 663 nz = bi[row + 1] - diag_offset[row] - 1; 664 pv += 36; 665 for (j = 0; j < nz; j++) { 666 x1 = pv[0]; 667 x2 = pv[1]; 668 x3 = pv[2]; 669 x4 = pv[3]; 670 x5 = pv[4]; 671 x6 = pv[5]; 672 x7 = pv[6]; 673 x8 = pv[7]; 674 x9 = pv[8]; 675 x10 = pv[9]; 676 x11 = pv[10]; 677 x12 = pv[11]; 678 x13 = pv[12]; 679 x14 = pv[13]; 680 x15 = pv[14]; 681 x16 = pv[15]; 682 x17 = pv[16]; 683 x18 = pv[17]; 684 x19 = pv[18]; 685 x20 = pv[19]; 686 x21 = pv[20]; 687 x22 = pv[21]; 688 x23 = pv[22]; 689 x24 = pv[23]; 690 x25 = pv[24]; 691 x26 = pv[25]; 692 x27 = pv[26]; 693 x28 = pv[27]; 694 x29 = pv[28]; 695 x30 = pv[29]; 696 x31 = pv[30]; 697 x32 = pv[31]; 698 x33 = pv[32]; 699 x34 = pv[33]; 700 x35 = pv[34]; 701 x36 = pv[35]; 702 x = rtmp + 36 * pj[j]; 703 x[0] -= m1 * x1 + m7 * x2 + m13 * x3 + m19 * x4 + m25 * x5 + m31 * x6; 704 x[1] -= m2 * x1 + m8 * x2 + m14 * x3 + m20 * x4 + m26 * x5 + m32 * x6; 705 x[2] -= m3 * x1 + m9 * x2 + m15 * x3 + m21 * x4 + m27 * x5 + m33 * x6; 706 x[3] -= m4 * x1 + m10 * x2 + m16 * x3 + m22 * x4 + m28 * x5 + m34 * x6; 707 x[4] -= m5 * x1 + m11 * x2 + m17 * x3 + m23 * x4 + m29 * x5 + m35 * x6; 708 x[5] -= m6 * x1 + m12 * x2 + m18 * x3 + m24 * x4 + m30 * x5 + m36 * x6; 709 710 x[6] -= m1 * x7 + m7 * x8 + m13 * x9 + m19 * x10 + m25 * x11 + m31 * x12; 711 x[7] -= m2 * x7 + m8 * x8 + m14 * x9 + m20 * x10 + m26 * x11 + m32 * x12; 712 x[8] -= m3 * x7 + m9 * x8 + m15 * x9 + m21 * x10 + m27 * x11 + m33 * x12; 713 x[9] -= m4 * x7 + m10 * x8 + m16 * x9 + m22 * x10 + m28 * x11 + m34 * x12; 714 x[10] -= m5 * x7 + m11 * x8 + m17 * x9 + m23 * x10 + m29 * x11 + m35 * x12; 715 x[11] -= m6 * x7 + m12 * x8 + m18 * x9 + m24 * x10 + m30 * x11 + m36 * x12; 716 717 x[12] -= m1 * x13 + m7 * x14 + m13 * x15 + m19 * x16 + m25 * x17 + m31 * x18; 718 x[13] -= m2 * x13 + m8 * x14 + m14 * x15 + m20 * x16 + m26 * x17 + m32 * x18; 719 x[14] -= m3 * x13 + m9 * x14 + m15 * x15 + m21 * x16 + m27 * x17 + m33 * x18; 720 x[15] -= m4 * x13 + m10 * x14 + m16 * x15 + m22 * x16 + m28 * x17 + m34 * x18; 721 x[16] -= m5 * x13 + m11 * x14 + m17 * x15 + m23 * x16 + m29 * x17 + m35 * x18; 722 x[17] -= m6 * x13 + m12 * x14 + m18 * x15 + m24 * x16 + m30 * x17 + m36 * x18; 723 724 x[18] -= m1 * x19 + m7 * x20 + m13 * x21 + m19 * x22 + m25 * x23 + m31 * x24; 725 x[19] -= m2 * x19 + m8 * x20 + m14 * x21 + m20 * x22 + m26 * x23 + m32 * x24; 726 x[20] -= m3 * x19 + m9 * x20 + m15 * x21 + m21 * x22 + m27 * x23 + m33 * x24; 727 x[21] -= m4 * x19 + m10 * x20 + m16 * x21 + m22 * x22 + m28 * x23 + m34 * x24; 728 x[22] -= m5 * x19 + m11 * x20 + m17 * x21 + m23 * x22 + m29 * x23 + m35 * x24; 729 x[23] -= m6 * x19 + m12 * x20 + m18 * x21 + m24 * x22 + m30 * x23 + m36 * x24; 730 731 x[24] -= m1 * x25 + m7 * x26 + m13 * x27 + m19 * x28 + m25 * x29 + m31 * x30; 732 x[25] -= m2 * x25 + m8 * x26 + m14 * x27 + m20 * x28 + m26 * x29 + m32 * x30; 733 x[26] -= m3 * x25 + m9 * x26 + m15 * x27 + m21 * x28 + m27 * x29 + m33 * x30; 734 x[27] -= m4 * x25 + m10 * x26 + m16 * x27 + m22 * x28 + m28 * x29 + m34 * x30; 735 x[28] -= m5 * x25 + m11 * x26 + m17 * x27 + m23 * x28 + m29 * x29 + m35 * x30; 736 x[29] -= m6 * x25 + m12 * x26 + m18 * x27 + m24 * x28 + m30 * x29 + m36 * x30; 737 738 x[30] -= m1 * x31 + m7 * x32 + m13 * x33 + m19 * x34 + m25 * x35 + m31 * x36; 739 x[31] -= m2 * x31 + m8 * x32 + m14 * x33 + m20 * x34 + m26 * x35 + m32 * x36; 740 x[32] -= m3 * x31 + m9 * x32 + m15 * x33 + m21 * x34 + m27 * x35 + m33 * x36; 741 x[33] -= m4 * x31 + m10 * x32 + m16 * x33 + m22 * x34 + m28 * x35 + m34 * x36; 742 x[34] -= m5 * x31 + m11 * x32 + m17 * x33 + m23 * x34 + m29 * x35 + m35 * x36; 743 x[35] -= m6 * x31 + m12 * x32 + m18 * x33 + m24 * x34 + m30 * x35 + m36 * x36; 744 745 pv += 36; 746 } 747 PetscCall(PetscLogFlops(432.0 * nz + 396.0)); 748 } 749 row = *ajtmp++; 750 } 751 /* finished row so stick it into b->a */ 752 pv = ba + 36 * bi[i]; 753 pj = bj + bi[i]; 754 nz = bi[i + 1] - bi[i]; 755 for (j = 0; j < nz; j++) { 756 x = rtmp + 36 * pj[j]; 757 pv[0] = x[0]; 758 pv[1] = x[1]; 759 pv[2] = x[2]; 760 pv[3] = x[3]; 761 pv[4] = x[4]; 762 pv[5] = x[5]; 763 pv[6] = x[6]; 764 pv[7] = x[7]; 765 pv[8] = x[8]; 766 pv[9] = x[9]; 767 pv[10] = x[10]; 768 pv[11] = x[11]; 769 pv[12] = x[12]; 770 pv[13] = x[13]; 771 pv[14] = x[14]; 772 pv[15] = x[15]; 773 pv[16] = x[16]; 774 pv[17] = x[17]; 775 pv[18] = x[18]; 776 pv[19] = x[19]; 777 pv[20] = x[20]; 778 pv[21] = x[21]; 779 pv[22] = x[22]; 780 pv[23] = x[23]; 781 pv[24] = x[24]; 782 pv[25] = x[25]; 783 pv[26] = x[26]; 784 pv[27] = x[27]; 785 pv[28] = x[28]; 786 pv[29] = x[29]; 787 pv[30] = x[30]; 788 pv[31] = x[31]; 789 pv[32] = x[32]; 790 pv[33] = x[33]; 791 pv[34] = x[34]; 792 pv[35] = x[35]; 793 pv += 36; 794 } 795 /* invert diagonal block */ 796 w = ba + 36 * diag_offset[i]; 797 PetscCall(PetscKernel_A_gets_inverse_A_6(w, shift, allowzeropivot, &zeropivotdetected)); 798 if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 799 } 800 801 PetscCall(PetscFree(rtmp)); 802 803 C->ops->solve = MatSolve_SeqBAIJ_6_NaturalOrdering_inplace; 804 C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_6_NaturalOrdering_inplace; 805 C->assembled = PETSC_TRUE; 806 807 PetscCall(PetscLogFlops(1.333333333333 * 6 * 6 * 6 * b->mbs)); /* from inverting diagonal blocks */ 808 PetscFunctionReturn(0); 809 } 810 811 PetscErrorCode MatLUFactorNumeric_SeqBAIJ_6_NaturalOrdering(Mat B, Mat A, const MatFactorInfo *info) { 812 Mat C = B; 813 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)C->data; 814 PetscInt i, j, k, nz, nzL, row; 815 const PetscInt n = a->mbs, *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j; 816 const PetscInt *ajtmp, *bjtmp, *bdiag = b->diag, *pj, bs2 = a->bs2; 817 MatScalar *rtmp, *pc, *mwork, *v, *pv, *aa = a->a; 818 PetscInt flg; 819 PetscReal shift = info->shiftamount; 820 PetscBool allowzeropivot, zeropivotdetected; 821 822 PetscFunctionBegin; 823 allowzeropivot = PetscNot(A->erroriffailure); 824 825 /* generate work space needed by the factorization */ 826 PetscCall(PetscMalloc2(bs2 * n, &rtmp, bs2, &mwork)); 827 PetscCall(PetscArrayzero(rtmp, bs2 * n)); 828 829 for (i = 0; i < n; i++) { 830 /* zero rtmp */ 831 /* L part */ 832 nz = bi[i + 1] - bi[i]; 833 bjtmp = bj + bi[i]; 834 for (j = 0; j < nz; j++) { PetscCall(PetscArrayzero(rtmp + bs2 * bjtmp[j], bs2)); } 835 836 /* U part */ 837 nz = bdiag[i] - bdiag[i + 1]; 838 bjtmp = bj + bdiag[i + 1] + 1; 839 for (j = 0; j < nz; j++) { PetscCall(PetscArrayzero(rtmp + bs2 * bjtmp[j], bs2)); } 840 841 /* load in initial (unfactored row) */ 842 nz = ai[i + 1] - ai[i]; 843 ajtmp = aj + ai[i]; 844 v = aa + bs2 * ai[i]; 845 for (j = 0; j < nz; j++) { PetscCall(PetscArraycpy(rtmp + bs2 * ajtmp[j], v + bs2 * j, bs2)); } 846 847 /* elimination */ 848 bjtmp = bj + bi[i]; 849 nzL = bi[i + 1] - bi[i]; 850 for (k = 0; k < nzL; k++) { 851 row = bjtmp[k]; 852 pc = rtmp + bs2 * row; 853 for (flg = 0, j = 0; j < bs2; j++) { 854 if (pc[j] != 0.0) { 855 flg = 1; 856 break; 857 } 858 } 859 if (flg) { 860 pv = b->a + bs2 * bdiag[row]; 861 /* PetscKernel_A_gets_A_times_B(bs,pc,pv,mwork); *pc = *pc * (*pv); */ 862 PetscCall(PetscKernel_A_gets_A_times_B_6(pc, pv, mwork)); 863 864 pj = b->j + bdiag[row + 1] + 1; /* beginning of U(row,:) */ 865 pv = b->a + bs2 * (bdiag[row + 1] + 1); 866 nz = bdiag[row] - bdiag[row + 1] - 1; /* num of entries inU(row,:), excluding diag */ 867 for (j = 0; j < nz; j++) { 868 /* PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j); */ 869 /* rtmp+bs2*pj[j] = rtmp+bs2*pj[j] - (*pc)*(pv+bs2*j) */ 870 v = rtmp + bs2 * pj[j]; 871 PetscCall(PetscKernel_A_gets_A_minus_B_times_C_6(v, pc, pv)); 872 pv += bs2; 873 } 874 PetscCall(PetscLogFlops(432.0 * nz + 396)); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */ 875 } 876 } 877 878 /* finished row so stick it into b->a */ 879 /* L part */ 880 pv = b->a + bs2 * bi[i]; 881 pj = b->j + bi[i]; 882 nz = bi[i + 1] - bi[i]; 883 for (j = 0; j < nz; j++) { PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2)); } 884 885 /* Mark diagonal and invert diagonal for simpler triangular solves */ 886 pv = b->a + bs2 * bdiag[i]; 887 pj = b->j + bdiag[i]; 888 PetscCall(PetscArraycpy(pv, rtmp + bs2 * pj[0], bs2)); 889 PetscCall(PetscKernel_A_gets_inverse_A_6(pv, shift, allowzeropivot, &zeropivotdetected)); 890 if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 891 892 /* U part */ 893 pv = b->a + bs2 * (bdiag[i + 1] + 1); 894 pj = b->j + bdiag[i + 1] + 1; 895 nz = bdiag[i] - bdiag[i + 1] - 1; 896 for (j = 0; j < nz; j++) { PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2)); } 897 } 898 PetscCall(PetscFree2(rtmp, mwork)); 899 900 C->ops->solve = MatSolve_SeqBAIJ_6_NaturalOrdering; 901 C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_6_NaturalOrdering; 902 C->assembled = PETSC_TRUE; 903 904 PetscCall(PetscLogFlops(1.333333333333 * 6 * 6 * 6 * n)); /* from inverting diagonal blocks */ 905 PetscFunctionReturn(0); 906 } 907