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