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