Actual source code: baijsolvnat1.c
1: #include <../src/mat/impls/baij/seq/baij.h>
2: #include <petsc/private/kernels/blockinvert.h>
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: */
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;
18: PetscFunctionBegin;
19: PetscCall(VecGetArrayRead(bb, &b));
20: PetscCall(VecGetArray(xx, &x));
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: }
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;
70: PetscFunctionBegin;
71: if (!n) PetscFunctionReturn(PETSC_SUCCESS);
73: PetscCall(VecGetArrayRead(bb, &b));
74: PetscCall(VecGetArray(xx, &x));
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: }
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;
103: PetscFunctionBegin;
104: if (!n) PetscFunctionReturn(PETSC_SUCCESS);
106: PetscCall(VecGetArrayRead(bb, &b));
107: PetscCall(VecGetArray(xx, &x));
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: }
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: }
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;
134: PetscFunctionBegin;
135: if (!n) PetscFunctionReturn(PETSC_SUCCESS);
137: PetscCall(VecGetArrayRead(bb, &b));
138: PetscCall(VecGetArray(xx, &x));
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: }
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: }
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: }