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