1 #include <../src/mat/impls/baij/seq/baij.h>
2 #include <petsc/private/kernels/blockinvert.h>
3
4 /*
5 Special case where the matrix was ILU(0) factored in the natural
6 ordering. This eliminates the need for the column and row permutation.
7 */
MatSolve_SeqBAIJ_1_NaturalOrdering_inplace(Mat A,Vec bb,Vec xx)8 PetscErrorCode MatSolve_SeqBAIJ_1_NaturalOrdering_inplace(Mat A, Vec bb, Vec xx)
9 {
10 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
11 const PetscInt n = a->mbs, *vi, *ai = a->i, *aj = a->j, *diag = a->diag;
12 const MatScalar *aa = a->a, *v;
13 PetscScalar *x;
14 const PetscScalar *b;
15 PetscScalar s1, x1;
16 PetscInt jdx, idt, idx, nz, i;
17
18 PetscFunctionBegin;
19 PetscCall(VecGetArrayRead(bb, &b));
20 PetscCall(VecGetArray(xx, &x));
21
22 /* forward solve the lower triangular */
23 idx = 0;
24 x[0] = b[0];
25 for (i = 1; i < n; i++) {
26 v = aa + ai[i];
27 vi = aj + ai[i];
28 nz = diag[i] - ai[i];
29 idx += 1;
30 s1 = b[idx];
31 while (nz--) {
32 jdx = *vi++;
33 x1 = x[jdx];
34 s1 -= v[0] * x1;
35 v += 1;
36 }
37 x[idx] = s1;
38 }
39 /* backward solve the upper triangular */
40 for (i = n - 1; i >= 0; i--) {
41 v = aa + diag[i] + 1;
42 vi = aj + diag[i] + 1;
43 nz = ai[i + 1] - diag[i] - 1;
44 idt = i;
45 s1 = x[idt];
46 while (nz--) {
47 idx = *vi++;
48 x1 = x[idx];
49 s1 -= v[0] * x1;
50 v += 1;
51 }
52 v = aa + diag[i];
53 x[idt] = v[0] * s1;
54 }
55 PetscCall(VecRestoreArrayRead(bb, &b));
56 PetscCall(VecRestoreArray(xx, &x));
57 PetscCall(PetscLogFlops(2.0 * (a->nz) - A->cmap->n));
58 PetscFunctionReturn(PETSC_SUCCESS);
59 }
60
MatForwardSolve_SeqBAIJ_1_NaturalOrdering(Mat A,Vec bb,Vec xx)61 PetscErrorCode MatForwardSolve_SeqBAIJ_1_NaturalOrdering(Mat A, Vec bb, Vec xx)
62 {
63 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
64 const PetscInt n = a->mbs, *ai = a->i, *aj = a->j, *vi;
65 PetscScalar *x, sum;
66 const PetscScalar *b;
67 const MatScalar *aa = a->a, *v;
68 PetscInt i, nz;
69
70 PetscFunctionBegin;
71 if (!n) PetscFunctionReturn(PETSC_SUCCESS);
72
73 PetscCall(VecGetArrayRead(bb, &b));
74 PetscCall(VecGetArray(xx, &x));
75
76 /* forward solve the lower triangular */
77 x[0] = b[0];
78 v = aa;
79 vi = aj;
80 for (i = 1; i < n; i++) {
81 nz = ai[i + 1] - ai[i];
82 sum = b[i];
83 PetscSparseDenseMinusDot(sum, x, v, vi, nz);
84 v += nz;
85 vi += nz;
86 x[i] = sum;
87 }
88 PetscCall(PetscLogFlops(a->nz - A->cmap->n));
89 PetscCall(VecRestoreArrayRead(bb, &b));
90 PetscCall(VecRestoreArray(xx, &x));
91 PetscFunctionReturn(PETSC_SUCCESS);
92 }
93
MatBackwardSolve_SeqBAIJ_1_NaturalOrdering(Mat A,Vec bb,Vec xx)94 PetscErrorCode MatBackwardSolve_SeqBAIJ_1_NaturalOrdering(Mat A, Vec bb, Vec xx)
95 {
96 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
97 const PetscInt n = a->mbs, *aj = a->j, *adiag = a->diag, *vi;
98 PetscScalar *x, sum;
99 const PetscScalar *b;
100 const MatScalar *aa = a->a, *v;
101 PetscInt i, nz;
102
103 PetscFunctionBegin;
104 if (!n) PetscFunctionReturn(PETSC_SUCCESS);
105
106 PetscCall(VecGetArrayRead(bb, &b));
107 PetscCall(VecGetArray(xx, &x));
108
109 /* backward solve the upper triangular */
110 for (i = n - 1; i >= 0; i--) {
111 v = aa + adiag[i + 1] + 1;
112 vi = aj + adiag[i + 1] + 1;
113 nz = adiag[i] - adiag[i + 1] - 1;
114 sum = b[i];
115 PetscSparseDenseMinusDot(sum, x, v, vi, nz);
116 x[i] = sum * v[nz]; /* x[i]=aa[adiag[i]]*sum; v++; */
117 }
118
119 PetscCall(PetscLogFlops(2.0 * a->nz - A->cmap->n));
120 PetscCall(VecRestoreArrayRead(bb, &b));
121 PetscCall(VecRestoreArray(xx, &x));
122 PetscFunctionReturn(PETSC_SUCCESS);
123 }
124
MatSolve_SeqBAIJ_1_NaturalOrdering(Mat A,Vec bb,Vec xx)125 PetscErrorCode MatSolve_SeqBAIJ_1_NaturalOrdering(Mat A, Vec bb, Vec xx)
126 {
127 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
128 const PetscInt n = a->mbs, *ai = a->i, *aj = a->j, *adiag = a->diag, *vi;
129 PetscScalar *x, sum;
130 const PetscScalar *b;
131 const MatScalar *aa = a->a, *v;
132 PetscInt i, nz;
133
134 PetscFunctionBegin;
135 if (!n) PetscFunctionReturn(PETSC_SUCCESS);
136
137 PetscCall(VecGetArrayRead(bb, &b));
138 PetscCall(VecGetArray(xx, &x));
139
140 /* forward solve the lower triangular */
141 x[0] = b[0];
142 v = aa;
143 vi = aj;
144 for (i = 1; i < n; i++) {
145 nz = ai[i + 1] - ai[i];
146 sum = b[i];
147 PetscSparseDenseMinusDot(sum, x, v, vi, nz);
148 v += nz;
149 vi += nz;
150 x[i] = sum;
151 }
152
153 /* backward solve the upper triangular */
154 for (i = n - 1; i >= 0; i--) {
155 v = aa + adiag[i + 1] + 1;
156 vi = aj + adiag[i + 1] + 1;
157 nz = adiag[i] - adiag[i + 1] - 1;
158 sum = x[i];
159 PetscSparseDenseMinusDot(sum, x, v, vi, nz);
160 x[i] = sum * v[nz]; /* x[i]=aa[adiag[i]]*sum; v++; */
161 }
162
163 PetscCall(PetscLogFlops(2.0 * a->nz - A->cmap->n));
164 PetscCall(VecRestoreArrayRead(bb, &b));
165 PetscCall(VecRestoreArray(xx, &x));
166 PetscFunctionReturn(PETSC_SUCCESS);
167 }
168