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