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