1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: baijfact2.c,v 1.14 1998/10/13 02:11:28 bsmith Exp bsmith $"; 3 #endif 4 /* 5 Factorization code for BAIJ format. 6 */ 7 8 #include "src/mat/impls/baij/seq/baij.h" 9 #include "src/vec/vecimpl.h" 10 #include "src/inline/ilu.h" 11 #include "src/inline/dot.h" 12 13 /* ----------------------------------------------------------- */ 14 #undef __FUNC__ 15 #define __FUNC__ "MatSolve_SeqBAIJ_N" 16 int MatSolve_SeqBAIJ_N(Mat A,Vec bb,Vec xx) 17 { 18 Mat_SeqBAIJ *a=(Mat_SeqBAIJ *)A->data; 19 IS iscol=a->col,isrow=a->row; 20 int *r,*c,ierr,i,n=a->mbs,*vi,*ai=a->i,*aj=a->j; 21 int nz,bs=a->bs,bs2=a->bs2,*rout,*cout; 22 Scalar *aa=a->a,*sum; 23 Scalar *x,*b,*lsum,*tmp,*v; 24 25 PetscFunctionBegin; 26 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 27 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 28 tmp = a->solve_work; 29 30 ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout; 31 ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout + (n-1); 32 33 /* forward solve the lower triangular */ 34 PetscMemcpy(tmp,b + bs*(*r++), bs*sizeof(Scalar)); 35 for ( i=1; i<n; i++ ) { 36 v = aa + bs2*ai[i]; 37 vi = aj + ai[i]; 38 nz = a->diag[i] - ai[i]; 39 sum = tmp + bs*i; 40 PetscMemcpy(sum,b+bs*(*r++),bs*sizeof(Scalar)); 41 while (nz--) { 42 Kernel_v_gets_v_minus_A_times_w(bs,sum,v,tmp+bs*(*vi++)); 43 v += bs2; 44 } 45 } 46 /* backward solve the upper triangular */ 47 lsum = a->solve_work + a->n; 48 for ( i=n-1; i>=0; i-- ){ 49 v = aa + bs2*(a->diag[i] + 1); 50 vi = aj + a->diag[i] + 1; 51 nz = ai[i+1] - a->diag[i] - 1; 52 PetscMemcpy(lsum,tmp+i*bs,bs*sizeof(Scalar)); 53 while (nz--) { 54 Kernel_v_gets_v_minus_A_times_w(bs,lsum,v,tmp+bs*(*vi++)); 55 v += bs2; 56 } 57 Kernel_w_gets_A_times_v(bs,lsum,aa+bs2*a->diag[i],tmp+i*bs); 58 PetscMemcpy(x + bs*(*c--),tmp+i*bs,bs*sizeof(Scalar)); 59 } 60 61 ierr = ISRestoreIndices(isrow,&rout); CHKERRQ(ierr); 62 ierr = ISRestoreIndices(iscol,&cout); CHKERRQ(ierr); 63 ierr = VecRestoreArray(bb,&b); CHKERRQ(ierr); 64 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 65 PLogFlops(2*(a->bs2)*(a->nz) - a->n); 66 PetscFunctionReturn(0); 67 } 68 69 #undef __FUNC__ 70 #define __FUNC__ "MatSolve_SeqBAIJ_7" 71 int MatSolve_SeqBAIJ_7(Mat A,Vec bb,Vec xx) 72 { 73 Mat_SeqBAIJ *a=(Mat_SeqBAIJ *)A->data; 74 IS iscol=a->col,isrow=a->row; 75 int *r,*c,ierr,i,n=a->mbs,*vi,*ai=a->i,*aj=a->j,nz,idx,idt,idc,*rout,*cout; 76 int *diag = a->diag; 77 Scalar *aa=a->a,sum1,sum2,sum3,sum4,sum5,sum6,sum7,x1,x2,x3,x4,x5,x6,x7; 78 Scalar *x,*b,*tmp,*v; 79 80 PetscFunctionBegin; 81 ierr = VecGetArray(bb,&b); CHKERRQ(ierr); 82 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 83 tmp = a->solve_work; 84 85 ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout; 86 ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout + (n-1); 87 88 /* forward solve the lower triangular */ 89 idx = 7*(*r++); 90 tmp[0] = b[idx]; tmp[1] = b[1+idx]; 91 tmp[2] = b[2+idx]; tmp[3] = b[3+idx]; tmp[4] = b[4+idx]; 92 tmp[5] = b[5+idx]; tmp[6] = b[6+idx]; 93 94 for ( i=1; i<n; i++ ) { 95 v = aa + 49*ai[i]; 96 vi = aj + ai[i]; 97 nz = diag[i] - ai[i]; 98 idx = 7*(*r++); 99 sum1 = b[idx];sum2 = b[1+idx];sum3 = b[2+idx];sum4 = b[3+idx]; 100 sum5 = b[4+idx];sum6 = b[5+idx];sum7 = b[6+idx]; 101 while (nz--) { 102 idx = 7*(*vi++); 103 x1 = tmp[idx]; x2 = tmp[1+idx];x3 = tmp[2+idx]; 104 x4 = tmp[3+idx];x5 = tmp[4+idx]; 105 x6 = tmp[5+idx];x7 = tmp[6+idx]; 106 sum1 -= v[0]*x1 + v[7]*x2 + v[14]*x3 + v[21]*x4 + v[28]*x5 + v[35]*x6 + v[42]*x7; 107 sum2 -= v[1]*x1 + v[8]*x2 + v[15]*x3 + v[22]*x4 + v[29]*x5 + v[36]*x6 + v[43]*x7; 108 sum3 -= v[2]*x1 + v[9]*x2 + v[16]*x3 + v[23]*x4 + v[30]*x5 + v[37]*x6 + v[44]*x7; 109 sum4 -= v[3]*x1 + v[10]*x2 + v[17]*x3 + v[24]*x4 + v[31]*x5 + v[38]*x6 + v[45]*x7; 110 sum5 -= v[4]*x1 + v[11]*x2 + v[18]*x3 + v[25]*x4 + v[32]*x5 + v[39]*x6 + v[46]*x7; 111 sum6 -= v[5]*x1 + v[12]*x2 + v[19]*x3 + v[26]*x4 + v[33]*x5 + v[40]*x6 + v[47]*x7; 112 sum7 -= v[6]*x1 + v[13]*x2 + v[20]*x3 + v[27]*x4 + v[34]*x5 + v[41]*x6 + v[48]*x7; 113 v += 49; 114 } 115 idx = 7*i; 116 tmp[idx] = sum1;tmp[1+idx] = sum2; 117 tmp[2+idx] = sum3;tmp[3+idx] = sum4; tmp[4+idx] = sum5; 118 tmp[5+idx] = sum6;tmp[6+idx] = sum7; 119 } 120 /* backward solve the upper triangular */ 121 for ( i=n-1; i>=0; i-- ){ 122 v = aa + 49*diag[i] + 49; 123 vi = aj + diag[i] + 1; 124 nz = ai[i+1] - diag[i] - 1; 125 idt = 7*i; 126 sum1 = tmp[idt]; sum2 = tmp[1+idt]; 127 sum3 = tmp[2+idt];sum4 = tmp[3+idt]; sum5 = tmp[4+idt]; 128 sum6 = tmp[5+idt];sum7 = tmp[6+idt]; 129 while (nz--) { 130 idx = 7*(*vi++); 131 x1 = tmp[idx]; x2 = tmp[1+idx]; 132 x3 = tmp[2+idx]; x4 = tmp[3+idx]; x5 = tmp[4+idx]; 133 x6 = tmp[5+idx]; x7 = tmp[6+idx]; 134 sum1 -= v[0]*x1 + v[7]*x2 + v[14]*x3 + v[21]*x4 + v[28]*x5 + v[35]*x6 + v[42]*x7; 135 sum2 -= v[1]*x1 + v[8]*x2 + v[15]*x3 + v[22]*x4 + v[29]*x5 + v[36]*x6 + v[43]*x7; 136 sum3 -= v[2]*x1 + v[9]*x2 + v[16]*x3 + v[23]*x4 + v[30]*x5 + v[37]*x6 + v[44]*x7; 137 sum4 -= v[3]*x1 + v[10]*x2 + v[17]*x3 + v[24]*x4 + v[31]*x5 + v[38]*x6 + v[45]*x7; 138 sum5 -= v[4]*x1 + v[11]*x2 + v[18]*x3 + v[25]*x4 + v[32]*x5 + v[39]*x6 + v[46]*x7; 139 sum6 -= v[5]*x1 + v[12]*x2 + v[19]*x3 + v[26]*x4 + v[33]*x5 + v[40]*x6 + v[47]*x7; 140 sum7 -= v[6]*x1 + v[13]*x2 + v[20]*x3 + v[27]*x4 + v[34]*x5 + v[41]*x6 + v[48]*x7; 141 v += 49; 142 } 143 idc = 7*(*c--); 144 v = aa + 49*diag[i]; 145 x[idc] = tmp[idt] = v[0]*sum1+v[7]*sum2+v[14]*sum3+ 146 v[21]*sum4+v[28]*sum5+v[35]*sum6+v[42]*sum7; 147 x[1+idc] = tmp[1+idt] = v[1]*sum1+v[8]*sum2+v[15]*sum3+ 148 v[22]*sum4+v[29]*sum5+v[36]*sum6+v[43]*sum7; 149 x[2+idc] = tmp[2+idt] = v[2]*sum1+v[9]*sum2+v[16]*sum3+ 150 v[23]*sum4+v[30]*sum5+v[37]*sum6+v[44]*sum7; 151 x[3+idc] = tmp[3+idt] = v[3]*sum1+v[10]*sum2+v[17]*sum3+ 152 v[24]*sum4+v[31]*sum5+v[38]*sum6+v[45]*sum7; 153 x[4+idc] = tmp[4+idt] = v[4]*sum1+v[11]*sum2+v[18]*sum3+ 154 v[25]*sum4+v[32]*sum5+v[39]*sum6+v[46]*sum7; 155 x[5+idc] = tmp[5+idt] = v[5]*sum1+v[12]*sum2+v[19]*sum3+ 156 v[26]*sum4+v[33]*sum5+v[40]*sum6+v[47]*sum7; 157 x[6+idc] = tmp[6+idt] = v[6]*sum1+v[13]*sum2+v[20]*sum3+ 158 v[27]*sum4+v[34]*sum5+v[41]*sum6+v[48]*sum7; 159 } 160 161 ierr = ISRestoreIndices(isrow,&rout); CHKERRQ(ierr); 162 ierr = ISRestoreIndices(iscol,&cout); CHKERRQ(ierr); 163 ierr = VecRestoreArray(bb,&b); CHKERRQ(ierr); 164 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 165 PLogFlops(2*49*(a->nz) - a->n); 166 PetscFunctionReturn(0); 167 } 168 169 #undef __FUNC__ 170 #define __FUNC__ "MatSolve_SeqBAIJ_5" 171 int MatSolve_SeqBAIJ_5(Mat A,Vec bb,Vec xx) 172 { 173 Mat_SeqBAIJ *a=(Mat_SeqBAIJ *)A->data; 174 IS iscol=a->col,isrow=a->row; 175 int *r,*c,ierr,i,n=a->mbs,*vi,*ai=a->i,*aj=a->j,nz,idx,idt,idc,*rout,*cout; 176 int *diag = a->diag; 177 Scalar *aa=a->a,sum1,sum2,sum3,sum4,sum5,x1,x2,x3,x4,x5; 178 Scalar *x,*b,*tmp,*v; 179 180 PetscFunctionBegin; 181 ierr = VecGetArray(bb,&b); CHKERRQ(ierr); 182 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 183 tmp = a->solve_work; 184 185 ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout; 186 ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout + (n-1); 187 188 /* forward solve the lower triangular */ 189 idx = 5*(*r++); 190 tmp[0] = b[idx]; tmp[1] = b[1+idx]; 191 tmp[2] = b[2+idx]; tmp[3] = b[3+idx]; tmp[4] = b[4+idx]; 192 for ( i=1; i<n; i++ ) { 193 v = aa + 25*ai[i]; 194 vi = aj + ai[i]; 195 nz = diag[i] - ai[i]; 196 idx = 5*(*r++); 197 sum1 = b[idx];sum2 = b[1+idx];sum3 = b[2+idx];sum4 = b[3+idx]; 198 sum5 = b[4+idx]; 199 while (nz--) { 200 idx = 5*(*vi++); 201 x1 = tmp[idx]; x2 = tmp[1+idx];x3 = tmp[2+idx]; 202 x4 = tmp[3+idx];x5 = tmp[4+idx]; 203 sum1 -= v[0]*x1 + v[5]*x2 + v[10]*x3 + v[15]*x4 + v[20]*x5; 204 sum2 -= v[1]*x1 + v[6]*x2 + v[11]*x3 + v[16]*x4 + v[21]*x5; 205 sum3 -= v[2]*x1 + v[7]*x2 + v[12]*x3 + v[17]*x4 + v[22]*x5; 206 sum4 -= v[3]*x1 + v[8]*x2 + v[13]*x3 + v[18]*x4 + v[23]*x5; 207 sum5 -= v[4]*x1 + v[9]*x2 + v[14]*x3 + v[19]*x4 + v[24]*x5; 208 v += 25; 209 } 210 idx = 5*i; 211 tmp[idx] = sum1;tmp[1+idx] = sum2; 212 tmp[2+idx] = sum3;tmp[3+idx] = sum4; tmp[4+idx] = sum5; 213 } 214 /* backward solve the upper triangular */ 215 for ( i=n-1; i>=0; i-- ){ 216 v = aa + 25*diag[i] + 25; 217 vi = aj + diag[i] + 1; 218 nz = ai[i+1] - diag[i] - 1; 219 idt = 5*i; 220 sum1 = tmp[idt]; sum2 = tmp[1+idt]; 221 sum3 = tmp[2+idt];sum4 = tmp[3+idt]; sum5 = tmp[4+idt]; 222 while (nz--) { 223 idx = 5*(*vi++); 224 x1 = tmp[idx]; x2 = tmp[1+idx]; 225 x3 = tmp[2+idx]; x4 = tmp[3+idx]; x5 = tmp[4+idx]; 226 sum1 -= v[0]*x1 + v[5]*x2 + v[10]*x3 + v[15]*x4 + v[20]*x5; 227 sum2 -= v[1]*x1 + v[6]*x2 + v[11]*x3 + v[16]*x4 + v[21]*x5; 228 sum3 -= v[2]*x1 + v[7]*x2 + v[12]*x3 + v[17]*x4 + v[22]*x5; 229 sum4 -= v[3]*x1 + v[8]*x2 + v[13]*x3 + v[18]*x4 + v[23]*x5; 230 sum5 -= v[4]*x1 + v[9]*x2 + v[14]*x3 + v[19]*x4 + v[24]*x5; 231 v += 25; 232 } 233 idc = 5*(*c--); 234 v = aa + 25*diag[i]; 235 x[idc] = tmp[idt] = v[0]*sum1+v[5]*sum2+v[10]*sum3+ 236 v[15]*sum4+v[20]*sum5; 237 x[1+idc] = tmp[1+idt] = v[1]*sum1+v[6]*sum2+v[11]*sum3+ 238 v[16]*sum4+v[21]*sum5; 239 x[2+idc] = tmp[2+idt] = v[2]*sum1+v[7]*sum2+v[12]*sum3+ 240 v[17]*sum4+v[22]*sum5; 241 x[3+idc] = tmp[3+idt] = v[3]*sum1+v[8]*sum2+v[13]*sum3+ 242 v[18]*sum4+v[23]*sum5; 243 x[4+idc] = tmp[4+idt] = v[4]*sum1+v[9]*sum2+v[14]*sum3+ 244 v[19]*sum4+v[24]*sum5; 245 } 246 247 ierr = ISRestoreIndices(isrow,&rout); CHKERRQ(ierr); 248 ierr = ISRestoreIndices(iscol,&cout); CHKERRQ(ierr); 249 ierr = VecRestoreArray(bb,&b); CHKERRQ(ierr); 250 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 251 PLogFlops(2*25*(a->nz) - a->n); 252 PetscFunctionReturn(0); 253 } 254 255 #undef __FUNC__ 256 #define __FUNC__ "MatSolve_SeqBAIJ_4" 257 int MatSolve_SeqBAIJ_4(Mat A,Vec bb,Vec xx) 258 { 259 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data; 260 IS iscol=a->col,isrow=a->row; 261 int *r,*c,ierr,i,n=a->mbs,*vi,*ai=a->i,*aj=a->j,nz,idx,idt,idc,*rout,*cout; 262 int *diag = a->diag; 263 Scalar *aa=a->a,sum1,sum2,sum3,sum4,x1,x2,x3,x4; 264 Scalar *x,*b,*tmp,*v; 265 266 PetscFunctionBegin; 267 ierr = VecGetArray(bb,&b); CHKERRQ(ierr); 268 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 269 tmp = a->solve_work; 270 271 ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout; 272 ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout + (n-1); 273 274 /* forward solve the lower triangular */ 275 idx = 4*(*r++); 276 tmp[0] = b[idx]; tmp[1] = b[1+idx]; 277 tmp[2] = b[2+idx]; tmp[3] = b[3+idx]; 278 for ( i=1; i<n; i++ ) { 279 v = aa + 16*ai[i]; 280 vi = aj + ai[i]; 281 nz = diag[i] - ai[i]; 282 idx = 4*(*r++); 283 sum1 = b[idx];sum2 = b[1+idx];sum3 = b[2+idx];sum4 = b[3+idx]; 284 while (nz--) { 285 idx = 4*(*vi++); 286 x1 = tmp[idx];x2 = tmp[1+idx];x3 = tmp[2+idx];x4 = tmp[3+idx]; 287 sum1 -= v[0]*x1 + v[4]*x2 + v[8]*x3 + v[12]*x4; 288 sum2 -= v[1]*x1 + v[5]*x2 + v[9]*x3 + v[13]*x4; 289 sum3 -= v[2]*x1 + v[6]*x2 + v[10]*x3 + v[14]*x4; 290 sum4 -= v[3]*x1 + v[7]*x2 + v[11]*x3 + v[15]*x4; 291 v += 16; 292 } 293 idx = 4*i; 294 tmp[idx] = sum1;tmp[1+idx] = sum2; 295 tmp[2+idx] = sum3;tmp[3+idx] = sum4; 296 } 297 /* backward solve the upper triangular */ 298 for ( i=n-1; i>=0; i-- ){ 299 v = aa + 16*diag[i] + 16; 300 vi = aj + diag[i] + 1; 301 nz = ai[i+1] - diag[i] - 1; 302 idt = 4*i; 303 sum1 = tmp[idt]; sum2 = tmp[1+idt]; 304 sum3 = tmp[2+idt];sum4 = tmp[3+idt]; 305 while (nz--) { 306 idx = 4*(*vi++); 307 x1 = tmp[idx]; x2 = tmp[1+idx]; 308 x3 = tmp[2+idx]; x4 = tmp[3+idx]; 309 sum1 -= v[0]*x1 + v[4]*x2 + v[8]*x3 + v[12]*x4; 310 sum2 -= v[1]*x1 + v[5]*x2 + v[9]*x3 + v[13]*x4; 311 sum3 -= v[2]*x1 + v[6]*x2 + v[10]*x3 + v[14]*x4; 312 sum4 -= v[3]*x1 + v[7]*x2 + v[11]*x3 + v[15]*x4; 313 v += 16; 314 } 315 idc = 4*(*c--); 316 v = aa + 16*diag[i]; 317 x[idc] = tmp[idt] = v[0]*sum1+v[4]*sum2+v[8]*sum3+v[12]*sum4; 318 x[1+idc] = tmp[1+idt] = v[1]*sum1+v[5]*sum2+v[9]*sum3+v[13]*sum4; 319 x[2+idc] = tmp[2+idt] = v[2]*sum1+v[6]*sum2+v[10]*sum3+v[14]*sum4; 320 x[3+idc] = tmp[3+idt] = v[3]*sum1+v[7]*sum2+v[11]*sum3+v[15]*sum4; 321 } 322 323 ierr = ISRestoreIndices(isrow,&rout); CHKERRQ(ierr); 324 ierr = ISRestoreIndices(iscol,&cout); CHKERRQ(ierr); 325 ierr = VecRestoreArray(bb,&b); CHKERRQ(ierr); 326 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 327 PLogFlops(2*16*(a->nz) - a->n); 328 PetscFunctionReturn(0); 329 } 330 331 /* 332 Special case where the matrix was ILU(0) factored in the natural 333 ordering. This eliminates the need for the column and row permutation. 334 */ 335 #undef __FUNC__ 336 #define __FUNC__ "MatSolve_SeqBAIJ_4_NaturalOrdering" 337 int MatSolve_SeqBAIJ_4_NaturalOrdering(Mat A,Vec bb,Vec xx) 338 { 339 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data; 340 int n=a->mbs,*ai=a->i,*aj=a->j; 341 int ierr,*diag = a->diag; 342 Scalar *aa=a->a; 343 Scalar *x,*b; 344 345 PetscFunctionBegin; 346 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 347 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 348 349 #if defined(USE_FORTRAN_KERNEL_SOLVEBAIJBLAS) 350 { 351 static Scalar w[2000]; /* very BAD need to fix */ 352 fortransolvebaij4blas_(&n,x,ai,aj,diag,aa,b,w); 353 } 354 #elif defined(USE_FORTRAN_KERNEL_SOLVEBAIJ) 355 { 356 static Scalar w[2000]; /* very BAD need to fix */ 357 fortransolvebaij4_(&n,x,ai,aj,diag,aa,b,w); 358 } 359 #elif defined(USE_FORTRAN_KERNEL_SOLVEBAIJUNROLL) 360 fortransolvebaij4unroll_(&n,x,ai,aj,diag,aa,b); 361 #else 362 { 363 Scalar sum1,sum2,sum3,sum4,x1,x2,x3,x4,*v; 364 int jdx,idt,idx,nz,*vi,i; 365 366 /* forward solve the lower triangular */ 367 idx = 0; 368 x[0] = b[0]; x[1] = b[1]; x[2] = b[2]; x[3] = b[3]; 369 for ( i=1; i<n; i++ ) { 370 v = aa + 16*ai[i]; 371 vi = aj + ai[i]; 372 nz = diag[i] - ai[i]; 373 idx += 4; 374 sum1 = b[idx];sum2 = b[1+idx];sum3 = b[2+idx];sum4 = b[3+idx]; 375 while (nz--) { 376 jdx = 4*(*vi++); 377 x1 = x[jdx];x2 = x[1+jdx];x3 = x[2+jdx];x4 = x[3+jdx]; 378 sum1 -= v[0]*x1 + v[4]*x2 + v[8]*x3 + v[12]*x4; 379 sum2 -= v[1]*x1 + v[5]*x2 + v[9]*x3 + v[13]*x4; 380 sum3 -= v[2]*x1 + v[6]*x2 + v[10]*x3 + v[14]*x4; 381 sum4 -= v[3]*x1 + v[7]*x2 + v[11]*x3 + v[15]*x4; 382 v += 16; 383 } 384 x[idx] = sum1; 385 x[1+idx] = sum2; 386 x[2+idx] = sum3; 387 x[3+idx] = sum4; 388 } 389 /* backward solve the upper triangular */ 390 for ( i=n-1; i>=0; i-- ){ 391 v = aa + 16*diag[i] + 16; 392 vi = aj + diag[i] + 1; 393 nz = ai[i+1] - diag[i] - 1; 394 idt = 4*i; 395 sum1 = x[idt]; sum2 = x[1+idt]; 396 sum3 = x[2+idt];sum4 = x[3+idt]; 397 while (nz--) { 398 idx = 4*(*vi++); 399 x1 = x[idx]; x2 = x[1+idx];x3 = x[2+idx]; x4 = x[3+idx]; 400 sum1 -= v[0]*x1 + v[4]*x2 + v[8]*x3 + v[12]*x4; 401 sum2 -= v[1]*x1 + v[5]*x2 + v[9]*x3 + v[13]*x4; 402 sum3 -= v[2]*x1 + v[6]*x2 + v[10]*x3 + v[14]*x4; 403 sum4 -= v[3]*x1 + v[7]*x2 + v[11]*x3 + v[15]*x4; 404 v += 16; 405 } 406 v = aa + 16*diag[i]; 407 x[idt] = v[0]*sum1 + v[4]*sum2 + v[8]*sum3 + v[12]*sum4; 408 x[1+idt] = v[1]*sum1 + v[5]*sum2 + v[9]*sum3 + v[13]*sum4; 409 x[2+idt] = v[2]*sum1 + v[6]*sum2 + v[10]*sum3 + v[14]*sum4; 410 x[3+idt] = v[3]*sum1 + v[7]*sum2 + v[11]*sum3 + v[15]*sum4; 411 } 412 } 413 #endif 414 415 ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr); 416 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 417 PLogFlops(2*16*(a->nz) - a->n); 418 PetscFunctionReturn(0); 419 } 420 421 422 #undef __FUNC__ 423 #define __FUNC__ "MatSolve_SeqBAIJ_3" 424 int MatSolve_SeqBAIJ_3(Mat A,Vec bb,Vec xx) 425 { 426 Mat_SeqBAIJ *a=(Mat_SeqBAIJ *)A->data; 427 IS iscol=a->col,isrow=a->row; 428 int *r,*c,ierr,i,n=a->mbs,*vi,*ai=a->i,*aj=a->j,nz,idx,idt,idc,*rout,*cout; 429 int *diag = a->diag; 430 Scalar *aa=a->a,sum1,sum2,sum3,x1,x2,x3; 431 Scalar *x,*b,*tmp,*v; 432 433 PetscFunctionBegin; 434 ierr = VecGetArray(bb,&b); CHKERRQ(ierr); 435 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 436 tmp = a->solve_work; 437 438 ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout; 439 ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout + (n-1); 440 441 /* forward solve the lower triangular */ 442 idx = 3*(*r++); 443 tmp[0] = b[idx]; tmp[1] = b[1+idx]; tmp[2] = b[2+idx]; 444 for ( i=1; i<n; i++ ) { 445 v = aa + 9*ai[i]; 446 vi = aj + ai[i]; 447 nz = diag[i] - ai[i]; 448 idx = 3*(*r++); 449 sum1 = b[idx]; sum2 = b[1+idx]; sum3 = b[2+idx]; 450 while (nz--) { 451 idx = 3*(*vi++); 452 x1 = tmp[idx]; x2 = tmp[1+idx]; x3 = tmp[2+idx]; 453 sum1 -= v[0]*x1 + v[3]*x2 + v[6]*x3; 454 sum2 -= v[1]*x1 + v[4]*x2 + v[7]*x3; 455 sum3 -= v[2]*x1 + v[5]*x2 + v[8]*x3; 456 v += 9; 457 } 458 idx = 3*i; 459 tmp[idx] = sum1; tmp[1+idx] = sum2; tmp[2+idx] = sum3; 460 } 461 /* backward solve the upper triangular */ 462 for ( i=n-1; i>=0; i-- ){ 463 v = aa + 9*diag[i] + 9; 464 vi = aj + diag[i] + 1; 465 nz = ai[i+1] - diag[i] - 1; 466 idt = 3*i; 467 sum1 = tmp[idt]; sum2 = tmp[1+idt]; sum3 = tmp[2+idt]; 468 while (nz--) { 469 idx = 3*(*vi++); 470 x1 = tmp[idx]; x2 = tmp[1+idx]; x3 = tmp[2+idx]; 471 sum1 -= v[0]*x1 + v[3]*x2 + v[6]*x3; 472 sum2 -= v[1]*x1 + v[4]*x2 + v[7]*x3; 473 sum3 -= v[2]*x1 + v[5]*x2 + v[8]*x3; 474 v += 9; 475 } 476 idc = 3*(*c--); 477 v = aa + 9*diag[i]; 478 x[idc] = tmp[idt] = v[0]*sum1 + v[3]*sum2 + v[6]*sum3; 479 x[1+idc] = tmp[1+idt] = v[1]*sum1 + v[4]*sum2 + v[7]*sum3; 480 x[2+idc] = tmp[2+idt] = v[2]*sum1 + v[5]*sum2 + v[8]*sum3; 481 } 482 ierr = ISRestoreIndices(isrow,&rout); CHKERRQ(ierr); 483 ierr = ISRestoreIndices(iscol,&cout); CHKERRQ(ierr); 484 ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr); 485 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 486 PLogFlops(2*9*(a->nz) - a->n); 487 PetscFunctionReturn(0); 488 } 489 490 #undef __FUNC__ 491 #define __FUNC__ "MatSolve_SeqBAIJ_2" 492 int MatSolve_SeqBAIJ_2(Mat A,Vec bb,Vec xx) 493 { 494 Mat_SeqBAIJ *a=(Mat_SeqBAIJ *)A->data; 495 IS iscol=a->col,isrow=a->row; 496 int *r,*c,ierr,i,n=a->mbs,*vi,*ai=a->i,*aj=a->j,nz,idx,idt,idc,*rout,*cout; 497 int *diag = a->diag; 498 Scalar *aa=a->a,sum1,sum2,x1,x2; 499 Scalar *x,*b,*tmp,*v; 500 501 PetscFunctionBegin; 502 ierr = VecGetArray(bb,&b); CHKERRQ(ierr); 503 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 504 tmp = a->solve_work; 505 506 ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout; 507 ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout + (n-1); 508 509 /* forward solve the lower triangular */ 510 idx = 2*(*r++); 511 tmp[0] = b[idx]; tmp[1] = b[1+idx]; 512 for ( i=1; i<n; i++ ) { 513 v = aa + 4*ai[i]; 514 vi = aj + ai[i]; 515 nz = diag[i] - ai[i]; 516 idx = 2*(*r++); 517 sum1 = b[idx]; sum2 = b[1+idx]; 518 while (nz--) { 519 idx = 2*(*vi++); 520 x1 = tmp[idx]; x2 = tmp[1+idx]; 521 sum1 -= v[0]*x1 + v[2]*x2; 522 sum2 -= v[1]*x1 + v[3]*x2; 523 v += 4; 524 } 525 idx = 2*i; 526 tmp[idx] = sum1; tmp[1+idx] = sum2; 527 } 528 /* backward solve the upper triangular */ 529 for ( i=n-1; i>=0; i-- ){ 530 v = aa + 4*diag[i] + 4; 531 vi = aj + diag[i] + 1; 532 nz = ai[i+1] - diag[i] - 1; 533 idt = 2*i; 534 sum1 = tmp[idt]; sum2 = tmp[1+idt]; 535 while (nz--) { 536 idx = 2*(*vi++); 537 x1 = tmp[idx]; x2 = tmp[1+idx]; 538 sum1 -= v[0]*x1 + v[2]*x2; 539 sum2 -= v[1]*x1 + v[3]*x2; 540 v += 4; 541 } 542 idc = 2*(*c--); 543 v = aa + 4*diag[i]; 544 x[idc] = tmp[idt] = v[0]*sum1 + v[2]*sum2; 545 x[1+idc] = tmp[1+idt] = v[1]*sum1 + v[3]*sum2; 546 } 547 ierr = ISRestoreIndices(isrow,&rout); CHKERRQ(ierr); 548 ierr = ISRestoreIndices(iscol,&cout); CHKERRQ(ierr); 549 ierr = VecRestoreArray(bb,&b); CHKERRQ(ierr); 550 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 551 PLogFlops(2*4*(a->nz) - a->n); 552 PetscFunctionReturn(0); 553 } 554 555 556 #undef __FUNC__ 557 #define __FUNC__ "MatSolve_SeqBAIJ_1" 558 int MatSolve_SeqBAIJ_1(Mat A,Vec bb,Vec xx) 559 { 560 Mat_SeqBAIJ *a=(Mat_SeqBAIJ *)A->data; 561 IS iscol=a->col,isrow=a->row; 562 int *r,*c,ierr,i,n=a->mbs,*vi,*ai=a->i,*aj=a->j,nz,*rout,*cout; 563 int *diag = a->diag; 564 Scalar *aa=a->a,sum1; 565 Scalar *x,*b,*tmp,*v; 566 567 PetscFunctionBegin; 568 if (!n) PetscFunctionReturn(0); 569 570 ierr = VecGetArray(bb,&b); CHKERRQ(ierr); 571 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 572 tmp = a->solve_work; 573 574 ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout; 575 ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout + (n-1); 576 577 /* forward solve the lower triangular */ 578 tmp[0] = b[*r++]; 579 for ( i=1; i<n; i++ ) { 580 v = aa + ai[i]; 581 vi = aj + ai[i]; 582 nz = diag[i] - ai[i]; 583 sum1 = b[*r++]; 584 while (nz--) { 585 sum1 -= (*v++)*tmp[*vi++]; 586 } 587 tmp[i] = sum1; 588 } 589 /* backward solve the upper triangular */ 590 for ( i=n-1; i>=0; i-- ){ 591 v = aa + diag[i] + 1; 592 vi = aj + diag[i] + 1; 593 nz = ai[i+1] - diag[i] - 1; 594 sum1 = tmp[i]; 595 while (nz--) { 596 sum1 -= (*v++)*tmp[*vi++]; 597 } 598 x[*c--] = tmp[i] = aa[diag[i]]*sum1; 599 } 600 601 ierr = ISRestoreIndices(isrow,&rout); CHKERRQ(ierr); 602 ierr = ISRestoreIndices(iscol,&cout); CHKERRQ(ierr); 603 ierr = VecRestoreArray(bb,&b); CHKERRQ(ierr); 604 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 605 PLogFlops(2*1*(a->nz) - a->n); 606 PetscFunctionReturn(0); 607 } 608 609 extern int MatLUFactorNumeric_SeqBAIJ_4_NaturalOrdering(Mat,Mat*); 610 extern int MatSolve_SeqBAIJ_4_NaturalOrdering(Mat,Vec,Vec); 611 /* ----------------------------------------------------------------*/ 612 /* 613 This code is virtually identical to MatILUFactorSymbolic_SeqAIJ 614 except that the data structure of Mat_SeqAIJ is slightly different. 615 Not a good example of code reuse. 616 */ 617 #undef __FUNC__ 618 #define __FUNC__ "MatILUFactorSymbolic_SeqBAIJ" 619 int MatILUFactorSymbolic_SeqBAIJ(Mat A,IS isrow,IS iscol,double f,int levels, 620 Mat *fact) 621 { 622 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *) A->data, *b; 623 IS isicol; 624 int *r,*ic, ierr, prow, n = a->mbs, *ai = a->i, *aj = a->j; 625 int *ainew,*ajnew, jmax,*fill, *xi, nz, *im,*ajfill,*flev; 626 int *dloc, idx, row,m,fm, nzf, nzi,len, realloc = 0; 627 int incrlev,nnz,i,bs = a->bs,bs2 = a->bs2; 628 PetscTruth col_identity, row_identity; 629 630 PetscFunctionBegin; 631 ierr = ISInvertPermutation(iscol,&isicol); CHKERRQ(ierr); 632 633 /* special case that simply copies fill pattern */ 634 PetscValidHeaderSpecific(isrow,IS_COOKIE); 635 PetscValidHeaderSpecific(iscol,IS_COOKIE); 636 ISIdentity(isrow,&row_identity); ISIdentity(iscol,&col_identity); 637 if (levels == 0 && row_identity && col_identity) { 638 ierr = MatDuplicate_SeqBAIJ(A,MAT_DO_NOT_COPY_VALUES,fact); CHKERRQ(ierr); 639 (*fact)->factor = FACTOR_LU; 640 b = (Mat_SeqBAIJ *) (*fact)->data; 641 if (!b->diag) { 642 ierr = MatMarkDiag_SeqBAIJ(*fact); CHKERRQ(ierr); 643 } 644 b->row = isrow; 645 b->col = iscol; 646 b->icol = isicol; 647 b->solve_work = (Scalar *) PetscMalloc((b->m+1+b->bs)*sizeof(Scalar));CHKPTRQ(b->solve_work); 648 /* 649 Blocksize 4 has a special faster solver for ILU(0) factorization 650 with natural ordering 651 */ 652 if (b->bs == 4) { 653 (*fact)->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_4_NaturalOrdering; 654 (*fact)->ops->solve = MatSolve_SeqBAIJ_4_NaturalOrdering; 655 } 656 PetscFunctionReturn(0); 657 } 658 659 ierr = ISGetIndices(isrow,&r); CHKERRQ(ierr); 660 ierr = ISGetIndices(isicol,&ic); CHKERRQ(ierr); 661 662 /* get new row pointers */ 663 ainew = (int *) PetscMalloc( (n+1)*sizeof(int) ); CHKPTRQ(ainew); 664 ainew[0] = 0; 665 /* don't know how many column pointers are needed so estimate */ 666 jmax = (int) (f*ai[n] + 1); 667 ajnew = (int *) PetscMalloc( (jmax)*sizeof(int) ); CHKPTRQ(ajnew); 668 /* ajfill is level of fill for each fill entry */ 669 ajfill = (int *) PetscMalloc( (jmax)*sizeof(int) ); CHKPTRQ(ajfill); 670 /* fill is a linked list of nonzeros in active row */ 671 fill = (int *) PetscMalloc( (n+1)*sizeof(int)); CHKPTRQ(fill); 672 /* im is level for each filled value */ 673 im = (int *) PetscMalloc( (n+1)*sizeof(int)); CHKPTRQ(im); 674 /* dloc is location of diagonal in factor */ 675 dloc = (int *) PetscMalloc( (n+1)*sizeof(int)); CHKPTRQ(dloc); 676 dloc[0] = 0; 677 for ( prow=0; prow<n; prow++ ) { 678 /* first copy previous fill into linked list */ 679 nzf = nz = ai[r[prow]+1] - ai[r[prow]]; 680 if (!nz) SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,1,"Empty row in matrix"); 681 xi = aj + ai[r[prow]]; 682 fill[n] = n; 683 while (nz--) { 684 fm = n; 685 idx = ic[*xi++]; 686 do { 687 m = fm; 688 fm = fill[m]; 689 } while (fm < idx); 690 fill[m] = idx; 691 fill[idx] = fm; 692 im[idx] = 0; 693 } 694 nzi = 0; 695 row = fill[n]; 696 while ( row < prow ) { 697 incrlev = im[row] + 1; 698 nz = dloc[row]; 699 xi = ajnew + ainew[row] + nz; 700 flev = ajfill + ainew[row] + nz + 1; 701 nnz = ainew[row+1] - ainew[row] - nz - 1; 702 if (*xi++ != row) { 703 SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,0,"Zero pivot: try running with -pc_ilu_nonzeros_along_diagonal"); 704 } 705 fm = row; 706 while (nnz-- > 0) { 707 idx = *xi++; 708 if (*flev + incrlev > levels) { 709 flev++; 710 continue; 711 } 712 do { 713 m = fm; 714 fm = fill[m]; 715 } while (fm < idx); 716 if (fm != idx) { 717 im[idx] = *flev + incrlev; 718 fill[m] = idx; 719 fill[idx] = fm; 720 fm = idx; 721 nzf++; 722 } else { 723 if (im[idx] > *flev + incrlev) im[idx] = *flev+incrlev; 724 } 725 flev++; 726 } 727 row = fill[row]; 728 nzi++; 729 } 730 /* copy new filled row into permanent storage */ 731 ainew[prow+1] = ainew[prow] + nzf; 732 if (ainew[prow+1] > jmax) { 733 734 /* estimate how much additional space we will need */ 735 /* use the strategy suggested by David Hysom <hysom@perch-t.icase.edu> */ 736 /* just double the memory each time */ 737 int maxadd = jmax; 738 /* maxadd = (int) (((f*ai[n]+1)*(n-prow+5))/n); */ 739 if (maxadd < nzf) maxadd = (n-prow)*(nzf+1); 740 jmax += maxadd; 741 742 /* allocate a longer ajnew and ajfill */ 743 xi = (int *) PetscMalloc( jmax*sizeof(int) );CHKPTRQ(xi); 744 PetscMemcpy(xi,ajnew,ainew[prow]*sizeof(int)); 745 PetscFree(ajnew); 746 ajnew = xi; 747 xi = (int *) PetscMalloc( jmax*sizeof(int) );CHKPTRQ(xi); 748 PetscMemcpy(xi,ajfill,ainew[prow]*sizeof(int)); 749 PetscFree(ajfill); 750 ajfill = xi; 751 realloc++; /* count how many reallocations are needed */ 752 } 753 xi = ajnew + ainew[prow]; 754 flev = ajfill + ainew[prow]; 755 dloc[prow] = nzi; 756 fm = fill[n]; 757 while (nzf--) { 758 *xi++ = fm; 759 *flev++ = im[fm]; 760 fm = fill[fm]; 761 } 762 } 763 PetscFree(ajfill); 764 ierr = ISRestoreIndices(isrow,&r); CHKERRQ(ierr); 765 ierr = ISRestoreIndices(isicol,&ic); CHKERRQ(ierr); 766 PetscFree(fill); PetscFree(im); 767 768 { 769 double af = ((double)ainew[n])/((double)ai[n]); 770 PLogInfo(A,"MatILUFactorSymbolic_SeqBAIJ:Reallocs %d Fill ratio:given %g needed %g\n", 771 realloc,f,af); 772 PLogInfo(A,"MatILUFactorSymbolic_SeqBAIJ:Run with -pc_ilu_fill %g or use \n",af); 773 PLogInfo(A,"MatILUFactorSymbolic_SeqBAIJ:PCILUSetFill(pc,%g);\n",af); 774 PLogInfo(A,"MatILUFactorSymbolic_SeqBAIJ:for best performance.\n"); 775 } 776 777 /* put together the new matrix */ 778 ierr = MatCreateSeqBAIJ(A->comm,bs,bs*n,bs*n,0,PETSC_NULL,fact);CHKERRQ(ierr); 779 PLogObjectParent(*fact,isicol); 780 b = (Mat_SeqBAIJ *) (*fact)->data; 781 PetscFree(b->imax); 782 b->singlemalloc = 0; 783 len = bs2*ainew[n]*sizeof(Scalar); 784 /* the next line frees the default space generated by the Create() */ 785 PetscFree(b->a); PetscFree(b->ilen); 786 b->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(b->a); 787 b->j = ajnew; 788 b->i = ainew; 789 for ( i=0; i<n; i++ ) dloc[i] += ainew[i]; 790 b->diag = dloc; 791 b->ilen = 0; 792 b->imax = 0; 793 b->row = isrow; 794 b->col = iscol; 795 b->icol = isicol; 796 b->solve_work = (Scalar *) PetscMalloc( (bs*n+bs)*sizeof(Scalar));CHKPTRQ(b->solve_work); 797 /* In b structure: Free imax, ilen, old a, old j. 798 Allocate dloc, solve_work, new a, new j */ 799 PLogObjectMemory(*fact,(ainew[n]-n)*(sizeof(int))+bs2*ainew[n]*sizeof(Scalar)); 800 b->maxnz = b->nz = ainew[n]; 801 (*fact)->factor = FACTOR_LU; 802 803 (*fact)->info.factor_mallocs = realloc; 804 (*fact)->info.fill_ratio_given = f; 805 (*fact)->info.fill_ratio_needed = ((double)ainew[n])/((double)ai[prow]); 806 807 PetscFunctionReturn(0); 808 } 809 810 811 812 813