1 #include <../src/mat/impls/baij/seq/baij.h> 2 3 PetscErrorCode MatSolveTranspose_SeqBAIJ_1_NaturalOrdering(Mat A, Vec bb, Vec xx) { 4 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data; 5 const PetscInt *adiag = a->diag, *ai = a->i, *aj = a->j, *vi; 6 PetscInt i, n = a->mbs, j; 7 PetscInt nz; 8 PetscScalar *x, *tmp, s1; 9 const MatScalar *aa = a->a, *v; 10 const PetscScalar *b; 11 12 PetscFunctionBegin; 13 PetscCall(VecGetArrayRead(bb, &b)); 14 PetscCall(VecGetArray(xx, &x)); 15 tmp = a->solve_work; 16 17 /* copy the b into temp work space according to permutation */ 18 for (i = 0; i < n; i++) tmp[i] = b[i]; 19 20 /* forward solve the U^T */ 21 for (i = 0; i < n; i++) { 22 v = aa + adiag[i + 1] + 1; 23 vi = aj + adiag[i + 1] + 1; 24 nz = adiag[i] - adiag[i + 1] - 1; 25 s1 = tmp[i]; 26 s1 *= v[nz]; /* multiply by inverse of diagonal entry */ 27 for (j = 0; j < nz; j++) tmp[vi[j]] -= s1 * v[j]; 28 tmp[i] = s1; 29 } 30 31 /* backward solve the L^T */ 32 for (i = n - 1; i >= 0; i--) { 33 v = aa + ai[i]; 34 vi = aj + ai[i]; 35 nz = ai[i + 1] - ai[i]; 36 s1 = tmp[i]; 37 for (j = 0; j < nz; j++) tmp[vi[j]] -= s1 * v[j]; 38 } 39 40 /* copy tmp into x according to permutation */ 41 for (i = 0; i < n; i++) x[i] = tmp[i]; 42 43 PetscCall(VecRestoreArrayRead(bb, &b)); 44 PetscCall(VecRestoreArray(xx, &x)); 45 46 PetscCall(PetscLogFlops(2.0 * a->nz - A->cmap->n)); 47 PetscFunctionReturn(0); 48 } 49 50 PetscErrorCode MatSolveTranspose_SeqBAIJ_1_NaturalOrdering_inplace(Mat A, Vec bb, Vec xx) { 51 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data; 52 PetscInt i, nz; 53 const PetscInt *diag = a->diag, n = a->mbs, *vi, *ai = a->i, *aj = a->j; 54 const MatScalar *aa = a->a, *v; 55 PetscScalar s1, *x; 56 57 PetscFunctionBegin; 58 PetscCall(VecCopy(bb, xx)); 59 PetscCall(VecGetArray(xx, &x)); 60 61 /* forward solve the U^T */ 62 for (i = 0; i < n; i++) { 63 v = aa + diag[i]; 64 /* multiply by the inverse of the block diagonal */ 65 s1 = (*v++) * x[i]; 66 vi = aj + diag[i] + 1; 67 nz = ai[i + 1] - diag[i] - 1; 68 while (nz--) { x[*vi++] -= (*v++) * s1; } 69 x[i] = s1; 70 } 71 /* backward solve the L^T */ 72 for (i = n - 1; i >= 0; i--) { 73 v = aa + diag[i] - 1; 74 vi = aj + diag[i] - 1; 75 nz = diag[i] - ai[i]; 76 s1 = x[i]; 77 while (nz--) { x[*vi--] -= (*v--) * s1; } 78 } 79 PetscCall(VecRestoreArray(xx, &x)); 80 PetscCall(PetscLogFlops(2.0 * (a->nz) - A->cmap->n)); 81 PetscFunctionReturn(0); 82 } 83