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(0); 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(0); 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(0); 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(0); 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(0); 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(0); 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(0); 167 } 168