1 #include <../src/mat/impls/baij/seq/baij.h> 2 #include <petsc/private/kernels/blockinvert.h> 3 4 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; 19 20 PetscFunctionBegin; 21 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 22 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 23 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 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 58 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 59 ierr = PetscLogFlops(2.0*(a->nz) - A->cmap->n);CHKERRQ(ierr); 60 PetscFunctionReturn(0); 61 } 62 63 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; 73 74 PetscFunctionBegin; 75 if (!n) PetscFunctionReturn(0); 76 77 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 78 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 79 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 ierr = PetscLogFlops(a->nz - A->cmap->n);CHKERRQ(ierr); 93 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 94 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 95 PetscFunctionReturn(0); 96 } 97 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; 107 108 PetscFunctionBegin; 109 if (!n) PetscFunctionReturn(0); 110 111 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 112 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 113 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 } 123 124 ierr = PetscLogFlops(2.0*a->nz - A->cmap->n);CHKERRQ(ierr); 125 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 126 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 127 PetscFunctionReturn(0); 128 } 129 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; 139 140 PetscFunctionBegin; 141 if (!n) PetscFunctionReturn(0); 142 143 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 144 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 145 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 } 158 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 } 168 169 ierr = PetscLogFlops(2.0*a->nz - A->cmap->n);CHKERRQ(ierr); 170 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 171 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 172 PetscFunctionReturn(0); 173 } 174 175