1 2 /* 3 Factorization code for BAIJ format. 4 */ 5 6 #include <../src/mat/impls/baij/seq/baij.h> 7 #include <petsc/private/kernels/blockinvert.h> 8 #include <petscbt.h> 9 #include <../src/mat/utils/freespace.h> 10 11 /* ----------------------------------------------------------------*/ 12 extern PetscErrorCode MatDuplicateNoCreate_SeqBAIJ(Mat,Mat,MatDuplicateOption,PetscBool); 13 14 /* 15 This is not much faster than MatLUFactorNumeric_SeqBAIJ_N() but the solve is faster at least sometimes 16 */ 17 PetscErrorCode MatLUFactorNumeric_SeqBAIJ_15_NaturalOrdering(Mat B,Mat A,const MatFactorInfo *info) 18 { 19 Mat C =B; 20 Mat_SeqBAIJ *a=(Mat_SeqBAIJ*)A->data,*b=(Mat_SeqBAIJ*)C->data; 21 PetscInt i,j,k,ipvt[15]; 22 const PetscInt n=a->mbs,*ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j,*ajtmp,*bjtmp,*bdiag=b->diag,*pj; 23 PetscInt nz,nzL,row; 24 MatScalar *rtmp,*pc,*mwork,*pv,*vv,work[225]; 25 const MatScalar *v,*aa=a->a; 26 PetscInt bs2 = a->bs2,bs=A->rmap->bs,flg; 27 PetscInt sol_ver; 28 PetscBool allowzeropivot,zeropivotdetected; 29 30 PetscFunctionBegin; 31 allowzeropivot = PetscNot(A->erroriffailure); 32 PetscCall(PetscOptionsGetInt(NULL,((PetscObject)A)->prefix,"-sol_ver",&sol_ver,NULL)); 33 34 /* generate work space needed by the factorization */ 35 PetscCall(PetscMalloc2(bs2*n,&rtmp,bs2,&mwork)); 36 PetscCall(PetscArrayzero(rtmp,bs2*n)); 37 38 for (i=0; i<n; i++) { 39 /* zero rtmp */ 40 /* L part */ 41 nz = bi[i+1] - bi[i]; 42 bjtmp = bj + bi[i]; 43 for (j=0; j<nz; j++) { 44 PetscCall(PetscArrayzero(rtmp+bs2*bjtmp[j],bs2)); 45 } 46 47 /* U part */ 48 nz = bdiag[i] - bdiag[i+1]; 49 bjtmp = bj + bdiag[i+1]+1; 50 for (j=0; j<nz; j++) { 51 PetscCall(PetscArrayzero(rtmp+bs2*bjtmp[j],bs2)); 52 } 53 54 /* load in initial (unfactored row) */ 55 nz = ai[i+1] - ai[i]; 56 ajtmp = aj + ai[i]; 57 v = aa + bs2*ai[i]; 58 for (j=0; j<nz; j++) { 59 PetscCall(PetscArraycpy(rtmp+bs2*ajtmp[j],v+bs2*j,bs2)); 60 } 61 62 /* elimination */ 63 bjtmp = bj + bi[i]; 64 nzL = bi[i+1] - bi[i]; 65 for (k=0; k < nzL; k++) { 66 row = bjtmp[k]; 67 pc = rtmp + bs2*row; 68 for (flg=0,j=0; j<bs2; j++) { 69 if (pc[j]!=0.0) { 70 flg = 1; 71 break; 72 } 73 } 74 if (flg) { 75 pv = b->a + bs2*bdiag[row]; 76 PetscKernel_A_gets_A_times_B(bs,pc,pv,mwork); 77 /* PetscCall(PetscKernel_A_gets_A_times_B_15(pc,pv,mwork)); */ 78 pj = b->j + bdiag[row+1]+1; /* beginning of U(row,:) */ 79 pv = b->a + bs2*(bdiag[row+1]+1); 80 nz = bdiag[row] - bdiag[row+1] - 1; /* num of entries inU(row,:), excluding diag */ 81 for (j=0; j<nz; j++) { 82 vv = rtmp + bs2*pj[j]; 83 PetscKernel_A_gets_A_minus_B_times_C(bs,vv,pc,pv); 84 /* PetscCall(PetscKernel_A_gets_A_minus_B_times_C_15(vv,pc,pv)); */ 85 pv += bs2; 86 } 87 PetscCall(PetscLogFlops(2.0*bs2*bs*(nz+1)-bs2)); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */ 88 } 89 } 90 91 /* finished row so stick it into b->a */ 92 /* L part */ 93 pv = b->a + bs2*bi[i]; 94 pj = b->j + bi[i]; 95 nz = bi[i+1] - bi[i]; 96 for (j=0; j<nz; j++) { 97 PetscCall(PetscArraycpy(pv+bs2*j,rtmp+bs2*pj[j],bs2)); 98 } 99 100 /* Mark diagonal and invert diagonal for simpler triangular solves */ 101 pv = b->a + bs2*bdiag[i]; 102 pj = b->j + bdiag[i]; 103 PetscCall(PetscArraycpy(pv,rtmp+bs2*pj[0],bs2)); 104 PetscCall(PetscKernel_A_gets_inverse_A_15(pv,ipvt,work,info->shiftamount,allowzeropivot,&zeropivotdetected)); 105 if (zeropivotdetected) C->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 106 107 /* U part */ 108 pv = b->a + bs2*(bdiag[i+1]+1); 109 pj = b->j + bdiag[i+1]+1; 110 nz = bdiag[i] - bdiag[i+1] - 1; 111 for (j=0; j<nz; j++) { 112 PetscCall(PetscArraycpy(pv+bs2*j,rtmp+bs2*pj[j],bs2)); 113 } 114 } 115 116 PetscCall(PetscFree2(rtmp,mwork)); 117 118 C->ops->solve = MatSolve_SeqBAIJ_15_NaturalOrdering_ver1; 119 C->ops->solvetranspose = MatSolve_SeqBAIJ_N_NaturalOrdering; 120 C->assembled = PETSC_TRUE; 121 122 PetscCall(PetscLogFlops(1.333333333333*bs*bs2*b->mbs)); /* from inverting diagonal blocks */ 123 PetscFunctionReturn(0); 124 } 125 126 PetscErrorCode MatLUFactorNumeric_SeqBAIJ_N(Mat B,Mat A,const MatFactorInfo *info) 127 { 128 Mat C =B; 129 Mat_SeqBAIJ *a =(Mat_SeqBAIJ*)A->data,*b=(Mat_SeqBAIJ*)C->data; 130 IS isrow = b->row,isicol = b->icol; 131 const PetscInt *r,*ic; 132 PetscInt i,j,k,n=a->mbs,*ai=a->i,*aj=a->j,*bi=b->i,*bj=b->j; 133 PetscInt *ajtmp,*bjtmp,nz,nzL,row,*bdiag=b->diag,*pj; 134 MatScalar *rtmp,*pc,*mwork,*v,*pv,*aa=a->a; 135 PetscInt bs=A->rmap->bs,bs2 = a->bs2,*v_pivots,flg; 136 MatScalar *v_work; 137 PetscBool col_identity,row_identity,both_identity; 138 PetscBool allowzeropivot,zeropivotdetected; 139 140 PetscFunctionBegin; 141 PetscCall(ISGetIndices(isrow,&r)); 142 PetscCall(ISGetIndices(isicol,&ic)); 143 allowzeropivot = PetscNot(A->erroriffailure); 144 145 PetscCall(PetscCalloc1(bs2*n,&rtmp)); 146 147 /* generate work space needed by dense LU factorization */ 148 PetscCall(PetscMalloc3(bs,&v_work,bs2,&mwork,bs,&v_pivots)); 149 150 for (i=0; i<n; i++) { 151 /* zero rtmp */ 152 /* L part */ 153 nz = bi[i+1] - bi[i]; 154 bjtmp = bj + bi[i]; 155 for (j=0; j<nz; j++) { 156 PetscCall(PetscArrayzero(rtmp+bs2*bjtmp[j],bs2)); 157 } 158 159 /* U part */ 160 nz = bdiag[i] - bdiag[i+1]; 161 bjtmp = bj + bdiag[i+1]+1; 162 for (j=0; j<nz; j++) { 163 PetscCall(PetscArrayzero(rtmp+bs2*bjtmp[j],bs2)); 164 } 165 166 /* load in initial (unfactored row) */ 167 nz = ai[r[i]+1] - ai[r[i]]; 168 ajtmp = aj + ai[r[i]]; 169 v = aa + bs2*ai[r[i]]; 170 for (j=0; j<nz; j++) { 171 PetscCall(PetscArraycpy(rtmp+bs2*ic[ajtmp[j]],v+bs2*j,bs2)); 172 } 173 174 /* elimination */ 175 bjtmp = bj + bi[i]; 176 nzL = bi[i+1] - bi[i]; 177 for (k=0; k < nzL; k++) { 178 row = bjtmp[k]; 179 pc = rtmp + bs2*row; 180 for (flg=0,j=0; j<bs2; j++) { 181 if (pc[j]!=0.0) { 182 flg = 1; 183 break; 184 } 185 } 186 if (flg) { 187 pv = b->a + bs2*bdiag[row]; 188 PetscKernel_A_gets_A_times_B(bs,pc,pv,mwork); /* *pc = *pc * (*pv); */ 189 pj = b->j + bdiag[row+1]+1; /* beginning of U(row,:) */ 190 pv = b->a + bs2*(bdiag[row+1]+1); 191 nz = bdiag[row] - bdiag[row+1] - 1; /* num of entries inU(row,:), excluding diag */ 192 for (j=0; j<nz; j++) { 193 PetscKernel_A_gets_A_minus_B_times_C(bs,rtmp+bs2*pj[j],pc,pv+bs2*j); 194 } 195 PetscCall(PetscLogFlops(2.0*bs2*bs*(nz+1)-bs2)); /* flops = 2*bs^3*nz + 2*bs^3 - bs2) */ 196 } 197 } 198 199 /* finished row so stick it into b->a */ 200 /* L part */ 201 pv = b->a + bs2*bi[i]; 202 pj = b->j + bi[i]; 203 nz = bi[i+1] - bi[i]; 204 for (j=0; j<nz; j++) { 205 PetscCall(PetscArraycpy(pv+bs2*j,rtmp+bs2*pj[j],bs2)); 206 } 207 208 /* Mark diagonal and invert diagonal for simpler triangular solves */ 209 pv = b->a + bs2*bdiag[i]; 210 pj = b->j + bdiag[i]; 211 PetscCall(PetscArraycpy(pv,rtmp+bs2*pj[0],bs2)); 212 213 PetscCall(PetscKernel_A_gets_inverse_A(bs,pv,v_pivots,v_work,allowzeropivot,&zeropivotdetected)); 214 if (zeropivotdetected) B->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 215 216 /* U part */ 217 pv = b->a + bs2*(bdiag[i+1]+1); 218 pj = b->j + bdiag[i+1]+1; 219 nz = bdiag[i] - bdiag[i+1] - 1; 220 for (j=0; j<nz; j++) { 221 PetscCall(PetscArraycpy(pv+bs2*j,rtmp+bs2*pj[j],bs2)); 222 } 223 } 224 225 PetscCall(PetscFree(rtmp)); 226 PetscCall(PetscFree3(v_work,mwork,v_pivots)); 227 PetscCall(ISRestoreIndices(isicol,&ic)); 228 PetscCall(ISRestoreIndices(isrow,&r)); 229 230 PetscCall(ISIdentity(isrow,&row_identity)); 231 PetscCall(ISIdentity(isicol,&col_identity)); 232 233 both_identity = (PetscBool) (row_identity && col_identity); 234 if (both_identity) { 235 switch (bs) { 236 case 9: 237 #if defined(PETSC_HAVE_IMMINTRIN_H) && defined(__AVX2__) && defined(__FMA__) && defined(PETSC_USE_REAL_DOUBLE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_64BIT_INDICES) 238 C->ops->solve = MatSolve_SeqBAIJ_9_NaturalOrdering; 239 #else 240 C->ops->solve = MatSolve_SeqBAIJ_N_NaturalOrdering; 241 #endif 242 break; 243 case 11: 244 C->ops->solve = MatSolve_SeqBAIJ_11_NaturalOrdering; 245 break; 246 case 12: 247 C->ops->solve = MatSolve_SeqBAIJ_12_NaturalOrdering; 248 break; 249 case 13: 250 C->ops->solve = MatSolve_SeqBAIJ_13_NaturalOrdering; 251 break; 252 case 14: 253 C->ops->solve = MatSolve_SeqBAIJ_14_NaturalOrdering; 254 break; 255 default: 256 C->ops->solve = MatSolve_SeqBAIJ_N_NaturalOrdering; 257 break; 258 } 259 } else { 260 C->ops->solve = MatSolve_SeqBAIJ_N; 261 } 262 C->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_N; 263 264 C->assembled = PETSC_TRUE; 265 266 PetscCall(PetscLogFlops(1.333333333333*bs*bs2*b->mbs)); /* from inverting diagonal blocks */ 267 PetscFunctionReturn(0); 268 } 269 270 /* 271 ilu(0) with natural ordering under new data structure. 272 See MatILUFactorSymbolic_SeqAIJ_ilu0() for detailed description 273 because this code is almost identical to MatILUFactorSymbolic_SeqAIJ_ilu0_inplace(). 274 */ 275 276 PetscErrorCode MatILUFactorSymbolic_SeqBAIJ_ilu0(Mat fact,Mat A,IS isrow,IS iscol,const MatFactorInfo *info) 277 { 278 279 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b; 280 PetscInt n=a->mbs,*ai=a->i,*aj,*adiag=a->diag,bs2 = a->bs2; 281 PetscInt i,j,nz,*bi,*bj,*bdiag,bi_temp; 282 283 PetscFunctionBegin; 284 PetscCall(MatDuplicateNoCreate_SeqBAIJ(fact,A,MAT_DO_NOT_COPY_VALUES,PETSC_FALSE)); 285 b = (Mat_SeqBAIJ*)(fact)->data; 286 287 /* allocate matrix arrays for new data structure */ 288 PetscCall(PetscMalloc3(bs2*ai[n]+1,&b->a,ai[n]+1,&b->j,n+1,&b->i)); 289 PetscCall(PetscLogObjectMemory((PetscObject)fact,ai[n]*(bs2*sizeof(PetscScalar)+sizeof(PetscInt))+(n+1)*sizeof(PetscInt))); 290 291 b->singlemalloc = PETSC_TRUE; 292 b->free_a = PETSC_TRUE; 293 b->free_ij = PETSC_TRUE; 294 fact->preallocated = PETSC_TRUE; 295 fact->assembled = PETSC_TRUE; 296 if (!b->diag) { 297 PetscCall(PetscMalloc1(n+1,&b->diag)); 298 PetscCall(PetscLogObjectMemory((PetscObject)fact,(n+1)*sizeof(PetscInt))); 299 } 300 bdiag = b->diag; 301 302 if (n > 0) { 303 PetscCall(PetscArrayzero(b->a,bs2*ai[n])); 304 } 305 306 /* set bi and bj with new data structure */ 307 bi = b->i; 308 bj = b->j; 309 310 /* L part */ 311 bi[0] = 0; 312 for (i=0; i<n; i++) { 313 nz = adiag[i] - ai[i]; 314 bi[i+1] = bi[i] + nz; 315 aj = a->j + ai[i]; 316 for (j=0; j<nz; j++) { 317 *bj = aj[j]; bj++; 318 } 319 } 320 321 /* U part */ 322 bi_temp = bi[n]; 323 bdiag[n] = bi[n]-1; 324 for (i=n-1; i>=0; i--) { 325 nz = ai[i+1] - adiag[i] - 1; 326 bi_temp = bi_temp + nz + 1; 327 aj = a->j + adiag[i] + 1; 328 for (j=0; j<nz; j++) { 329 *bj = aj[j]; bj++; 330 } 331 /* diag[i] */ 332 *bj = i; bj++; 333 bdiag[i] = bi_temp - 1; 334 } 335 PetscFunctionReturn(0); 336 } 337 338 PetscErrorCode MatILUFactorSymbolic_SeqBAIJ(Mat fact,Mat A,IS isrow,IS iscol,const MatFactorInfo *info) 339 { 340 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b; 341 IS isicol; 342 const PetscInt *r,*ic; 343 PetscInt n=a->mbs,*ai=a->i,*aj=a->j,d; 344 PetscInt *bi,*cols,nnz,*cols_lvl; 345 PetscInt *bdiag,prow,fm,nzbd,reallocs=0,dcount=0; 346 PetscInt i,levels,diagonal_fill; 347 PetscBool col_identity,row_identity,both_identity; 348 PetscReal f; 349 PetscInt nlnk,*lnk,*lnk_lvl=NULL; 350 PetscBT lnkbt; 351 PetscInt nzi,*bj,**bj_ptr,**bjlvl_ptr; 352 PetscFreeSpaceList free_space =NULL,current_space=NULL; 353 PetscFreeSpaceList free_space_lvl=NULL,current_space_lvl=NULL; 354 PetscBool missing; 355 PetscInt bs=A->rmap->bs,bs2=a->bs2; 356 357 PetscFunctionBegin; 358 PetscCheckFalse(A->rmap->n != A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Must be square matrix, rows %" PetscInt_FMT " columns %" PetscInt_FMT,A->rmap->n,A->cmap->n); 359 if (bs>1) { /* check shifttype */ 360 PetscCheckFalse(info->shifttype == MAT_SHIFT_NONZERO || info->shifttype == MAT_SHIFT_POSITIVE_DEFINITE,PETSC_COMM_SELF,PETSC_ERR_SUP,"Only MAT_SHIFT_NONE and MAT_SHIFT_INBLOCKS are supported for BAIJ matrix"); 361 } 362 363 PetscCall(MatMissingDiagonal(A,&missing,&d)); 364 PetscCheck(!missing,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %" PetscInt_FMT,d); 365 366 f = info->fill; 367 levels = (PetscInt)info->levels; 368 diagonal_fill = (PetscInt)info->diagonal_fill; 369 370 PetscCall(ISInvertPermutation(iscol,PETSC_DECIDE,&isicol)); 371 372 PetscCall(ISIdentity(isrow,&row_identity)); 373 PetscCall(ISIdentity(iscol,&col_identity)); 374 375 both_identity = (PetscBool) (row_identity && col_identity); 376 377 if (!levels && both_identity) { 378 /* special case: ilu(0) with natural ordering */ 379 PetscCall(MatILUFactorSymbolic_SeqBAIJ_ilu0(fact,A,isrow,iscol,info)); 380 PetscCall(MatSeqBAIJSetNumericFactorization(fact,both_identity)); 381 382 fact->factortype = MAT_FACTOR_ILU; 383 (fact)->info.factor_mallocs = 0; 384 (fact)->info.fill_ratio_given = info->fill; 385 (fact)->info.fill_ratio_needed = 1.0; 386 387 b = (Mat_SeqBAIJ*)(fact)->data; 388 b->row = isrow; 389 b->col = iscol; 390 b->icol = isicol; 391 PetscCall(PetscObjectReference((PetscObject)isrow)); 392 PetscCall(PetscObjectReference((PetscObject)iscol)); 393 b->pivotinblocks = (info->pivotinblocks) ? PETSC_TRUE : PETSC_FALSE; 394 395 PetscCall(PetscMalloc1((n+1)*bs,&b->solve_work)); 396 PetscFunctionReturn(0); 397 } 398 399 PetscCall(ISGetIndices(isrow,&r)); 400 PetscCall(ISGetIndices(isicol,&ic)); 401 402 /* get new row pointers */ 403 PetscCall(PetscMalloc1(n+1,&bi)); 404 bi[0] = 0; 405 /* bdiag is location of diagonal in factor */ 406 PetscCall(PetscMalloc1(n+1,&bdiag)); 407 bdiag[0] = 0; 408 409 PetscCall(PetscMalloc2(n,&bj_ptr,n,&bjlvl_ptr)); 410 411 /* create a linked list for storing column indices of the active row */ 412 nlnk = n + 1; 413 PetscCall(PetscIncompleteLLCreate(n,n,nlnk,lnk,lnk_lvl,lnkbt)); 414 415 /* initial FreeSpace size is f*(ai[n]+1) */ 416 PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(f,ai[n]+1),&free_space)); 417 current_space = free_space; 418 PetscCall(PetscFreeSpaceGet(PetscRealIntMultTruncate(f,ai[n]+1),&free_space_lvl)); 419 current_space_lvl = free_space_lvl; 420 421 for (i=0; i<n; i++) { 422 nzi = 0; 423 /* copy current row into linked list */ 424 nnz = ai[r[i]+1] - ai[r[i]]; 425 PetscCheck(nnz,PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Empty row in matrix: row in original ordering %" PetscInt_FMT " in permuted ordering %" PetscInt_FMT,r[i],i); 426 cols = aj + ai[r[i]]; 427 lnk[i] = -1; /* marker to indicate if diagonal exists */ 428 PetscCall(PetscIncompleteLLInit(nnz,cols,n,ic,&nlnk,lnk,lnk_lvl,lnkbt)); 429 nzi += nlnk; 430 431 /* make sure diagonal entry is included */ 432 if (diagonal_fill && lnk[i] == -1) { 433 fm = n; 434 while (lnk[fm] < i) fm = lnk[fm]; 435 lnk[i] = lnk[fm]; /* insert diagonal into linked list */ 436 lnk[fm] = i; 437 lnk_lvl[i] = 0; 438 nzi++; dcount++; 439 } 440 441 /* add pivot rows into the active row */ 442 nzbd = 0; 443 prow = lnk[n]; 444 while (prow < i) { 445 nnz = bdiag[prow]; 446 cols = bj_ptr[prow] + nnz + 1; 447 cols_lvl = bjlvl_ptr[prow] + nnz + 1; 448 nnz = bi[prow+1] - bi[prow] - nnz - 1; 449 450 PetscCall(PetscILULLAddSorted(nnz,cols,levels,cols_lvl,prow,&nlnk,lnk,lnk_lvl,lnkbt,prow)); 451 nzi += nlnk; 452 prow = lnk[prow]; 453 nzbd++; 454 } 455 bdiag[i] = nzbd; 456 bi[i+1] = bi[i] + nzi; 457 458 /* if free space is not available, make more free space */ 459 if (current_space->local_remaining<nzi) { 460 nnz = PetscIntMultTruncate(2,PetscIntMultTruncate(nzi,(n - i))); /* estimated and max additional space needed */ 461 PetscCall(PetscFreeSpaceGet(nnz,¤t_space)); 462 PetscCall(PetscFreeSpaceGet(nnz,¤t_space_lvl)); 463 reallocs++; 464 } 465 466 /* copy data into free_space and free_space_lvl, then initialize lnk */ 467 PetscCall(PetscIncompleteLLClean(n,n,nzi,lnk,lnk_lvl,current_space->array,current_space_lvl->array,lnkbt)); 468 469 bj_ptr[i] = current_space->array; 470 bjlvl_ptr[i] = current_space_lvl->array; 471 472 /* make sure the active row i has diagonal entry */ 473 PetscCheck(*(bj_ptr[i]+bdiag[i]) == i,PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Row %" PetscInt_FMT " has missing diagonal in factored matrix\ntry running with -pc_factor_nonzeros_along_diagonal or -pc_factor_diagonal_fill",i); 474 475 current_space->array += nzi; 476 current_space->local_used += nzi; 477 current_space->local_remaining -= nzi; 478 479 current_space_lvl->array += nzi; 480 current_space_lvl->local_used += nzi; 481 current_space_lvl->local_remaining -= nzi; 482 } 483 484 PetscCall(ISRestoreIndices(isrow,&r)); 485 PetscCall(ISRestoreIndices(isicol,&ic)); 486 487 /* copy free_space into bj and free free_space; set bi, bj, bdiag in new datastructure; */ 488 PetscCall(PetscMalloc1(bi[n]+1,&bj)); 489 PetscCall(PetscFreeSpaceContiguous_LU(&free_space,bj,n,bi,bdiag)); 490 491 PetscCall(PetscIncompleteLLDestroy(lnk,lnkbt)); 492 PetscCall(PetscFreeSpaceDestroy(free_space_lvl)); 493 PetscCall(PetscFree2(bj_ptr,bjlvl_ptr)); 494 495 #if defined(PETSC_USE_INFO) 496 { 497 PetscReal af = ((PetscReal)(bdiag[0]+1))/((PetscReal)ai[n]); 498 PetscCall(PetscInfo(A,"Reallocs %" PetscInt_FMT " Fill ratio:given %g needed %g\n",reallocs,(double)f,(double)af)); 499 PetscCall(PetscInfo(A,"Run with -[sub_]pc_factor_fill %g or use \n",(double)af)); 500 PetscCall(PetscInfo(A,"PCFactorSetFill([sub]pc,%g);\n",(double)af)); 501 PetscCall(PetscInfo(A,"for best performance.\n")); 502 if (diagonal_fill) { 503 PetscCall(PetscInfo(A,"Detected and replaced %" PetscInt_FMT " missing diagonals\n",dcount)); 504 } 505 } 506 #endif 507 508 /* put together the new matrix */ 509 PetscCall(MatSeqBAIJSetPreallocation(fact,bs,MAT_SKIP_ALLOCATION,NULL)); 510 PetscCall(PetscLogObjectParent((PetscObject)fact,(PetscObject)isicol)); 511 512 b = (Mat_SeqBAIJ*)(fact)->data; 513 b->free_a = PETSC_TRUE; 514 b->free_ij = PETSC_TRUE; 515 b->singlemalloc = PETSC_FALSE; 516 517 PetscCall(PetscMalloc1(bs2*(bdiag[0]+1),&b->a)); 518 519 b->j = bj; 520 b->i = bi; 521 b->diag = bdiag; 522 b->free_diag = PETSC_TRUE; 523 b->ilen = NULL; 524 b->imax = NULL; 525 b->row = isrow; 526 b->col = iscol; 527 PetscCall(PetscObjectReference((PetscObject)isrow)); 528 PetscCall(PetscObjectReference((PetscObject)iscol)); 529 b->icol = isicol; 530 531 PetscCall(PetscMalloc1(bs*n+bs,&b->solve_work)); 532 /* In b structure: Free imax, ilen, old a, old j. 533 Allocate bdiag, solve_work, new a, new j */ 534 PetscCall(PetscLogObjectMemory((PetscObject)fact,(bdiag[0]+1) * (sizeof(PetscInt)+bs2*sizeof(PetscScalar)))); 535 b->maxnz = b->nz = bdiag[0]+1; 536 537 fact->info.factor_mallocs = reallocs; 538 fact->info.fill_ratio_given = f; 539 fact->info.fill_ratio_needed = ((PetscReal)(bdiag[0]+1))/((PetscReal)ai[n]); 540 541 PetscCall(MatSeqBAIJSetNumericFactorization(fact,both_identity)); 542 PetscFunctionReturn(0); 543 } 544 545 /* 546 This code is virtually identical to MatILUFactorSymbolic_SeqAIJ 547 except that the data structure of Mat_SeqAIJ is slightly different. 548 Not a good example of code reuse. 549 */ 550 PetscErrorCode MatILUFactorSymbolic_SeqBAIJ_inplace(Mat fact,Mat A,IS isrow,IS iscol,const MatFactorInfo *info) 551 { 552 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b; 553 IS isicol; 554 const PetscInt *r,*ic,*ai = a->i,*aj = a->j,*xi; 555 PetscInt prow,n = a->mbs,*ainew,*ajnew,jmax,*fill,nz,*im,*ajfill,*flev,*xitmp; 556 PetscInt *dloc,idx,row,m,fm,nzf,nzi,reallocate = 0,dcount = 0; 557 PetscInt incrlev,nnz,i,bs = A->rmap->bs,bs2 = a->bs2,levels,diagonal_fill,dd; 558 PetscBool col_identity,row_identity,both_identity,flg; 559 PetscReal f; 560 561 PetscFunctionBegin; 562 PetscCall(MatMissingDiagonal_SeqBAIJ(A,&flg,&dd)); 563 PetscCheck(!flg,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix A is missing diagonal entry in row %" PetscInt_FMT,dd); 564 565 f = info->fill; 566 levels = (PetscInt)info->levels; 567 diagonal_fill = (PetscInt)info->diagonal_fill; 568 569 PetscCall(ISInvertPermutation(iscol,PETSC_DECIDE,&isicol)); 570 571 PetscCall(ISIdentity(isrow,&row_identity)); 572 PetscCall(ISIdentity(iscol,&col_identity)); 573 both_identity = (PetscBool) (row_identity && col_identity); 574 575 if (!levels && both_identity) { /* special case copy the nonzero structure */ 576 PetscCall(MatDuplicateNoCreate_SeqBAIJ(fact,A,MAT_DO_NOT_COPY_VALUES,PETSC_TRUE)); 577 PetscCall(MatSeqBAIJSetNumericFactorization_inplace(fact,both_identity)); 578 579 fact->factortype = MAT_FACTOR_ILU; 580 b = (Mat_SeqBAIJ*)fact->data; 581 b->row = isrow; 582 b->col = iscol; 583 PetscCall(PetscObjectReference((PetscObject)isrow)); 584 PetscCall(PetscObjectReference((PetscObject)iscol)); 585 b->icol = isicol; 586 b->pivotinblocks = (info->pivotinblocks) ? PETSC_TRUE : PETSC_FALSE; 587 588 PetscCall(PetscMalloc1((n+1)*bs,&b->solve_work)); 589 PetscFunctionReturn(0); 590 } 591 592 /* general case perform the symbolic factorization */ 593 PetscCall(ISGetIndices(isrow,&r)); 594 PetscCall(ISGetIndices(isicol,&ic)); 595 596 /* get new row pointers */ 597 PetscCall(PetscMalloc1(n+1,&ainew)); 598 ainew[0] = 0; 599 /* don't know how many column pointers are needed so estimate */ 600 jmax = (PetscInt)(f*ai[n] + 1); 601 PetscCall(PetscMalloc1(jmax,&ajnew)); 602 /* ajfill is level of fill for each fill entry */ 603 PetscCall(PetscMalloc1(jmax,&ajfill)); 604 /* fill is a linked list of nonzeros in active row */ 605 PetscCall(PetscMalloc1(n+1,&fill)); 606 /* im is level for each filled value */ 607 PetscCall(PetscMalloc1(n+1,&im)); 608 /* dloc is location of diagonal in factor */ 609 PetscCall(PetscMalloc1(n+1,&dloc)); 610 dloc[0] = 0; 611 for (prow=0; prow<n; prow++) { 612 613 /* copy prow into linked list */ 614 nzf = nz = ai[r[prow]+1] - ai[r[prow]]; 615 PetscCheck(nz,PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Empty row in matrix: row in original ordering %" PetscInt_FMT " in permuted ordering %" PetscInt_FMT,r[prow],prow); 616 xi = aj + ai[r[prow]]; 617 fill[n] = n; 618 fill[prow] = -1; /* marker for diagonal entry */ 619 while (nz--) { 620 fm = n; 621 idx = ic[*xi++]; 622 do { 623 m = fm; 624 fm = fill[m]; 625 } while (fm < idx); 626 fill[m] = idx; 627 fill[idx] = fm; 628 im[idx] = 0; 629 } 630 631 /* make sure diagonal entry is included */ 632 if (diagonal_fill && fill[prow] == -1) { 633 fm = n; 634 while (fill[fm] < prow) fm = fill[fm]; 635 fill[prow] = fill[fm]; /* insert diagonal into linked list */ 636 fill[fm] = prow; 637 im[prow] = 0; 638 nzf++; 639 dcount++; 640 } 641 642 nzi = 0; 643 row = fill[n]; 644 while (row < prow) { 645 incrlev = im[row] + 1; 646 nz = dloc[row]; 647 xi = ajnew + ainew[row] + nz + 1; 648 flev = ajfill + ainew[row] + nz + 1; 649 nnz = ainew[row+1] - ainew[row] - nz - 1; 650 fm = row; 651 while (nnz-- > 0) { 652 idx = *xi++; 653 if (*flev + incrlev > levels) { 654 flev++; 655 continue; 656 } 657 do { 658 m = fm; 659 fm = fill[m]; 660 } while (fm < idx); 661 if (fm != idx) { 662 im[idx] = *flev + incrlev; 663 fill[m] = idx; 664 fill[idx] = fm; 665 fm = idx; 666 nzf++; 667 } else if (im[idx] > *flev + incrlev) im[idx] = *flev+incrlev; 668 flev++; 669 } 670 row = fill[row]; 671 nzi++; 672 } 673 /* copy new filled row into permanent storage */ 674 ainew[prow+1] = ainew[prow] + nzf; 675 if (ainew[prow+1] > jmax) { 676 677 /* estimate how much additional space we will need */ 678 /* use the strategy suggested by David Hysom <hysom@perch-t.icase.edu> */ 679 /* just double the memory each time */ 680 PetscInt maxadd = jmax; 681 /* maxadd = (int)(((f*ai[n]+1)*(n-prow+5))/n); */ 682 if (maxadd < nzf) maxadd = (n-prow)*(nzf+1); 683 jmax += maxadd; 684 685 /* allocate a longer ajnew and ajfill */ 686 PetscCall(PetscMalloc1(jmax,&xitmp)); 687 PetscCall(PetscArraycpy(xitmp,ajnew,ainew[prow])); 688 PetscCall(PetscFree(ajnew)); 689 ajnew = xitmp; 690 PetscCall(PetscMalloc1(jmax,&xitmp)); 691 PetscCall(PetscArraycpy(xitmp,ajfill,ainew[prow])); 692 PetscCall(PetscFree(ajfill)); 693 ajfill = xitmp; 694 reallocate++; /* count how many reallocations are needed */ 695 } 696 xitmp = ajnew + ainew[prow]; 697 flev = ajfill + ainew[prow]; 698 dloc[prow] = nzi; 699 fm = fill[n]; 700 while (nzf--) { 701 *xitmp++ = fm; 702 *flev++ = im[fm]; 703 fm = fill[fm]; 704 } 705 /* make sure row has diagonal entry */ 706 PetscCheckFalse(ajnew[ainew[prow]+dloc[prow]] != prow,PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Row %" PetscInt_FMT " has missing diagonal in factored matrix\n\ 707 try running with -pc_factor_nonzeros_along_diagonal or -pc_factor_diagonal_fill",prow); 708 } 709 PetscCall(PetscFree(ajfill)); 710 PetscCall(ISRestoreIndices(isrow,&r)); 711 PetscCall(ISRestoreIndices(isicol,&ic)); 712 PetscCall(PetscFree(fill)); 713 PetscCall(PetscFree(im)); 714 715 #if defined(PETSC_USE_INFO) 716 { 717 PetscReal af = ((PetscReal)ainew[n])/((PetscReal)ai[n]); 718 PetscCall(PetscInfo(A,"Reallocs %" PetscInt_FMT " Fill ratio:given %g needed %g\n",reallocate,(double)f,(double)af)); 719 PetscCall(PetscInfo(A,"Run with -pc_factor_fill %g or use \n",(double)af)); 720 PetscCall(PetscInfo(A,"PCFactorSetFill(pc,%g);\n",(double)af)); 721 PetscCall(PetscInfo(A,"for best performance.\n")); 722 if (diagonal_fill) { 723 PetscCall(PetscInfo(A,"Detected and replaced %" PetscInt_FMT " missing diagonals\n",dcount)); 724 } 725 } 726 #endif 727 728 /* put together the new matrix */ 729 PetscCall(MatSeqBAIJSetPreallocation(fact,bs,MAT_SKIP_ALLOCATION,NULL)); 730 PetscCall(PetscLogObjectParent((PetscObject)fact,(PetscObject)isicol)); 731 b = (Mat_SeqBAIJ*)fact->data; 732 733 b->free_a = PETSC_TRUE; 734 b->free_ij = PETSC_TRUE; 735 b->singlemalloc = PETSC_FALSE; 736 737 PetscCall(PetscMalloc1(bs2*ainew[n],&b->a)); 738 739 b->j = ajnew; 740 b->i = ainew; 741 for (i=0; i<n; i++) dloc[i] += ainew[i]; 742 b->diag = dloc; 743 b->free_diag = PETSC_TRUE; 744 b->ilen = NULL; 745 b->imax = NULL; 746 b->row = isrow; 747 b->col = iscol; 748 b->pivotinblocks = (info->pivotinblocks) ? PETSC_TRUE : PETSC_FALSE; 749 750 PetscCall(PetscObjectReference((PetscObject)isrow)); 751 PetscCall(PetscObjectReference((PetscObject)iscol)); 752 b->icol = isicol; 753 PetscCall(PetscMalloc1(bs*n+bs,&b->solve_work)); 754 /* In b structure: Free imax, ilen, old a, old j. 755 Allocate dloc, solve_work, new a, new j */ 756 PetscCall(PetscLogObjectMemory((PetscObject)fact,(ainew[n]-n)*(sizeof(PetscInt))+bs2*ainew[n]*sizeof(PetscScalar))); 757 b->maxnz = b->nz = ainew[n]; 758 759 fact->info.factor_mallocs = reallocate; 760 fact->info.fill_ratio_given = f; 761 fact->info.fill_ratio_needed = ((PetscReal)ainew[n])/((PetscReal)ai[prow]); 762 763 PetscCall(MatSeqBAIJSetNumericFactorization_inplace(fact,both_identity)); 764 PetscFunctionReturn(0); 765 } 766 767 PetscErrorCode MatSetUnfactored_SeqBAIJ_4_NaturalOrdering_SSE(Mat A) 768 { 769 /* Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; */ 770 /* int i,*AJ=a->j,nz=a->nz; */ 771 772 PetscFunctionBegin; 773 /* Undo Column scaling */ 774 /* while (nz--) { */ 775 /* AJ[i] = AJ[i]/4; */ 776 /* } */ 777 /* This should really invoke a push/pop logic, but we don't have that yet. */ 778 A->ops->setunfactored = NULL; 779 PetscFunctionReturn(0); 780 } 781 782 PetscErrorCode MatSetUnfactored_SeqBAIJ_4_NaturalOrdering_SSE_usj(Mat A) 783 { 784 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 785 PetscInt *AJ=a->j,nz=a->nz; 786 unsigned short *aj=(unsigned short*)AJ; 787 788 PetscFunctionBegin; 789 /* Is this really necessary? */ 790 while (nz--) { 791 AJ[nz] = (int)((unsigned int)aj[nz]); /* First extend, then convert to signed. */ 792 } 793 A->ops->setunfactored = NULL; 794 PetscFunctionReturn(0); 795 } 796