1 #include <../src/mat/impls/baij/seq/baij.h> 2 #include <petsc/private/kernels/blockinvert.h> 3 4 PetscErrorCode MatSolveTranspose_SeqBAIJ_7_inplace(Mat A,Vec bb,Vec xx) 5 { 6 Mat_SeqBAIJ *a =(Mat_SeqBAIJ*)A->data; 7 IS iscol=a->col,isrow=a->row; 8 PetscErrorCode ierr; 9 const PetscInt *r,*c,*rout,*cout; 10 const PetscInt *diag=a->diag,n=a->mbs,*vi,*ai=a->i,*aj=a->j; 11 PetscInt i,nz,idx,idt,ii,ic,ir,oidx; 12 const MatScalar *aa=a->a,*v; 13 PetscScalar s1,s2,s3,s4,s5,s6,s7,x1,x2,x3,x4,x5,x6,x7,*x,*t; 14 const PetscScalar *b; 15 16 PetscFunctionBegin; 17 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 18 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 19 t = a->solve_work; 20 21 ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout; 22 ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout; 23 24 /* copy the b into temp work space according to permutation */ 25 ii = 0; 26 for (i=0; i<n; i++) { 27 ic = 7*c[i]; 28 t[ii] = b[ic]; 29 t[ii+1] = b[ic+1]; 30 t[ii+2] = b[ic+2]; 31 t[ii+3] = b[ic+3]; 32 t[ii+4] = b[ic+4]; 33 t[ii+5] = b[ic+5]; 34 t[ii+6] = b[ic+6]; 35 ii += 7; 36 } 37 38 /* forward solve the U^T */ 39 idx = 0; 40 for (i=0; i<n; i++) { 41 42 v = aa + 49*diag[i]; 43 /* multiply by the inverse of the block diagonal */ 44 x1 = t[idx]; x2 = t[1+idx]; x3 = t[2+idx]; x4 = t[3+idx]; x5 = t[4+idx]; 45 x6 = t[5+idx]; x7 = t[6+idx]; 46 s1 = v[0]*x1 + v[1]*x2 + v[2]*x3 + v[3]*x4 + v[4]*x5 + v[5]*x6 + v[6]*x7; 47 s2 = v[7]*x1 + v[8]*x2 + v[9]*x3 + v[10]*x4 + v[11]*x5 + v[12]*x6 + v[13]*x7; 48 s3 = v[14]*x1 + v[15]*x2 + v[16]*x3 + v[17]*x4 + v[18]*x5 + v[19]*x6 + v[20]*x7; 49 s4 = v[21]*x1 + v[22]*x2 + v[23]*x3 + v[24]*x4 + v[25]*x5 + v[26]*x6 + v[27]*x7; 50 s5 = v[28]*x1 + v[29]*x2 + v[30]*x3 + v[31]*x4 + v[32]*x5 + v[33]*x6 + v[34]*x7; 51 s6 = v[35]*x1 + v[36]*x2 + v[37]*x3 + v[38]*x4 + v[39]*x5 + v[40]*x6 + v[41]*x7; 52 s7 = v[42]*x1 + v[43]*x2 + v[44]*x3 + v[45]*x4 + v[46]*x5 + v[47]*x6 + v[48]*x7; 53 v += 49; 54 55 vi = aj + diag[i] + 1; 56 nz = ai[i+1] - diag[i] - 1; 57 while (nz--) { 58 oidx = 7*(*vi++); 59 t[oidx] -= v[0]*s1 + v[1]*s2 + v[2]*s3 + v[3]*s4 + v[4]*s5 + v[5]*s6 + v[6]*s7; 60 t[oidx+1] -= v[7]*s1 + v[8]*s2 + v[9]*s3 + v[10]*s4 + v[11]*s5 + v[12]*s6 + v[13]*s7; 61 t[oidx+2] -= v[14]*s1 + v[15]*s2 + v[16]*s3 + v[17]*s4 + v[18]*s5 + v[19]*s6 + v[20]*s7; 62 t[oidx+3] -= v[21]*s1 + v[22]*s2 + v[23]*s3 + v[24]*s4 + v[25]*s5 + v[26]*s6 + v[27]*s7; 63 t[oidx+4] -= v[28]*s1 + v[29]*s2 + v[30]*s3 + v[31]*s4 + v[32]*s5 + v[33]*s6 + v[34]*s7; 64 t[oidx+5] -= v[35]*s1 + v[36]*s2 + v[37]*s3 + v[38]*s4 + v[39]*s5 + v[40]*s6 + v[41]*s7; 65 t[oidx+6] -= v[42]*s1 + v[43]*s2 + v[44]*s3 + v[45]*s4 + v[46]*s5 + v[47]*s6 + v[48]*s7; 66 v += 49; 67 } 68 t[idx] = s1;t[1+idx] = s2; t[2+idx] = s3;t[3+idx] = s4; t[4+idx] = s5; 69 t[5+idx] = s6;t[6+idx] = s7; 70 idx += 7; 71 } 72 /* backward solve the L^T */ 73 for (i=n-1; i>=0; i--) { 74 v = aa + 49*diag[i] - 49; 75 vi = aj + diag[i] - 1; 76 nz = diag[i] - ai[i]; 77 idt = 7*i; 78 s1 = t[idt]; s2 = t[1+idt]; s3 = t[2+idt];s4 = t[3+idt]; s5 = t[4+idt]; 79 s6 = t[5+idt];s7 = t[6+idt]; 80 while (nz--) { 81 idx = 7*(*vi--); 82 t[idx] -= v[0]*s1 + v[1]*s2 + v[2]*s3 + v[3]*s4 + v[4]*s5 + v[5]*s6 + v[6]*s7; 83 t[idx+1] -= v[7]*s1 + v[8]*s2 + v[9]*s3 + v[10]*s4 + v[11]*s5 + v[12]*s6 + v[13]*s7; 84 t[idx+2] -= v[14]*s1 + v[15]*s2 + v[16]*s3 + v[17]*s4 + v[18]*s5 + v[19]*s6 + v[20]*s7; 85 t[idx+3] -= v[21]*s1 + v[22]*s2 + v[23]*s3 + v[24]*s4 + v[25]*s5 + v[26]*s6 + v[27]*s7; 86 t[idx+4] -= v[28]*s1 + v[29]*s2 + v[30]*s3 + v[31]*s4 + v[32]*s5 + v[33]*s6 + v[34]*s7; 87 t[idx+5] -= v[35]*s1 + v[36]*s2 + v[37]*s3 + v[38]*s4 + v[39]*s5 + v[40]*s6 + v[41]*s7; 88 t[idx+6] -= v[42]*s1 + v[43]*s2 + v[44]*s3 + v[45]*s4 + v[46]*s5 + v[47]*s6 + v[48]*s7; 89 v -= 49; 90 } 91 } 92 93 /* copy t into x according to permutation */ 94 ii = 0; 95 for (i=0; i<n; i++) { 96 ir = 7*r[i]; 97 x[ir] = t[ii]; 98 x[ir+1] = t[ii+1]; 99 x[ir+2] = t[ii+2]; 100 x[ir+3] = t[ii+3]; 101 x[ir+4] = t[ii+4]; 102 x[ir+5] = t[ii+5]; 103 x[ir+6] = t[ii+6]; 104 ii += 7; 105 } 106 107 ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr); 108 ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr); 109 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 110 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 111 ierr = PetscLogFlops(2.0*49*(a->nz) - 7.0*A->cmap->n);CHKERRQ(ierr); 112 PetscFunctionReturn(0); 113 } 114 PetscErrorCode MatSolveTranspose_SeqBAIJ_7(Mat A,Vec bb,Vec xx) 115 { 116 Mat_SeqBAIJ *a=(Mat_SeqBAIJ*)A->data; 117 PetscErrorCode ierr; 118 IS iscol=a->col,isrow=a->row; 119 const PetscInt n =a->mbs,*vi,*ai=a->i,*aj=a->j,*diag=a->diag; 120 const PetscInt *r,*c,*rout,*cout; 121 PetscInt nz,idx,idt,j,i,oidx,ii,ic,ir; 122 const PetscInt bs =A->rmap->bs,bs2=a->bs2; 123 const MatScalar *aa=a->a,*v; 124 PetscScalar s1,s2,s3,s4,s5,s6,s7,x1,x2,x3,x4,x5,x6,x7,*x,*t; 125 const PetscScalar *b; 126 127 PetscFunctionBegin; 128 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 129 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 130 t = a->solve_work; 131 132 ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout; 133 ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout; 134 135 /* copy b into temp work space according to permutation */ 136 for (i=0; i<n; i++) { 137 ii = bs*i; ic = bs*c[i]; 138 t[ii] = b[ic]; t[ii+1] = b[ic+1]; t[ii+2] = b[ic+2]; t[ii+3] = b[ic+3]; 139 t[ii+4] = b[ic+4]; t[ii+5] = b[ic+5]; t[ii+6] = b[ic+6]; 140 } 141 142 /* forward solve the U^T */ 143 idx = 0; 144 for (i=0; i<n; i++) { 145 v = aa + bs2*diag[i]; 146 /* multiply by the inverse of the block diagonal */ 147 x1 = t[idx]; x2 = t[1+idx]; x3 = t[2+idx]; x4 = t[3+idx]; x5 = t[4+idx]; 148 x6 = t[5+idx]; x7 = t[6+idx]; 149 s1 = v[0]*x1 + v[1]*x2 + v[2]*x3 + v[3]*x4 + v[4]*x5 + v[5]*x6 + v[6]*x7; 150 s2 = v[7]*x1 + v[8]*x2 + v[9]*x3 + v[10]*x4 + v[11]*x5 + v[12]*x6 + v[13]*x7; 151 s3 = v[14]*x1 + v[15]*x2 + v[16]*x3 + v[17]*x4 + v[18]*x5 + v[19]*x6 + v[20]*x7; 152 s4 = v[21]*x1 + v[22]*x2 + v[23]*x3 + v[24]*x4 + v[25]*x5 + v[26]*x6 + v[27]*x7; 153 s5 = v[28]*x1 + v[29]*x2 + v[30]*x3 + v[31]*x4 + v[32]*x5 + v[33]*x6 + v[34]*x7; 154 s6 = v[35]*x1 + v[36]*x2 + v[37]*x3 + v[38]*x4 + v[39]*x5 + v[40]*x6 + v[41]*x7; 155 s7 = v[42]*x1 + v[43]*x2 + v[44]*x3 + v[45]*x4 + v[46]*x5 + v[47]*x6 + v[48]*x7; 156 v -= bs2; 157 158 vi = aj + diag[i] - 1; 159 nz = diag[i] - diag[i+1] - 1; 160 for (j=0; j>-nz; j--) { 161 oidx = bs*vi[j]; 162 t[oidx] -= v[0]*s1 + v[1]*s2 + v[2]*s3 + v[3]*s4 + v[4]*s5 + v[5]*s6 + v[6]*s7; 163 t[oidx+1] -= v[7]*s1 + v[8]*s2 + v[9]*s3 + v[10]*s4 + v[11]*s5 + v[12]*s6 + v[13]*s7; 164 t[oidx+2] -= v[14]*s1 + v[15]*s2 + v[16]*s3 + v[17]*s4 + v[18]*s5 + v[19]*s6 + v[20]*s7; 165 t[oidx+3] -= v[21]*s1 + v[22]*s2 + v[23]*s3 + v[24]*s4 + v[25]*s5 + v[26]*s6 + v[27]*s7; 166 t[oidx+4] -= v[28]*s1 + v[29]*s2 + v[30]*s3 + v[31]*s4 + v[32]*s5 + v[33]*s6 + v[34]*s7; 167 t[oidx+5] -= v[35]*s1 + v[36]*s2 + v[37]*s3 + v[38]*s4 + v[39]*s5 + v[40]*s6 + v[41]*s7; 168 t[oidx+6] -= v[42]*s1 + v[43]*s2 + v[44]*s3 + v[45]*s4 + v[46]*s5 + v[47]*s6 + v[48]*s7; 169 v -= bs2; 170 } 171 t[idx] = s1;t[1+idx] = s2; t[2+idx] = s3; t[3+idx] = s4; t[4+idx] =s5; 172 t[5+idx] = s6; t[6+idx] = s7; 173 idx += bs; 174 } 175 /* backward solve the L^T */ 176 for (i=n-1; i>=0; i--) { 177 v = aa + bs2*ai[i]; 178 vi = aj + ai[i]; 179 nz = ai[i+1] - ai[i]; 180 idt = bs*i; 181 s1 = t[idt]; s2 = t[1+idt]; s3 = t[2+idt]; s4 = t[3+idt]; s5 = t[4+idt]; 182 s6 = t[5+idt]; s7 = t[6+idt]; 183 for (j=0; j<nz; j++) { 184 idx = bs*vi[j]; 185 t[idx] -= v[0]*s1 + v[1]*s2 + v[2]*s3 + v[3]*s4 + v[4]*s5 + v[5]*s6 + v[6]*s7; 186 t[idx+1] -= v[7]*s1 + v[8]*s2 + v[9]*s3 + v[10]*s4 + v[11]*s5 + v[12]*s6 + v[13]*s7; 187 t[idx+2] -= v[14]*s1 + v[15]*s2 + v[16]*s3 + v[17]*s4 + v[18]*s5 + v[19]*s6 + v[20]*s7; 188 t[idx+3] -= v[21]*s1 + v[22]*s2 + v[23]*s3 + v[24]*s4 + v[25]*s5 + v[26]*s6 + v[27]*s7; 189 t[idx+4] -= v[28]*s1 + v[29]*s2 + v[30]*s3 + v[31]*s4 + v[32]*s5 + v[33]*s6 + v[34]*s7; 190 t[idx+5] -= v[35]*s1 + v[36]*s2 + v[37]*s3 + v[38]*s4 + v[39]*s5 + v[40]*s6 + v[41]*s7; 191 t[idx+6] -= v[42]*s1 + v[43]*s2 + v[44]*s3 + v[45]*s4 + v[46]*s5 + v[47]*s6 + v[48]*s7; 192 v += bs2; 193 } 194 } 195 196 /* copy t into x according to permutation */ 197 for (i=0; i<n; i++) { 198 ii = bs*i; ir = bs*r[i]; 199 x[ir] = t[ii]; x[ir+1] = t[ii+1]; x[ir+2] = t[ii+2]; x[ir+3] = t[ii+3]; 200 x[ir+4] = t[ii+4]; x[ir+5] = t[ii+5]; x[ir+6] = t[ii+6]; 201 } 202 203 ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr); 204 ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr); 205 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 206 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 207 ierr = PetscLogFlops(2.0*bs2*(a->nz) - bs*A->cmap->n);CHKERRQ(ierr); 208 PetscFunctionReturn(0); 209 } 210