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: }