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