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