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