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