1 2 /* 3 Factorization code for BAIJ format. 4 */ 5 #include <../src/mat/impls/baij/seq/baij.h> 6 #include <petsc/private/kernels/blockinvert.h> 7 8 PetscErrorCode MatLUFactorNumeric_SeqBAIJ_2(Mat B,Mat A,const MatFactorInfo *info) 9 { 10 Mat C =B; 11 Mat_SeqBAIJ *a =(Mat_SeqBAIJ*)A->data,*b=(Mat_SeqBAIJ*)C->data; 12 IS isrow = b->row,isicol = b->icol; 13 PetscErrorCode ierr; 14 const PetscInt *r,*ic; 15 PetscInt i,j,k,nz,nzL,row,*pj; 16 const PetscInt n=a->mbs,*ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,bs2=a->bs2; 17 const PetscInt *ajtmp,*bjtmp,*bdiag=b->diag; 18 MatScalar *rtmp,*pc,*mwork,*pv; 19 MatScalar *aa=a->a,*v; 20 PetscInt flg; 21 PetscReal shift = info->shiftamount; 22 PetscBool allowzeropivot,zeropivotdetected; 23 24 PetscFunctionBegin; 25 ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr); 26 ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr); 27 allowzeropivot = PetscNot(A->erroriffailure); 28 29 /* generate work space needed by the factorization */ 30 ierr = PetscMalloc2(bs2*n,&rtmp,bs2,&mwork);CHKERRQ(ierr); 31 ierr = PetscArrayzero(rtmp,bs2*n);CHKERRQ(ierr); 32 33 for (i=0; i<n; i++) { 34 /* zero rtmp */ 35 /* L part */ 36 nz = bi[i+1] - bi[i]; 37 bjtmp = bj + bi[i]; 38 for (j=0; j<nz; j++) { 39 ierr = PetscArrayzero(rtmp+bs2*bjtmp[j],bs2);CHKERRQ(ierr); 40 } 41 42 /* U part */ 43 nz = bdiag[i] - bdiag[i+1]; 44 bjtmp = bj + bdiag[i+1]+1; 45 for (j=0; j<nz; j++) { 46 ierr = PetscArrayzero(rtmp+bs2*bjtmp[j],bs2);CHKERRQ(ierr); 47 } 48 49 /* load in initial (unfactored row) */ 50 nz = ai[r[i]+1] - ai[r[i]]; 51 ajtmp = aj + ai[r[i]]; 52 v = aa + bs2*ai[r[i]]; 53 for (j=0; j<nz; j++) { 54 ierr = PetscArraycpy(rtmp+bs2*ic[ajtmp[j]],v+bs2*j,bs2);CHKERRQ(ierr); 55 } 56 57 /* elimination */ 58 bjtmp = bj + bi[i]; 59 nzL = bi[i+1] - bi[i]; 60 for (k=0; k < nzL; k++) { 61 row = bjtmp[k]; 62 pc = rtmp + bs2*row; 63 for (flg=0,j=0; j<bs2; j++) { 64 if (pc[j] != (PetscScalar)0.0) { 65 flg = 1; 66 break; 67 } 68 } 69 if (flg) { 70 pv = b->a + bs2*bdiag[row]; 71 /* PetscKernel_A_gets_A_times_B(bs,pc,pv,mwork); *pc = *pc * (*pv); */ 72 ierr = PetscKernel_A_gets_A_times_B_2(pc,pv,mwork);CHKERRQ(ierr); 73 74 pj = b->j + bdiag[row+1]+1; /* begining of U(row,:) */ 75 pv = b->a + bs2*(bdiag[row+1]+1); 76 nz = bdiag[row] - bdiag[row+1] - 1; /* num of entries inU(row,:), excluding diag */ 77 for (j=0; j<nz; j++) { 78 /* PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j); */ 79 /* rtmp+bs2*pj[j] = rtmp+bs2*pj[j] - (*pc)*(pv+bs2*j) */ 80 v = rtmp + 4*pj[j]; 81 ierr = PetscKernel_A_gets_A_minus_B_times_C_2(v,pc,pv);CHKERRQ(ierr); 82 pv += 4; 83 } 84 ierr = PetscLogFlops(16.0*nz+12);CHKERRQ(ierr); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */ 85 } 86 } 87 88 /* finished row so stick it into b->a */ 89 /* L part */ 90 pv = b->a + bs2*bi[i]; 91 pj = b->j + bi[i]; 92 nz = bi[i+1] - bi[i]; 93 for (j=0; j<nz; j++) { 94 ierr = PetscArraycpy(pv+bs2*j,rtmp+bs2*pj[j],bs2);CHKERRQ(ierr); 95 } 96 97 /* Mark diagonal and invert diagonal for simplier triangular solves */ 98 pv = b->a + bs2*bdiag[i]; 99 pj = b->j + bdiag[i]; 100 ierr = PetscArraycpy(pv,rtmp+bs2*pj[0],bs2);CHKERRQ(ierr); 101 ierr = PetscKernel_A_gets_inverse_A_2(pv,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 102 if (zeropivotdetected) B->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 103 104 /* U part */ 105 pv = b->a + bs2*(bdiag[i+1]+1); 106 pj = b->j + bdiag[i+1]+1; 107 nz = bdiag[i] - bdiag[i+1] - 1; 108 for (j=0; j<nz; j++) { 109 ierr = PetscArraycpy(pv+bs2*j,rtmp+bs2*pj[j],bs2);CHKERRQ(ierr); 110 } 111 } 112 113 ierr = PetscFree2(rtmp,mwork);CHKERRQ(ierr); 114 ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr); 115 ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr); 116 117 C->ops->solve = MatSolve_SeqBAIJ_2; 118 C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2; 119 C->assembled = PETSC_TRUE; 120 121 ierr = PetscLogFlops(1.333333333333*2*2*2*n);CHKERRQ(ierr); /* from inverting diagonal blocks */ 122 PetscFunctionReturn(0); 123 } 124 125 PetscErrorCode MatLUFactorNumeric_SeqBAIJ_2_NaturalOrdering(Mat B,Mat A,const MatFactorInfo *info) 126 { 127 Mat C =B; 128 Mat_SeqBAIJ *a=(Mat_SeqBAIJ*)A->data,*b=(Mat_SeqBAIJ*)C->data; 129 PetscErrorCode ierr; 130 PetscInt i,j,k,nz,nzL,row,*pj; 131 const PetscInt n=a->mbs,*ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,bs2=a->bs2; 132 const PetscInt *ajtmp,*bjtmp,*bdiag=b->diag; 133 MatScalar *rtmp,*pc,*mwork,*pv; 134 MatScalar *aa=a->a,*v; 135 PetscInt flg; 136 PetscReal shift = info->shiftamount; 137 PetscBool allowzeropivot,zeropivotdetected; 138 139 PetscFunctionBegin; 140 allowzeropivot = PetscNot(A->erroriffailure); 141 142 /* generate work space needed by the factorization */ 143 ierr = PetscMalloc2(bs2*n,&rtmp,bs2,&mwork);CHKERRQ(ierr); 144 ierr = PetscArrayzero(rtmp,bs2*n);CHKERRQ(ierr); 145 146 for (i=0; i<n; i++) { 147 /* zero rtmp */ 148 /* L part */ 149 nz = bi[i+1] - bi[i]; 150 bjtmp = bj + bi[i]; 151 for (j=0; j<nz; j++) { 152 ierr = PetscArrayzero(rtmp+bs2*bjtmp[j],bs2);CHKERRQ(ierr); 153 } 154 155 /* U part */ 156 nz = bdiag[i] - bdiag[i+1]; 157 bjtmp = bj + bdiag[i+1]+1; 158 for (j=0; j<nz; j++) { 159 ierr = PetscArrayzero(rtmp+bs2*bjtmp[j],bs2);CHKERRQ(ierr); 160 } 161 162 /* load in initial (unfactored row) */ 163 nz = ai[i+1] - ai[i]; 164 ajtmp = aj + ai[i]; 165 v = aa + bs2*ai[i]; 166 for (j=0; j<nz; j++) { 167 ierr = PetscArraycpy(rtmp+bs2*ajtmp[j],v+bs2*j,bs2);CHKERRQ(ierr); 168 } 169 170 /* elimination */ 171 bjtmp = bj + bi[i]; 172 nzL = bi[i+1] - bi[i]; 173 for (k=0; k < nzL; k++) { 174 row = bjtmp[k]; 175 pc = rtmp + bs2*row; 176 for (flg=0,j=0; j<bs2; j++) { 177 if (pc[j]!=(PetscScalar)0.0) { 178 flg = 1; 179 break; 180 } 181 } 182 if (flg) { 183 pv = b->a + bs2*bdiag[row]; 184 /* PetscKernel_A_gets_A_times_B(bs,pc,pv,mwork); *pc = *pc * (*pv); */ 185 ierr = PetscKernel_A_gets_A_times_B_2(pc,pv,mwork);CHKERRQ(ierr); 186 187 pj = b->j + bdiag[row+1]+1; /* beginning of U(row,:) */ 188 pv = b->a + bs2*(bdiag[row+1]+1); 189 nz = bdiag[row]-bdiag[row+1] - 1; /* num of entries in U(row,:) excluding diag */ 190 for (j=0; j<nz; j++) { 191 /* PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j); */ 192 /* rtmp+bs2*pj[j] = rtmp+bs2*pj[j] - (*pc)*(pv+bs2*j) */ 193 v = rtmp + 4*pj[j]; 194 ierr = PetscKernel_A_gets_A_minus_B_times_C_2(v,pc,pv);CHKERRQ(ierr); 195 pv += 4; 196 } 197 ierr = PetscLogFlops(16.0*nz+12);CHKERRQ(ierr); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */ 198 } 199 } 200 201 /* finished row so stick it into b->a */ 202 /* L part */ 203 pv = b->a + bs2*bi[i]; 204 pj = b->j + bi[i]; 205 nz = bi[i+1] - bi[i]; 206 for (j=0; j<nz; j++) { 207 ierr = PetscArraycpy(pv+bs2*j,rtmp+bs2*pj[j],bs2);CHKERRQ(ierr); 208 } 209 210 /* Mark diagonal and invert diagonal for simplier triangular solves */ 211 pv = b->a + bs2*bdiag[i]; 212 pj = b->j + bdiag[i]; 213 ierr = PetscArraycpy(pv,rtmp+bs2*pj[0],bs2);CHKERRQ(ierr); 214 ierr = PetscKernel_A_gets_inverse_A_2(pv,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 215 if (zeropivotdetected) B->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 216 217 /* U part */ 218 /* 219 pv = b->a + bs2*bi[2*n-i]; 220 pj = b->j + bi[2*n-i]; 221 nz = bi[2*n-i+1] - bi[2*n-i] - 1; 222 */ 223 pv = b->a + bs2*(bdiag[i+1]+1); 224 pj = b->j + bdiag[i+1]+1; 225 nz = bdiag[i] - bdiag[i+1] - 1; 226 for (j=0; j<nz; j++) { 227 ierr = PetscArraycpy(pv+bs2*j,rtmp+bs2*pj[j],bs2);CHKERRQ(ierr); 228 } 229 } 230 ierr = PetscFree2(rtmp,mwork);CHKERRQ(ierr); 231 232 C->ops->solve = MatSolve_SeqBAIJ_2_NaturalOrdering; 233 C->ops->forwardsolve = MatForwardSolve_SeqBAIJ_2_NaturalOrdering; 234 C->ops->backwardsolve = MatBackwardSolve_SeqBAIJ_2_NaturalOrdering; 235 C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2_NaturalOrdering; 236 C->assembled = PETSC_TRUE; 237 238 ierr = PetscLogFlops(1.333333333333*2*2*2*n);CHKERRQ(ierr); /* from inverting diagonal blocks */ 239 PetscFunctionReturn(0); 240 } 241 242 PetscErrorCode MatLUFactorNumeric_SeqBAIJ_2_inplace(Mat B,Mat A,const MatFactorInfo *info) 243 { 244 Mat C = B; 245 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b = (Mat_SeqBAIJ*)C->data; 246 IS isrow = b->row,isicol = b->icol; 247 PetscErrorCode ierr; 248 const PetscInt *r,*ic; 249 PetscInt i,j,n = a->mbs,*bi = b->i,*bj = b->j; 250 PetscInt *ajtmpold,*ajtmp,nz,row; 251 PetscInt *diag_offset=b->diag,idx,*ai=a->i,*aj=a->j,*pj; 252 MatScalar *pv,*v,*rtmp,m1,m2,m3,m4,*pc,*w,*x,x1,x2,x3,x4; 253 MatScalar p1,p2,p3,p4; 254 MatScalar *ba = b->a,*aa = a->a; 255 PetscReal shift = info->shiftamount; 256 PetscBool allowzeropivot,zeropivotdetected; 257 258 PetscFunctionBegin; 259 allowzeropivot = PetscNot(A->erroriffailure); 260 ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr); 261 ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr); 262 ierr = PetscMalloc1(4*(n+1),&rtmp);CHKERRQ(ierr); 263 264 for (i=0; i<n; i++) { 265 nz = bi[i+1] - bi[i]; 266 ajtmp = bj + bi[i]; 267 for (j=0; j<nz; j++) { 268 x = rtmp+4*ajtmp[j]; x[0] = x[1] = x[2] = x[3] = 0.0; 269 } 270 /* load in initial (unfactored row) */ 271 idx = r[i]; 272 nz = ai[idx+1] - ai[idx]; 273 ajtmpold = aj + ai[idx]; 274 v = aa + 4*ai[idx]; 275 for (j=0; j<nz; j++) { 276 x = rtmp+4*ic[ajtmpold[j]]; 277 x[0] = v[0]; x[1] = v[1]; x[2] = v[2]; x[3] = v[3]; 278 v += 4; 279 } 280 row = *ajtmp++; 281 while (row < i) { 282 pc = rtmp + 4*row; 283 p1 = pc[0]; p2 = pc[1]; p3 = pc[2]; p4 = pc[3]; 284 if (p1 != (PetscScalar)0.0 || p2 != (PetscScalar)0.0 || p3 != (PetscScalar)0.0 || p4 != (PetscScalar)0.0) { 285 pv = ba + 4*diag_offset[row]; 286 pj = bj + diag_offset[row] + 1; 287 x1 = pv[0]; x2 = pv[1]; x3 = pv[2]; x4 = pv[3]; 288 pc[0] = m1 = p1*x1 + p3*x2; 289 pc[1] = m2 = p2*x1 + p4*x2; 290 pc[2] = m3 = p1*x3 + p3*x4; 291 pc[3] = m4 = p2*x3 + p4*x4; 292 nz = bi[row+1] - diag_offset[row] - 1; 293 pv += 4; 294 for (j=0; j<nz; j++) { 295 x1 = pv[0]; x2 = pv[1]; x3 = pv[2]; x4 = pv[3]; 296 x = rtmp + 4*pj[j]; 297 x[0] -= m1*x1 + m3*x2; 298 x[1] -= m2*x1 + m4*x2; 299 x[2] -= m1*x3 + m3*x4; 300 x[3] -= m2*x3 + m4*x4; 301 pv += 4; 302 } 303 ierr = PetscLogFlops(16.0*nz+12.0);CHKERRQ(ierr); 304 } 305 row = *ajtmp++; 306 } 307 /* finished row so stick it into b->a */ 308 pv = ba + 4*bi[i]; 309 pj = bj + bi[i]; 310 nz = bi[i+1] - bi[i]; 311 for (j=0; j<nz; j++) { 312 x = rtmp+4*pj[j]; 313 pv[0] = x[0]; pv[1] = x[1]; pv[2] = x[2]; pv[3] = x[3]; 314 pv += 4; 315 } 316 /* invert diagonal block */ 317 w = ba + 4*diag_offset[i]; 318 ierr = PetscKernel_A_gets_inverse_A_2(w,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 319 if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 320 } 321 322 ierr = PetscFree(rtmp);CHKERRQ(ierr); 323 ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr); 324 ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr); 325 326 C->ops->solve = MatSolve_SeqBAIJ_2_inplace; 327 C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2_inplace; 328 C->assembled = PETSC_TRUE; 329 330 ierr = PetscLogFlops(1.333333333333*8*b->mbs);CHKERRQ(ierr); /* from inverting diagonal blocks */ 331 PetscFunctionReturn(0); 332 } 333 /* 334 Version for when blocks are 2 by 2 Using natural ordering 335 */ 336 PetscErrorCode MatLUFactorNumeric_SeqBAIJ_2_NaturalOrdering_inplace(Mat C,Mat A,const MatFactorInfo *info) 337 { 338 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b = (Mat_SeqBAIJ*)C->data; 339 PetscErrorCode ierr; 340 PetscInt i,j,n = a->mbs,*bi = b->i,*bj = b->j; 341 PetscInt *ajtmpold,*ajtmp,nz,row; 342 PetscInt *diag_offset = b->diag,*ai=a->i,*aj=a->j,*pj; 343 MatScalar *pv,*v,*rtmp,*pc,*w,*x; 344 MatScalar p1,p2,p3,p4,m1,m2,m3,m4,x1,x2,x3,x4; 345 MatScalar *ba = b->a,*aa = a->a; 346 PetscReal shift = info->shiftamount; 347 PetscBool allowzeropivot,zeropivotdetected; 348 349 PetscFunctionBegin; 350 allowzeropivot = PetscNot(A->erroriffailure); 351 ierr = PetscMalloc1(4*(n+1),&rtmp);CHKERRQ(ierr); 352 for (i=0; i<n; i++) { 353 nz = bi[i+1] - bi[i]; 354 ajtmp = bj + bi[i]; 355 for (j=0; j<nz; j++) { 356 x = rtmp+4*ajtmp[j]; 357 x[0] = x[1] = x[2] = x[3] = 0.0; 358 } 359 /* load in initial (unfactored row) */ 360 nz = ai[i+1] - ai[i]; 361 ajtmpold = aj + ai[i]; 362 v = aa + 4*ai[i]; 363 for (j=0; j<nz; j++) { 364 x = rtmp+4*ajtmpold[j]; 365 x[0] = v[0]; x[1] = v[1]; x[2] = v[2]; x[3] = v[3]; 366 v += 4; 367 } 368 row = *ajtmp++; 369 while (row < i) { 370 pc = rtmp + 4*row; 371 p1 = pc[0]; p2 = pc[1]; p3 = pc[2]; p4 = pc[3]; 372 if (p1 != (PetscScalar)0.0 || p2 != (PetscScalar)0.0 || p3 != (PetscScalar)0.0 || p4 != (PetscScalar)0.0) { 373 pv = ba + 4*diag_offset[row]; 374 pj = bj + diag_offset[row] + 1; 375 x1 = pv[0]; x2 = pv[1]; x3 = pv[2]; x4 = pv[3]; 376 pc[0] = m1 = p1*x1 + p3*x2; 377 pc[1] = m2 = p2*x1 + p4*x2; 378 pc[2] = m3 = p1*x3 + p3*x4; 379 pc[3] = m4 = p2*x3 + p4*x4; 380 nz = bi[row+1] - diag_offset[row] - 1; 381 pv += 4; 382 for (j=0; j<nz; j++) { 383 x1 = pv[0]; x2 = pv[1]; x3 = pv[2]; x4 = pv[3]; 384 x = rtmp + 4*pj[j]; 385 x[0] -= m1*x1 + m3*x2; 386 x[1] -= m2*x1 + m4*x2; 387 x[2] -= m1*x3 + m3*x4; 388 x[3] -= m2*x3 + m4*x4; 389 pv += 4; 390 } 391 ierr = PetscLogFlops(16.0*nz+12.0);CHKERRQ(ierr); 392 } 393 row = *ajtmp++; 394 } 395 /* finished row so stick it into b->a */ 396 pv = ba + 4*bi[i]; 397 pj = bj + bi[i]; 398 nz = bi[i+1] - bi[i]; 399 for (j=0; j<nz; j++) { 400 x = rtmp+4*pj[j]; 401 pv[0] = x[0]; pv[1] = x[1]; pv[2] = x[2]; pv[3] = x[3]; 402 /* 403 printf(" col %d:",pj[j]); 404 PetscInt j1; 405 for (j1=0; j1<4; j1++) printf(" %g,",*(pv+j1)); 406 printf("\n"); 407 */ 408 pv += 4; 409 } 410 /* invert diagonal block */ 411 w = ba + 4*diag_offset[i]; 412 ierr = PetscKernel_A_gets_inverse_A_2(w,shift, allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 413 if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 414 } 415 416 ierr = PetscFree(rtmp);CHKERRQ(ierr); 417 418 C->ops->solve = MatSolve_SeqBAIJ_2_NaturalOrdering_inplace; 419 C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_2_NaturalOrdering_inplace; 420 C->assembled = PETSC_TRUE; 421 422 ierr = PetscLogFlops(1.333333333333*8*b->mbs);CHKERRQ(ierr); /* from inverting diagonal blocks */ 423 PetscFunctionReturn(0); 424 } 425 426 /* ----------------------------------------------------------- */ 427 /* 428 Version for when blocks are 1 by 1. 429 */ 430 PetscErrorCode MatLUFactorNumeric_SeqBAIJ_1(Mat B,Mat A,const MatFactorInfo *info) 431 { 432 Mat C =B; 433 Mat_SeqBAIJ *a =(Mat_SeqBAIJ*)A->data,*b=(Mat_SeqBAIJ*)C->data; 434 IS isrow = b->row,isicol = b->icol; 435 PetscErrorCode ierr; 436 const PetscInt *r,*ic,*ics; 437 const PetscInt n=a->mbs,*ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,*bdiag=b->diag; 438 PetscInt i,j,k,nz,nzL,row,*pj; 439 const PetscInt *ajtmp,*bjtmp; 440 MatScalar *rtmp,*pc,multiplier,*pv; 441 const MatScalar *aa=a->a,*v; 442 PetscBool row_identity,col_identity; 443 FactorShiftCtx sctx; 444 const PetscInt *ddiag; 445 PetscReal rs; 446 MatScalar d; 447 448 PetscFunctionBegin; 449 /* MatPivotSetUp(): initialize shift context sctx */ 450 ierr = PetscMemzero(&sctx,sizeof(FactorShiftCtx));CHKERRQ(ierr); 451 452 if (info->shifttype == (PetscReal) MAT_SHIFT_POSITIVE_DEFINITE) { /* set sctx.shift_top=max{rs} */ 453 ddiag = a->diag; 454 sctx.shift_top = info->zeropivot; 455 for (i=0; i<n; i++) { 456 /* calculate sum(|aij|)-RealPart(aii), amt of shift needed for this row */ 457 d = (aa)[ddiag[i]]; 458 rs = -PetscAbsScalar(d) - PetscRealPart(d); 459 v = aa+ai[i]; 460 nz = ai[i+1] - ai[i]; 461 for (j=0; j<nz; j++) rs += PetscAbsScalar(v[j]); 462 if (rs>sctx.shift_top) sctx.shift_top = rs; 463 } 464 sctx.shift_top *= 1.1; 465 sctx.nshift_max = 5; 466 sctx.shift_lo = 0.; 467 sctx.shift_hi = 1.; 468 } 469 470 ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr); 471 ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr); 472 ierr = PetscMalloc1(n+1,&rtmp);CHKERRQ(ierr); 473 ics = ic; 474 475 do { 476 sctx.newshift = PETSC_FALSE; 477 for (i=0; i<n; i++) { 478 /* zero rtmp */ 479 /* L part */ 480 nz = bi[i+1] - bi[i]; 481 bjtmp = bj + bi[i]; 482 for (j=0; j<nz; j++) rtmp[bjtmp[j]] = 0.0; 483 484 /* U part */ 485 nz = bdiag[i]-bdiag[i+1]; 486 bjtmp = bj + bdiag[i+1]+1; 487 for (j=0; j<nz; j++) rtmp[bjtmp[j]] = 0.0; 488 489 /* load in initial (unfactored row) */ 490 nz = ai[r[i]+1] - ai[r[i]]; 491 ajtmp = aj + ai[r[i]]; 492 v = aa + ai[r[i]]; 493 for (j=0; j<nz; j++) rtmp[ics[ajtmp[j]]] = v[j]; 494 495 /* ZeropivotApply() */ 496 rtmp[i] += sctx.shift_amount; /* shift the diagonal of the matrix */ 497 498 /* elimination */ 499 bjtmp = bj + bi[i]; 500 row = *bjtmp++; 501 nzL = bi[i+1] - bi[i]; 502 for (k=0; k < nzL; k++) { 503 pc = rtmp + row; 504 if (*pc != (PetscScalar)0.0) { 505 pv = b->a + bdiag[row]; 506 multiplier = *pc * (*pv); 507 *pc = multiplier; 508 509 pj = b->j + bdiag[row+1]+1; /* beginning of U(row,:) */ 510 pv = b->a + bdiag[row+1]+1; 511 nz = bdiag[row]-bdiag[row+1]-1; /* num of entries in U(row,:) excluding diag */ 512 for (j=0; j<nz; j++) rtmp[pj[j]] -= multiplier * pv[j]; 513 ierr = PetscLogFlops(2.0*nz);CHKERRQ(ierr); 514 } 515 row = *bjtmp++; 516 } 517 518 /* finished row so stick it into b->a */ 519 rs = 0.0; 520 /* L part */ 521 pv = b->a + bi[i]; 522 pj = b->j + bi[i]; 523 nz = bi[i+1] - bi[i]; 524 for (j=0; j<nz; j++) { 525 pv[j] = rtmp[pj[j]]; rs += PetscAbsScalar(pv[j]); 526 } 527 528 /* U part */ 529 pv = b->a + bdiag[i+1]+1; 530 pj = b->j + bdiag[i+1]+1; 531 nz = bdiag[i] - bdiag[i+1]-1; 532 for (j=0; j<nz; j++) { 533 pv[j] = rtmp[pj[j]]; rs += PetscAbsScalar(pv[j]); 534 } 535 536 sctx.rs = rs; 537 sctx.pv = rtmp[i]; 538 ierr = MatPivotCheck(B,A,info,&sctx,i);CHKERRQ(ierr); 539 if (sctx.newshift) break; /* break for-loop */ 540 rtmp[i] = sctx.pv; /* sctx.pv might be updated in the case of MAT_SHIFT_INBLOCKS */ 541 542 /* Mark diagonal and invert diagonal for simplier triangular solves */ 543 pv = b->a + bdiag[i]; 544 *pv = (PetscScalar)1.0/rtmp[i]; 545 546 } /* endof for (i=0; i<n; i++) { */ 547 548 /* MatPivotRefine() */ 549 if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE && !sctx.newshift && sctx.shift_fraction>0 && sctx.nshift<sctx.nshift_max) { 550 /* 551 * if no shift in this attempt & shifting & started shifting & can refine, 552 * then try lower shift 553 */ 554 sctx.shift_hi = sctx.shift_fraction; 555 sctx.shift_fraction = (sctx.shift_hi+sctx.shift_lo)/2.; 556 sctx.shift_amount = sctx.shift_fraction * sctx.shift_top; 557 sctx.newshift = PETSC_TRUE; 558 sctx.nshift++; 559 } 560 } while (sctx.newshift); 561 562 ierr = PetscFree(rtmp);CHKERRQ(ierr); 563 ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr); 564 ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr); 565 566 ierr = ISIdentity(isrow,&row_identity);CHKERRQ(ierr); 567 ierr = ISIdentity(isicol,&col_identity);CHKERRQ(ierr); 568 if (row_identity && col_identity) { 569 C->ops->solve = MatSolve_SeqBAIJ_1_NaturalOrdering; 570 C->ops->forwardsolve = MatForwardSolve_SeqBAIJ_1_NaturalOrdering; 571 C->ops->backwardsolve = MatBackwardSolve_SeqBAIJ_1_NaturalOrdering; 572 C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1_NaturalOrdering; 573 } else { 574 C->ops->solve = MatSolve_SeqBAIJ_1; 575 C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1; 576 } 577 C->assembled = PETSC_TRUE; 578 ierr = PetscLogFlops(C->cmap->n);CHKERRQ(ierr); 579 580 /* MatShiftView(A,info,&sctx) */ 581 if (sctx.nshift) { 582 if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) { 583 ierr = PetscInfo4(A,"number of shift_pd tries %D, shift_amount %g, diagonal shifted up by %e fraction top_value %e\n",sctx.nshift,(double)sctx.shift_amount,(double)sctx.shift_fraction,(double)sctx.shift_top);CHKERRQ(ierr); 584 } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) { 585 ierr = PetscInfo2(A,"number of shift_nz tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);CHKERRQ(ierr); 586 } else if (info->shifttype == (PetscReal)MAT_SHIFT_INBLOCKS) { 587 ierr = PetscInfo2(A,"number of shift_inblocks applied %D, each shift_amount %g\n",sctx.nshift,(double)info->shiftamount);CHKERRQ(ierr); 588 } 589 } 590 PetscFunctionReturn(0); 591 } 592 593 PetscErrorCode MatLUFactorNumeric_SeqBAIJ_1_inplace(Mat C,Mat A,const MatFactorInfo *info) 594 { 595 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b = (Mat_SeqBAIJ*)C->data; 596 IS isrow = b->row,isicol = b->icol; 597 PetscErrorCode ierr; 598 const PetscInt *r,*ic; 599 PetscInt i,j,n = a->mbs,*bi = b->i,*bj = b->j; 600 PetscInt *ajtmpold,*ajtmp,nz,row,*ai = a->i,*aj = a->j; 601 PetscInt *diag_offset = b->diag,diag,*pj; 602 MatScalar *pv,*v,*rtmp,multiplier,*pc; 603 MatScalar *ba = b->a,*aa = a->a; 604 PetscBool row_identity, col_identity; 605 606 PetscFunctionBegin; 607 ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr); 608 ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr); 609 ierr = PetscMalloc1(n+1,&rtmp);CHKERRQ(ierr); 610 611 for (i=0; i<n; i++) { 612 nz = bi[i+1] - bi[i]; 613 ajtmp = bj + bi[i]; 614 for (j=0; j<nz; j++) rtmp[ajtmp[j]] = 0.0; 615 616 /* load in initial (unfactored row) */ 617 nz = ai[r[i]+1] - ai[r[i]]; 618 ajtmpold = aj + ai[r[i]]; 619 v = aa + ai[r[i]]; 620 for (j=0; j<nz; j++) rtmp[ic[ajtmpold[j]]] = v[j]; 621 622 row = *ajtmp++; 623 while (row < i) { 624 pc = rtmp + row; 625 if (*pc != 0.0) { 626 pv = ba + diag_offset[row]; 627 pj = bj + diag_offset[row] + 1; 628 multiplier = *pc * *pv++; 629 *pc = multiplier; 630 nz = bi[row+1] - diag_offset[row] - 1; 631 for (j=0; j<nz; j++) rtmp[pj[j]] -= multiplier * pv[j]; 632 ierr = PetscLogFlops(1.0+2.0*nz);CHKERRQ(ierr); 633 } 634 row = *ajtmp++; 635 } 636 /* finished row so stick it into b->a */ 637 pv = ba + bi[i]; 638 pj = bj + bi[i]; 639 nz = bi[i+1] - bi[i]; 640 for (j=0; j<nz; j++) pv[j] = rtmp[pj[j]]; 641 diag = diag_offset[i] - bi[i]; 642 /* check pivot entry for current row */ 643 if (pv[diag] == 0.0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot: row in original ordering %D in permuted ordering %D",r[i],i); 644 pv[diag] = 1.0/pv[diag]; 645 } 646 647 ierr = PetscFree(rtmp);CHKERRQ(ierr); 648 ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr); 649 ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr); 650 ierr = ISIdentity(isrow,&row_identity);CHKERRQ(ierr); 651 ierr = ISIdentity(isicol,&col_identity);CHKERRQ(ierr); 652 if (row_identity && col_identity) { 653 C->ops->solve = MatSolve_SeqBAIJ_1_NaturalOrdering_inplace; 654 C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1_NaturalOrdering_inplace; 655 } else { 656 C->ops->solve = MatSolve_SeqBAIJ_1_inplace; 657 C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_1_inplace; 658 } 659 C->assembled = PETSC_TRUE; 660 ierr = PetscLogFlops(C->cmap->n);CHKERRQ(ierr); 661 PetscFunctionReturn(0); 662 } 663 664 static PetscErrorCode MatFactorGetSolverType_petsc(Mat A,MatSolverType *type) 665 { 666 PetscFunctionBegin; 667 *type = MATSOLVERPETSC; 668 PetscFunctionReturn(0); 669 } 670 671 PETSC_INTERN PetscErrorCode MatGetFactor_seqbaij_petsc(Mat A,MatFactorType ftype,Mat *B) 672 { 673 PetscInt n = A->rmap->n; 674 PetscErrorCode ierr; 675 676 PetscFunctionBegin; 677 #if defined(PETSC_USE_COMPLEX) 678 if (A->hermitian && (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Hermitian Factor is not supported"); 679 #endif 680 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 681 ierr = MatSetSizes(*B,n,n,n,n);CHKERRQ(ierr); 682 if (ftype == MAT_FACTOR_LU || ftype == MAT_FACTOR_ILU || ftype == MAT_FACTOR_ILUDT) { 683 ierr = MatSetType(*B,MATSEQBAIJ);CHKERRQ(ierr); 684 685 (*B)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqBAIJ; 686 (*B)->ops->ilufactorsymbolic = MatILUFactorSymbolic_SeqBAIJ; 687 ierr = PetscStrallocpy(MATORDERINGND,(char**)&(*B)->preferredordering[MAT_FACTOR_LU]);CHKERRQ(ierr); 688 ierr = PetscStrallocpy(MATORDERINGNATURAL,(char**)&(*B)->preferredordering[MAT_FACTOR_ILU]);CHKERRQ(ierr); 689 ierr = PetscStrallocpy(MATORDERINGNATURAL,(char**)&(*B)->preferredordering[MAT_FACTOR_ILUDT]);CHKERRQ(ierr); 690 } else if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 691 ierr = MatSetType(*B,MATSEQSBAIJ);CHKERRQ(ierr); 692 ierr = MatSeqSBAIJSetPreallocation(*B,A->rmap->bs,MAT_SKIP_ALLOCATION,NULL);CHKERRQ(ierr); 693 694 (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqBAIJ; 695 (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqBAIJ; 696 /* Future optimization would be direct symbolic and numerical factorization for BAIJ to support orderings and Cholesky, instead of first converting to SBAIJ */ 697 ierr = PetscStrallocpy(MATORDERINGNATURAL,(char**)&(*B)->preferredordering[MAT_FACTOR_CHOLESKY]);CHKERRQ(ierr); 698 ierr = PetscStrallocpy(MATORDERINGNATURAL,(char**)&(*B)->preferredordering[MAT_FACTOR_ICC]);CHKERRQ(ierr); 699 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported"); 700 (*B)->factortype = ftype; 701 (*B)->canuseordering = PETSC_TRUE; 702 703 ierr = PetscFree((*B)->solvertype);CHKERRQ(ierr); 704 ierr = PetscStrallocpy(MATSOLVERPETSC,&(*B)->solvertype);CHKERRQ(ierr); 705 ierr = PetscObjectComposeFunction((PetscObject)*B,"MatFactorGetSolverType_C",MatFactorGetSolverType_petsc);CHKERRQ(ierr); 706 PetscFunctionReturn(0); 707 } 708 709 /* ----------------------------------------------------------- */ 710 PetscErrorCode MatLUFactor_SeqBAIJ(Mat A,IS row,IS col,const MatFactorInfo *info) 711 { 712 PetscErrorCode ierr; 713 Mat C; 714 715 PetscFunctionBegin; 716 ierr = MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_LU,&C);CHKERRQ(ierr); 717 ierr = MatLUFactorSymbolic(C,A,row,col,info);CHKERRQ(ierr); 718 ierr = MatLUFactorNumeric(C,A,info);CHKERRQ(ierr); 719 720 A->ops->solve = C->ops->solve; 721 A->ops->solvetranspose = C->ops->solvetranspose; 722 723 ierr = MatHeaderMerge(A,&C);CHKERRQ(ierr); 724 ierr = PetscLogObjectParent((PetscObject)A,(PetscObject)((Mat_SeqBAIJ*)(A->data))->icol);CHKERRQ(ierr); 725 PetscFunctionReturn(0); 726 } 727 728 #include <../src/mat/impls/sbaij/seq/sbaij.h> 729 PetscErrorCode MatCholeskyFactorNumeric_SeqBAIJ_N(Mat C,Mat A,const MatFactorInfo *info) 730 { 731 PetscErrorCode ierr; 732 Mat_SeqBAIJ *a=(Mat_SeqBAIJ*)A->data; 733 Mat_SeqSBAIJ *b=(Mat_SeqSBAIJ*)C->data; 734 IS ip=b->row; 735 const PetscInt *rip; 736 PetscInt i,j,mbs=a->mbs,bs=A->rmap->bs,*bi=b->i,*bj=b->j,*bcol; 737 PetscInt *ai=a->i,*aj=a->j; 738 PetscInt k,jmin,jmax,*jl,*il,col,nexti,ili,nz; 739 MatScalar *rtmp,*ba=b->a,*bval,*aa=a->a,dk,uikdi; 740 PetscReal rs; 741 FactorShiftCtx sctx; 742 743 PetscFunctionBegin; 744 if (bs > 1) { /* convert A to a SBAIJ matrix and apply Cholesky factorization from it */ 745 if (!a->sbaijMat) { 746 ierr = MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&a->sbaijMat);CHKERRQ(ierr); 747 } 748 ierr = (a->sbaijMat)->ops->choleskyfactornumeric(C,a->sbaijMat,info);CHKERRQ(ierr); 749 ierr = MatDestroy(&a->sbaijMat);CHKERRQ(ierr); 750 PetscFunctionReturn(0); 751 } 752 753 /* MatPivotSetUp(): initialize shift context sctx */ 754 ierr = PetscMemzero(&sctx,sizeof(FactorShiftCtx));CHKERRQ(ierr); 755 756 ierr = ISGetIndices(ip,&rip);CHKERRQ(ierr); 757 ierr = PetscMalloc3(mbs,&rtmp,mbs,&il,mbs,&jl);CHKERRQ(ierr); 758 759 sctx.shift_amount = 0.; 760 sctx.nshift = 0; 761 do { 762 sctx.newshift = PETSC_FALSE; 763 for (i=0; i<mbs; i++) { 764 rtmp[i] = 0.0; jl[i] = mbs; il[0] = 0; 765 } 766 767 for (k = 0; k<mbs; k++) { 768 bval = ba + bi[k]; 769 /* initialize k-th row by the perm[k]-th row of A */ 770 jmin = ai[rip[k]]; jmax = ai[rip[k]+1]; 771 for (j = jmin; j < jmax; j++) { 772 col = rip[aj[j]]; 773 if (col >= k) { /* only take upper triangular entry */ 774 rtmp[col] = aa[j]; 775 *bval++ = 0.0; /* for in-place factorization */ 776 } 777 } 778 779 /* shift the diagonal of the matrix */ 780 if (sctx.nshift) rtmp[k] += sctx.shift_amount; 781 782 /* modify k-th row by adding in those rows i with U(i,k)!=0 */ 783 dk = rtmp[k]; 784 i = jl[k]; /* first row to be added to k_th row */ 785 786 while (i < k) { 787 nexti = jl[i]; /* next row to be added to k_th row */ 788 789 /* compute multiplier, update diag(k) and U(i,k) */ 790 ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */ 791 uikdi = -ba[ili]*ba[bi[i]]; /* diagonal(k) */ 792 dk += uikdi*ba[ili]; 793 ba[ili] = uikdi; /* -U(i,k) */ 794 795 /* add multiple of row i to k-th row */ 796 jmin = ili + 1; jmax = bi[i+1]; 797 if (jmin < jmax) { 798 for (j=jmin; j<jmax; j++) rtmp[bj[j]] += uikdi*ba[j]; 799 /* update il and jl for row i */ 800 il[i] = jmin; 801 j = bj[jmin]; jl[i] = jl[j]; jl[j] = i; 802 } 803 i = nexti; 804 } 805 806 /* shift the diagonals when zero pivot is detected */ 807 /* compute rs=sum of abs(off-diagonal) */ 808 rs = 0.0; 809 jmin = bi[k]+1; 810 nz = bi[k+1] - jmin; 811 if (nz) { 812 bcol = bj + jmin; 813 while (nz--) { 814 rs += PetscAbsScalar(rtmp[*bcol]); 815 bcol++; 816 } 817 } 818 819 sctx.rs = rs; 820 sctx.pv = dk; 821 ierr = MatPivotCheck(C,A,info,&sctx,k);CHKERRQ(ierr); 822 if (sctx.newshift) break; 823 dk = sctx.pv; 824 825 /* copy data into U(k,:) */ 826 ba[bi[k]] = 1.0/dk; /* U(k,k) */ 827 jmin = bi[k]+1; jmax = bi[k+1]; 828 if (jmin < jmax) { 829 for (j=jmin; j<jmax; j++) { 830 col = bj[j]; ba[j] = rtmp[col]; rtmp[col] = 0.0; 831 } 832 /* add the k-th row into il and jl */ 833 il[k] = jmin; 834 i = bj[jmin]; jl[k] = jl[i]; jl[i] = k; 835 } 836 } 837 } while (sctx.newshift); 838 ierr = PetscFree3(rtmp,il,jl);CHKERRQ(ierr); 839 840 ierr = ISRestoreIndices(ip,&rip);CHKERRQ(ierr); 841 842 C->assembled = PETSC_TRUE; 843 C->preallocated = PETSC_TRUE; 844 845 ierr = PetscLogFlops(C->rmap->N);CHKERRQ(ierr); 846 if (sctx.nshift) { 847 if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) { 848 ierr = PetscInfo2(A,"number of shiftpd tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);CHKERRQ(ierr); 849 } else if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) { 850 ierr = PetscInfo2(A,"number of shiftnz tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);CHKERRQ(ierr); 851 } 852 } 853 PetscFunctionReturn(0); 854 } 855 856 PetscErrorCode MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering(Mat C,Mat A,const MatFactorInfo *info) 857 { 858 Mat_SeqBAIJ *a=(Mat_SeqBAIJ*)A->data; 859 Mat_SeqSBAIJ *b=(Mat_SeqSBAIJ*)C->data; 860 PetscErrorCode ierr; 861 PetscInt i,j,am=a->mbs; 862 PetscInt *ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j; 863 PetscInt k,jmin,*jl,*il,nexti,ili,*acol,*bcol,nz; 864 MatScalar *rtmp,*ba=b->a,*aa=a->a,dk,uikdi,*aval,*bval; 865 PetscReal rs; 866 FactorShiftCtx sctx; 867 868 PetscFunctionBegin; 869 /* MatPivotSetUp(): initialize shift context sctx */ 870 ierr = PetscMemzero(&sctx,sizeof(FactorShiftCtx));CHKERRQ(ierr); 871 872 ierr = PetscMalloc3(am,&rtmp,am,&il,am,&jl);CHKERRQ(ierr); 873 874 do { 875 sctx.newshift = PETSC_FALSE; 876 for (i=0; i<am; i++) { 877 rtmp[i] = 0.0; jl[i] = am; il[0] = 0; 878 } 879 880 for (k = 0; k<am; k++) { 881 /* initialize k-th row with elements nonzero in row perm(k) of A */ 882 nz = ai[k+1] - ai[k]; 883 acol = aj + ai[k]; 884 aval = aa + ai[k]; 885 bval = ba + bi[k]; 886 while (nz--) { 887 if (*acol < k) { /* skip lower triangular entries */ 888 acol++; aval++; 889 } else { 890 rtmp[*acol++] = *aval++; 891 *bval++ = 0.0; /* for in-place factorization */ 892 } 893 } 894 895 /* shift the diagonal of the matrix */ 896 if (sctx.nshift) rtmp[k] += sctx.shift_amount; 897 898 /* modify k-th row by adding in those rows i with U(i,k)!=0 */ 899 dk = rtmp[k]; 900 i = jl[k]; /* first row to be added to k_th row */ 901 902 while (i < k) { 903 nexti = jl[i]; /* next row to be added to k_th row */ 904 /* compute multiplier, update D(k) and U(i,k) */ 905 ili = il[i]; /* index of first nonzero element in U(i,k:bms-1) */ 906 uikdi = -ba[ili]*ba[bi[i]]; 907 dk += uikdi*ba[ili]; 908 ba[ili] = uikdi; /* -U(i,k) */ 909 910 /* add multiple of row i to k-th row ... */ 911 jmin = ili + 1; 912 nz = bi[i+1] - jmin; 913 if (nz > 0) { 914 bcol = bj + jmin; 915 bval = ba + jmin; 916 while (nz--) rtmp[*bcol++] += uikdi*(*bval++); 917 /* update il and jl for i-th row */ 918 il[i] = jmin; 919 j = bj[jmin]; jl[i] = jl[j]; jl[j] = i; 920 } 921 i = nexti; 922 } 923 924 /* shift the diagonals when zero pivot is detected */ 925 /* compute rs=sum of abs(off-diagonal) */ 926 rs = 0.0; 927 jmin = bi[k]+1; 928 nz = bi[k+1] - jmin; 929 if (nz) { 930 bcol = bj + jmin; 931 while (nz--) { 932 rs += PetscAbsScalar(rtmp[*bcol]); 933 bcol++; 934 } 935 } 936 937 sctx.rs = rs; 938 sctx.pv = dk; 939 ierr = MatPivotCheck(C,A,info,&sctx,k);CHKERRQ(ierr); 940 if (sctx.newshift) break; /* sctx.shift_amount is updated */ 941 dk = sctx.pv; 942 943 /* copy data into U(k,:) */ 944 ba[bi[k]] = 1.0/dk; 945 jmin = bi[k]+1; 946 nz = bi[k+1] - jmin; 947 if (nz) { 948 bcol = bj + jmin; 949 bval = ba + jmin; 950 while (nz--) { 951 *bval++ = rtmp[*bcol]; 952 rtmp[*bcol++] = 0.0; 953 } 954 /* add k-th row into il and jl */ 955 il[k] = jmin; 956 i = bj[jmin]; jl[k] = jl[i]; jl[i] = k; 957 } 958 } 959 } while (sctx.newshift); 960 ierr = PetscFree3(rtmp,il,jl);CHKERRQ(ierr); 961 962 C->ops->solve = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace; 963 C->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace; 964 C->assembled = PETSC_TRUE; 965 C->preallocated = PETSC_TRUE; 966 967 ierr = PetscLogFlops(C->rmap->N);CHKERRQ(ierr); 968 if (sctx.nshift) { 969 if (info->shifttype == (PetscReal)MAT_SHIFT_NONZERO) { 970 ierr = PetscInfo2(A,"number of shiftnz tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);CHKERRQ(ierr); 971 } else if (info->shifttype == (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE) { 972 ierr = PetscInfo2(A,"number of shiftpd tries %D, shift_amount %g\n",sctx.nshift,(double)sctx.shift_amount);CHKERRQ(ierr); 973 } 974 } 975 PetscFunctionReturn(0); 976 } 977 978 #include <petscbt.h> 979 #include <../src/mat/utils/freespace.h> 980 PetscErrorCode MatICCFactorSymbolic_SeqBAIJ(Mat fact,Mat A,IS perm,const MatFactorInfo *info) 981 { 982 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 983 Mat_SeqSBAIJ *b; 984 Mat B; 985 PetscErrorCode ierr; 986 PetscBool perm_identity,missing; 987 PetscInt reallocs=0,i,*ai=a->i,*aj=a->j,am=a->mbs,bs=A->rmap->bs,*ui; 988 const PetscInt *rip; 989 PetscInt jmin,jmax,nzk,k,j,*jl,prow,*il,nextprow; 990 PetscInt nlnk,*lnk,*lnk_lvl=NULL,ncols,ncols_upper,*cols,*cols_lvl,*uj,**uj_ptr,**uj_lvl_ptr; 991 PetscReal fill =info->fill,levels=info->levels; 992 PetscFreeSpaceList free_space =NULL,current_space=NULL; 993 PetscFreeSpaceList free_space_lvl=NULL,current_space_lvl=NULL; 994 PetscBT lnkbt; 995 996 PetscFunctionBegin; 997 ierr = MatMissingDiagonal(A,&missing,&i);CHKERRQ(ierr); 998 if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",i); 999 1000 if (bs > 1) { 1001 if (!a->sbaijMat) { 1002 ierr = MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&a->sbaijMat);CHKERRQ(ierr); 1003 } 1004 (fact)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ; /* undue the change made in MatGetFactor_seqbaij_petsc */ 1005 1006 ierr = MatICCFactorSymbolic(fact,a->sbaijMat,perm,info);CHKERRQ(ierr); 1007 PetscFunctionReturn(0); 1008 } 1009 1010 ierr = ISIdentity(perm,&perm_identity);CHKERRQ(ierr); 1011 ierr = ISGetIndices(perm,&rip);CHKERRQ(ierr); 1012 1013 /* special case that simply copies fill pattern */ 1014 if (!levels && perm_identity) { 1015 ierr = PetscMalloc1(am+1,&ui);CHKERRQ(ierr); 1016 for (i=0; i<am; i++) ui[i] = ai[i+1] - a->diag[i]; /* ui: rowlengths - changes when !perm_identity */ 1017 B = fact; 1018 ierr = MatSeqSBAIJSetPreallocation(B,1,0,ui);CHKERRQ(ierr); 1019 1020 1021 b = (Mat_SeqSBAIJ*)B->data; 1022 uj = b->j; 1023 for (i=0; i<am; i++) { 1024 aj = a->j + a->diag[i]; 1025 for (j=0; j<ui[i]; j++) *uj++ = *aj++; 1026 b->ilen[i] = ui[i]; 1027 } 1028 ierr = PetscFree(ui);CHKERRQ(ierr); 1029 1030 B->factortype = MAT_FACTOR_NONE; 1031 1032 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1033 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1034 B->factortype = MAT_FACTOR_ICC; 1035 1036 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering; 1037 PetscFunctionReturn(0); 1038 } 1039 1040 /* initialization */ 1041 ierr = PetscMalloc1(am+1,&ui);CHKERRQ(ierr); 1042 ui[0] = 0; 1043 ierr = PetscMalloc1(2*am+1,&cols_lvl);CHKERRQ(ierr); 1044 1045 /* jl: linked list for storing indices of the pivot rows 1046 il: il[i] points to the 1st nonzero entry of U(i,k:am-1) */ 1047 ierr = PetscMalloc4(am,&uj_ptr,am,&uj_lvl_ptr,am,&il,am,&jl);CHKERRQ(ierr); 1048 for (i=0; i<am; i++) { 1049 jl[i] = am; il[i] = 0; 1050 } 1051 1052 /* create and initialize a linked list for storing column indices of the active row k */ 1053 nlnk = am + 1; 1054 ierr = PetscIncompleteLLCreate(am,am,nlnk,lnk,lnk_lvl,lnkbt);CHKERRQ(ierr); 1055 1056 /* initial FreeSpace size is fill*(ai[am]+am)/2 */ 1057 ierr = PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ai[am]/2,am/2)),&free_space);CHKERRQ(ierr); 1058 1059 current_space = free_space; 1060 1061 ierr = PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ai[am]/2,am/2)),&free_space_lvl);CHKERRQ(ierr); 1062 current_space_lvl = free_space_lvl; 1063 1064 for (k=0; k<am; k++) { /* for each active row k */ 1065 /* initialize lnk by the column indices of row rip[k] of A */ 1066 nzk = 0; 1067 ncols = ai[rip[k]+1] - ai[rip[k]]; 1068 ncols_upper = 0; 1069 cols = cols_lvl + am; 1070 for (j=0; j<ncols; j++) { 1071 i = rip[*(aj + ai[rip[k]] + j)]; 1072 if (i >= k) { /* only take upper triangular entry */ 1073 cols[ncols_upper] = i; 1074 cols_lvl[ncols_upper] = -1; /* initialize level for nonzero entries */ 1075 ncols_upper++; 1076 } 1077 } 1078 ierr = PetscIncompleteLLAdd(ncols_upper,cols,levels,cols_lvl,am,nlnk,lnk,lnk_lvl,lnkbt);CHKERRQ(ierr); 1079 nzk += nlnk; 1080 1081 /* update lnk by computing fill-in for each pivot row to be merged in */ 1082 prow = jl[k]; /* 1st pivot row */ 1083 1084 while (prow < k) { 1085 nextprow = jl[prow]; 1086 1087 /* merge prow into k-th row */ 1088 jmin = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:am-1) */ 1089 jmax = ui[prow+1]; 1090 ncols = jmax-jmin; 1091 i = jmin - ui[prow]; 1092 cols = uj_ptr[prow] + i; /* points to the 2nd nzero entry in U(prow,k:am-1) */ 1093 for (j=0; j<ncols; j++) cols_lvl[j] = *(uj_lvl_ptr[prow] + i + j); 1094 ierr = PetscIncompleteLLAddSorted(ncols,cols,levels,cols_lvl,am,nlnk,lnk,lnk_lvl,lnkbt);CHKERRQ(ierr); 1095 nzk += nlnk; 1096 1097 /* update il and jl for prow */ 1098 if (jmin < jmax) { 1099 il[prow] = jmin; 1100 1101 j = *cols; jl[prow] = jl[j]; jl[j] = prow; 1102 } 1103 prow = nextprow; 1104 } 1105 1106 /* if free space is not available, make more free space */ 1107 if (current_space->local_remaining<nzk) { 1108 i = am - k + 1; /* num of unfactored rows */ 1109 i = PetscMin(PetscIntMultTruncate(i,nzk), PetscIntMultTruncate(i,i-1)); /* i*nzk, i*(i-1): estimated and max additional space needed */ 1110 ierr = PetscFreeSpaceGet(i,¤t_space);CHKERRQ(ierr); 1111 ierr = PetscFreeSpaceGet(i,¤t_space_lvl);CHKERRQ(ierr); 1112 reallocs++; 1113 } 1114 1115 /* copy data into free_space and free_space_lvl, then initialize lnk */ 1116 ierr = PetscIncompleteLLClean(am,am,nzk,lnk,lnk_lvl,current_space->array,current_space_lvl->array,lnkbt);CHKERRQ(ierr); 1117 1118 /* add the k-th row into il and jl */ 1119 if (nzk-1 > 0) { 1120 i = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:am-1) */ 1121 jl[k] = jl[i]; jl[i] = k; 1122 il[k] = ui[k] + 1; 1123 } 1124 uj_ptr[k] = current_space->array; 1125 uj_lvl_ptr[k] = current_space_lvl->array; 1126 1127 current_space->array += nzk; 1128 current_space->local_used += nzk; 1129 current_space->local_remaining -= nzk; 1130 1131 current_space_lvl->array += nzk; 1132 current_space_lvl->local_used += nzk; 1133 current_space_lvl->local_remaining -= nzk; 1134 1135 ui[k+1] = ui[k] + nzk; 1136 } 1137 1138 ierr = ISRestoreIndices(perm,&rip);CHKERRQ(ierr); 1139 ierr = PetscFree4(uj_ptr,uj_lvl_ptr,il,jl);CHKERRQ(ierr); 1140 ierr = PetscFree(cols_lvl);CHKERRQ(ierr); 1141 1142 /* copy free_space into uj and free free_space; set uj in new datastructure; */ 1143 ierr = PetscMalloc1(ui[am]+1,&uj);CHKERRQ(ierr); 1144 ierr = PetscFreeSpaceContiguous(&free_space,uj);CHKERRQ(ierr); 1145 ierr = PetscIncompleteLLDestroy(lnk,lnkbt);CHKERRQ(ierr); 1146 ierr = PetscFreeSpaceDestroy(free_space_lvl);CHKERRQ(ierr); 1147 1148 /* put together the new matrix in MATSEQSBAIJ format */ 1149 B = fact; 1150 ierr = MatSeqSBAIJSetPreallocation(B,1,MAT_SKIP_ALLOCATION,NULL);CHKERRQ(ierr); 1151 1152 b = (Mat_SeqSBAIJ*)B->data; 1153 b->singlemalloc = PETSC_FALSE; 1154 b->free_a = PETSC_TRUE; 1155 b->free_ij = PETSC_TRUE; 1156 1157 ierr = PetscMalloc1(ui[am]+1,&b->a);CHKERRQ(ierr); 1158 1159 b->j = uj; 1160 b->i = ui; 1161 b->diag = NULL; 1162 b->ilen = NULL; 1163 b->imax = NULL; 1164 b->row = perm; 1165 b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */ 1166 1167 ierr = PetscObjectReference((PetscObject)perm);CHKERRQ(ierr); 1168 1169 b->icol = perm; 1170 1171 ierr = PetscObjectReference((PetscObject)perm);CHKERRQ(ierr); 1172 ierr = PetscMalloc1(am+1,&b->solve_work);CHKERRQ(ierr); 1173 ierr = PetscLogObjectMemory((PetscObject)B,(ui[am]-am)*(sizeof(PetscInt)+sizeof(MatScalar)));CHKERRQ(ierr); 1174 1175 b->maxnz = b->nz = ui[am]; 1176 1177 B->info.factor_mallocs = reallocs; 1178 B->info.fill_ratio_given = fill; 1179 if (ai[am] != 0.) { 1180 /* nonzeros in lower triangular part of A (includign diagonals)= (ai[am]+am)/2 */ 1181 B->info.fill_ratio_needed = ((PetscReal)2*ui[am])/(ai[am]+am); 1182 } else { 1183 B->info.fill_ratio_needed = 0.0; 1184 } 1185 #if defined(PETSC_USE_INFO) 1186 if (ai[am] != 0) { 1187 PetscReal af = B->info.fill_ratio_needed; 1188 ierr = PetscInfo3(A,"Reallocs %D Fill ratio:given %g needed %g\n",reallocs,(double)fill,(double)af);CHKERRQ(ierr); 1189 ierr = PetscInfo1(A,"Run with -pc_factor_fill %g or use \n",(double)af);CHKERRQ(ierr); 1190 ierr = PetscInfo1(A,"PCFactorSetFill(pc,%g) for best performance.\n",(double)af);CHKERRQ(ierr); 1191 } else { 1192 ierr = PetscInfo(A,"Empty matrix\n");CHKERRQ(ierr); 1193 } 1194 #endif 1195 if (perm_identity) { 1196 B->ops->solve = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace; 1197 B->ops->solvetranspose = MatSolve_SeqSBAIJ_1_NaturalOrdering_inplace; 1198 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering; 1199 } else { 1200 (fact)->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N; 1201 } 1202 PetscFunctionReturn(0); 1203 } 1204 1205 PetscErrorCode MatCholeskyFactorSymbolic_SeqBAIJ(Mat fact,Mat A,IS perm,const MatFactorInfo *info) 1206 { 1207 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 1208 Mat_SeqSBAIJ *b; 1209 Mat B; 1210 PetscErrorCode ierr; 1211 PetscBool perm_identity,missing; 1212 PetscReal fill = info->fill; 1213 const PetscInt *rip; 1214 PetscInt i,mbs=a->mbs,bs=A->rmap->bs,*ai=a->i,*aj=a->j,reallocs=0,prow; 1215 PetscInt *jl,jmin,jmax,nzk,*ui,k,j,*il,nextprow; 1216 PetscInt nlnk,*lnk,ncols,ncols_upper,*cols,*uj,**ui_ptr,*uj_ptr; 1217 PetscFreeSpaceList free_space=NULL,current_space=NULL; 1218 PetscBT lnkbt; 1219 1220 PetscFunctionBegin; 1221 if (bs > 1) { /* convert to seqsbaij */ 1222 if (!a->sbaijMat) { 1223 ierr = MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&a->sbaijMat);CHKERRQ(ierr); 1224 } 1225 (fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ; /* undue the change made in MatGetFactor_seqbaij_petsc */ 1226 1227 ierr = MatCholeskyFactorSymbolic(fact,a->sbaijMat,perm,info);CHKERRQ(ierr); 1228 PetscFunctionReturn(0); 1229 } 1230 1231 ierr = MatMissingDiagonal(A,&missing,&i);CHKERRQ(ierr); 1232 if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",i); 1233 1234 /* check whether perm is the identity mapping */ 1235 ierr = ISIdentity(perm,&perm_identity);CHKERRQ(ierr); 1236 if (!perm_identity) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Matrix reordering is not supported"); 1237 ierr = ISGetIndices(perm,&rip);CHKERRQ(ierr); 1238 1239 /* initialization */ 1240 ierr = PetscMalloc1(mbs+1,&ui);CHKERRQ(ierr); 1241 ui[0] = 0; 1242 1243 /* jl: linked list for storing indices of the pivot rows 1244 il: il[i] points to the 1st nonzero entry of U(i,k:mbs-1) */ 1245 ierr = PetscMalloc4(mbs,&ui_ptr,mbs,&il,mbs,&jl,mbs,&cols);CHKERRQ(ierr); 1246 for (i=0; i<mbs; i++) { 1247 jl[i] = mbs; il[i] = 0; 1248 } 1249 1250 /* create and initialize a linked list for storing column indices of the active row k */ 1251 nlnk = mbs + 1; 1252 ierr = PetscLLCreate(mbs,mbs,nlnk,lnk,lnkbt);CHKERRQ(ierr); 1253 1254 /* initial FreeSpace size is fill* (ai[mbs]+mbs)/2 */ 1255 ierr = PetscFreeSpaceGet(PetscRealIntMultTruncate(fill,PetscIntSumTruncate(ai[mbs]/2,mbs/2)),&free_space);CHKERRQ(ierr); 1256 1257 current_space = free_space; 1258 1259 for (k=0; k<mbs; k++) { /* for each active row k */ 1260 /* initialize lnk by the column indices of row rip[k] of A */ 1261 nzk = 0; 1262 ncols = ai[rip[k]+1] - ai[rip[k]]; 1263 ncols_upper = 0; 1264 for (j=0; j<ncols; j++) { 1265 i = rip[*(aj + ai[rip[k]] + j)]; 1266 if (i >= k) { /* only take upper triangular entry */ 1267 cols[ncols_upper] = i; 1268 ncols_upper++; 1269 } 1270 } 1271 ierr = PetscLLAdd(ncols_upper,cols,mbs,nlnk,lnk,lnkbt);CHKERRQ(ierr); 1272 nzk += nlnk; 1273 1274 /* update lnk by computing fill-in for each pivot row to be merged in */ 1275 prow = jl[k]; /* 1st pivot row */ 1276 1277 while (prow < k) { 1278 nextprow = jl[prow]; 1279 /* merge prow into k-th row */ 1280 jmin = il[prow] + 1; /* index of the 2nd nzero entry in U(prow,k:mbs-1) */ 1281 jmax = ui[prow+1]; 1282 ncols = jmax-jmin; 1283 uj_ptr = ui_ptr[prow] + jmin - ui[prow]; /* points to the 2nd nzero entry in U(prow,k:mbs-1) */ 1284 ierr = PetscLLAddSorted(ncols,uj_ptr,mbs,nlnk,lnk,lnkbt);CHKERRQ(ierr); 1285 nzk += nlnk; 1286 1287 /* update il and jl for prow */ 1288 if (jmin < jmax) { 1289 il[prow] = jmin; 1290 j = *uj_ptr; 1291 jl[prow] = jl[j]; 1292 jl[j] = prow; 1293 } 1294 prow = nextprow; 1295 } 1296 1297 /* if free space is not available, make more free space */ 1298 if (current_space->local_remaining<nzk) { 1299 i = mbs - k + 1; /* num of unfactored rows */ 1300 i = PetscMin(PetscIntMultTruncate(i,nzk), PetscIntMultTruncate(i,i-1)); /* i*nzk, i*(i-1): estimated and max additional space needed */ 1301 ierr = PetscFreeSpaceGet(i,¤t_space);CHKERRQ(ierr); 1302 reallocs++; 1303 } 1304 1305 /* copy data into free space, then initialize lnk */ 1306 ierr = PetscLLClean(mbs,mbs,nzk,lnk,current_space->array,lnkbt);CHKERRQ(ierr); 1307 1308 /* add the k-th row into il and jl */ 1309 if (nzk-1 > 0) { 1310 i = current_space->array[1]; /* col value of the first nonzero element in U(k, k+1:mbs-1) */ 1311 jl[k] = jl[i]; jl[i] = k; 1312 il[k] = ui[k] + 1; 1313 } 1314 ui_ptr[k] = current_space->array; 1315 current_space->array += nzk; 1316 current_space->local_used += nzk; 1317 current_space->local_remaining -= nzk; 1318 1319 ui[k+1] = ui[k] + nzk; 1320 } 1321 1322 ierr = ISRestoreIndices(perm,&rip);CHKERRQ(ierr); 1323 ierr = PetscFree4(ui_ptr,il,jl,cols);CHKERRQ(ierr); 1324 1325 /* copy free_space into uj and free free_space; set uj in new datastructure; */ 1326 ierr = PetscMalloc1(ui[mbs]+1,&uj);CHKERRQ(ierr); 1327 ierr = PetscFreeSpaceContiguous(&free_space,uj);CHKERRQ(ierr); 1328 ierr = PetscLLDestroy(lnk,lnkbt);CHKERRQ(ierr); 1329 1330 /* put together the new matrix in MATSEQSBAIJ format */ 1331 B = fact; 1332 ierr = MatSeqSBAIJSetPreallocation(B,bs,MAT_SKIP_ALLOCATION,NULL);CHKERRQ(ierr); 1333 1334 b = (Mat_SeqSBAIJ*)B->data; 1335 b->singlemalloc = PETSC_FALSE; 1336 b->free_a = PETSC_TRUE; 1337 b->free_ij = PETSC_TRUE; 1338 1339 ierr = PetscMalloc1(ui[mbs]+1,&b->a);CHKERRQ(ierr); 1340 1341 b->j = uj; 1342 b->i = ui; 1343 b->diag = NULL; 1344 b->ilen = NULL; 1345 b->imax = NULL; 1346 b->row = perm; 1347 b->pivotinblocks = PETSC_FALSE; /* need to get from MatFactorInfo */ 1348 1349 ierr = PetscObjectReference((PetscObject)perm);CHKERRQ(ierr); 1350 b->icol = perm; 1351 ierr = PetscObjectReference((PetscObject)perm);CHKERRQ(ierr); 1352 ierr = PetscMalloc1(mbs+1,&b->solve_work);CHKERRQ(ierr); 1353 ierr = PetscLogObjectMemory((PetscObject)B,(ui[mbs]-mbs)*(sizeof(PetscInt)+sizeof(MatScalar)));CHKERRQ(ierr); 1354 b->maxnz = b->nz = ui[mbs]; 1355 1356 B->info.factor_mallocs = reallocs; 1357 B->info.fill_ratio_given = fill; 1358 if (ai[mbs] != 0.) { 1359 /* nonzeros in lower triangular part of A = (ai[mbs]+mbs)/2 */ 1360 B->info.fill_ratio_needed = ((PetscReal)2*ui[mbs])/(ai[mbs]+mbs); 1361 } else { 1362 B->info.fill_ratio_needed = 0.0; 1363 } 1364 #if defined(PETSC_USE_INFO) 1365 if (ai[mbs] != 0.) { 1366 PetscReal af = B->info.fill_ratio_needed; 1367 ierr = PetscInfo3(A,"Reallocs %D Fill ratio:given %g needed %g\n",reallocs,(double)fill,(double)af);CHKERRQ(ierr); 1368 ierr = PetscInfo1(A,"Run with -pc_factor_fill %g or use \n",(double)af);CHKERRQ(ierr); 1369 ierr = PetscInfo1(A,"PCFactorSetFill(pc,%g) for best performance.\n",(double)af);CHKERRQ(ierr); 1370 } else { 1371 ierr = PetscInfo(A,"Empty matrix\n");CHKERRQ(ierr); 1372 } 1373 #endif 1374 if (perm_identity) { 1375 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N_NaturalOrdering; 1376 } else { 1377 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqBAIJ_N; 1378 } 1379 PetscFunctionReturn(0); 1380 } 1381 1382 PetscErrorCode MatSolve_SeqBAIJ_N_NaturalOrdering(Mat A,Vec bb,Vec xx) 1383 { 1384 Mat_SeqBAIJ *a=(Mat_SeqBAIJ*)A->data; 1385 PetscErrorCode ierr; 1386 const PetscInt *ai=a->i,*aj=a->j,*adiag=a->diag,*vi; 1387 PetscInt i,k,n=a->mbs; 1388 PetscInt nz,bs=A->rmap->bs,bs2=a->bs2; 1389 const MatScalar *aa=a->a,*v; 1390 PetscScalar *x,*s,*t,*ls; 1391 const PetscScalar *b; 1392 1393 PetscFunctionBegin; 1394 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 1395 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1396 t = a->solve_work; 1397 1398 /* forward solve the lower triangular */ 1399 ierr = PetscArraycpy(t,b,bs);CHKERRQ(ierr); /* copy 1st block of b to t */ 1400 1401 for (i=1; i<n; i++) { 1402 v = aa + bs2*ai[i]; 1403 vi = aj + ai[i]; 1404 nz = ai[i+1] - ai[i]; 1405 s = t + bs*i; 1406 ierr = PetscArraycpy(s,b+bs*i,bs);CHKERRQ(ierr); /* copy i_th block of b to t */ 1407 for (k=0;k<nz;k++) { 1408 PetscKernel_v_gets_v_minus_A_times_w(bs,s,v,t+bs*vi[k]); 1409 v += bs2; 1410 } 1411 } 1412 1413 /* backward solve the upper triangular */ 1414 ls = a->solve_work + A->cmap->n; 1415 for (i=n-1; i>=0; i--) { 1416 v = aa + bs2*(adiag[i+1]+1); 1417 vi = aj + adiag[i+1]+1; 1418 nz = adiag[i] - adiag[i+1]-1; 1419 ierr = PetscArraycpy(ls,t+i*bs,bs);CHKERRQ(ierr); 1420 for (k=0; k<nz; k++) { 1421 PetscKernel_v_gets_v_minus_A_times_w(bs,ls,v,t+bs*vi[k]); 1422 v += bs2; 1423 } 1424 PetscKernel_w_gets_A_times_v(bs,ls,aa+bs2*adiag[i],t+i*bs); /* *inv(diagonal[i]) */ 1425 ierr = PetscArraycpy(x+i*bs,t+i*bs,bs);CHKERRQ(ierr); 1426 } 1427 1428 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1429 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1430 ierr = PetscLogFlops(2.0*(a->bs2)*(a->nz) - A->rmap->bs*A->cmap->n);CHKERRQ(ierr); 1431 PetscFunctionReturn(0); 1432 } 1433 1434 PetscErrorCode MatSolve_SeqBAIJ_N(Mat A,Vec bb,Vec xx) 1435 { 1436 Mat_SeqBAIJ *a =(Mat_SeqBAIJ*)A->data; 1437 IS iscol=a->col,isrow=a->row; 1438 PetscErrorCode ierr; 1439 const PetscInt *r,*c,*rout,*cout,*ai=a->i,*aj=a->j,*adiag=a->diag,*vi; 1440 PetscInt i,m,n=a->mbs; 1441 PetscInt nz,bs=A->rmap->bs,bs2=a->bs2; 1442 const MatScalar *aa=a->a,*v; 1443 PetscScalar *x,*s,*t,*ls; 1444 const PetscScalar *b; 1445 1446 PetscFunctionBegin; 1447 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 1448 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1449 t = a->solve_work; 1450 1451 ierr = ISGetIndices(isrow,&rout);CHKERRQ(ierr); r = rout; 1452 ierr = ISGetIndices(iscol,&cout);CHKERRQ(ierr); c = cout; 1453 1454 /* forward solve the lower triangular */ 1455 ierr = PetscArraycpy(t,b+bs*r[0],bs);CHKERRQ(ierr); 1456 for (i=1; i<n; i++) { 1457 v = aa + bs2*ai[i]; 1458 vi = aj + ai[i]; 1459 nz = ai[i+1] - ai[i]; 1460 s = t + bs*i; 1461 ierr = PetscArraycpy(s,b+bs*r[i],bs);CHKERRQ(ierr); 1462 for (m=0; m<nz; m++) { 1463 PetscKernel_v_gets_v_minus_A_times_w(bs,s,v,t+bs*vi[m]); 1464 v += bs2; 1465 } 1466 } 1467 1468 /* backward solve the upper triangular */ 1469 ls = a->solve_work + A->cmap->n; 1470 for (i=n-1; i>=0; i--) { 1471 v = aa + bs2*(adiag[i+1]+1); 1472 vi = aj + adiag[i+1]+1; 1473 nz = adiag[i] - adiag[i+1] - 1; 1474 ierr = PetscArraycpy(ls,t+i*bs,bs);CHKERRQ(ierr); 1475 for (m=0; m<nz; m++) { 1476 PetscKernel_v_gets_v_minus_A_times_w(bs,ls,v,t+bs*vi[m]); 1477 v += bs2; 1478 } 1479 PetscKernel_w_gets_A_times_v(bs,ls,v,t+i*bs); /* *inv(diagonal[i]) */ 1480 ierr = PetscArraycpy(x + bs*c[i],t+i*bs,bs);CHKERRQ(ierr); 1481 } 1482 ierr = ISRestoreIndices(isrow,&rout);CHKERRQ(ierr); 1483 ierr = ISRestoreIndices(iscol,&cout);CHKERRQ(ierr); 1484 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1485 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1486 ierr = PetscLogFlops(2.0*(a->bs2)*(a->nz) - A->rmap->bs*A->cmap->n);CHKERRQ(ierr); 1487 PetscFunctionReturn(0); 1488 } 1489 1490 /* 1491 For each block in an block array saves the largest absolute value in the block into another array 1492 */ 1493 static PetscErrorCode MatBlockAbs_private(PetscInt nbs,PetscInt bs2,PetscScalar *blockarray,PetscReal *absarray) 1494 { 1495 PetscErrorCode ierr; 1496 PetscInt i,j; 1497 1498 PetscFunctionBegin; 1499 ierr = PetscArrayzero(absarray,nbs+1);CHKERRQ(ierr); 1500 for (i=0; i<nbs; i++) { 1501 for (j=0; j<bs2; j++) { 1502 if (absarray[i] < PetscAbsScalar(blockarray[i*nbs+j])) absarray[i] = PetscAbsScalar(blockarray[i*nbs+j]); 1503 } 1504 } 1505 PetscFunctionReturn(0); 1506 } 1507 1508 /* 1509 This needs to be renamed and called by the regular MatILUFactor_SeqBAIJ when drop tolerance is used 1510 */ 1511 PetscErrorCode MatILUDTFactor_SeqBAIJ(Mat A,IS isrow,IS iscol,const MatFactorInfo *info,Mat *fact) 1512 { 1513 Mat B = *fact; 1514 Mat_SeqBAIJ *a=(Mat_SeqBAIJ*)A->data,*b; 1515 IS isicol; 1516 PetscErrorCode ierr; 1517 const PetscInt *r,*ic; 1518 PetscInt i,mbs=a->mbs,bs=A->rmap->bs,bs2=a->bs2,*ai=a->i,*aj=a->j,*ajtmp,*adiag; 1519 PetscInt *bi,*bj,*bdiag; 1520 1521 PetscInt row,nzi,nzi_bl,nzi_bu,*im,dtcount,nzi_al,nzi_au; 1522 PetscInt nlnk,*lnk; 1523 PetscBT lnkbt; 1524 PetscBool row_identity,icol_identity; 1525 MatScalar *aatmp,*pv,*batmp,*ba,*rtmp,*pc,*multiplier,*vtmp; 1526 PetscInt j,nz,*pj,*bjtmp,k,ncut,*jtmp; 1527 1528 PetscReal dt=info->dt; /* shift=info->shiftamount; */ 1529 PetscInt nnz_max; 1530 PetscBool missing; 1531 PetscReal *vtmp_abs; 1532 MatScalar *v_work; 1533 PetscInt *v_pivots; 1534 PetscBool allowzeropivot,zeropivotdetected=PETSC_FALSE; 1535 1536 PetscFunctionBegin; 1537 /* ------- symbolic factorization, can be reused ---------*/ 1538 ierr = MatMissingDiagonal(A,&missing,&i);CHKERRQ(ierr); 1539 if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",i); 1540 adiag=a->diag; 1541 1542 ierr = ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);CHKERRQ(ierr); 1543 1544 /* bdiag is location of diagonal in factor */ 1545 ierr = PetscMalloc1(mbs+1,&bdiag);CHKERRQ(ierr); 1546 1547 /* allocate row pointers bi */ 1548 ierr = PetscMalloc1(2*mbs+2,&bi);CHKERRQ(ierr); 1549 1550 /* allocate bj and ba; max num of nonzero entries is (ai[n]+2*n*dtcount+2) */ 1551 dtcount = (PetscInt)info->dtcount; 1552 if (dtcount > mbs-1) dtcount = mbs-1; 1553 nnz_max = ai[mbs]+2*mbs*dtcount +2; 1554 /* printf("MatILUDTFactor_SeqBAIJ, bs %d, ai[mbs] %d, nnz_max %d, dtcount %d\n",bs,ai[mbs],nnz_max,dtcount); */ 1555 ierr = PetscMalloc1(nnz_max,&bj);CHKERRQ(ierr); 1556 nnz_max = nnz_max*bs2; 1557 ierr = PetscMalloc1(nnz_max,&ba);CHKERRQ(ierr); 1558 1559 /* put together the new matrix */ 1560 ierr = MatSeqBAIJSetPreallocation(B,bs,MAT_SKIP_ALLOCATION,NULL);CHKERRQ(ierr); 1561 ierr = PetscLogObjectParent((PetscObject)B,(PetscObject)isicol);CHKERRQ(ierr); 1562 1563 b = (Mat_SeqBAIJ*)(B)->data; 1564 b->free_a = PETSC_TRUE; 1565 b->free_ij = PETSC_TRUE; 1566 b->singlemalloc = PETSC_FALSE; 1567 1568 b->a = ba; 1569 b->j = bj; 1570 b->i = bi; 1571 b->diag = bdiag; 1572 b->ilen = NULL; 1573 b->imax = NULL; 1574 b->row = isrow; 1575 b->col = iscol; 1576 1577 ierr = PetscObjectReference((PetscObject)isrow);CHKERRQ(ierr); 1578 ierr = PetscObjectReference((PetscObject)iscol);CHKERRQ(ierr); 1579 1580 b->icol = isicol; 1581 ierr = PetscMalloc1(bs*(mbs+1),&b->solve_work);CHKERRQ(ierr); 1582 ierr = PetscLogObjectMemory((PetscObject)B,nnz_max*(sizeof(PetscInt)+sizeof(MatScalar)));CHKERRQ(ierr); 1583 b->maxnz = nnz_max/bs2; 1584 1585 (B)->factortype = MAT_FACTOR_ILUDT; 1586 (B)->info.factor_mallocs = 0; 1587 (B)->info.fill_ratio_given = ((PetscReal)nnz_max)/((PetscReal)(ai[mbs]*bs2)); 1588 /* ------- end of symbolic factorization ---------*/ 1589 ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr); 1590 ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr); 1591 1592 /* linked list for storing column indices of the active row */ 1593 nlnk = mbs + 1; 1594 ierr = PetscLLCreate(mbs,mbs,nlnk,lnk,lnkbt);CHKERRQ(ierr); 1595 1596 /* im: used by PetscLLAddSortedLU(); jtmp: working array for column indices of active row */ 1597 ierr = PetscMalloc2(mbs,&im,mbs,&jtmp);CHKERRQ(ierr); 1598 /* rtmp, vtmp: working arrays for sparse and contiguous row entries of active row */ 1599 ierr = PetscMalloc2(mbs*bs2,&rtmp,mbs*bs2,&vtmp);CHKERRQ(ierr); 1600 ierr = PetscMalloc1(mbs+1,&vtmp_abs);CHKERRQ(ierr); 1601 ierr = PetscMalloc3(bs,&v_work,bs2,&multiplier,bs,&v_pivots);CHKERRQ(ierr); 1602 1603 allowzeropivot = PetscNot(A->erroriffailure); 1604 bi[0] = 0; 1605 bdiag[0] = (nnz_max/bs2)-1; /* location of diagonal in factor B */ 1606 bi[2*mbs+1] = bdiag[0]+1; /* endof bj and ba array */ 1607 for (i=0; i<mbs; i++) { 1608 /* copy initial fill into linked list */ 1609 nzi = ai[r[i]+1] - ai[r[i]]; 1610 if (!nzi) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Empty row in matrix: row in original ordering %D in permuted ordering %D",r[i],i); 1611 nzi_al = adiag[r[i]] - ai[r[i]]; 1612 nzi_au = ai[r[i]+1] - adiag[r[i]] -1; 1613 1614 /* load in initial unfactored row */ 1615 ajtmp = aj + ai[r[i]]; 1616 ierr = PetscLLAddPerm(nzi,ajtmp,ic,mbs,nlnk,lnk,lnkbt);CHKERRQ(ierr); 1617 ierr = PetscArrayzero(rtmp,mbs*bs2);CHKERRQ(ierr); 1618 aatmp = a->a + bs2*ai[r[i]]; 1619 for (j=0; j<nzi; j++) { 1620 ierr = PetscArraycpy(rtmp+bs2*ic[ajtmp[j]],aatmp+bs2*j,bs2);CHKERRQ(ierr); 1621 } 1622 1623 /* add pivot rows into linked list */ 1624 row = lnk[mbs]; 1625 while (row < i) { 1626 nzi_bl = bi[row+1] - bi[row] + 1; 1627 bjtmp = bj + bdiag[row+1]+1; /* points to 1st column next to the diagonal in U */ 1628 ierr = PetscLLAddSortedLU(bjtmp,row,nlnk,lnk,lnkbt,i,nzi_bl,im);CHKERRQ(ierr); 1629 nzi += nlnk; 1630 row = lnk[row]; 1631 } 1632 1633 /* copy data from lnk into jtmp, then initialize lnk */ 1634 ierr = PetscLLClean(mbs,mbs,nzi,lnk,jtmp,lnkbt);CHKERRQ(ierr); 1635 1636 /* numerical factorization */ 1637 bjtmp = jtmp; 1638 row = *bjtmp++; /* 1st pivot row */ 1639 1640 while (row < i) { 1641 pc = rtmp + bs2*row; 1642 pv = ba + bs2*bdiag[row]; /* inv(diag) of the pivot row */ 1643 PetscKernel_A_gets_A_times_B(bs,pc,pv,multiplier); /* pc= multiplier = pc*inv(diag[row]) */ 1644 ierr = MatBlockAbs_private(1,bs2,pc,vtmp_abs);CHKERRQ(ierr); 1645 if (vtmp_abs[0] > dt) { /* apply tolerance dropping rule */ 1646 pj = bj + bdiag[row+1] + 1; /* point to 1st entry of U(row,:) */ 1647 pv = ba + bs2*(bdiag[row+1] + 1); 1648 nz = bdiag[row] - bdiag[row+1] - 1; /* num of entries in U(row,:), excluding diagonal */ 1649 for (j=0; j<nz; j++) { 1650 PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j); 1651 } 1652 /* ierr = PetscLogFlops(bslog*(nz+1.0)-bs);CHKERRQ(ierr); */ 1653 } 1654 row = *bjtmp++; 1655 } 1656 1657 /* copy sparse rtmp into contiguous vtmp; separate L and U part */ 1658 nzi_bl = 0; j = 0; 1659 while (jtmp[j] < i) { /* L-part. Note: jtmp is sorted */ 1660 ierr = PetscArraycpy(vtmp+bs2*j,rtmp+bs2*jtmp[j],bs2);CHKERRQ(ierr); 1661 nzi_bl++; j++; 1662 } 1663 nzi_bu = nzi - nzi_bl -1; 1664 1665 while (j < nzi) { /* U-part */ 1666 ierr = PetscArraycpy(vtmp+bs2*j,rtmp+bs2*jtmp[j],bs2);CHKERRQ(ierr); 1667 j++; 1668 } 1669 1670 ierr = MatBlockAbs_private(nzi,bs2,vtmp,vtmp_abs);CHKERRQ(ierr); 1671 1672 bjtmp = bj + bi[i]; 1673 batmp = ba + bs2*bi[i]; 1674 /* apply level dropping rule to L part */ 1675 ncut = nzi_al + dtcount; 1676 if (ncut < nzi_bl) { 1677 ierr = PetscSortSplitReal(ncut,nzi_bl,vtmp_abs,jtmp);CHKERRQ(ierr); 1678 ierr = PetscSortIntWithScalarArray(ncut,jtmp,vtmp);CHKERRQ(ierr); 1679 } else { 1680 ncut = nzi_bl; 1681 } 1682 for (j=0; j<ncut; j++) { 1683 bjtmp[j] = jtmp[j]; 1684 ierr = PetscArraycpy(batmp+bs2*j,rtmp+bs2*bjtmp[j],bs2);CHKERRQ(ierr); 1685 } 1686 bi[i+1] = bi[i] + ncut; 1687 nzi = ncut + 1; 1688 1689 /* apply level dropping rule to U part */ 1690 ncut = nzi_au + dtcount; 1691 if (ncut < nzi_bu) { 1692 ierr = PetscSortSplitReal(ncut,nzi_bu,vtmp_abs+nzi_bl+1,jtmp+nzi_bl+1);CHKERRQ(ierr); 1693 ierr = PetscSortIntWithScalarArray(ncut,jtmp+nzi_bl+1,vtmp+nzi_bl+1);CHKERRQ(ierr); 1694 } else { 1695 ncut = nzi_bu; 1696 } 1697 nzi += ncut; 1698 1699 /* mark bdiagonal */ 1700 bdiag[i+1] = bdiag[i] - (ncut + 1); 1701 bi[2*mbs - i] = bi[2*mbs - i +1] - (ncut + 1); 1702 1703 bjtmp = bj + bdiag[i]; 1704 batmp = ba + bs2*bdiag[i]; 1705 ierr = PetscArraycpy(batmp,rtmp+bs2*i,bs2);CHKERRQ(ierr); 1706 *bjtmp = i; 1707 1708 bjtmp = bj + bdiag[i+1]+1; 1709 batmp = ba + (bdiag[i+1]+1)*bs2; 1710 1711 for (k=0; k<ncut; k++) { 1712 bjtmp[k] = jtmp[nzi_bl+1+k]; 1713 ierr = PetscArraycpy(batmp+bs2*k,rtmp+bs2*bjtmp[k],bs2);CHKERRQ(ierr); 1714 } 1715 1716 im[i] = nzi; /* used by PetscLLAddSortedLU() */ 1717 1718 /* invert diagonal block for simplier triangular solves - add shift??? */ 1719 batmp = ba + bs2*bdiag[i]; 1720 1721 ierr = PetscKernel_A_gets_inverse_A(bs,batmp,v_pivots,v_work,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 1722 if (zeropivotdetected) B->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1723 } /* for (i=0; i<mbs; i++) */ 1724 ierr = PetscFree3(v_work,multiplier,v_pivots);CHKERRQ(ierr); 1725 1726 /* printf("end of L %d, beginning of U %d\n",bi[mbs],bdiag[mbs]); */ 1727 if (bi[mbs] >= bdiag[mbs]) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"end of L array %d cannot >= the beginning of U array %d",bi[mbs],bdiag[mbs]); 1728 1729 ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr); 1730 ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr); 1731 1732 ierr = PetscLLDestroy(lnk,lnkbt);CHKERRQ(ierr); 1733 1734 ierr = PetscFree2(im,jtmp);CHKERRQ(ierr); 1735 ierr = PetscFree2(rtmp,vtmp);CHKERRQ(ierr); 1736 1737 ierr = PetscLogFlops(bs2*B->cmap->n);CHKERRQ(ierr); 1738 b->maxnz = b->nz = bi[mbs] + bdiag[0] - bdiag[mbs]; 1739 1740 ierr = ISIdentity(isrow,&row_identity);CHKERRQ(ierr); 1741 ierr = ISIdentity(isicol,&icol_identity);CHKERRQ(ierr); 1742 if (row_identity && icol_identity) { 1743 B->ops->solve = MatSolve_SeqBAIJ_N_NaturalOrdering; 1744 } else { 1745 B->ops->solve = MatSolve_SeqBAIJ_N; 1746 } 1747 1748 B->ops->solveadd = NULL; 1749 B->ops->solvetranspose = NULL; 1750 B->ops->solvetransposeadd = NULL; 1751 B->ops->matsolve = NULL; 1752 B->assembled = PETSC_TRUE; 1753 B->preallocated = PETSC_TRUE; 1754 PetscFunctionReturn(0); 1755 } 1756