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