1 /* 2 Factorization code for BAIJ format. 3 */ 4 #include <../src/mat/impls/baij/seq/baij.h> 5 #include <petsc/private/kernels/blockinvert.h> 6 /* 7 Version for when blocks are 7 by 7 8 */ 9 PetscErrorCode MatLUFactorNumeric_SeqBAIJ_7_inplace(Mat C, Mat A, const MatFactorInfo *info) 10 { 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(PETSC_SUCCESS); 454 } 455 456 PetscErrorCode MatLUFactorNumeric_SeqBAIJ_7(Mat B, Mat A, const MatFactorInfo *info) 457 { 458 Mat C = B; 459 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)C->data; 460 IS isrow = b->row, isicol = b->icol; 461 const PetscInt *r, *ic; 462 PetscInt i, j, k, nz, nzL, row; 463 const PetscInt n = a->mbs, *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j; 464 const PetscInt *ajtmp, *bjtmp, *bdiag = b->diag, *pj, bs2 = a->bs2; 465 MatScalar *rtmp, *pc, *mwork, *v, *pv, *aa = a->a; 466 PetscInt flg; 467 PetscReal shift = info->shiftamount; 468 PetscBool allowzeropivot, zeropivotdetected; 469 470 PetscFunctionBegin; 471 allowzeropivot = PetscNot(A->erroriffailure); 472 PetscCall(ISGetIndices(isrow, &r)); 473 PetscCall(ISGetIndices(isicol, &ic)); 474 475 /* generate work space needed by the factorization */ 476 PetscCall(PetscMalloc2(bs2 * n, &rtmp, bs2, &mwork)); 477 PetscCall(PetscArrayzero(rtmp, bs2 * n)); 478 479 for (i = 0; i < n; i++) { 480 /* zero rtmp */ 481 /* L part */ 482 nz = bi[i + 1] - bi[i]; 483 bjtmp = bj + bi[i]; 484 for (j = 0; j < nz; j++) PetscCall(PetscArrayzero(rtmp + bs2 * bjtmp[j], bs2)); 485 486 /* U part */ 487 nz = bdiag[i] - bdiag[i + 1]; 488 bjtmp = bj + bdiag[i + 1] + 1; 489 for (j = 0; j < nz; j++) PetscCall(PetscArrayzero(rtmp + bs2 * bjtmp[j], bs2)); 490 491 /* load in initial (unfactored row) */ 492 nz = ai[r[i] + 1] - ai[r[i]]; 493 ajtmp = aj + ai[r[i]]; 494 v = aa + bs2 * ai[r[i]]; 495 for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(rtmp + bs2 * ic[ajtmp[j]], v + bs2 * j, bs2)); 496 497 /* elimination */ 498 bjtmp = bj + bi[i]; 499 nzL = bi[i + 1] - bi[i]; 500 for (k = 0; k < nzL; k++) { 501 row = bjtmp[k]; 502 pc = rtmp + bs2 * row; 503 for (flg = 0, j = 0; j < bs2; j++) { 504 if (pc[j] != 0.0) { 505 flg = 1; 506 break; 507 } 508 } 509 if (flg) { 510 pv = b->a + bs2 * bdiag[row]; 511 /* PetscKernel_A_gets_A_times_B(bs,pc,pv,mwork); *pc = *pc * (*pv); */ 512 PetscCall(PetscKernel_A_gets_A_times_B_7(pc, pv, mwork)); 513 514 pj = b->j + bdiag[row + 1] + 1; /* beginning of U(row,:) */ 515 pv = b->a + bs2 * (bdiag[row + 1] + 1); 516 nz = bdiag[row] - bdiag[row + 1] - 1; /* num of entries inU(row,:), excluding diag */ 517 for (j = 0; j < nz; j++) { 518 /* PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j); */ 519 /* rtmp+bs2*pj[j] = rtmp+bs2*pj[j] - (*pc)*(pv+bs2*j) */ 520 v = rtmp + bs2 * pj[j]; 521 PetscCall(PetscKernel_A_gets_A_minus_B_times_C_7(v, pc, pv)); 522 pv += bs2; 523 } 524 PetscCall(PetscLogFlops(686.0 * nz + 637)); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */ 525 } 526 } 527 528 /* finished row so stick it into b->a */ 529 /* L part */ 530 pv = b->a + bs2 * bi[i]; 531 pj = b->j + bi[i]; 532 nz = bi[i + 1] - bi[i]; 533 for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2)); 534 535 /* Mark diagonal and invert diagonal for simpler triangular solves */ 536 pv = b->a + bs2 * bdiag[i]; 537 pj = b->j + bdiag[i]; 538 PetscCall(PetscArraycpy(pv, rtmp + bs2 * pj[0], bs2)); 539 PetscCall(PetscKernel_A_gets_inverse_A_7(pv, shift, allowzeropivot, &zeropivotdetected)); 540 if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 541 542 /* U part */ 543 pv = b->a + bs2 * (bdiag[i + 1] + 1); 544 pj = b->j + bdiag[i + 1] + 1; 545 nz = bdiag[i] - bdiag[i + 1] - 1; 546 for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2)); 547 } 548 549 PetscCall(PetscFree2(rtmp, mwork)); 550 PetscCall(ISRestoreIndices(isicol, &ic)); 551 PetscCall(ISRestoreIndices(isrow, &r)); 552 553 C->ops->solve = MatSolve_SeqBAIJ_7; 554 C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_7; 555 C->assembled = PETSC_TRUE; 556 557 PetscCall(PetscLogFlops(1.333333333333 * 7 * 7 * 7 * n)); /* from inverting diagonal blocks */ 558 PetscFunctionReturn(PETSC_SUCCESS); 559 } 560 561 PetscErrorCode MatLUFactorNumeric_SeqBAIJ_7_NaturalOrdering_inplace(Mat C, Mat A, const MatFactorInfo *info) 562 { 563 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)C->data; 564 PetscInt i, j, n = a->mbs, *bi = b->i, *bj = b->j; 565 PetscInt *ajtmpold, *ajtmp, nz, row; 566 PetscInt *diag_offset = b->diag, *ai = a->i, *aj = a->j, *pj; 567 MatScalar *pv, *v, *rtmp, *pc, *w, *x; 568 MatScalar x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15; 569 MatScalar x16, x17, x18, x19, x20, x21, x22, x23, x24, x25; 570 MatScalar p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12, p13, p14, p15; 571 MatScalar p16, p17, p18, p19, p20, p21, p22, p23, p24, p25; 572 MatScalar m1, m2, m3, m4, m5, m6, m7, m8, m9, m10, m11, m12, m13, m14, m15; 573 MatScalar m16, m17, m18, m19, m20, m21, m22, m23, m24, m25; 574 MatScalar p26, p27, p28, p29, p30, p31, p32, p33, p34, p35, p36; 575 MatScalar p37, p38, p39, p40, p41, p42, p43, p44, p45, p46, p47, p48, p49; 576 MatScalar x26, x27, x28, x29, x30, x31, x32, x33, x34, x35, x36; 577 MatScalar x37, x38, x39, x40, x41, x42, x43, x44, x45, x46, x47, x48, x49; 578 MatScalar m26, m27, m28, m29, m30, m31, m32, m33, m34, m35, m36; 579 MatScalar m37, m38, m39, m40, m41, m42, m43, m44, m45, m46, m47, m48, m49; 580 MatScalar *ba = b->a, *aa = a->a; 581 PetscReal shift = info->shiftamount; 582 PetscBool allowzeropivot, zeropivotdetected; 583 584 PetscFunctionBegin; 585 allowzeropivot = PetscNot(A->erroriffailure); 586 PetscCall(PetscMalloc1(49 * (n + 1), &rtmp)); 587 for (i = 0; i < n; i++) { 588 nz = bi[i + 1] - bi[i]; 589 ajtmp = bj + bi[i]; 590 for (j = 0; j < nz; j++) { 591 x = rtmp + 49 * ajtmp[j]; 592 x[0] = x[1] = x[2] = x[3] = x[4] = x[5] = x[6] = x[7] = x[8] = x[9] = 0.0; 593 x[10] = x[11] = x[12] = x[13] = x[14] = x[15] = x[16] = x[17] = 0.0; 594 x[18] = x[19] = x[20] = x[21] = x[22] = x[23] = x[24] = x[25] = 0.0; 595 x[26] = x[27] = x[28] = x[29] = x[30] = x[31] = x[32] = x[33] = 0.0; 596 x[34] = x[35] = x[36] = x[37] = x[38] = x[39] = x[40] = x[41] = 0.0; 597 x[42] = x[43] = x[44] = x[45] = x[46] = x[47] = x[48] = 0.0; 598 } 599 /* load in initial (unfactored row) */ 600 nz = ai[i + 1] - ai[i]; 601 ajtmpold = aj + ai[i]; 602 v = aa + 49 * ai[i]; 603 for (j = 0; j < nz; j++) { 604 x = rtmp + 49 * ajtmpold[j]; 605 x[0] = v[0]; 606 x[1] = v[1]; 607 x[2] = v[2]; 608 x[3] = v[3]; 609 x[4] = v[4]; 610 x[5] = v[5]; 611 x[6] = v[6]; 612 x[7] = v[7]; 613 x[8] = v[8]; 614 x[9] = v[9]; 615 x[10] = v[10]; 616 x[11] = v[11]; 617 x[12] = v[12]; 618 x[13] = v[13]; 619 x[14] = v[14]; 620 x[15] = v[15]; 621 x[16] = v[16]; 622 x[17] = v[17]; 623 x[18] = v[18]; 624 x[19] = v[19]; 625 x[20] = v[20]; 626 x[21] = v[21]; 627 x[22] = v[22]; 628 x[23] = v[23]; 629 x[24] = v[24]; 630 x[25] = v[25]; 631 x[26] = v[26]; 632 x[27] = v[27]; 633 x[28] = v[28]; 634 x[29] = v[29]; 635 x[30] = v[30]; 636 x[31] = v[31]; 637 x[32] = v[32]; 638 x[33] = v[33]; 639 x[34] = v[34]; 640 x[35] = v[35]; 641 x[36] = v[36]; 642 x[37] = v[37]; 643 x[38] = v[38]; 644 x[39] = v[39]; 645 x[40] = v[40]; 646 x[41] = v[41]; 647 x[42] = v[42]; 648 x[43] = v[43]; 649 x[44] = v[44]; 650 x[45] = v[45]; 651 x[46] = v[46]; 652 x[47] = v[47]; 653 x[48] = v[48]; 654 v += 49; 655 } 656 row = *ajtmp++; 657 while (row < i) { 658 pc = rtmp + 49 * row; 659 p1 = pc[0]; 660 p2 = pc[1]; 661 p3 = pc[2]; 662 p4 = pc[3]; 663 p5 = pc[4]; 664 p6 = pc[5]; 665 p7 = pc[6]; 666 p8 = pc[7]; 667 p9 = pc[8]; 668 p10 = pc[9]; 669 p11 = pc[10]; 670 p12 = pc[11]; 671 p13 = pc[12]; 672 p14 = pc[13]; 673 p15 = pc[14]; 674 p16 = pc[15]; 675 p17 = pc[16]; 676 p18 = pc[17]; 677 p19 = pc[18]; 678 p20 = pc[19]; 679 p21 = pc[20]; 680 p22 = pc[21]; 681 p23 = pc[22]; 682 p24 = pc[23]; 683 p25 = pc[24]; 684 p26 = pc[25]; 685 p27 = pc[26]; 686 p28 = pc[27]; 687 p29 = pc[28]; 688 p30 = pc[29]; 689 p31 = pc[30]; 690 p32 = pc[31]; 691 p33 = pc[32]; 692 p34 = pc[33]; 693 p35 = pc[34]; 694 p36 = pc[35]; 695 p37 = pc[36]; 696 p38 = pc[37]; 697 p39 = pc[38]; 698 p40 = pc[39]; 699 p41 = pc[40]; 700 p42 = pc[41]; 701 p43 = pc[42]; 702 p44 = pc[43]; 703 p45 = pc[44]; 704 p46 = pc[45]; 705 p47 = pc[46]; 706 p48 = pc[47]; 707 p49 = pc[48]; 708 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) { 709 pv = ba + 49 * diag_offset[row]; 710 pj = bj + diag_offset[row] + 1; 711 x1 = pv[0]; 712 x2 = pv[1]; 713 x3 = pv[2]; 714 x4 = pv[3]; 715 x5 = pv[4]; 716 x6 = pv[5]; 717 x7 = pv[6]; 718 x8 = pv[7]; 719 x9 = pv[8]; 720 x10 = pv[9]; 721 x11 = pv[10]; 722 x12 = pv[11]; 723 x13 = pv[12]; 724 x14 = pv[13]; 725 x15 = pv[14]; 726 x16 = pv[15]; 727 x17 = pv[16]; 728 x18 = pv[17]; 729 x19 = pv[18]; 730 x20 = pv[19]; 731 x21 = pv[20]; 732 x22 = pv[21]; 733 x23 = pv[22]; 734 x24 = pv[23]; 735 x25 = pv[24]; 736 x26 = pv[25]; 737 x27 = pv[26]; 738 x28 = pv[27]; 739 x29 = pv[28]; 740 x30 = pv[29]; 741 x31 = pv[30]; 742 x32 = pv[31]; 743 x33 = pv[32]; 744 x34 = pv[33]; 745 x35 = pv[34]; 746 x36 = pv[35]; 747 x37 = pv[36]; 748 x38 = pv[37]; 749 x39 = pv[38]; 750 x40 = pv[39]; 751 x41 = pv[40]; 752 x42 = pv[41]; 753 x43 = pv[42]; 754 x44 = pv[43]; 755 x45 = pv[44]; 756 x46 = pv[45]; 757 x47 = pv[46]; 758 x48 = pv[47]; 759 x49 = pv[48]; 760 pc[0] = m1 = p1 * x1 + p8 * x2 + p15 * x3 + p22 * x4 + p29 * x5 + p36 * x6 + p43 * x7; 761 pc[1] = m2 = p2 * x1 + p9 * x2 + p16 * x3 + p23 * x4 + p30 * x5 + p37 * x6 + p44 * x7; 762 pc[2] = m3 = p3 * x1 + p10 * x2 + p17 * x3 + p24 * x4 + p31 * x5 + p38 * x6 + p45 * x7; 763 pc[3] = m4 = p4 * x1 + p11 * x2 + p18 * x3 + p25 * x4 + p32 * x5 + p39 * x6 + p46 * x7; 764 pc[4] = m5 = p5 * x1 + p12 * x2 + p19 * x3 + p26 * x4 + p33 * x5 + p40 * x6 + p47 * x7; 765 pc[5] = m6 = p6 * x1 + p13 * x2 + p20 * x3 + p27 * x4 + p34 * x5 + p41 * x6 + p48 * x7; 766 pc[6] = m7 = p7 * x1 + p14 * x2 + p21 * x3 + p28 * x4 + p35 * x5 + p42 * x6 + p49 * x7; 767 768 pc[7] = m8 = p1 * x8 + p8 * x9 + p15 * x10 + p22 * x11 + p29 * x12 + p36 * x13 + p43 * x14; 769 pc[8] = m9 = p2 * x8 + p9 * x9 + p16 * x10 + p23 * x11 + p30 * x12 + p37 * x13 + p44 * x14; 770 pc[9] = m10 = p3 * x8 + p10 * x9 + p17 * x10 + p24 * x11 + p31 * x12 + p38 * x13 + p45 * x14; 771 pc[10] = m11 = p4 * x8 + p11 * x9 + p18 * x10 + p25 * x11 + p32 * x12 + p39 * x13 + p46 * x14; 772 pc[11] = m12 = p5 * x8 + p12 * x9 + p19 * x10 + p26 * x11 + p33 * x12 + p40 * x13 + p47 * x14; 773 pc[12] = m13 = p6 * x8 + p13 * x9 + p20 * x10 + p27 * x11 + p34 * x12 + p41 * x13 + p48 * x14; 774 pc[13] = m14 = p7 * x8 + p14 * x9 + p21 * x10 + p28 * x11 + p35 * x12 + p42 * x13 + p49 * x14; 775 776 pc[14] = m15 = p1 * x15 + p8 * x16 + p15 * x17 + p22 * x18 + p29 * x19 + p36 * x20 + p43 * x21; 777 pc[15] = m16 = p2 * x15 + p9 * x16 + p16 * x17 + p23 * x18 + p30 * x19 + p37 * x20 + p44 * x21; 778 pc[16] = m17 = p3 * x15 + p10 * x16 + p17 * x17 + p24 * x18 + p31 * x19 + p38 * x20 + p45 * x21; 779 pc[17] = m18 = p4 * x15 + p11 * x16 + p18 * x17 + p25 * x18 + p32 * x19 + p39 * x20 + p46 * x21; 780 pc[18] = m19 = p5 * x15 + p12 * x16 + p19 * x17 + p26 * x18 + p33 * x19 + p40 * x20 + p47 * x21; 781 pc[19] = m20 = p6 * x15 + p13 * x16 + p20 * x17 + p27 * x18 + p34 * x19 + p41 * x20 + p48 * x21; 782 pc[20] = m21 = p7 * x15 + p14 * x16 + p21 * x17 + p28 * x18 + p35 * x19 + p42 * x20 + p49 * x21; 783 784 pc[21] = m22 = p1 * x22 + p8 * x23 + p15 * x24 + p22 * x25 + p29 * x26 + p36 * x27 + p43 * x28; 785 pc[22] = m23 = p2 * x22 + p9 * x23 + p16 * x24 + p23 * x25 + p30 * x26 + p37 * x27 + p44 * x28; 786 pc[23] = m24 = p3 * x22 + p10 * x23 + p17 * x24 + p24 * x25 + p31 * x26 + p38 * x27 + p45 * x28; 787 pc[24] = m25 = p4 * x22 + p11 * x23 + p18 * x24 + p25 * x25 + p32 * x26 + p39 * x27 + p46 * x28; 788 pc[25] = m26 = p5 * x22 + p12 * x23 + p19 * x24 + p26 * x25 + p33 * x26 + p40 * x27 + p47 * x28; 789 pc[26] = m27 = p6 * x22 + p13 * x23 + p20 * x24 + p27 * x25 + p34 * x26 + p41 * x27 + p48 * x28; 790 pc[27] = m28 = p7 * x22 + p14 * x23 + p21 * x24 + p28 * x25 + p35 * x26 + p42 * x27 + p49 * x28; 791 792 pc[28] = m29 = p1 * x29 + p8 * x30 + p15 * x31 + p22 * x32 + p29 * x33 + p36 * x34 + p43 * x35; 793 pc[29] = m30 = p2 * x29 + p9 * x30 + p16 * x31 + p23 * x32 + p30 * x33 + p37 * x34 + p44 * x35; 794 pc[30] = m31 = p3 * x29 + p10 * x30 + p17 * x31 + p24 * x32 + p31 * x33 + p38 * x34 + p45 * x35; 795 pc[31] = m32 = p4 * x29 + p11 * x30 + p18 * x31 + p25 * x32 + p32 * x33 + p39 * x34 + p46 * x35; 796 pc[32] = m33 = p5 * x29 + p12 * x30 + p19 * x31 + p26 * x32 + p33 * x33 + p40 * x34 + p47 * x35; 797 pc[33] = m34 = p6 * x29 + p13 * x30 + p20 * x31 + p27 * x32 + p34 * x33 + p41 * x34 + p48 * x35; 798 pc[34] = m35 = p7 * x29 + p14 * x30 + p21 * x31 + p28 * x32 + p35 * x33 + p42 * x34 + p49 * x35; 799 800 pc[35] = m36 = p1 * x36 + p8 * x37 + p15 * x38 + p22 * x39 + p29 * x40 + p36 * x41 + p43 * x42; 801 pc[36] = m37 = p2 * x36 + p9 * x37 + p16 * x38 + p23 * x39 + p30 * x40 + p37 * x41 + p44 * x42; 802 pc[37] = m38 = p3 * x36 + p10 * x37 + p17 * x38 + p24 * x39 + p31 * x40 + p38 * x41 + p45 * x42; 803 pc[38] = m39 = p4 * x36 + p11 * x37 + p18 * x38 + p25 * x39 + p32 * x40 + p39 * x41 + p46 * x42; 804 pc[39] = m40 = p5 * x36 + p12 * x37 + p19 * x38 + p26 * x39 + p33 * x40 + p40 * x41 + p47 * x42; 805 pc[40] = m41 = p6 * x36 + p13 * x37 + p20 * x38 + p27 * x39 + p34 * x40 + p41 * x41 + p48 * x42; 806 pc[41] = m42 = p7 * x36 + p14 * x37 + p21 * x38 + p28 * x39 + p35 * x40 + p42 * x41 + p49 * x42; 807 808 pc[42] = m43 = p1 * x43 + p8 * x44 + p15 * x45 + p22 * x46 + p29 * x47 + p36 * x48 + p43 * x49; 809 pc[43] = m44 = p2 * x43 + p9 * x44 + p16 * x45 + p23 * x46 + p30 * x47 + p37 * x48 + p44 * x49; 810 pc[44] = m45 = p3 * x43 + p10 * x44 + p17 * x45 + p24 * x46 + p31 * x47 + p38 * x48 + p45 * x49; 811 pc[45] = m46 = p4 * x43 + p11 * x44 + p18 * x45 + p25 * x46 + p32 * x47 + p39 * x48 + p46 * x49; 812 pc[46] = m47 = p5 * x43 + p12 * x44 + p19 * x45 + p26 * x46 + p33 * x47 + p40 * x48 + p47 * x49; 813 pc[47] = m48 = p6 * x43 + p13 * x44 + p20 * x45 + p27 * x46 + p34 * x47 + p41 * x48 + p48 * x49; 814 pc[48] = m49 = p7 * x43 + p14 * x44 + p21 * x45 + p28 * x46 + p35 * x47 + p42 * x48 + p49 * x49; 815 816 nz = bi[row + 1] - diag_offset[row] - 1; 817 pv += 49; 818 for (j = 0; j < nz; j++) { 819 x1 = pv[0]; 820 x2 = pv[1]; 821 x3 = pv[2]; 822 x4 = pv[3]; 823 x5 = pv[4]; 824 x6 = pv[5]; 825 x7 = pv[6]; 826 x8 = pv[7]; 827 x9 = pv[8]; 828 x10 = pv[9]; 829 x11 = pv[10]; 830 x12 = pv[11]; 831 x13 = pv[12]; 832 x14 = pv[13]; 833 x15 = pv[14]; 834 x16 = pv[15]; 835 x17 = pv[16]; 836 x18 = pv[17]; 837 x19 = pv[18]; 838 x20 = pv[19]; 839 x21 = pv[20]; 840 x22 = pv[21]; 841 x23 = pv[22]; 842 x24 = pv[23]; 843 x25 = pv[24]; 844 x26 = pv[25]; 845 x27 = pv[26]; 846 x28 = pv[27]; 847 x29 = pv[28]; 848 x30 = pv[29]; 849 x31 = pv[30]; 850 x32 = pv[31]; 851 x33 = pv[32]; 852 x34 = pv[33]; 853 x35 = pv[34]; 854 x36 = pv[35]; 855 x37 = pv[36]; 856 x38 = pv[37]; 857 x39 = pv[38]; 858 x40 = pv[39]; 859 x41 = pv[40]; 860 x42 = pv[41]; 861 x43 = pv[42]; 862 x44 = pv[43]; 863 x45 = pv[44]; 864 x46 = pv[45]; 865 x47 = pv[46]; 866 x48 = pv[47]; 867 x49 = pv[48]; 868 x = rtmp + 49 * pj[j]; 869 x[0] -= m1 * x1 + m8 * x2 + m15 * x3 + m22 * x4 + m29 * x5 + m36 * x6 + m43 * x7; 870 x[1] -= m2 * x1 + m9 * x2 + m16 * x3 + m23 * x4 + m30 * x5 + m37 * x6 + m44 * x7; 871 x[2] -= m3 * x1 + m10 * x2 + m17 * x3 + m24 * x4 + m31 * x5 + m38 * x6 + m45 * x7; 872 x[3] -= m4 * x1 + m11 * x2 + m18 * x3 + m25 * x4 + m32 * x5 + m39 * x6 + m46 * x7; 873 x[4] -= m5 * x1 + m12 * x2 + m19 * x3 + m26 * x4 + m33 * x5 + m40 * x6 + m47 * x7; 874 x[5] -= m6 * x1 + m13 * x2 + m20 * x3 + m27 * x4 + m34 * x5 + m41 * x6 + m48 * x7; 875 x[6] -= m7 * x1 + m14 * x2 + m21 * x3 + m28 * x4 + m35 * x5 + m42 * x6 + m49 * x7; 876 877 x[7] -= m1 * x8 + m8 * x9 + m15 * x10 + m22 * x11 + m29 * x12 + m36 * x13 + m43 * x14; 878 x[8] -= m2 * x8 + m9 * x9 + m16 * x10 + m23 * x11 + m30 * x12 + m37 * x13 + m44 * x14; 879 x[9] -= m3 * x8 + m10 * x9 + m17 * x10 + m24 * x11 + m31 * x12 + m38 * x13 + m45 * x14; 880 x[10] -= m4 * x8 + m11 * x9 + m18 * x10 + m25 * x11 + m32 * x12 + m39 * x13 + m46 * x14; 881 x[11] -= m5 * x8 + m12 * x9 + m19 * x10 + m26 * x11 + m33 * x12 + m40 * x13 + m47 * x14; 882 x[12] -= m6 * x8 + m13 * x9 + m20 * x10 + m27 * x11 + m34 * x12 + m41 * x13 + m48 * x14; 883 x[13] -= m7 * x8 + m14 * x9 + m21 * x10 + m28 * x11 + m35 * x12 + m42 * x13 + m49 * x14; 884 885 x[14] -= m1 * x15 + m8 * x16 + m15 * x17 + m22 * x18 + m29 * x19 + m36 * x20 + m43 * x21; 886 x[15] -= m2 * x15 + m9 * x16 + m16 * x17 + m23 * x18 + m30 * x19 + m37 * x20 + m44 * x21; 887 x[16] -= m3 * x15 + m10 * x16 + m17 * x17 + m24 * x18 + m31 * x19 + m38 * x20 + m45 * x21; 888 x[17] -= m4 * x15 + m11 * x16 + m18 * x17 + m25 * x18 + m32 * x19 + m39 * x20 + m46 * x21; 889 x[18] -= m5 * x15 + m12 * x16 + m19 * x17 + m26 * x18 + m33 * x19 + m40 * x20 + m47 * x21; 890 x[19] -= m6 * x15 + m13 * x16 + m20 * x17 + m27 * x18 + m34 * x19 + m41 * x20 + m48 * x21; 891 x[20] -= m7 * x15 + m14 * x16 + m21 * x17 + m28 * x18 + m35 * x19 + m42 * x20 + m49 * x21; 892 893 x[21] -= m1 * x22 + m8 * x23 + m15 * x24 + m22 * x25 + m29 * x26 + m36 * x27 + m43 * x28; 894 x[22] -= m2 * x22 + m9 * x23 + m16 * x24 + m23 * x25 + m30 * x26 + m37 * x27 + m44 * x28; 895 x[23] -= m3 * x22 + m10 * x23 + m17 * x24 + m24 * x25 + m31 * x26 + m38 * x27 + m45 * x28; 896 x[24] -= m4 * x22 + m11 * x23 + m18 * x24 + m25 * x25 + m32 * x26 + m39 * x27 + m46 * x28; 897 x[25] -= m5 * x22 + m12 * x23 + m19 * x24 + m26 * x25 + m33 * x26 + m40 * x27 + m47 * x28; 898 x[26] -= m6 * x22 + m13 * x23 + m20 * x24 + m27 * x25 + m34 * x26 + m41 * x27 + m48 * x28; 899 x[27] -= m7 * x22 + m14 * x23 + m21 * x24 + m28 * x25 + m35 * x26 + m42 * x27 + m49 * x28; 900 901 x[28] -= m1 * x29 + m8 * x30 + m15 * x31 + m22 * x32 + m29 * x33 + m36 * x34 + m43 * x35; 902 x[29] -= m2 * x29 + m9 * x30 + m16 * x31 + m23 * x32 + m30 * x33 + m37 * x34 + m44 * x35; 903 x[30] -= m3 * x29 + m10 * x30 + m17 * x31 + m24 * x32 + m31 * x33 + m38 * x34 + m45 * x35; 904 x[31] -= m4 * x29 + m11 * x30 + m18 * x31 + m25 * x32 + m32 * x33 + m39 * x34 + m46 * x35; 905 x[32] -= m5 * x29 + m12 * x30 + m19 * x31 + m26 * x32 + m33 * x33 + m40 * x34 + m47 * x35; 906 x[33] -= m6 * x29 + m13 * x30 + m20 * x31 + m27 * x32 + m34 * x33 + m41 * x34 + m48 * x35; 907 x[34] -= m7 * x29 + m14 * x30 + m21 * x31 + m28 * x32 + m35 * x33 + m42 * x34 + m49 * x35; 908 909 x[35] -= m1 * x36 + m8 * x37 + m15 * x38 + m22 * x39 + m29 * x40 + m36 * x41 + m43 * x42; 910 x[36] -= m2 * x36 + m9 * x37 + m16 * x38 + m23 * x39 + m30 * x40 + m37 * x41 + m44 * x42; 911 x[37] -= m3 * x36 + m10 * x37 + m17 * x38 + m24 * x39 + m31 * x40 + m38 * x41 + m45 * x42; 912 x[38] -= m4 * x36 + m11 * x37 + m18 * x38 + m25 * x39 + m32 * x40 + m39 * x41 + m46 * x42; 913 x[39] -= m5 * x36 + m12 * x37 + m19 * x38 + m26 * x39 + m33 * x40 + m40 * x41 + m47 * x42; 914 x[40] -= m6 * x36 + m13 * x37 + m20 * x38 + m27 * x39 + m34 * x40 + m41 * x41 + m48 * x42; 915 x[41] -= m7 * x36 + m14 * x37 + m21 * x38 + m28 * x39 + m35 * x40 + m42 * x41 + m49 * x42; 916 917 x[42] -= m1 * x43 + m8 * x44 + m15 * x45 + m22 * x46 + m29 * x47 + m36 * x48 + m43 * x49; 918 x[43] -= m2 * x43 + m9 * x44 + m16 * x45 + m23 * x46 + m30 * x47 + m37 * x48 + m44 * x49; 919 x[44] -= m3 * x43 + m10 * x44 + m17 * x45 + m24 * x46 + m31 * x47 + m38 * x48 + m45 * x49; 920 x[45] -= m4 * x43 + m11 * x44 + m18 * x45 + m25 * x46 + m32 * x47 + m39 * x48 + m46 * x49; 921 x[46] -= m5 * x43 + m12 * x44 + m19 * x45 + m26 * x46 + m33 * x47 + m40 * x48 + m47 * x49; 922 x[47] -= m6 * x43 + m13 * x44 + m20 * x45 + m27 * x46 + m34 * x47 + m41 * x48 + m48 * x49; 923 x[48] -= m7 * x43 + m14 * x44 + m21 * x45 + m28 * x46 + m35 * x47 + m42 * x48 + m49 * x49; 924 pv += 49; 925 } 926 PetscCall(PetscLogFlops(686.0 * nz + 637.0)); 927 } 928 row = *ajtmp++; 929 } 930 /* finished row so stick it into b->a */ 931 pv = ba + 49 * bi[i]; 932 pj = bj + bi[i]; 933 nz = bi[i + 1] - bi[i]; 934 for (j = 0; j < nz; j++) { 935 x = rtmp + 49 * pj[j]; 936 pv[0] = x[0]; 937 pv[1] = x[1]; 938 pv[2] = x[2]; 939 pv[3] = x[3]; 940 pv[4] = x[4]; 941 pv[5] = x[5]; 942 pv[6] = x[6]; 943 pv[7] = x[7]; 944 pv[8] = x[8]; 945 pv[9] = x[9]; 946 pv[10] = x[10]; 947 pv[11] = x[11]; 948 pv[12] = x[12]; 949 pv[13] = x[13]; 950 pv[14] = x[14]; 951 pv[15] = x[15]; 952 pv[16] = x[16]; 953 pv[17] = x[17]; 954 pv[18] = x[18]; 955 pv[19] = x[19]; 956 pv[20] = x[20]; 957 pv[21] = x[21]; 958 pv[22] = x[22]; 959 pv[23] = x[23]; 960 pv[24] = x[24]; 961 pv[25] = x[25]; 962 pv[26] = x[26]; 963 pv[27] = x[27]; 964 pv[28] = x[28]; 965 pv[29] = x[29]; 966 pv[30] = x[30]; 967 pv[31] = x[31]; 968 pv[32] = x[32]; 969 pv[33] = x[33]; 970 pv[34] = x[34]; 971 pv[35] = x[35]; 972 pv[36] = x[36]; 973 pv[37] = x[37]; 974 pv[38] = x[38]; 975 pv[39] = x[39]; 976 pv[40] = x[40]; 977 pv[41] = x[41]; 978 pv[42] = x[42]; 979 pv[43] = x[43]; 980 pv[44] = x[44]; 981 pv[45] = x[45]; 982 pv[46] = x[46]; 983 pv[47] = x[47]; 984 pv[48] = x[48]; 985 pv += 49; 986 } 987 /* invert diagonal block */ 988 w = ba + 49 * diag_offset[i]; 989 PetscCall(PetscKernel_A_gets_inverse_A_7(w, shift, allowzeropivot, &zeropivotdetected)); 990 if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 991 } 992 993 PetscCall(PetscFree(rtmp)); 994 995 C->ops->solve = MatSolve_SeqBAIJ_7_NaturalOrdering_inplace; 996 C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_7_NaturalOrdering_inplace; 997 C->assembled = PETSC_TRUE; 998 999 PetscCall(PetscLogFlops(1.333333333333 * 7 * 7 * 7 * b->mbs)); /* from inverting diagonal blocks */ 1000 PetscFunctionReturn(PETSC_SUCCESS); 1001 } 1002 1003 PetscErrorCode MatLUFactorNumeric_SeqBAIJ_7_NaturalOrdering(Mat B, Mat A, const MatFactorInfo *info) 1004 { 1005 Mat C = B; 1006 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data, *b = (Mat_SeqBAIJ *)C->data; 1007 PetscInt i, j, k, nz, nzL, row; 1008 const PetscInt n = a->mbs, *ai = a->i, *aj = a->j, *bi = b->i, *bj = b->j; 1009 const PetscInt *ajtmp, *bjtmp, *bdiag = b->diag, *pj, bs2 = a->bs2; 1010 MatScalar *rtmp, *pc, *mwork, *v, *pv, *aa = a->a; 1011 PetscInt flg; 1012 PetscReal shift = info->shiftamount; 1013 PetscBool allowzeropivot, zeropivotdetected; 1014 1015 PetscFunctionBegin; 1016 allowzeropivot = PetscNot(A->erroriffailure); 1017 1018 /* generate work space needed by the factorization */ 1019 PetscCall(PetscMalloc2(bs2 * n, &rtmp, bs2, &mwork)); 1020 PetscCall(PetscArrayzero(rtmp, bs2 * n)); 1021 1022 for (i = 0; i < n; i++) { 1023 /* zero rtmp */ 1024 /* L part */ 1025 nz = bi[i + 1] - bi[i]; 1026 bjtmp = bj + bi[i]; 1027 for (j = 0; j < nz; j++) PetscCall(PetscArrayzero(rtmp + bs2 * bjtmp[j], bs2)); 1028 1029 /* U part */ 1030 nz = bdiag[i] - bdiag[i + 1]; 1031 bjtmp = bj + bdiag[i + 1] + 1; 1032 for (j = 0; j < nz; j++) PetscCall(PetscArrayzero(rtmp + bs2 * bjtmp[j], bs2)); 1033 1034 /* load in initial (unfactored row) */ 1035 nz = ai[i + 1] - ai[i]; 1036 ajtmp = aj + ai[i]; 1037 v = aa + bs2 * ai[i]; 1038 for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(rtmp + bs2 * ajtmp[j], v + bs2 * j, bs2)); 1039 1040 /* elimination */ 1041 bjtmp = bj + bi[i]; 1042 nzL = bi[i + 1] - bi[i]; 1043 for (k = 0; k < nzL; k++) { 1044 row = bjtmp[k]; 1045 pc = rtmp + bs2 * row; 1046 for (flg = 0, j = 0; j < bs2; j++) { 1047 if (pc[j] != 0.0) { 1048 flg = 1; 1049 break; 1050 } 1051 } 1052 if (flg) { 1053 pv = b->a + bs2 * bdiag[row]; 1054 /* PetscKernel_A_gets_A_times_B(bs,pc,pv,mwork); *pc = *pc * (*pv); */ 1055 PetscCall(PetscKernel_A_gets_A_times_B_7(pc, pv, mwork)); 1056 1057 pj = b->j + bdiag[row + 1] + 1; /* beginning of U(row,:) */ 1058 pv = b->a + bs2 * (bdiag[row + 1] + 1); 1059 nz = bdiag[row] - bdiag[row + 1] - 1; /* num of entries inU(row,:), excluding diag */ 1060 for (j = 0; j < nz; j++) { 1061 /* PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j); */ 1062 /* rtmp+bs2*pj[j] = rtmp+bs2*pj[j] - (*pc)*(pv+bs2*j) */ 1063 v = rtmp + bs2 * pj[j]; 1064 PetscCall(PetscKernel_A_gets_A_minus_B_times_C_7(v, pc, pv)); 1065 pv += bs2; 1066 } 1067 PetscCall(PetscLogFlops(686.0 * nz + 637)); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */ 1068 } 1069 } 1070 1071 /* finished row so stick it into b->a */ 1072 /* L part */ 1073 pv = b->a + bs2 * bi[i]; 1074 pj = b->j + bi[i]; 1075 nz = bi[i + 1] - bi[i]; 1076 for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2)); 1077 1078 /* Mark diagonal and invert diagonal for simpler triangular solves */ 1079 pv = b->a + bs2 * bdiag[i]; 1080 pj = b->j + bdiag[i]; 1081 PetscCall(PetscArraycpy(pv, rtmp + bs2 * pj[0], bs2)); 1082 PetscCall(PetscKernel_A_gets_inverse_A_7(pv, shift, allowzeropivot, &zeropivotdetected)); 1083 if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1084 1085 /* U part */ 1086 pv = b->a + bs2 * (bdiag[i + 1] + 1); 1087 pj = b->j + bdiag[i + 1] + 1; 1088 nz = bdiag[i] - bdiag[i + 1] - 1; 1089 for (j = 0; j < nz; j++) PetscCall(PetscArraycpy(pv + bs2 * j, rtmp + bs2 * pj[j], bs2)); 1090 } 1091 PetscCall(PetscFree2(rtmp, mwork)); 1092 1093 C->ops->solve = MatSolve_SeqBAIJ_7_NaturalOrdering; 1094 C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_7_NaturalOrdering; 1095 C->assembled = PETSC_TRUE; 1096 1097 PetscCall(PetscLogFlops(1.333333333333 * 7 * 7 * 7 * n)); /* from inverting diagonal blocks */ 1098 PetscFunctionReturn(PETSC_SUCCESS); 1099 } 1100