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