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