1 #include <../src/mat/impls/baij/seq/baij.h> 2 #include <petsc/private/kernels/blockinvert.h> 3 4 PetscErrorCode MatSolveTranspose_SeqBAIJ_7_inplace(Mat A, Vec bb, Vec xx) 5 { 6 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data; 7 IS iscol = a->col, isrow = a->row; 8 const PetscInt *r, *c, *rout, *cout; 9 const PetscInt *diag = a->diag, n = a->mbs, *vi, *ai = a->i, *aj = a->j; 10 PetscInt i, nz, idx, idt, ii, ic, ir, oidx; 11 const MatScalar *aa = a->a, *v; 12 PetscScalar s1, s2, s3, s4, s5, s6, s7, x1, x2, x3, x4, x5, x6, x7, *x, *t; 13 const PetscScalar *b; 14 15 PetscFunctionBegin; 16 PetscCall(VecGetArrayRead(bb, &b)); 17 PetscCall(VecGetArray(xx, &x)); 18 t = a->solve_work; 19 20 PetscCall(ISGetIndices(isrow, &rout)); 21 r = rout; 22 PetscCall(ISGetIndices(iscol, &cout)); 23 c = cout; 24 25 /* copy the b into temp work space according to permutation */ 26 ii = 0; 27 for (i = 0; i < n; i++) { 28 ic = 7 * c[i]; 29 t[ii] = b[ic]; 30 t[ii + 1] = b[ic + 1]; 31 t[ii + 2] = b[ic + 2]; 32 t[ii + 3] = b[ic + 3]; 33 t[ii + 4] = b[ic + 4]; 34 t[ii + 5] = b[ic + 5]; 35 t[ii + 6] = b[ic + 6]; 36 ii += 7; 37 } 38 39 /* forward solve the U^T */ 40 idx = 0; 41 for (i = 0; i < n; i++) { 42 v = aa + 49 * diag[i]; 43 /* multiply by the inverse of the block diagonal */ 44 x1 = t[idx]; 45 x2 = t[1 + idx]; 46 x3 = t[2 + idx]; 47 x4 = t[3 + idx]; 48 x5 = t[4 + idx]; 49 x6 = t[5 + idx]; 50 x7 = t[6 + idx]; 51 s1 = v[0] * x1 + v[1] * x2 + v[2] * x3 + v[3] * x4 + v[4] * x5 + v[5] * x6 + v[6] * x7; 52 s2 = v[7] * x1 + v[8] * x2 + v[9] * x3 + v[10] * x4 + v[11] * x5 + v[12] * x6 + v[13] * x7; 53 s3 = v[14] * x1 + v[15] * x2 + v[16] * x3 + v[17] * x4 + v[18] * x5 + v[19] * x6 + v[20] * x7; 54 s4 = v[21] * x1 + v[22] * x2 + v[23] * x3 + v[24] * x4 + v[25] * x5 + v[26] * x6 + v[27] * x7; 55 s5 = v[28] * x1 + v[29] * x2 + v[30] * x3 + v[31] * x4 + v[32] * x5 + v[33] * x6 + v[34] * x7; 56 s6 = v[35] * x1 + v[36] * x2 + v[37] * x3 + v[38] * x4 + v[39] * x5 + v[40] * x6 + v[41] * x7; 57 s7 = v[42] * x1 + v[43] * x2 + v[44] * x3 + v[45] * x4 + v[46] * x5 + v[47] * x6 + v[48] * x7; 58 v += 49; 59 60 vi = aj + diag[i] + 1; 61 nz = ai[i + 1] - diag[i] - 1; 62 while (nz--) { 63 oidx = 7 * (*vi++); 64 t[oidx] -= v[0] * s1 + v[1] * s2 + v[2] * s3 + v[3] * s4 + v[4] * s5 + v[5] * s6 + v[6] * s7; 65 t[oidx + 1] -= v[7] * s1 + v[8] * s2 + v[9] * s3 + v[10] * s4 + v[11] * s5 + v[12] * s6 + v[13] * s7; 66 t[oidx + 2] -= v[14] * s1 + v[15] * s2 + v[16] * s3 + v[17] * s4 + v[18] * s5 + v[19] * s6 + v[20] * s7; 67 t[oidx + 3] -= v[21] * s1 + v[22] * s2 + v[23] * s3 + v[24] * s4 + v[25] * s5 + v[26] * s6 + v[27] * s7; 68 t[oidx + 4] -= v[28] * s1 + v[29] * s2 + v[30] * s3 + v[31] * s4 + v[32] * s5 + v[33] * s6 + v[34] * s7; 69 t[oidx + 5] -= v[35] * s1 + v[36] * s2 + v[37] * s3 + v[38] * s4 + v[39] * s5 + v[40] * s6 + v[41] * s7; 70 t[oidx + 6] -= v[42] * s1 + v[43] * s2 + v[44] * s3 + v[45] * s4 + v[46] * s5 + v[47] * s6 + v[48] * s7; 71 v += 49; 72 } 73 t[idx] = s1; 74 t[1 + idx] = s2; 75 t[2 + idx] = s3; 76 t[3 + idx] = s4; 77 t[4 + idx] = s5; 78 t[5 + idx] = s6; 79 t[6 + idx] = s7; 80 idx += 7; 81 } 82 /* backward solve the L^T */ 83 for (i = n - 1; i >= 0; i--) { 84 v = aa + 49 * diag[i] - 49; 85 vi = aj + diag[i] - 1; 86 nz = diag[i] - ai[i]; 87 idt = 7 * i; 88 s1 = t[idt]; 89 s2 = t[1 + idt]; 90 s3 = t[2 + idt]; 91 s4 = t[3 + idt]; 92 s5 = t[4 + idt]; 93 s6 = t[5 + idt]; 94 s7 = t[6 + idt]; 95 while (nz--) { 96 idx = 7 * (*vi--); 97 t[idx] -= v[0] * s1 + v[1] * s2 + v[2] * s3 + v[3] * s4 + v[4] * s5 + v[5] * s6 + v[6] * s7; 98 t[idx + 1] -= v[7] * s1 + v[8] * s2 + v[9] * s3 + v[10] * s4 + v[11] * s5 + v[12] * s6 + v[13] * s7; 99 t[idx + 2] -= v[14] * s1 + v[15] * s2 + v[16] * s3 + v[17] * s4 + v[18] * s5 + v[19] * s6 + v[20] * s7; 100 t[idx + 3] -= v[21] * s1 + v[22] * s2 + v[23] * s3 + v[24] * s4 + v[25] * s5 + v[26] * s6 + v[27] * s7; 101 t[idx + 4] -= v[28] * s1 + v[29] * s2 + v[30] * s3 + v[31] * s4 + v[32] * s5 + v[33] * s6 + v[34] * s7; 102 t[idx + 5] -= v[35] * s1 + v[36] * s2 + v[37] * s3 + v[38] * s4 + v[39] * s5 + v[40] * s6 + v[41] * s7; 103 t[idx + 6] -= v[42] * s1 + v[43] * s2 + v[44] * s3 + v[45] * s4 + v[46] * s5 + v[47] * s6 + v[48] * s7; 104 v -= 49; 105 } 106 } 107 108 /* copy t into x according to permutation */ 109 ii = 0; 110 for (i = 0; i < n; i++) { 111 ir = 7 * r[i]; 112 x[ir] = t[ii]; 113 x[ir + 1] = t[ii + 1]; 114 x[ir + 2] = t[ii + 2]; 115 x[ir + 3] = t[ii + 3]; 116 x[ir + 4] = t[ii + 4]; 117 x[ir + 5] = t[ii + 5]; 118 x[ir + 6] = t[ii + 6]; 119 ii += 7; 120 } 121 122 PetscCall(ISRestoreIndices(isrow, &rout)); 123 PetscCall(ISRestoreIndices(iscol, &cout)); 124 PetscCall(VecRestoreArrayRead(bb, &b)); 125 PetscCall(VecRestoreArray(xx, &x)); 126 PetscCall(PetscLogFlops(2.0 * 49 * (a->nz) - 7.0 * A->cmap->n)); 127 PetscFunctionReturn(PETSC_SUCCESS); 128 } 129 PetscErrorCode MatSolveTranspose_SeqBAIJ_7(Mat A, Vec bb, Vec xx) 130 { 131 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data; 132 IS iscol = a->col, isrow = a->row; 133 const PetscInt n = a->mbs, *vi, *ai = a->i, *aj = a->j, *diag = a->diag; 134 const PetscInt *r, *c, *rout, *cout; 135 PetscInt nz, idx, idt, j, i, oidx, ii, ic, ir; 136 const PetscInt bs = A->rmap->bs, bs2 = a->bs2; 137 const MatScalar *aa = a->a, *v; 138 PetscScalar s1, s2, s3, s4, s5, s6, s7, x1, x2, x3, x4, x5, x6, x7, *x, *t; 139 const PetscScalar *b; 140 141 PetscFunctionBegin; 142 PetscCall(VecGetArrayRead(bb, &b)); 143 PetscCall(VecGetArray(xx, &x)); 144 t = a->solve_work; 145 146 PetscCall(ISGetIndices(isrow, &rout)); 147 r = rout; 148 PetscCall(ISGetIndices(iscol, &cout)); 149 c = cout; 150 151 /* copy b into temp work space according to permutation */ 152 for (i = 0; i < n; i++) { 153 ii = bs * i; 154 ic = bs * c[i]; 155 t[ii] = b[ic]; 156 t[ii + 1] = b[ic + 1]; 157 t[ii + 2] = b[ic + 2]; 158 t[ii + 3] = b[ic + 3]; 159 t[ii + 4] = b[ic + 4]; 160 t[ii + 5] = b[ic + 5]; 161 t[ii + 6] = b[ic + 6]; 162 } 163 164 /* forward solve the U^T */ 165 idx = 0; 166 for (i = 0; i < n; i++) { 167 v = aa + bs2 * diag[i]; 168 /* multiply by the inverse of the block diagonal */ 169 x1 = t[idx]; 170 x2 = t[1 + idx]; 171 x3 = t[2 + idx]; 172 x4 = t[3 + idx]; 173 x5 = t[4 + idx]; 174 x6 = t[5 + idx]; 175 x7 = t[6 + idx]; 176 s1 = v[0] * x1 + v[1] * x2 + v[2] * x3 + v[3] * x4 + v[4] * x5 + v[5] * x6 + v[6] * x7; 177 s2 = v[7] * x1 + v[8] * x2 + v[9] * x3 + v[10] * x4 + v[11] * x5 + v[12] * x6 + v[13] * x7; 178 s3 = v[14] * x1 + v[15] * x2 + v[16] * x3 + v[17] * x4 + v[18] * x5 + v[19] * x6 + v[20] * x7; 179 s4 = v[21] * x1 + v[22] * x2 + v[23] * x3 + v[24] * x4 + v[25] * x5 + v[26] * x6 + v[27] * x7; 180 s5 = v[28] * x1 + v[29] * x2 + v[30] * x3 + v[31] * x4 + v[32] * x5 + v[33] * x6 + v[34] * x7; 181 s6 = v[35] * x1 + v[36] * x2 + v[37] * x3 + v[38] * x4 + v[39] * x5 + v[40] * x6 + v[41] * x7; 182 s7 = v[42] * x1 + v[43] * x2 + v[44] * x3 + v[45] * x4 + v[46] * x5 + v[47] * x6 + v[48] * x7; 183 v -= bs2; 184 185 vi = aj + diag[i] - 1; 186 nz = diag[i] - diag[i + 1] - 1; 187 for (j = 0; j > -nz; j--) { 188 oidx = bs * vi[j]; 189 t[oidx] -= v[0] * s1 + v[1] * s2 + v[2] * s3 + v[3] * s4 + v[4] * s5 + v[5] * s6 + v[6] * s7; 190 t[oidx + 1] -= v[7] * s1 + v[8] * s2 + v[9] * s3 + v[10] * s4 + v[11] * s5 + v[12] * s6 + v[13] * s7; 191 t[oidx + 2] -= v[14] * s1 + v[15] * s2 + v[16] * s3 + v[17] * s4 + v[18] * s5 + v[19] * s6 + v[20] * s7; 192 t[oidx + 3] -= v[21] * s1 + v[22] * s2 + v[23] * s3 + v[24] * s4 + v[25] * s5 + v[26] * s6 + v[27] * s7; 193 t[oidx + 4] -= v[28] * s1 + v[29] * s2 + v[30] * s3 + v[31] * s4 + v[32] * s5 + v[33] * s6 + v[34] * s7; 194 t[oidx + 5] -= v[35] * s1 + v[36] * s2 + v[37] * s3 + v[38] * s4 + v[39] * s5 + v[40] * s6 + v[41] * s7; 195 t[oidx + 6] -= v[42] * s1 + v[43] * s2 + v[44] * s3 + v[45] * s4 + v[46] * s5 + v[47] * s6 + v[48] * s7; 196 v -= bs2; 197 } 198 t[idx] = s1; 199 t[1 + idx] = s2; 200 t[2 + idx] = s3; 201 t[3 + idx] = s4; 202 t[4 + idx] = s5; 203 t[5 + idx] = s6; 204 t[6 + idx] = s7; 205 idx += bs; 206 } 207 /* backward solve the L^T */ 208 for (i = n - 1; i >= 0; i--) { 209 v = aa + bs2 * ai[i]; 210 vi = aj + ai[i]; 211 nz = ai[i + 1] - ai[i]; 212 idt = bs * i; 213 s1 = t[idt]; 214 s2 = t[1 + idt]; 215 s3 = t[2 + idt]; 216 s4 = t[3 + idt]; 217 s5 = t[4 + idt]; 218 s6 = t[5 + idt]; 219 s7 = t[6 + idt]; 220 for (j = 0; j < nz; j++) { 221 idx = bs * vi[j]; 222 t[idx] -= v[0] * s1 + v[1] * s2 + v[2] * s3 + v[3] * s4 + v[4] * s5 + v[5] * s6 + v[6] * s7; 223 t[idx + 1] -= v[7] * s1 + v[8] * s2 + v[9] * s3 + v[10] * s4 + v[11] * s5 + v[12] * s6 + v[13] * s7; 224 t[idx + 2] -= v[14] * s1 + v[15] * s2 + v[16] * s3 + v[17] * s4 + v[18] * s5 + v[19] * s6 + v[20] * s7; 225 t[idx + 3] -= v[21] * s1 + v[22] * s2 + v[23] * s3 + v[24] * s4 + v[25] * s5 + v[26] * s6 + v[27] * s7; 226 t[idx + 4] -= v[28] * s1 + v[29] * s2 + v[30] * s3 + v[31] * s4 + v[32] * s5 + v[33] * s6 + v[34] * s7; 227 t[idx + 5] -= v[35] * s1 + v[36] * s2 + v[37] * s3 + v[38] * s4 + v[39] * s5 + v[40] * s6 + v[41] * s7; 228 t[idx + 6] -= v[42] * s1 + v[43] * s2 + v[44] * s3 + v[45] * s4 + v[46] * s5 + v[47] * s6 + v[48] * s7; 229 v += bs2; 230 } 231 } 232 233 /* copy t into x according to permutation */ 234 for (i = 0; i < n; i++) { 235 ii = bs * i; 236 ir = bs * r[i]; 237 x[ir] = t[ii]; 238 x[ir + 1] = t[ii + 1]; 239 x[ir + 2] = t[ii + 2]; 240 x[ir + 3] = t[ii + 3]; 241 x[ir + 4] = t[ii + 4]; 242 x[ir + 5] = t[ii + 5]; 243 x[ir + 6] = t[ii + 6]; 244 } 245 246 PetscCall(ISRestoreIndices(isrow, &rout)); 247 PetscCall(ISRestoreIndices(iscol, &cout)); 248 PetscCall(VecRestoreArrayRead(bb, &b)); 249 PetscCall(VecRestoreArray(xx, &x)); 250 PetscCall(PetscLogFlops(2.0 * bs2 * (a->nz) - bs * A->cmap->n)); 251 PetscFunctionReturn(PETSC_SUCCESS); 252 } 253