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;

 21:   VecGetArrayRead(bb,&b);
 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:   }
 57:   VecRestoreArrayRead(bb,&b);
 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);

 77:   VecGetArrayRead(bb,&b);
 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);
 93:   VecRestoreArrayRead(bb,&b);
 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);

111:   VecGetArrayRead(bb,&b);
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;
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:   }

124:   PetscLogFlops(2.0*a->nz - A->cmap->n);
125:   VecRestoreArrayRead(bb,&b);
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);

143:   VecGetArrayRead(bb,&b);
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;
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:   }

169:   PetscLogFlops(2.0*a->nz - A->cmap->n);
170:   VecRestoreArrayRead(bb,&b);
171:   VecRestoreArray(xx,&x);
172:   return(0);
173: }