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: VecGetArrayRead(bb, &b);
19: VecGetArray(xx, &x);
21: /* forward solve the lower triangular */
22: idx = 0;
23: x[0] = b[0];
24: for (i = 1; i < n; i++) {
25: v = aa + ai[i];
26: vi = aj + ai[i];
27: nz = diag[i] - ai[i];
28: idx += 1;
29: s1 = b[idx];
30: while (nz--) {
31: jdx = *vi++;
32: x1 = x[jdx];
33: s1 -= v[0] * x1;
34: v += 1;
35: }
36: x[idx] = s1;
37: }
38: /* backward solve the upper triangular */
39: for (i = n - 1; i >= 0; i--) {
40: v = aa + diag[i] + 1;
41: vi = aj + diag[i] + 1;
42: nz = ai[i + 1] - diag[i] - 1;
43: idt = i;
44: s1 = x[idt];
45: while (nz--) {
46: idx = *vi++;
47: x1 = x[idx];
48: s1 -= v[0] * x1;
49: v += 1;
50: }
51: v = aa + diag[i];
52: x[idt] = v[0] * s1;
53: }
54: VecRestoreArrayRead(bb, &b);
55: VecRestoreArray(xx, &x);
56: PetscLogFlops(2.0 * (a->nz) - A->cmap->n);
57: return 0;
58: }
60: PetscErrorCode MatForwardSolve_SeqBAIJ_1_NaturalOrdering(Mat A, Vec bb, Vec xx)
61: {
62: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
63: const PetscInt n = a->mbs, *ai = a->i, *aj = a->j, *vi;
64: PetscScalar *x, sum;
65: const PetscScalar *b;
66: const MatScalar *aa = a->a, *v;
67: PetscInt i, nz;
69: if (!n) return 0;
71: VecGetArrayRead(bb, &b);
72: VecGetArray(xx, &x);
74: /* forward solve the lower triangular */
75: x[0] = b[0];
76: v = aa;
77: vi = aj;
78: for (i = 1; i < n; i++) {
79: nz = ai[i + 1] - ai[i];
80: sum = b[i];
81: PetscSparseDenseMinusDot(sum, x, v, vi, nz);
82: v += nz;
83: vi += nz;
84: x[i] = sum;
85: }
86: PetscLogFlops(a->nz - A->cmap->n);
87: VecRestoreArrayRead(bb, &b);
88: VecRestoreArray(xx, &x);
89: return 0;
90: }
92: PetscErrorCode MatBackwardSolve_SeqBAIJ_1_NaturalOrdering(Mat A, Vec bb, Vec xx)
93: {
94: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
95: const PetscInt n = a->mbs, *aj = a->j, *adiag = a->diag, *vi;
96: PetscScalar *x, sum;
97: const PetscScalar *b;
98: const MatScalar *aa = a->a, *v;
99: PetscInt i, nz;
101: if (!n) return 0;
103: VecGetArrayRead(bb, &b);
104: VecGetArray(xx, &x);
106: /* backward solve the upper triangular */
107: for (i = n - 1; i >= 0; i--) {
108: v = aa + adiag[i + 1] + 1;
109: vi = aj + adiag[i + 1] + 1;
110: nz = adiag[i] - adiag[i + 1] - 1;
111: sum = b[i];
112: PetscSparseDenseMinusDot(sum, x, v, vi, nz);
113: x[i] = sum * v[nz]; /* x[i]=aa[adiag[i]]*sum; v++; */
114: }
116: PetscLogFlops(2.0 * a->nz - A->cmap->n);
117: VecRestoreArrayRead(bb, &b);
118: VecRestoreArray(xx, &x);
119: return 0;
120: }
122: PetscErrorCode MatSolve_SeqBAIJ_1_NaturalOrdering(Mat A, Vec bb, Vec xx)
123: {
124: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
125: const PetscInt n = a->mbs, *ai = a->i, *aj = a->j, *adiag = a->diag, *vi;
126: PetscScalar *x, sum;
127: const PetscScalar *b;
128: const MatScalar *aa = a->a, *v;
129: PetscInt i, nz;
131: if (!n) return 0;
133: VecGetArrayRead(bb, &b);
134: VecGetArray(xx, &x);
136: /* forward solve the lower triangular */
137: x[0] = b[0];
138: v = aa;
139: vi = aj;
140: for (i = 1; i < n; i++) {
141: nz = ai[i + 1] - ai[i];
142: sum = b[i];
143: PetscSparseDenseMinusDot(sum, x, v, vi, nz);
144: v += nz;
145: vi += nz;
146: x[i] = sum;
147: }
149: /* backward solve the upper triangular */
150: for (i = n - 1; i >= 0; i--) {
151: v = aa + adiag[i + 1] + 1;
152: vi = aj + adiag[i + 1] + 1;
153: nz = adiag[i] - adiag[i + 1] - 1;
154: sum = x[i];
155: PetscSparseDenseMinusDot(sum, x, v, vi, nz);
156: x[i] = sum * v[nz]; /* x[i]=aa[adiag[i]]*sum; v++; */
157: }
159: PetscLogFlops(2.0 * a->nz - A->cmap->n);
160: VecRestoreArrayRead(bb, &b);
161: VecRestoreArray(xx, &x);
162: return 0;
163: }