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