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