1 #include <../src/mat/impls/baij/seq/baij.h> 2 3 PetscErrorCode MatSolveTranspose_SeqBAIJ_4_NaturalOrdering_inplace(Mat A, Vec bb, Vec xx) { 4 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data; 5 const PetscInt *diag = a->diag, n = a->mbs, *vi, *ai = a->i, *aj = a->j; 6 PetscInt i, nz, idx, idt, oidx; 7 const MatScalar *aa = a->a, *v; 8 PetscScalar s1, s2, s3, s4, x1, x2, x3, x4, *x; 9 10 PetscFunctionBegin; 11 PetscCall(VecCopy(bb, xx)); 12 PetscCall(VecGetArray(xx, &x)); 13 14 /* forward solve the U^T */ 15 idx = 0; 16 for (i = 0; i < n; i++) { 17 v = aa + 16 * diag[i]; 18 /* multiply by the inverse of the block diagonal */ 19 x1 = x[idx]; 20 x2 = x[1 + idx]; 21 x3 = x[2 + idx]; 22 x4 = x[3 + idx]; 23 s1 = v[0] * x1 + v[1] * x2 + v[2] * x3 + v[3] * x4; 24 s2 = v[4] * x1 + v[5] * x2 + v[6] * x3 + v[7] * x4; 25 s3 = v[8] * x1 + v[9] * x2 + v[10] * x3 + v[11] * x4; 26 s4 = v[12] * x1 + v[13] * x2 + v[14] * x3 + v[15] * x4; 27 v += 16; 28 29 vi = aj + diag[i] + 1; 30 nz = ai[i + 1] - diag[i] - 1; 31 while (nz--) { 32 oidx = 4 * (*vi++); 33 x[oidx] -= v[0] * s1 + v[1] * s2 + v[2] * s3 + v[3] * s4; 34 x[oidx + 1] -= v[4] * s1 + v[5] * s2 + v[6] * s3 + v[7] * s4; 35 x[oidx + 2] -= v[8] * s1 + v[9] * s2 + v[10] * s3 + v[11] * s4; 36 x[oidx + 3] -= v[12] * s1 + v[13] * s2 + v[14] * s3 + v[15] * s4; 37 v += 16; 38 } 39 x[idx] = s1; 40 x[1 + idx] = s2; 41 x[2 + idx] = s3; 42 x[3 + idx] = s4; 43 idx += 4; 44 } 45 /* backward solve the L^T */ 46 for (i = n - 1; i >= 0; i--) { 47 v = aa + 16 * diag[i] - 16; 48 vi = aj + diag[i] - 1; 49 nz = diag[i] - ai[i]; 50 idt = 4 * i; 51 s1 = x[idt]; 52 s2 = x[1 + idt]; 53 s3 = x[2 + idt]; 54 s4 = x[3 + idt]; 55 while (nz--) { 56 idx = 4 * (*vi--); 57 x[idx] -= v[0] * s1 + v[1] * s2 + v[2] * s3 + v[3] * s4; 58 x[idx + 1] -= v[4] * s1 + v[5] * s2 + v[6] * s3 + v[7] * s4; 59 x[idx + 2] -= v[8] * s1 + v[9] * s2 + v[10] * s3 + v[11] * s4; 60 x[idx + 3] -= v[12] * s1 + v[13] * s2 + v[14] * s3 + v[15] * s4; 61 v -= 16; 62 } 63 } 64 PetscCall(VecRestoreArray(xx, &x)); 65 PetscCall(PetscLogFlops(2.0 * 16 * (a->nz) - 4.0 * A->cmap->n)); 66 PetscFunctionReturn(0); 67 } 68 69 PetscErrorCode MatSolveTranspose_SeqBAIJ_4_NaturalOrdering(Mat A, Vec bb, Vec xx) { 70 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data; 71 const PetscInt n = a->mbs, *vi, *ai = a->i, *aj = a->j, *diag = a->diag; 72 PetscInt nz, idx, idt, j, i, oidx; 73 const PetscInt bs = A->rmap->bs, bs2 = a->bs2; 74 const MatScalar *aa = a->a, *v; 75 PetscScalar s1, s2, s3, s4, x1, x2, x3, x4, *x; 76 77 PetscFunctionBegin; 78 PetscCall(VecCopy(bb, xx)); 79 PetscCall(VecGetArray(xx, &x)); 80 81 /* forward solve the U^T */ 82 idx = 0; 83 for (i = 0; i < n; i++) { 84 v = aa + bs2 * diag[i]; 85 /* multiply by the inverse of the block diagonal */ 86 x1 = x[idx]; 87 x2 = x[1 + idx]; 88 x3 = x[2 + idx]; 89 x4 = x[3 + idx]; 90 s1 = v[0] * x1 + v[1] * x2 + v[2] * x3 + v[3] * x4; 91 s2 = v[4] * x1 + v[5] * x2 + v[6] * x3 + v[7] * x4; 92 s3 = v[8] * x1 + v[9] * x2 + v[10] * x3 + v[11] * x4; 93 s4 = v[12] * x1 + v[13] * x2 + v[14] * x3 + v[15] * x4; 94 v -= bs2; 95 96 vi = aj + diag[i] - 1; 97 nz = diag[i] - diag[i + 1] - 1; 98 for (j = 0; j > -nz; j--) { 99 oidx = bs * vi[j]; 100 x[oidx] -= v[0] * s1 + v[1] * s2 + v[2] * s3 + v[3] * s4; 101 x[oidx + 1] -= v[4] * s1 + v[5] * s2 + v[6] * s3 + v[7] * s4; 102 x[oidx + 2] -= v[8] * s1 + v[9] * s2 + v[10] * s3 + v[11] * s4; 103 x[oidx + 3] -= v[12] * s1 + v[13] * s2 + v[14] * s3 + v[15] * s4; 104 v -= bs2; 105 } 106 x[idx] = s1; 107 x[1 + idx] = s2; 108 x[2 + idx] = s3; 109 x[3 + idx] = s4; 110 idx += bs; 111 } 112 /* backward solve the L^T */ 113 for (i = n - 1; i >= 0; i--) { 114 v = aa + bs2 * ai[i]; 115 vi = aj + ai[i]; 116 nz = ai[i + 1] - ai[i]; 117 idt = bs * i; 118 s1 = x[idt]; 119 s2 = x[1 + idt]; 120 s3 = x[2 + idt]; 121 s4 = x[3 + idt]; 122 for (j = 0; j < nz; j++) { 123 idx = bs * vi[j]; 124 x[idx] -= v[0] * s1 + v[1] * s2 + v[2] * s3 + v[3] * s4; 125 x[idx + 1] -= v[4] * s1 + v[5] * s2 + v[6] * s3 + v[7] * s4; 126 x[idx + 2] -= v[8] * s1 + v[9] * s2 + v[10] * s3 + v[11] * s4; 127 x[idx + 3] -= v[12] * s1 + v[13] * s2 + v[14] * s3 + v[15] * s4; 128 v += bs2; 129 } 130 } 131 PetscCall(VecRestoreArray(xx, &x)); 132 PetscCall(PetscLogFlops(2.0 * bs2 * (a->nz) - bs * A->cmap->n)); 133 PetscFunctionReturn(0); 134 } 135