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