1 #ifndef lint 2 static char vcid[] = "$Id: aijfact.c,v 1.42 1995/10/12 15:01:41 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 441 int MatILUFactorSymbolic_SeqAIJ(Mat A,IS isrow,IS iscol,double f,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 /* special case that simply copies fill pattern */ 455 if (levels == 0 && ISIsIdentity(isrow) && ISIsIdentity(iscol)) { 456 ierr = MatCopyPrivate_SeqAIJ(A,fact,DO_NOT_COPY_VALUES); CHKERRQ(ierr); 457 (*fact)->factor = FACTOR_LU; 458 b = (Mat_SeqAIJ *) (*fact)->data; 459 if (!b->diag) { 460 ierr = MatMarkDiag_SeqAIJ(*fact); CHKERRQ(ierr); 461 } 462 b->row = isrow; 463 b->col = iscol; 464 b->solve_work = (Scalar *) PETSCMALLOC((b->m+1)*sizeof(Scalar));CHKPTRQ(b->solve_work); 465 return 0; 466 } 467 468 ierr = ISInvertPermutation(iscol,&isicol); CHKERRQ(ierr); 469 ierr = ISGetIndices(isrow,&r); CHKERRQ(ierr); 470 ierr = ISGetIndices(isicol,&ic); CHKERRQ(ierr); 471 472 /* get new row pointers */ 473 ainew = (int *) PETSCMALLOC( (n+1)*sizeof(int) ); CHKPTRQ(ainew); 474 ainew[0] = -shift; 475 /* don't know how many column pointers are needed so estimate */ 476 jmax = (int) (f*(ai[n]+!shift)); 477 ajnew = (int *) PETSCMALLOC( (jmax)*sizeof(int) ); CHKPTRQ(ajnew); 478 /* ajfill is level of fill for each fill entry */ 479 ajfill = (int *) PETSCMALLOC( (jmax)*sizeof(int) ); CHKPTRQ(ajfill); 480 /* fill is a linked list of nonzeros in active row */ 481 fill = (int *) PETSCMALLOC( (n+1)*sizeof(int)); CHKPTRQ(fill); 482 /* im is level for each filled value */ 483 im = (int *) PETSCMALLOC( (n+1)*sizeof(int)); CHKPTRQ(im); 484 /* dloc is location of diagonal in factor */ 485 dloc = (int *) PETSCMALLOC( (n+1)*sizeof(int)); CHKPTRQ(dloc); 486 dloc[0] = 0; 487 488 for ( prow=0; prow<n; prow++ ) { 489 /* first copy previous fill into linked list */ 490 nzf = nz = ai[r[prow]+1] - ai[r[prow]]; 491 xi = aj + ai[r[prow]] + shift; 492 fill[n] = n; 493 while (nz--) { 494 fm = n; 495 idx = ic[*xi++ + shift]; 496 do { 497 m = fm; 498 fm = fill[m]; 499 } while (fm < idx); 500 fill[m] = idx; 501 fill[idx] = fm; 502 im[idx] = 0; 503 } 504 nzi = 0; 505 row = fill[n]; 506 while ( row < prow ) { 507 incrlev = im[row] + 1; 508 nz = dloc[row]; 509 xi = ajnew + ainew[row] + shift + nz; 510 flev = ajfill + ainew[row] + shift + nz + 1; 511 nnz = ainew[row+1] - ainew[row] - nz - 1; 512 if (*xi++ + shift != row) { 513 SETERRQ(1,"MatILUFactorSymbolic_SeqAIJ:zero pivot"); 514 } 515 fm = row; 516 while (nnz-- > 0) { 517 idx = *xi++ + shift; 518 if (*flev + incrlev > levels) { 519 flev++; 520 continue; 521 } 522 do { 523 m = fm; 524 fm = fill[m]; 525 } while (fm < idx); 526 if (fm != idx) { 527 im[idx] = *flev + incrlev; 528 fill[m] = idx; 529 fill[idx] = fm; 530 fm = idx; 531 nzf++; 532 } 533 else { 534 if (im[idx] > *flev + incrlev) im[idx] = *flev+incrlev; 535 } 536 flev++; 537 } 538 row = fill[row]; 539 nzi++; 540 } 541 /* copy new filled row into permanent storage */ 542 ainew[prow+1] = ainew[prow] + nzf; 543 if (ainew[prow+1] > jmax+1) { 544 /* allocate a longer ajnew */ 545 int maxadd; 546 maxadd = (int) ((f*(ai[n]+!shift)*(n-prow+5))/n); 547 if (maxadd < nzf) maxadd = (n-prow)*(nzf+1); 548 jmax += maxadd; 549 xi = (int *) PETSCMALLOC( jmax*sizeof(int) );CHKPTRQ(xi); 550 PetscMemcpy(xi,ajnew,(ainew[prow]+shift)*sizeof(int)); 551 PETSCFREE(ajnew); 552 ajnew = xi; 553 /* allocate a longer ajfill */ 554 xi = (int *) PETSCMALLOC( jmax*sizeof(int) );CHKPTRQ(xi); 555 PetscMemcpy(xi,ajfill,(ainew[prow]+shift)*sizeof(int)); 556 PETSCFREE(ajfill); 557 ajfill = xi; 558 realloc++; 559 } 560 xi = ajnew + ainew[prow] + shift; 561 flev = ajfill + ainew[prow] + shift; 562 dloc[prow] = nzi; 563 fm = fill[n]; 564 while (nzf--) { 565 *xi++ = fm - shift; 566 *flev++ = im[fm]; 567 fm = fill[fm]; 568 } 569 } 570 PETSCFREE(ajfill); 571 ierr = ISRestoreIndices(isrow,&r); CHKERRQ(ierr); 572 ierr = ISRestoreIndices(isicol,&ic); CHKERRQ(ierr); 573 ierr = ISDestroy(isicol); CHKERRQ(ierr); 574 PETSCFREE(fill); PETSCFREE(im); 575 576 PLogInfo((PetscObject)A, 577 "Info:MatILUFactorSymbolic_SeqAIJ:Realloc %d Fill ratio:given %g needed %g\n", 578 realloc,f,((double)ainew[n])/((double)ai[prow])); 579 580 /* put together the new matrix */ 581 ierr = MatCreateSeqAIJ(A->comm,n, n, 0, 0, fact); CHKERRQ(ierr); 582 b = (Mat_SeqAIJ *) (*fact)->data; 583 PETSCFREE(b->imax); 584 b->singlemalloc = 0; 585 len = (ainew[n] + shift)*sizeof(Scalar); 586 /* the next line frees the default space generated by the Create() */ 587 PETSCFREE(b->a); PETSCFREE(b->ilen); 588 b->a = (Scalar *) PETSCMALLOC( len ); CHKPTRQ(b->a); 589 b->j = ajnew; 590 b->i = ainew; 591 for ( i=0; i<n; i++ ) dloc[i] += ainew[i]; 592 b->diag = dloc; 593 b->ilen = 0; 594 b->imax = 0; 595 b->row = isrow; 596 b->col = iscol; 597 b->solve_work = (Scalar *) PETSCMALLOC( (n+1)*sizeof(Scalar)); 598 CHKPTRQ(b->solve_work); 599 /* In b structure: Free imax, ilen, old a, old j. 600 Allocate dloc, solve_work, new a, new j */ 601 PLogObjectMemory(*fact,(ainew[n]+shift-n) * (sizeof(int)+sizeof(Scalar))); 602 b->maxnz = b->nz = ainew[n] + shift; 603 (*fact)->factor = FACTOR_LU; 604 return 0; 605 } 606