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