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