1 #ifndef lint 2 static char vcid[] = "$Id: aijfact.c,v 1.55 1996/01/24 05:45:53 bsmith Exp bsmith $"; 3 #endif 4 5 #include "aij.h" 6 /* 7 Factorization code for AIJ format. 8 */ 9 10 int MatLUFactorSymbolic_SeqAIJ(Mat A,IS isrow,IS iscol,double f,Mat *B) 11 { 12 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data, *b; 13 IS isicol; 14 int *r,*ic, ierr, i, n = a->m, *ai = a->i, *aj = a->j; 15 int *ainew,*ajnew, jmax,*fill, *ajtmp, nz,shift = a->indexshift; 16 int *idnew, idx, row,m,fm, nnz, nzi,len, realloc = 0,nzbd,*im; 17 18 if (n != a->n) SETERRQ(1,"MatLUFactorSymbolic_SeqAIJ:Must be square"); 19 if (!isrow) SETERRQ(1,"MatLUFactorSymbolic_SeqAIJ:Must have row permutation"); 20 if (!iscol) SETERRQ(1,"MatLUFactorSymbolic_SeqAIJ:Must have col. permutation"); 21 22 ierr = ISInvertPermutation(iscol,&isicol); CHKERRQ(ierr); 23 ISGetIndices(isrow,&r); ISGetIndices(isicol,&ic); 24 25 /* get new row pointers */ 26 ainew = (int *) PetscMalloc( (n+1)*sizeof(int) ); CHKPTRQ(ainew); 27 ainew[0] = -shift; 28 /* don't know how many column pointers are needed so estimate */ 29 jmax = (int) (f*ai[n]+(!shift)); 30 ajnew = (int *) PetscMalloc( (jmax)*sizeof(int) ); CHKPTRQ(ajnew); 31 /* fill is a linked list of nonzeros in active row */ 32 fill = (int *) PetscMalloc( (2*n+1)*sizeof(int)); CHKPTRQ(fill); 33 im = fill + n + 1; 34 /* idnew is location of diagonal in factor */ 35 idnew = (int *) PetscMalloc( (n+1)*sizeof(int)); CHKPTRQ(idnew); 36 idnew[0] = -shift; 37 38 for ( i=0; i<n; i++ ) { 39 /* first copy previous fill into linked list */ 40 nnz = nz = ai[r[i]+1] - ai[r[i]]; 41 ajtmp = aj + ai[r[i]] + shift; 42 fill[n] = n; 43 while (nz--) { 44 fm = n; 45 idx = ic[*ajtmp++ + shift]; 46 do { 47 m = fm; 48 fm = fill[m]; 49 } while (fm < idx); 50 fill[m] = idx; 51 fill[idx] = fm; 52 } 53 row = fill[n]; 54 while ( row < i ) { 55 ajtmp = ajnew + idnew[row] + (!shift); 56 nzbd = 1 + idnew[row] - ainew[row]; 57 nz = im[row] - nzbd; 58 fm = row; 59 while (nz-- > 0) { 60 idx = *ajtmp++ + shift; 61 nzbd++; 62 if (idx == i) im[row] = nzbd; 63 do { 64 m = fm; 65 fm = fill[m]; 66 } while (fm < idx); 67 if (fm != idx) { 68 fill[m] = idx; 69 fill[idx] = fm; 70 fm = idx; 71 nnz++; 72 } 73 } 74 row = fill[row]; 75 } 76 /* copy new filled row into permanent storage */ 77 ainew[i+1] = ainew[i] + nnz; 78 if (ainew[i+1] > jmax+1) { 79 /* allocate a longer ajnew */ 80 int maxadd; 81 maxadd = (int) ((f*(ai[n]+(!shift))*(n-i+5))/n); 82 if (maxadd < nnz) maxadd = (n-i)*(nnz+1); 83 jmax += maxadd; 84 ajtmp = (int *) PetscMalloc( jmax*sizeof(int) );CHKPTRQ(ajtmp); 85 PetscMemcpy(ajtmp,ajnew,(ainew[i]+shift)*sizeof(int)); 86 PetscFree(ajnew); 87 ajnew = ajtmp; 88 realloc++; /* count how many times we realloc */ 89 } 90 ajtmp = ajnew + ainew[i] + shift; 91 fm = fill[n]; 92 nzi = 0; 93 im[i] = nnz; 94 while (nnz--) { 95 if (fm < i) nzi++; 96 *ajtmp++ = fm - shift; 97 fm = fill[fm]; 98 } 99 idnew[i] = ainew[i] + nzi; 100 } 101 102 PLogInfo((PetscObject)A, 103 "Info:MatLUFactorSymbolic_SeqAIJ:Reallocs %d Fill ratio:given %g needed %g\n", 104 realloc,f,((double)ainew[n])/((double)ai[i])); 105 106 ierr = ISRestoreIndices(isrow,&r); CHKERRQ(ierr); 107 ierr = ISRestoreIndices(isicol,&ic); CHKERRQ(ierr); 108 109 PetscFree(fill); 110 111 /* put together the new matrix */ 112 ierr = MatCreateSeqAIJ(A->comm,n,n,0,PETSC_NULL,B); CHKERRQ(ierr); 113 PLogObjectParent(*B,isicol); 114 ierr = ISDestroy(isicol); CHKERRQ(ierr); 115 b = (Mat_SeqAIJ *) (*B)->data; 116 PetscFree(b->imax); 117 b->singlemalloc = 0; 118 len = (ainew[n] + shift)*sizeof(Scalar); 119 /* the next line frees the default space generated by the Create() */ 120 PetscFree(b->a); PetscFree(b->ilen); 121 b->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(b->a); 122 b->j = ajnew; 123 b->i = ainew; 124 b->diag = idnew; 125 b->ilen = 0; 126 b->imax = 0; 127 b->row = isrow; 128 b->col = iscol; 129 b->solve_work = (Scalar *) PetscMalloc( n*sizeof(Scalar)); 130 CHKPTRQ(b->solve_work); 131 /* In b structure: Free imax, ilen, old a, old j. 132 Allocate idnew, solve_work, new a, new j */ 133 PLogObjectMemory(*B,(ainew[n]+shift-n)*(sizeof(int)+sizeof(Scalar))); 134 b->maxnz = b->nz = ainew[n] + shift; 135 136 return 0; 137 } 138 /* ----------------------------------------------------------- */ 139 int Mat_AIJ_CheckInode(Mat); 140 141 int MatLUFactorNumeric_SeqAIJ(Mat A,Mat *B) 142 { 143 Mat C = *B; 144 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data, *b = (Mat_SeqAIJ *)C->data; 145 IS iscol = b->col, isrow = b->row, isicol; 146 int *r,*ic, ierr, i, j, n = a->m, *ai = b->i, *aj = b->j; 147 int *ajtmpold, *ajtmp, nz, row, *ics, shift = a->indexshift; 148 int *diag_offset = b->diag; 149 Scalar *rtmp,*v, *pc, multiplier; 150 /* These declarations are for optimizations. They reduce the number of 151 memory references that are made by locally storing information; the 152 word "register" used here with pointers can be viewed as "private" or 153 "known only to me" 154 */ 155 register Scalar *pv, *rtmps; 156 register int *pj; 157 158 ierr = ISInvertPermutation(iscol,&isicol); CHKERRQ(ierr); 159 PLogObjectParent(*B,isicol); 160 ierr = ISGetIndices(isrow,&r); CHKERRQ(ierr); 161 ierr = ISGetIndices(isicol,&ic); CHKERRQ(ierr); 162 rtmp = (Scalar *) PetscMalloc( (n+1)*sizeof(Scalar) ); CHKPTRQ(rtmp); 163 rtmps = rtmp + shift; ics = ic + shift; 164 165 for ( i=0; i<n; i++ ) { 166 nz = ai[i+1] - ai[i]; 167 ajtmp = aj + ai[i] + shift; 168 for ( j=0; j<nz; j++ ) rtmps[ajtmp[j]] = 0.0; 169 /* for(j = 0; j < n; ++j) rtmp[j] =0.0; */ 170 171 /* load in initial (unfactored row) */ 172 nz = a->i[r[i]+1] - a->i[r[i]]; 173 ajtmpold = a->j + a->i[r[i]] + shift; 174 v = a->a + a->i[r[i]] + shift; 175 for ( j=0; j<nz; j++ ) rtmp[ics[ajtmpold[j]]] = v[j]; 176 177 row = *ajtmp++ + shift; 178 while (row < i) { 179 pc = rtmp + row; 180 if (*pc != 0.0) { 181 pv = b->a + diag_offset[row] + shift; 182 pj = b->j + diag_offset[row] + (!shift); 183 multiplier = *pc * *pv++; 184 *pc = multiplier; 185 nz = ai[row+1] - diag_offset[row] - 1; 186 PLogFlops(2*nz); 187 /* The for-loop form can aid the compiler in overlapping 188 loads and stores */ 189 /*while (nz-->0) rtmps[*pj++] -= multiplier* *pv++; */ 190 { 191 int __i; 192 for (__i=0; __i<nz; __i++) rtmps[pj[__i]] -= multiplier * pv[__i]; 193 } 194 } 195 row = *ajtmp++ + shift; 196 } 197 /* finished row so stick it into b->a */ 198 pv = b->a + ai[i] + shift; 199 pj = b->j + ai[i] + shift; 200 nz = ai[i+1] - ai[i]; 201 if (rtmp[i] == 0.0) {SETERRQ(1,"MatLUFactorNumeric_SeqAIJ:Zero pivot");} 202 rtmp[i] = 1.0/rtmp[i]; 203 for ( j=0; j<nz; j++ ) {pv[j] = rtmps[pj[j]];} 204 } 205 206 PetscFree(rtmp); 207 ierr = ISRestoreIndices(isicol,&ic); CHKERRQ(ierr); 208 ierr = ISRestoreIndices(isrow,&r); CHKERRQ(ierr); 209 ierr = ISDestroy(isicol); CHKERRQ(ierr); 210 C->factor = FACTOR_LU; 211 ierr = Mat_AIJ_CheckInode(C); CHKERRQ(ierr); 212 C->assembled = PETSC_TRUE; 213 PLogFlops(b->n); 214 return 0; 215 } 216 /* ----------------------------------------------------------- */ 217 int MatLUFactor_SeqAIJ(Mat A,IS row,IS col,double f) 218 { 219 Mat_SeqAIJ *mat = (Mat_SeqAIJ *) A->data; 220 int ierr; 221 Mat C; 222 223 ierr = MatLUFactorSymbolic_SeqAIJ(A,row,col,f,&C); CHKERRQ(ierr); 224 ierr = MatLUFactorNumeric_SeqAIJ(A,&C); CHKERRQ(ierr); 225 226 /* free all the data structures from mat */ 227 PetscFree(mat->a); 228 if (!mat->singlemalloc) {PetscFree(mat->i); PetscFree(mat->j);} 229 if (mat->diag) PetscFree(mat->diag); 230 if (mat->ilen) PetscFree(mat->ilen); 231 if (mat->imax) PetscFree(mat->imax); 232 if (mat->solve_work) PetscFree(mat->solve_work); 233 if (mat->inode.size) PetscFree(mat->inode.size); 234 PetscFree(mat); 235 236 PetscMemcpy(A,C,sizeof(struct _Mat)); 237 PetscHeaderDestroy(C); 238 return 0; 239 } 240 /* ----------------------------------------------------------- */ 241 int MatSolve_SeqAIJ(Mat A,Vec bb, Vec xx) 242 { 243 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 244 IS iscol = a->col, isrow = a->row; 245 int *r,*c, ierr, i, n = a->m, *vi, *ai = a->i, *aj = a->j; 246 int nz,shift = a->indexshift; 247 Scalar *x,*b,*tmp, *tmps, *aa = a->a, sum, *v; 248 249 if (A->factor != FACTOR_LU) SETERRQ(1,"MatSolve_SeqAIJ:Not for unfactored matrix"); 250 251 ierr = VecGetArray(bb,&b); CHKERRQ(ierr); 252 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 253 tmp = a->solve_work; 254 255 ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr); 256 ierr = ISGetIndices(iscol,&c);CHKERRQ(ierr); c = c + (n-1); 257 258 /* forward solve the lower triangular */ 259 tmp[0] = b[*r++]; 260 tmps = tmp + shift; 261 for ( i=1; i<n; i++ ) { 262 v = aa + ai[i] + shift; 263 vi = aj + ai[i] + shift; 264 nz = a->diag[i] - ai[i]; 265 sum = b[*r++]; 266 while (nz--) sum -= *v++ * tmps[*vi++]; 267 tmp[i] = sum; 268 } 269 270 /* backward solve the upper triangular */ 271 for ( i=n-1; i>=0; i-- ){ 272 v = aa + a->diag[i] + (!shift); 273 vi = aj + a->diag[i] + (!shift); 274 nz = ai[i+1] - a->diag[i] - 1; 275 sum = tmp[i]; 276 while (nz--) sum -= *v++ * tmps[*vi++]; 277 x[*c--] = tmp[i] = sum*aa[a->diag[i]+shift]; 278 } 279 280 ierr = ISRestoreIndices(isrow,&r); CHKERRQ(ierr); 281 ierr = ISRestoreIndices(iscol,&c); CHKERRQ(ierr); 282 PLogFlops(2*a->nz - a->n); 283 return 0; 284 } 285 int MatSolveAdd_SeqAIJ(Mat A,Vec bb, Vec yy, Vec xx) 286 { 287 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 288 IS iscol = a->col, isrow = a->row; 289 int *r,*c, ierr, i, n = a->m, *vi, *ai = a->i, *aj = a->j; 290 int nz, shift = a->indexshift; 291 Scalar *x,*b,*tmp, *aa = a->a, sum, *v; 292 293 if (A->factor != FACTOR_LU) SETERRQ(1,"MatSolveAdd_SeqAIJ:Not for unfactored matrix"); 294 if (yy != xx) {ierr = VecCopy(yy,xx); CHKERRQ(ierr);} 295 296 ierr = VecGetArray(bb,&b); CHKERRQ(ierr); 297 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 298 tmp = a->solve_work; 299 300 ierr = ISGetIndices(isrow,&r); CHKERRQ(ierr); 301 ierr = ISGetIndices(iscol,&c); CHKERRQ(ierr); c = c + (n-1); 302 303 /* forward solve the lower triangular */ 304 tmp[0] = b[*r++]; 305 for ( i=1; i<n; i++ ) { 306 v = aa + ai[i] + shift; 307 vi = aj + ai[i] + shift; 308 nz = a->diag[i] - ai[i]; 309 sum = b[*r++]; 310 while (nz--) sum -= *v++ * tmp[*vi++ + shift]; 311 tmp[i] = sum; 312 } 313 314 /* backward solve the upper triangular */ 315 for ( i=n-1; i>=0; i-- ){ 316 v = aa + a->diag[i] + (!shift); 317 vi = aj + a->diag[i] + (!shift); 318 nz = ai[i+1] - a->diag[i] - 1; 319 sum = tmp[i]; 320 while (nz--) sum -= *v++ * tmp[*vi++ + shift]; 321 tmp[i] = sum*aa[a->diag[i]+shift]; 322 x[*c--] += tmp[i]; 323 } 324 325 ierr = ISRestoreIndices(isrow,&r); CHKERRQ(ierr); 326 ierr = ISRestoreIndices(iscol,&c); CHKERRQ(ierr); 327 PLogFlops(2*a->nz); 328 329 return 0; 330 } 331 /* -------------------------------------------------------------------*/ 332 int MatSolveTrans_SeqAIJ(Mat A,Vec bb, Vec xx) 333 { 334 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 335 IS iscol = a->col, isrow = a->row, invisrow,inviscol; 336 int *r,*c, ierr, i, n = a->m, *vi, *ai = a->i, *aj = a->j; 337 int nz,shift = a->indexshift; 338 Scalar *x,*b,*tmp, *aa = a->a, *v; 339 340 if (A->factor != FACTOR_LU) SETERRQ(1,"MatSolveTrans_SeqAIJ:Not unfactored matrix"); 341 ierr = VecGetArray(bb,&b); CHKERRQ(ierr); 342 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 343 tmp = a->solve_work; 344 345 /* invert the permutations */ 346 ierr = ISInvertPermutation(isrow,&invisrow); CHKERRQ(ierr); 347 ierr = ISInvertPermutation(iscol,&inviscol); CHKERRQ(ierr); 348 349 ierr = ISGetIndices(invisrow,&r); CHKERRQ(ierr); 350 ierr = ISGetIndices(inviscol,&c); CHKERRQ(ierr); 351 352 /* copy the b into temp work space according to permutation */ 353 for ( i=0; i<n; i++ ) tmp[c[i]] = b[i]; 354 355 /* forward solve the U^T */ 356 for ( i=0; i<n; i++ ) { 357 v = aa + a->diag[i] + shift; 358 vi = aj + a->diag[i] + (!shift); 359 nz = ai[i+1] - a->diag[i] - 1; 360 tmp[i] *= *v++; 361 while (nz--) { 362 tmp[*vi++ + shift] -= (*v++)*tmp[i]; 363 } 364 } 365 366 /* backward solve the L^T */ 367 for ( i=n-1; i>=0; i-- ){ 368 v = aa + a->diag[i] - 1 + shift; 369 vi = aj + a->diag[i] - 1 + shift; 370 nz = a->diag[i] - ai[i]; 371 while (nz--) { 372 tmp[*vi-- + shift] -= (*v--)*tmp[i]; 373 } 374 } 375 376 /* copy tmp into x according to permutation */ 377 for ( i=0; i<n; i++ ) x[r[i]] = tmp[i]; 378 379 ierr = ISRestoreIndices(invisrow,&r); CHKERRQ(ierr); 380 ierr = ISRestoreIndices(inviscol,&c); CHKERRQ(ierr); 381 ierr = ISDestroy(invisrow); CHKERRQ(ierr); 382 ierr = ISDestroy(inviscol); CHKERRQ(ierr); 383 384 PLogFlops(2*a->nz-a->n); 385 return 0; 386 } 387 388 int MatSolveTransAdd_SeqAIJ(Mat A,Vec bb, Vec zz,Vec xx) 389 { 390 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 391 IS iscol = a->col, isrow = a->row, invisrow,inviscol; 392 int *r,*c, ierr, i, n = a->m, *vi, *ai = a->i, *aj = a->j; 393 int nz,shift = a->indexshift; 394 Scalar *x,*b,*tmp, *aa = a->a, *v; 395 396 if (A->factor != FACTOR_LU)SETERRQ(1,"MatSolveTransAdd_SeqAIJ:Not unfactored matrix"); 397 if (zz != xx) VecCopy(zz,xx); 398 399 ierr = VecGetArray(bb,&b); CHKERRQ(ierr); 400 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 401 tmp = a->solve_work; 402 403 /* invert the permutations */ 404 ierr = ISInvertPermutation(isrow,&invisrow); CHKERRQ(ierr); 405 ierr = ISInvertPermutation(iscol,&inviscol); CHKERRQ(ierr); 406 ierr = ISGetIndices(invisrow,&r); CHKERRQ(ierr); 407 ierr = ISGetIndices(inviscol,&c); CHKERRQ(ierr); 408 409 /* copy the b into temp work space according to permutation */ 410 for ( i=0; i<n; i++ ) tmp[c[i]] = b[i]; 411 412 /* forward solve the U^T */ 413 for ( i=0; i<n; i++ ) { 414 v = aa + a->diag[i] + shift; 415 vi = aj + a->diag[i] + (!shift); 416 nz = ai[i+1] - a->diag[i] - 1; 417 tmp[i] *= *v++; 418 while (nz--) { 419 tmp[*vi++ + shift] -= (*v++)*tmp[i]; 420 } 421 } 422 423 /* backward solve the L^T */ 424 for ( i=n-1; i>=0; i-- ){ 425 v = aa + a->diag[i] - 1 + shift; 426 vi = aj + a->diag[i] - 1 + shift; 427 nz = a->diag[i] - ai[i]; 428 while (nz--) { 429 tmp[*vi-- + shift] -= (*v--)*tmp[i]; 430 } 431 } 432 433 /* copy tmp into x according to permutation */ 434 for ( i=0; i<n; i++ ) x[r[i]] += tmp[i]; 435 436 ierr = ISRestoreIndices(invisrow,&r); CHKERRQ(ierr); 437 ierr = ISRestoreIndices(inviscol,&c); CHKERRQ(ierr); 438 ierr = ISDestroy(invisrow); CHKERRQ(ierr); 439 ierr = ISDestroy(inviscol); CHKERRQ(ierr); 440 441 PLogFlops(2*a->nz); 442 return 0; 443 } 444 /* ----------------------------------------------------------------*/ 445 446 int MatILUFactorSymbolic_SeqAIJ(Mat A,IS isrow,IS iscol,double f,int levels,Mat *fact) 447 { 448 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data, *b; 449 IS isicol; 450 int *r,*ic, ierr, prow, n = a->m, *ai = a->i, *aj = a->j; 451 int *ainew,*ajnew, jmax,*fill, *xi, nz, *im,*ajfill,*flev; 452 int *dloc, idx, row,m,fm, nzf, nzi,len, realloc = 0; 453 int incrlev,nnz,i,shift = a->indexshift; 454 455 if (n != a->n) SETERRQ(1,"MatILUFactorSymbolic_SeqAIJ:Matrix must be square"); 456 if (!isrow) SETERRQ(1,"MatILUFactorSymbolic_SeqAIJ:Must have row permutation"); 457 if (!iscol) SETERRQ(1,"MatILUFactorSymbolic_SeqAIJ:Must have column permutation"); 458 459 /* special case that simply copies fill pattern */ 460 if (levels == 0 && ISIsIdentity(isrow) && ISIsIdentity(iscol)) { 461 ierr = MatConvertSameType_SeqAIJ(A,fact,DO_NOT_COPY_VALUES); CHKERRQ(ierr); 462 (*fact)->factor = FACTOR_LU; 463 b = (Mat_SeqAIJ *) (*fact)->data; 464 if (!b->diag) { 465 ierr = MatMarkDiag_SeqAIJ(*fact); CHKERRQ(ierr); 466 } 467 b->row = isrow; 468 b->col = iscol; 469 b->solve_work = (Scalar *) PetscMalloc((b->m+1)*sizeof(Scalar));CHKPTRQ(b->solve_work); 470 return 0; 471 } 472 473 ierr = ISInvertPermutation(iscol,&isicol); CHKERRQ(ierr); 474 ierr = ISGetIndices(isrow,&r); CHKERRQ(ierr); 475 ierr = ISGetIndices(isicol,&ic); CHKERRQ(ierr); 476 477 /* get new row pointers */ 478 ainew = (int *) PetscMalloc( (n+1)*sizeof(int) ); CHKPTRQ(ainew); 479 ainew[0] = -shift; 480 /* don't know how many column pointers are needed so estimate */ 481 jmax = (int) (f*(ai[n]+!shift)); 482 ajnew = (int *) PetscMalloc( (jmax)*sizeof(int) ); CHKPTRQ(ajnew); 483 /* ajfill is level of fill for each fill entry */ 484 ajfill = (int *) PetscMalloc( (jmax)*sizeof(int) ); CHKPTRQ(ajfill); 485 /* fill is a linked list of nonzeros in active row */ 486 fill = (int *) PetscMalloc( (n+1)*sizeof(int)); CHKPTRQ(fill); 487 /* im is level for each filled value */ 488 im = (int *) PetscMalloc( (n+1)*sizeof(int)); CHKPTRQ(im); 489 /* dloc is location of diagonal in factor */ 490 dloc = (int *) PetscMalloc( (n+1)*sizeof(int)); CHKPTRQ(dloc); 491 dloc[0] = 0; 492 for ( prow=0; prow<n; prow++ ) { 493 /* first copy previous fill into linked list */ 494 nzf = nz = ai[r[prow]+1] - ai[r[prow]]; 495 xi = aj + ai[r[prow]] + shift; 496 fill[n] = n; 497 while (nz--) { 498 fm = n; 499 idx = ic[*xi++ + shift]; 500 do { 501 m = fm; 502 fm = fill[m]; 503 } while (fm < idx); 504 fill[m] = idx; 505 fill[idx] = fm; 506 im[idx] = 0; 507 } 508 nzi = 0; 509 row = fill[n]; 510 while ( row < prow ) { 511 incrlev = im[row] + 1; 512 nz = dloc[row]; 513 xi = ajnew + ainew[row] + shift + nz; 514 flev = ajfill + ainew[row] + shift + nz + 1; 515 nnz = ainew[row+1] - ainew[row] - nz - 1; 516 if (*xi++ + shift != row) { 517 SETERRQ(1,"MatILUFactorSymbolic_SeqAIJ:zero pivot"); 518 } 519 fm = row; 520 while (nnz-- > 0) { 521 idx = *xi++ + shift; 522 if (*flev + incrlev > levels) { 523 flev++; 524 continue; 525 } 526 do { 527 m = fm; 528 fm = fill[m]; 529 } while (fm < idx); 530 if (fm != idx) { 531 im[idx] = *flev + incrlev; 532 fill[m] = idx; 533 fill[idx] = fm; 534 fm = idx; 535 nzf++; 536 } 537 else { 538 if (im[idx] > *flev + incrlev) im[idx] = *flev+incrlev; 539 } 540 flev++; 541 } 542 row = fill[row]; 543 nzi++; 544 } 545 /* copy new filled row into permanent storage */ 546 ainew[prow+1] = ainew[prow] + nzf; 547 if (ainew[prow+1] > jmax-shift) { 548 /* allocate a longer ajnew */ 549 int maxadd; 550 maxadd = (int) ((f*(ai[n]+!shift)*(n-prow+5))/n); 551 if (maxadd < nzf) maxadd = (n-prow)*(nzf+1); 552 jmax += maxadd; 553 xi = (int *) PetscMalloc( jmax*sizeof(int) );CHKPTRQ(xi); 554 PetscMemcpy(xi,ajnew,(ainew[prow]+shift)*sizeof(int)); 555 PetscFree(ajnew); 556 ajnew = xi; 557 /* allocate a longer ajfill */ 558 xi = (int *) PetscMalloc( jmax*sizeof(int) );CHKPTRQ(xi); 559 PetscMemcpy(xi,ajfill,(ainew[prow]+shift)*sizeof(int)); 560 PetscFree(ajfill); 561 ajfill = xi; 562 realloc++; 563 } 564 xi = ajnew + ainew[prow] + shift; 565 flev = ajfill + ainew[prow] + shift; 566 dloc[prow] = nzi; 567 fm = fill[n]; 568 while (nzf--) { 569 *xi++ = fm - shift; 570 *flev++ = im[fm]; 571 fm = fill[fm]; 572 } 573 } 574 PetscFree(ajfill); 575 ierr = ISRestoreIndices(isrow,&r); CHKERRQ(ierr); 576 ierr = ISRestoreIndices(isicol,&ic); CHKERRQ(ierr); 577 ierr = ISDestroy(isicol); CHKERRQ(ierr); 578 PetscFree(fill); PetscFree(im); 579 580 PLogInfo((PetscObject)A, 581 "Info:MatILUFactorSymbolic_SeqAIJ:Realloc %d Fill ratio:given %g needed %g\n", 582 realloc,f,((double)ainew[n])/((double)ai[prow])); 583 584 /* put together the new matrix */ 585 ierr = MatCreateSeqAIJ(A->comm,n,n,0,PETSC_NULL,fact); CHKERRQ(ierr); 586 b = (Mat_SeqAIJ *) (*fact)->data; 587 PetscFree(b->imax); 588 b->singlemalloc = 0; 589 len = (ainew[n] + shift)*sizeof(Scalar); 590 /* the next line frees the default space generated by the Create() */ 591 PetscFree(b->a); PetscFree(b->ilen); 592 b->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(b->a); 593 b->j = ajnew; 594 b->i = ainew; 595 for ( i=0; i<n; i++ ) dloc[i] += ainew[i]; 596 b->diag = dloc; 597 b->ilen = 0; 598 b->imax = 0; 599 b->row = isrow; 600 b->col = iscol; 601 b->solve_work = (Scalar *) PetscMalloc( (n+1)*sizeof(Scalar)); 602 CHKPTRQ(b->solve_work); 603 /* In b structure: Free imax, ilen, old a, old j. 604 Allocate dloc, solve_work, new a, new j */ 605 PLogObjectMemory(*fact,(ainew[n]+shift-n) * (sizeof(int)+sizeof(Scalar))); 606 b->maxnz = b->nz = ainew[n] + shift; 607 (*fact)->factor = FACTOR_LU; 608 return 0; 609 } 610 611 612 613 614