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