Actual source code: baijsolvnat1.c

petsc-master 2019-12-08
Report Typos and Errors
```  1:  #include <../src/mat/impls/baij/seq/baij.h>
2:  #include <petsc/private/kernels/blockinvert.h>

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: */
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;

22:   VecGetArray(xx,&x);

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:   }
58:   VecRestoreArray(xx,&x);
59:   PetscLogFlops(2.0*(a->nz) - A->cmap->n);
60:   return(0);
61: }

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;

75:   if (!n) return(0);

78:   VecGetArray(xx,&x);

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:   PetscLogFlops(a->nz - A->cmap->n);
94:   VecRestoreArray(xx,&x);
95:   return(0);
96: }

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;

109:   if (!n) return(0);

112:   VecGetArray(xx,&x);

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

124:   PetscLogFlops(2.0*a->nz - A->cmap->n);
126:   VecRestoreArray(xx,&x);
127:   return(0);
128: }

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;

141:   if (!n) return(0);

144:   VecGetArray(xx,&x);

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

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;