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 /* 9 This is used to set the numeric factorization for both LU and ILU symbolic factorization 10 */ 11 PetscErrorCode MatSeqBAIJSetNumericFactorization(Mat fact,PetscBool natural) 12 { 13 PetscFunctionBegin; 14 if (natural) { 15 switch (fact->rmap->bs) { 16 case 1: 17 fact->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_1; 18 break; 19 case 2: 20 fact->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_2_NaturalOrdering; 21 break; 22 case 3: 23 fact->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_3_NaturalOrdering; 24 break; 25 case 4: 26 fact->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_4_NaturalOrdering; 27 break; 28 case 5: 29 fact->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_5_NaturalOrdering; 30 break; 31 case 6: 32 fact->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_6_NaturalOrdering; 33 break; 34 case 7: 35 fact->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_7_NaturalOrdering; 36 break; 37 case 9: 38 #if defined(PETSC_HAVE_IMMINTRIN_H) && defined(__AVX2__) && defined(__FMA__) && defined(PETSC_USE_REAL_DOUBLE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_64BIT_INDICES) 39 fact->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_9_NaturalOrdering; 40 #else 41 fact->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_N; 42 #endif 43 break; 44 case 15: 45 fact->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_15_NaturalOrdering; 46 break; 47 default: 48 fact->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_N; 49 break; 50 } 51 } else { 52 switch (fact->rmap->bs) { 53 case 1: 54 fact->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_1; 55 break; 56 case 2: 57 fact->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_2; 58 break; 59 case 3: 60 fact->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_3; 61 break; 62 case 4: 63 fact->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_4; 64 break; 65 case 5: 66 fact->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_5; 67 break; 68 case 6: 69 fact->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_6; 70 break; 71 case 7: 72 fact->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_7; 73 break; 74 default: 75 fact->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_N; 76 break; 77 } 78 } 79 PetscFunctionReturn(0); 80 } 81 82 PetscErrorCode MatSeqBAIJSetNumericFactorization_inplace(Mat inA,PetscBool natural) 83 { 84 PetscFunctionBegin; 85 if (natural) { 86 switch (inA->rmap->bs) { 87 case 1: 88 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_1_inplace; 89 break; 90 case 2: 91 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_2_NaturalOrdering_inplace; 92 break; 93 case 3: 94 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_3_NaturalOrdering_inplace; 95 break; 96 case 4: 97 #if defined(PETSC_USE_REAL_MAT_SINGLE) 98 { 99 PetscBool sse_enabled_local; 100 PetscErrorCode ierr; 101 ierr = PetscSSEIsEnabled(inA->comm,&sse_enabled_local,NULL);CHKERRQ(ierr); 102 if (sse_enabled_local) { 103 # if defined(PETSC_HAVE_SSE) 104 int i,*AJ=a->j,nz=a->nz,n=a->mbs; 105 if (n==(unsigned short)n) { 106 unsigned short *aj=(unsigned short*)AJ; 107 for (i=0; i<nz; i++) aj[i] = (unsigned short)AJ[i]; 108 109 inA->ops->setunfactored = MatSetUnfactored_SeqBAIJ_4_NaturalOrdering_SSE_usj; 110 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_4_NaturalOrdering_SSE_usj; 111 112 ierr = PetscInfo(inA,"Using special SSE, in-place natural ordering, ushort j index factor BS=4\n");CHKERRQ(ierr); 113 } else { 114 /* Scale the column indices for easier indexing in MatSolve. */ 115 /* for (i=0;i<nz;i++) { */ 116 /* AJ[i] = AJ[i]*4; */ 117 /* } */ 118 inA->ops->setunfactored = MatSetUnfactored_SeqBAIJ_4_NaturalOrdering_SSE; 119 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_4_NaturalOrdering_SSE; 120 121 ierr = PetscInfo(inA,"Using special SSE, in-place natural ordering, int j index factor BS=4\n");CHKERRQ(ierr); 122 } 123 # else 124 /* This should never be reached. If so, problem in PetscSSEIsEnabled. */ 125 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SSE Hardware unavailable"); 126 # endif 127 } else { 128 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_4_NaturalOrdering_inplace; 129 } 130 } 131 #else 132 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_4_NaturalOrdering_inplace; 133 #endif 134 break; 135 case 5: 136 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_5_NaturalOrdering_inplace; 137 break; 138 case 6: 139 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_6_NaturalOrdering_inplace; 140 break; 141 case 7: 142 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_7_NaturalOrdering_inplace; 143 break; 144 default: 145 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_N_inplace; 146 break; 147 } 148 } else { 149 switch (inA->rmap->bs) { 150 case 1: 151 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_1_inplace; 152 break; 153 case 2: 154 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_2_inplace; 155 break; 156 case 3: 157 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_3_inplace; 158 break; 159 case 4: 160 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_4_inplace; 161 break; 162 case 5: 163 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_5_inplace; 164 break; 165 case 6: 166 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_6_inplace; 167 break; 168 case 7: 169 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_7_inplace; 170 break; 171 default: 172 inA->ops->lufactornumeric = MatLUFactorNumeric_SeqBAIJ_N_inplace; 173 break; 174 } 175 } 176 PetscFunctionReturn(0); 177 } 178 179 /* 180 The symbolic factorization code is identical to that for AIJ format, 181 except for very small changes since this is now a SeqBAIJ datastructure. 182 NOT good code reuse. 183 */ 184 #include <petscbt.h> 185 #include <../src/mat/utils/freespace.h> 186 187 PetscErrorCode MatLUFactorSymbolic_SeqBAIJ(Mat B,Mat A,IS isrow,IS iscol,const MatFactorInfo *info) 188 { 189 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b; 190 PetscInt n =a->mbs,bs = A->rmap->bs,bs2=a->bs2; 191 PetscBool row_identity,col_identity,both_identity; 192 IS isicol; 193 PetscErrorCode ierr; 194 const PetscInt *r,*ic; 195 PetscInt i,*ai=a->i,*aj=a->j; 196 PetscInt *bi,*bj,*ajtmp; 197 PetscInt *bdiag,row,nnz,nzi,reallocs=0,nzbd,*im; 198 PetscReal f; 199 PetscInt nlnk,*lnk,k,**bi_ptr; 200 PetscFreeSpaceList free_space=NULL,current_space=NULL; 201 PetscBT lnkbt; 202 PetscBool missing; 203 204 PetscFunctionBegin; 205 if (A->rmap->N != A->cmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"matrix must be square"); 206 ierr = MatMissingDiagonal(A,&missing,&i);CHKERRQ(ierr); 207 if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",i); 208 209 if (bs>1) { /* check shifttype */ 210 if (info->shifttype == MAT_SHIFT_NONZERO || info->shifttype == MAT_SHIFT_POSITIVE_DEFINITE) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only MAT_SHIFT_NONE and MAT_SHIFT_INBLOCKS are supported for BAIJ matrix"); 211 } 212 213 ierr = ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);CHKERRQ(ierr); 214 ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr); 215 ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr); 216 217 /* get new row and diagonal pointers, must be allocated separately because they will be given to the Mat_SeqAIJ and freed separately */ 218 ierr = PetscMalloc1(n+1,&bi);CHKERRQ(ierr); 219 ierr = PetscMalloc1(n+1,&bdiag);CHKERRQ(ierr); 220 bi[0] = bdiag[0] = 0; 221 222 /* linked list for storing column indices of the active row */ 223 nlnk = n + 1; 224 ierr = PetscLLCreate(n,n,nlnk,lnk,lnkbt);CHKERRQ(ierr); 225 226 ierr = PetscMalloc2(n+1,&bi_ptr,n+1,&im);CHKERRQ(ierr); 227 228 /* initial FreeSpace size is f*(ai[n]+1) */ 229 f = info->fill; 230 ierr = PetscFreeSpaceGet(PetscRealIntMultTruncate(f,ai[n]+1),&free_space);CHKERRQ(ierr); 231 232 current_space = free_space; 233 234 for (i=0; i<n; i++) { 235 /* copy previous fill into linked list */ 236 nzi = 0; 237 nnz = ai[r[i]+1] - ai[r[i]]; 238 ajtmp = aj + ai[r[i]]; 239 ierr = PetscLLAddPerm(nnz,ajtmp,ic,n,nlnk,lnk,lnkbt);CHKERRQ(ierr); 240 nzi += nlnk; 241 242 /* add pivot rows into linked list */ 243 row = lnk[n]; 244 while (row < i) { 245 nzbd = bdiag[row] + 1; /* num of entries in the row with column index <= row */ 246 ajtmp = bi_ptr[row] + nzbd; /* points to the entry next to the diagonal */ 247 ierr = PetscLLAddSortedLU(ajtmp,row,nlnk,lnk,lnkbt,i,nzbd,im);CHKERRQ(ierr); 248 nzi += nlnk; 249 row = lnk[row]; 250 } 251 bi[i+1] = bi[i] + nzi; 252 im[i] = nzi; 253 254 /* mark bdiag */ 255 nzbd = 0; 256 nnz = nzi; 257 k = lnk[n]; 258 while (nnz-- && k < i) { 259 nzbd++; 260 k = lnk[k]; 261 } 262 bdiag[i] = nzbd; /* note : bdaig[i] = nnzL as input for PetscFreeSpaceContiguous_LU() */ 263 264 /* if free space is not available, make more free space */ 265 if (current_space->local_remaining<nzi) { 266 nnz = PetscIntMultTruncate(2,PetscIntMultTruncate(n - i,nzi)); /* estimated and max additional space needed */ 267 ierr = PetscFreeSpaceGet(nnz,¤t_space);CHKERRQ(ierr); 268 reallocs++; 269 } 270 271 /* copy data into free space, then initialize lnk */ 272 ierr = PetscLLClean(n,n,nzi,lnk,current_space->array,lnkbt);CHKERRQ(ierr); 273 274 bi_ptr[i] = current_space->array; 275 current_space->array += nzi; 276 current_space->local_used += nzi; 277 current_space->local_remaining -= nzi; 278 } 279 280 ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr); 281 ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr); 282 283 /* copy free_space into bj and free free_space; set bi, bj, bdiag in new datastructure; */ 284 ierr = PetscMalloc1(bi[n]+1,&bj);CHKERRQ(ierr); 285 ierr = PetscFreeSpaceContiguous_LU(&free_space,bj,n,bi,bdiag);CHKERRQ(ierr); 286 ierr = PetscLLDestroy(lnk,lnkbt);CHKERRQ(ierr); 287 ierr = PetscFree2(bi_ptr,im);CHKERRQ(ierr); 288 289 /* put together the new matrix */ 290 ierr = MatSeqBAIJSetPreallocation(B,bs,MAT_SKIP_ALLOCATION,NULL);CHKERRQ(ierr); 291 ierr = PetscLogObjectParent((PetscObject)B,(PetscObject)isicol);CHKERRQ(ierr); 292 b = (Mat_SeqBAIJ*)(B)->data; 293 294 b->free_a = PETSC_TRUE; 295 b->free_ij = PETSC_TRUE; 296 b->singlemalloc = PETSC_FALSE; 297 298 ierr = PetscMalloc1((bdiag[0]+1)*bs2,&b->a);CHKERRQ(ierr); 299 b->j = bj; 300 b->i = bi; 301 b->diag = bdiag; 302 b->free_diag = PETSC_TRUE; 303 b->ilen = NULL; 304 b->imax = NULL; 305 b->row = isrow; 306 b->col = iscol; 307 b->pivotinblocks = (info->pivotinblocks) ? PETSC_TRUE : PETSC_FALSE; 308 309 ierr = PetscObjectReference((PetscObject)isrow);CHKERRQ(ierr); 310 ierr = PetscObjectReference((PetscObject)iscol);CHKERRQ(ierr); 311 b->icol = isicol; 312 ierr = PetscMalloc1(bs*n+bs,&b->solve_work);CHKERRQ(ierr); 313 ierr = PetscLogObjectMemory((PetscObject)B,(bdiag[0]+1)*(sizeof(PetscInt)+sizeof(PetscScalar)*bs2));CHKERRQ(ierr); 314 315 b->maxnz = b->nz = bdiag[0]+1; 316 317 B->factortype = MAT_FACTOR_LU; 318 B->info.factor_mallocs = reallocs; 319 B->info.fill_ratio_given = f; 320 321 if (ai[n] != 0) { 322 B->info.fill_ratio_needed = ((PetscReal)(bdiag[0]+1))/((PetscReal)ai[n]); 323 } else { 324 B->info.fill_ratio_needed = 0.0; 325 } 326 #if defined(PETSC_USE_INFO) 327 if (ai[n] != 0) { 328 PetscReal af = B->info.fill_ratio_needed; 329 ierr = PetscInfo3(A,"Reallocs %D Fill ratio:given %g needed %g\n",reallocs,(double)f,(double)af);CHKERRQ(ierr); 330 ierr = PetscInfo1(A,"Run with -pc_factor_fill %g or use \n",(double)af);CHKERRQ(ierr); 331 ierr = PetscInfo1(A,"PCFactorSetFill(pc,%g);\n",(double)af);CHKERRQ(ierr); 332 ierr = PetscInfo(A,"for best performance.\n");CHKERRQ(ierr); 333 } else { 334 ierr = PetscInfo(A,"Empty matrix\n");CHKERRQ(ierr); 335 } 336 #endif 337 338 ierr = ISIdentity(isrow,&row_identity);CHKERRQ(ierr); 339 ierr = ISIdentity(iscol,&col_identity);CHKERRQ(ierr); 340 341 both_identity = (PetscBool) (row_identity && col_identity); 342 343 ierr = MatSeqBAIJSetNumericFactorization(B,both_identity);CHKERRQ(ierr); 344 PetscFunctionReturn(0); 345 } 346 347 PetscErrorCode MatLUFactorSymbolic_SeqBAIJ_inplace(Mat B,Mat A,IS isrow,IS iscol,const MatFactorInfo *info) 348 { 349 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data,*b; 350 PetscInt n =a->mbs,bs = A->rmap->bs,bs2=a->bs2; 351 PetscBool row_identity,col_identity,both_identity; 352 IS isicol; 353 PetscErrorCode ierr; 354 const PetscInt *r,*ic; 355 PetscInt i,*ai=a->i,*aj=a->j; 356 PetscInt *bi,*bj,*ajtmp; 357 PetscInt *bdiag,row,nnz,nzi,reallocs=0,nzbd,*im; 358 PetscReal f; 359 PetscInt nlnk,*lnk,k,**bi_ptr; 360 PetscFreeSpaceList free_space=NULL,current_space=NULL; 361 PetscBT lnkbt; 362 PetscBool missing; 363 364 PetscFunctionBegin; 365 if (A->rmap->N != A->cmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"matrix must be square"); 366 ierr = MatMissingDiagonal(A,&missing,&i);CHKERRQ(ierr); 367 if (missing) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry %D",i); 368 369 ierr = ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);CHKERRQ(ierr); 370 ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr); 371 ierr = ISGetIndices(isicol,&ic);CHKERRQ(ierr); 372 373 /* get new row and diagonal pointers, must be allocated separately because they will be given to the Mat_SeqAIJ and freed separately */ 374 ierr = PetscMalloc1(n+1,&bi);CHKERRQ(ierr); 375 ierr = PetscMalloc1(n+1,&bdiag);CHKERRQ(ierr); 376 377 bi[0] = bdiag[0] = 0; 378 379 /* linked list for storing column indices of the active row */ 380 nlnk = n + 1; 381 ierr = PetscLLCreate(n,n,nlnk,lnk,lnkbt);CHKERRQ(ierr); 382 383 ierr = PetscMalloc2(n+1,&bi_ptr,n+1,&im);CHKERRQ(ierr); 384 385 /* initial FreeSpace size is f*(ai[n]+1) */ 386 f = info->fill; 387 ierr = PetscFreeSpaceGet(PetscRealIntMultTruncate(f,ai[n]+1),&free_space);CHKERRQ(ierr); 388 current_space = free_space; 389 390 for (i=0; i<n; i++) { 391 /* copy previous fill into linked list */ 392 nzi = 0; 393 nnz = ai[r[i]+1] - ai[r[i]]; 394 ajtmp = aj + ai[r[i]]; 395 ierr = PetscLLAddPerm(nnz,ajtmp,ic,n,nlnk,lnk,lnkbt);CHKERRQ(ierr); 396 nzi += nlnk; 397 398 /* add pivot rows into linked list */ 399 row = lnk[n]; 400 while (row < i) { 401 nzbd = bdiag[row] - bi[row] + 1; /* num of entries in the row with column index <= row */ 402 ajtmp = bi_ptr[row] + nzbd; /* points to the entry next to the diagonal */ 403 ierr = PetscLLAddSortedLU(ajtmp,row,nlnk,lnk,lnkbt,i,nzbd,im);CHKERRQ(ierr); 404 nzi += nlnk; 405 row = lnk[row]; 406 } 407 bi[i+1] = bi[i] + nzi; 408 im[i] = nzi; 409 410 /* mark bdiag */ 411 nzbd = 0; 412 nnz = nzi; 413 k = lnk[n]; 414 while (nnz-- && k < i) { 415 nzbd++; 416 k = lnk[k]; 417 } 418 bdiag[i] = bi[i] + nzbd; 419 420 /* if free space is not available, make more free space */ 421 if (current_space->local_remaining<nzi) { 422 nnz = PetscIntMultTruncate(n - i,nzi); /* estimated and max additional space needed */ 423 ierr = PetscFreeSpaceGet(nnz,¤t_space);CHKERRQ(ierr); 424 reallocs++; 425 } 426 427 /* copy data into free space, then initialize lnk */ 428 ierr = PetscLLClean(n,n,nzi,lnk,current_space->array,lnkbt);CHKERRQ(ierr); 429 430 bi_ptr[i] = current_space->array; 431 current_space->array += nzi; 432 current_space->local_used += nzi; 433 current_space->local_remaining -= nzi; 434 } 435 #if defined(PETSC_USE_INFO) 436 if (ai[n] != 0) { 437 PetscReal af = ((PetscReal)bi[n])/((PetscReal)ai[n]); 438 ierr = PetscInfo3(A,"Reallocs %D Fill ratio:given %g needed %g\n",reallocs,(double)f,(double)af);CHKERRQ(ierr); 439 ierr = PetscInfo1(A,"Run with -pc_factor_fill %g or use \n",(double)af);CHKERRQ(ierr); 440 ierr = PetscInfo1(A,"PCFactorSetFill(pc,%g);\n",(double)af);CHKERRQ(ierr); 441 ierr = PetscInfo(A,"for best performance.\n");CHKERRQ(ierr); 442 } else { 443 ierr = PetscInfo(A,"Empty matrix\n");CHKERRQ(ierr); 444 } 445 #endif 446 447 ierr = ISRestoreIndices(isrow,&r);CHKERRQ(ierr); 448 ierr = ISRestoreIndices(isicol,&ic);CHKERRQ(ierr); 449 450 /* destroy list of free space and other temporary array(s) */ 451 ierr = PetscMalloc1(bi[n]+1,&bj);CHKERRQ(ierr); 452 ierr = PetscFreeSpaceContiguous(&free_space,bj);CHKERRQ(ierr); 453 ierr = PetscLLDestroy(lnk,lnkbt);CHKERRQ(ierr); 454 ierr = PetscFree2(bi_ptr,im);CHKERRQ(ierr); 455 456 /* put together the new matrix */ 457 ierr = MatSeqBAIJSetPreallocation(B,bs,MAT_SKIP_ALLOCATION,NULL);CHKERRQ(ierr); 458 ierr = PetscLogObjectParent((PetscObject)B,(PetscObject)isicol);CHKERRQ(ierr); 459 b = (Mat_SeqBAIJ*)(B)->data; 460 461 b->free_a = PETSC_TRUE; 462 b->free_ij = PETSC_TRUE; 463 b->singlemalloc = PETSC_FALSE; 464 465 ierr = PetscMalloc1((bi[n]+1)*bs2,&b->a);CHKERRQ(ierr); 466 b->j = bj; 467 b->i = bi; 468 b->diag = bdiag; 469 b->free_diag = PETSC_TRUE; 470 b->ilen = NULL; 471 b->imax = NULL; 472 b->row = isrow; 473 b->col = iscol; 474 b->pivotinblocks = (info->pivotinblocks) ? PETSC_TRUE : PETSC_FALSE; 475 476 ierr = PetscObjectReference((PetscObject)isrow);CHKERRQ(ierr); 477 ierr = PetscObjectReference((PetscObject)iscol);CHKERRQ(ierr); 478 b->icol = isicol; 479 480 ierr = PetscMalloc1(bs*n+bs,&b->solve_work);CHKERRQ(ierr); 481 ierr = PetscLogObjectMemory((PetscObject)B,(bi[n]-n)*(sizeof(PetscInt)+sizeof(PetscScalar)*bs2));CHKERRQ(ierr); 482 483 b->maxnz = b->nz = bi[n]; 484 485 (B)->factortype = MAT_FACTOR_LU; 486 (B)->info.factor_mallocs = reallocs; 487 (B)->info.fill_ratio_given = f; 488 489 if (ai[n] != 0) { 490 (B)->info.fill_ratio_needed = ((PetscReal)bi[n])/((PetscReal)ai[n]); 491 } else { 492 (B)->info.fill_ratio_needed = 0.0; 493 } 494 495 ierr = ISIdentity(isrow,&row_identity);CHKERRQ(ierr); 496 ierr = ISIdentity(iscol,&col_identity);CHKERRQ(ierr); 497 498 both_identity = (PetscBool) (row_identity && col_identity); 499 500 ierr = MatSeqBAIJSetNumericFactorization_inplace(B,both_identity);CHKERRQ(ierr); 501 PetscFunctionReturn(0); 502 } 503 504