1 #include <../src/mat/impls/baij/seq/baij.h> 2 #include <petsc/private/kernels/blockinvert.h> 3 4 PetscErrorCode MatSolve_SeqBAIJ_7_NaturalOrdering_inplace(Mat A,Vec bb,Vec xx) 5 { 6 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 7 const PetscInt *diag=a->diag,n=a->mbs,*vi,*ai=a->i,*aj=a->j; 8 PetscErrorCode ierr; 9 PetscInt i,nz,idx,idt,jdx; 10 const MatScalar *aa=a->a,*v; 11 PetscScalar *x,s1,s2,s3,s4,s5,s6,s7,x1,x2,x3,x4,x5,x6,x7; 12 const PetscScalar *b; 13 14 PetscFunctionBegin; 15 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 16 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 17 /* forward solve the lower triangular */ 18 idx = 0; 19 x[0] = b[idx]; x[1] = b[1+idx]; x[2] = b[2+idx]; 20 x[3] = b[3+idx]; x[4] = b[4+idx]; x[5] = b[5+idx]; 21 x[6] = b[6+idx]; 22 for (i=1; i<n; i++) { 23 v = aa + 49*ai[i]; 24 vi = aj + ai[i]; 25 nz = diag[i] - ai[i]; 26 idx = 7*i; 27 s1 = b[idx]; s2 = b[1+idx]; s3 = b[2+idx]; 28 s4 = b[3+idx]; s5 = b[4+idx]; s6 = b[5+idx]; 29 s7 = b[6+idx]; 30 while (nz--) { 31 jdx = 7*(*vi++); 32 x1 = x[jdx]; x2 = x[1+jdx]; x3 = x[2+jdx]; 33 x4 = x[3+jdx]; x5 = x[4+jdx]; x6 = x[5+jdx]; 34 x7 = x[6+jdx]; 35 s1 -= v[0]*x1 + v[7]*x2 + v[14]*x3 + v[21]*x4 + v[28]*x5 + v[35]*x6 + v[42]*x7; 36 s2 -= v[1]*x1 + v[8]*x2 + v[15]*x3 + v[22]*x4 + v[29]*x5 + v[36]*x6 + v[43]*x7; 37 s3 -= v[2]*x1 + v[9]*x2 + v[16]*x3 + v[23]*x4 + v[30]*x5 + v[37]*x6 + v[44]*x7; 38 s4 -= v[3]*x1 + v[10]*x2 + v[17]*x3 + v[24]*x4 + v[31]*x5 + v[38]*x6 + v[45]*x7; 39 s5 -= v[4]*x1 + v[11]*x2 + v[18]*x3 + v[25]*x4 + v[32]*x5 + v[39]*x6 + v[46]*x7; 40 s6 -= v[5]*x1 + v[12]*x2 + v[19]*x3 + v[26]*x4 + v[33]*x5 + v[40]*x6 + v[47]*x7; 41 s7 -= v[6]*x1 + v[13]*x2 + v[20]*x3 + v[27]*x4 + v[34]*x5 + v[41]*x6 + v[48]*x7; 42 v += 49; 43 } 44 x[idx] = s1; 45 x[1+idx] = s2; 46 x[2+idx] = s3; 47 x[3+idx] = s4; 48 x[4+idx] = s5; 49 x[5+idx] = s6; 50 x[6+idx] = s7; 51 } 52 /* backward solve the upper triangular */ 53 for (i=n-1; i>=0; i--) { 54 v = aa + 49*diag[i] + 49; 55 vi = aj + diag[i] + 1; 56 nz = ai[i+1] - diag[i] - 1; 57 idt = 7*i; 58 s1 = x[idt]; s2 = x[1+idt]; 59 s3 = x[2+idt]; s4 = x[3+idt]; 60 s5 = x[4+idt]; s6 = x[5+idt]; 61 s7 = x[6+idt]; 62 while (nz--) { 63 idx = 7*(*vi++); 64 x1 = x[idx]; x2 = x[1+idx]; x3 = x[2+idx]; 65 x4 = x[3+idx]; x5 = x[4+idx]; x6 = x[5+idx]; 66 x7 = x[6+idx]; 67 s1 -= v[0]*x1 + v[7]*x2 + v[14]*x3 + v[21]*x4 + v[28]*x5 + v[35]*x6 + v[42]*x7; 68 s2 -= v[1]*x1 + v[8]*x2 + v[15]*x3 + v[22]*x4 + v[29]*x5 + v[36]*x6 + v[43]*x7; 69 s3 -= v[2]*x1 + v[9]*x2 + v[16]*x3 + v[23]*x4 + v[30]*x5 + v[37]*x6 + v[44]*x7; 70 s4 -= v[3]*x1 + v[10]*x2 + v[17]*x3 + v[24]*x4 + v[31]*x5 + v[38]*x6 + v[45]*x7; 71 s5 -= v[4]*x1 + v[11]*x2 + v[18]*x3 + v[25]*x4 + v[32]*x5 + v[39]*x6 + v[46]*x7; 72 s6 -= v[5]*x1 + v[12]*x2 + v[19]*x3 + v[26]*x4 + v[33]*x5 + v[40]*x6 + v[47]*x7; 73 s7 -= v[6]*x1 + v[13]*x2 + v[20]*x3 + v[27]*x4 + v[34]*x5 + v[41]*x6 + v[48]*x7; 74 v += 49; 75 } 76 v = aa + 49*diag[i]; 77 x[idt] = v[0]*s1 + v[7]*s2 + v[14]*s3 + v[21]*s4 78 + v[28]*s5 + v[35]*s6 + v[42]*s7; 79 x[1+idt] = v[1]*s1 + v[8]*s2 + v[15]*s3 + v[22]*s4 80 + v[29]*s5 + v[36]*s6 + v[43]*s7; 81 x[2+idt] = v[2]*s1 + v[9]*s2 + v[16]*s3 + v[23]*s4 82 + v[30]*s5 + v[37]*s6 + v[44]*s7; 83 x[3+idt] = v[3]*s1 + v[10]*s2 + v[17]*s3 + v[24]*s4 84 + v[31]*s5 + v[38]*s6 + v[45]*s7; 85 x[4+idt] = v[4]*s1 + v[11]*s2 + v[18]*s3 + v[25]*s4 86 + v[32]*s5 + v[39]*s6 + v[46]*s7; 87 x[5+idt] = v[5]*s1 + v[12]*s2 + v[19]*s3 + v[26]*s4 88 + v[33]*s5 + v[40]*s6 + v[47]*s7; 89 x[6+idt] = v[6]*s1 + v[13]*s2 + v[20]*s3 + v[27]*s4 90 + v[34]*s5 + v[41]*s6 + v[48]*s7; 91 } 92 93 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 94 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 95 ierr = PetscLogFlops(2.0*36*(a->nz) - 6.0*A->cmap->n);CHKERRQ(ierr); 96 PetscFunctionReturn(0); 97 } 98 99 PetscErrorCode MatSolve_SeqBAIJ_7_NaturalOrdering(Mat A,Vec bb,Vec xx) 100 { 101 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 102 const PetscInt n =a->mbs,*vi,*ai=a->i,*aj=a->j,*adiag=a->diag; 103 PetscErrorCode ierr; 104 PetscInt i,k,nz,idx,jdx,idt; 105 const PetscInt bs = A->rmap->bs,bs2 = a->bs2; 106 const MatScalar *aa=a->a,*v; 107 PetscScalar *x; 108 const PetscScalar *b; 109 PetscScalar s1,s2,s3,s4,s5,s6,s7,x1,x2,x3,x4,x5,x6,x7; 110 111 PetscFunctionBegin; 112 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 113 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 114 /* forward solve the lower triangular */ 115 idx = 0; 116 x[0] = b[idx]; x[1] = b[1+idx];x[2] = b[2+idx];x[3] = b[3+idx]; 117 x[4] = b[4+idx];x[5] = b[5+idx];x[6] = b[6+idx]; 118 for (i=1; i<n; i++) { 119 v = aa + bs2*ai[i]; 120 vi = aj + ai[i]; 121 nz = ai[i+1] - ai[i]; 122 idx = bs*i; 123 s1 = b[idx];s2 = b[1+idx];s3 = b[2+idx];s4 = b[3+idx]; 124 s5 = b[4+idx];s6 = b[5+idx];s7 = b[6+idx]; 125 for (k=0; k<nz; k++) { 126 jdx = bs*vi[k]; 127 x1 = x[jdx];x2 = x[1+jdx]; x3 =x[2+jdx];x4 =x[3+jdx]; 128 x5 = x[4+jdx]; x6 = x[5+jdx];x7 = x[6+jdx]; 129 s1 -= v[0]*x1 + v[7]*x2 + v[14]*x3 + v[21]*x4 + v[28]*x5 + v[35]*x6 + v[42]*x7; 130 s2 -= v[1]*x1 + v[8]*x2 + v[15]*x3 + v[22]*x4 + v[29]*x5 + v[36]*x6 + v[43]*x7; 131 s3 -= v[2]*x1 + v[9]*x2 + v[16]*x3 + v[23]*x4 + v[30]*x5 + v[37]*x6 + v[44]*x7; 132 s4 -= v[3]*x1 + v[10]*x2 + v[17]*x3 + v[24]*x4 + v[31]*x5 + v[38]*x6 + v[45]*x7; 133 s5 -= v[4]*x1 + v[11]*x2 + v[18]*x3 + v[25]*x4 + v[32]*x5 + v[39]*x6 + v[46]*x7; 134 s6 -= v[5]*x1 + v[12]*x2 + v[19]*x3 + v[26]*x4 + v[33]*x5 + v[40]*x6 + v[47]*x7; 135 s7 -= v[6]*x1 + v[13]*x2 + v[20]*x3 + v[27]*x4 + v[34]*x5 + v[41]*x6 + v[48]*x7; 136 v += bs2; 137 } 138 139 x[idx] = s1; 140 x[1+idx] = s2; 141 x[2+idx] = s3; 142 x[3+idx] = s4; 143 x[4+idx] = s5; 144 x[5+idx] = s6; 145 x[6+idx] = s7; 146 } 147 148 /* backward solve the upper triangular */ 149 for (i=n-1; i>=0; i--) { 150 v = aa + bs2*(adiag[i+1]+1); 151 vi = aj + adiag[i+1]+1; 152 nz = adiag[i] - adiag[i+1]-1; 153 idt = bs*i; 154 s1 = x[idt]; s2 = x[1+idt];s3 = x[2+idt];s4 = x[3+idt]; 155 s5 = x[4+idt];s6 = x[5+idt];s7 = x[6+idt]; 156 for (k=0; k<nz; k++) { 157 idx = bs*vi[k]; 158 x1 = x[idx]; x2 = x[1+idx]; x3 = x[2+idx];x4 = x[3+idx]; 159 x5 = x[4+idx];x6 = x[5+idx];x7 = x[6+idx]; 160 s1 -= v[0]*x1 + v[7]*x2 + v[14]*x3 + v[21]*x4 + v[28]*x5 + v[35]*x6 + v[42]*x7; 161 s2 -= v[1]*x1 + v[8]*x2 + v[15]*x3 + v[22]*x4 + v[29]*x5 + v[36]*x6 + v[43]*x7; 162 s3 -= v[2]*x1 + v[9]*x2 + v[16]*x3 + v[23]*x4 + v[30]*x5 + v[37]*x6 + v[44]*x7; 163 s4 -= v[3]*x1 + v[10]*x2 + v[17]*x3 + v[24]*x4 + v[31]*x5 + v[38]*x6 + v[45]*x7; 164 s5 -= v[4]*x1 + v[11]*x2 + v[18]*x3 + v[25]*x4 + v[32]*x5 + v[39]*x6 + v[46]*x7; 165 s6 -= v[5]*x1 + v[12]*x2 + v[19]*x3 + v[26]*x4 + v[33]*x5 + v[40]*x6 + v[47]*x7; 166 s7 -= v[6]*x1 + v[13]*x2 + v[20]*x3 + v[27]*x4 + v[34]*x5 + v[41]*x6 + v[48]*x7; 167 v += bs2; 168 } 169 /* x = inv_diagonal*x */ 170 x[idt] = v[0]*s1 + v[7]*s2 + v[14]*s3 + v[21]*s4 + v[28]*s5 + v[35]*s6 + v[42]*s7; 171 x[1+idt] = v[1]*s1 + v[8]*s2 + v[15]*s3 + v[22]*s4 + v[29]*s5 + v[36]*s6 + v[43]*s7; 172 x[2+idt] = v[2]*s1 + v[9]*s2 + v[16]*s3 + v[23]*s4 + v[30]*s5 + v[37]*s6 + v[44]*s7; 173 x[3+idt] = v[3]*s1 + v[10]*s2 + v[17]*s3 + v[24]*s4 + v[31]*s5 + v[38]*s6 + v[45]*s7; 174 x[4+idt] = v[4]*s1 + v[11]*s2 + v[18]*s3 + v[25]*s4 + v[32]*s5 + v[39]*s6 + v[46]*s7; 175 x[5+idt] = v[5]*s1 + v[12]*s2 + v[19]*s3 + v[26]*s4 + v[33]*s5 + v[40]*s6 + v[47]*s7; 176 x[6+idt] = v[6]*s1 + v[13]*s2 + v[20]*s3 + v[27]*s4 + v[34]*s5 + v[41]*s6 + v[48]*s7; 177 } 178 179 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 180 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 181 ierr = PetscLogFlops(2.0*bs2*(a->nz) - bs*A->cmap->n);CHKERRQ(ierr); 182 PetscFunctionReturn(0); 183 } 184 185