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