1 #ifndef lint 2 static char vcid[] = "$Id: aijfact.c,v 1.39 1995/09/22 20:17:46 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 ierr = ISDestroy(isicol); CHKERRQ(ierr); 111 PETSCFREE(fill); 112 113 /* put together the new matrix */ 114 ierr = MatCreateSeqAIJ(A->comm,n, n, 0, 0, B); 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 /* Cannot do this here because child is destroyed before parent created 137 PLogObjectParent(*B,isicol); */ 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,*pj, *ics, shift = a->indexshift; 148 Scalar *rtmp,*v, *pv, *pc, multiplier, *rtmps; 149 150 ierr = ISInvertPermutation(iscol,&isicol); CHKERRQ(ierr); 151 PLogObjectParent(*B,isicol); 152 ierr = ISGetIndices(isrow,&r); CHKERRQ(ierr); 153 ierr = ISGetIndices(isicol,&ic); CHKERRQ(ierr); 154 rtmp = (Scalar *) PETSCMALLOC( (n+1)*sizeof(Scalar) ); CHKPTRQ(rtmp); 155 rtmps = rtmp + shift; ics = ic + shift; 156 157 for ( i=0; i<n; i++ ) { 158 nz = ai[i+1] - ai[i]; 159 ajtmp = aj + ai[i] + shift; 160 for ( j=0; j<nz; j++ ) rtmps[ajtmp[j]] = 0.0; 161 162 /* load in initial (unfactored row) */ 163 nz = a->i[r[i]+1] - a->i[r[i]]; 164 ajtmpold = a->j + a->i[r[i]] + shift; 165 v = a->a + a->i[r[i]] + shift; 166 for ( j=0; j<nz; j++ ) rtmp[ics[ajtmpold[j]]] = v[j]; 167 168 row = *ajtmp++ + shift; 169 while (row < i) { 170 pc = rtmp + row; 171 if (*pc != 0.0) { 172 pv = b->a + b->diag[row] + shift; 173 pj = b->j + b->diag[row] + (!shift); 174 multiplier = *pc * *pv++; 175 *pc = multiplier; 176 nz = ai[row+1] - b->diag[row] - 1; 177 PLogFlops(2*nz); 178 while (nz-->0) rtmps[*pj++] -= multiplier* *pv++; 179 } 180 row = *ajtmp++ + shift; 181 } 182 /* finished row so stick it into b->a */ 183 pv = b->a + ai[i] + shift; 184 pj = b->j + ai[i] + shift; 185 nz = ai[i+1] - ai[i]; 186 if (rtmp[i] == 0.0) {SETERRQ(1,"MatLUFactorNumeric_SeqAIJ:Zero pivot");} 187 rtmp[i] = 1.0/rtmp[i]; 188 for ( j=0; j<nz; j++ ) {pv[j] = rtmps[pj[j]];} 189 } 190 PETSCFREE(rtmp); 191 ierr = ISRestoreIndices(isicol,&ic); CHKERRQ(ierr); 192 ierr = ISRestoreIndices(isrow,&r); CHKERRQ(ierr); 193 ierr = ISDestroy(isicol); CHKERRQ(ierr); 194 C->factor = FACTOR_LU; 195 b->assembled = 1; 196 PLogFlops(b->n); 197 return 0; 198 } 199 /* ----------------------------------------------------------- */ 200 int MatLUFactor_SeqAIJ(Mat A,IS row,IS col,double f) 201 { 202 Mat_SeqAIJ *mat = (Mat_SeqAIJ *) A->data; 203 int ierr; 204 Mat C; 205 206 ierr = MatLUFactorSymbolic_SeqAIJ(A,row,col,f,&C); CHKERRQ(ierr); 207 ierr = MatLUFactorNumeric_SeqAIJ(A,&C); CHKERRQ(ierr); 208 209 /* free all the data structures from mat */ 210 PETSCFREE(mat->a); 211 if (!mat->singlemalloc) {PETSCFREE(mat->i); PETSCFREE(mat->j);} 212 if (mat->diag) PETSCFREE(mat->diag); 213 if (mat->ilen) PETSCFREE(mat->ilen); 214 if (mat->imax) PETSCFREE(mat->imax); 215 if (mat->solve_work) PETSCFREE(mat->solve_work); 216 PETSCFREE(mat); 217 218 PetscMemcpy(A,C,sizeof(struct _Mat)); 219 PETSCHEADERDESTROY(C); 220 return 0; 221 } 222 /* ----------------------------------------------------------- */ 223 int MatSolve_SeqAIJ(Mat A,Vec bb, Vec xx) 224 { 225 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 226 IS iscol = a->col, isrow = a->row; 227 int *r,*c, ierr, i, n = a->m, *vi, *ai = a->i, *aj = a->j; 228 int nz,shift = a->indexshift; 229 Scalar *x,*b,*tmp, *tmps, *aa = a->a, sum, *v; 230 231 if (A->factor != FACTOR_LU) SETERRQ(1,"MatSolve_SeqAIJ:Not for unfactored matrix"); 232 233 ierr = VecGetArray(bb,&b); CHKERRQ(ierr); 234 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 235 tmp = a->solve_work; 236 237 ierr = ISGetIndices(isrow,&r);CHKERRQ(ierr); 238 ierr = ISGetIndices(iscol,&c);CHKERRQ(ierr); c = c + (n-1); 239 240 /* forward solve the lower triangular */ 241 tmp[0] = b[*r++]; 242 tmps = tmp + shift; 243 for ( i=1; i<n; i++ ) { 244 v = aa + ai[i] + shift; 245 vi = aj + ai[i] + shift; 246 nz = a->diag[i] - ai[i]; 247 sum = b[*r++]; 248 while (nz--) sum -= *v++ * tmps[*vi++]; 249 tmp[i] = sum; 250 } 251 252 /* backward solve the upper triangular */ 253 for ( i=n-1; i>=0; i-- ){ 254 v = aa + a->diag[i] + (!shift); 255 vi = aj + a->diag[i] + (!shift); 256 nz = ai[i+1] - a->diag[i] - 1; 257 sum = tmp[i]; 258 while (nz--) sum -= *v++ * tmps[*vi++]; 259 x[*c--] = tmp[i] = sum*aa[a->diag[i]+shift]; 260 } 261 262 ierr = ISRestoreIndices(isrow,&r); CHKERRQ(ierr); 263 ierr = ISRestoreIndices(iscol,&c); CHKERRQ(ierr); 264 PLogFlops(2*a->nz - a->n); 265 return 0; 266 } 267 int MatSolveAdd_SeqAIJ(Mat A,Vec bb, Vec yy, Vec xx) 268 { 269 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 270 IS iscol = a->col, isrow = a->row; 271 int *r,*c, ierr, i, n = a->m, *vi, *ai = a->i, *aj = a->j; 272 int nz, shift = a->indexshift; 273 Scalar *x,*b,*tmp, *aa = a->a, sum, *v; 274 275 if (A->factor != FACTOR_LU) SETERRQ(1,"MatSolveAdd_SeqAIJ:Not for unfactored matrix"); 276 if (yy != xx) {ierr = VecCopy(yy,xx); CHKERRQ(ierr);} 277 278 ierr = VecGetArray(bb,&b); CHKERRQ(ierr); 279 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 280 tmp = a->solve_work; 281 282 ierr = ISGetIndices(isrow,&r); CHKERRQ(ierr); 283 ierr = ISGetIndices(iscol,&c); CHKERRQ(ierr); c = c + (n-1); 284 285 /* forward solve the lower triangular */ 286 tmp[0] = b[*r++]; 287 for ( i=1; i<n; i++ ) { 288 v = aa + ai[i] + shift; 289 vi = aj + ai[i] + shift; 290 nz = a->diag[i] - ai[i]; 291 sum = b[*r++]; 292 while (nz--) sum -= *v++ * tmp[*vi++ + shift]; 293 tmp[i] = sum; 294 } 295 296 /* backward solve the upper triangular */ 297 for ( i=n-1; i>=0; i-- ){ 298 v = aa + a->diag[i] + (!shift); 299 vi = aj + a->diag[i] + (!shift); 300 nz = ai[i+1] - a->diag[i] - 1; 301 sum = tmp[i]; 302 while (nz--) sum -= *v++ * tmp[*vi++ + shift]; 303 tmp[i] = sum*aa[a->diag[i]+shift]; 304 x[*c--] += tmp[i]; 305 } 306 307 ierr = ISRestoreIndices(isrow,&r); CHKERRQ(ierr); 308 ierr = ISRestoreIndices(iscol,&c); CHKERRQ(ierr); 309 PLogFlops(2*a->nz); 310 311 return 0; 312 } 313 /* -------------------------------------------------------------------*/ 314 int MatSolveTrans_SeqAIJ(Mat A,Vec bb, Vec xx) 315 { 316 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 317 IS iscol = a->col, isrow = a->row, invisrow,inviscol; 318 int *r,*c, ierr, i, n = a->m, *vi, *ai = a->i, *aj = a->j; 319 int nz,shift = a->indexshift; 320 Scalar *x,*b,*tmp, *aa = a->a, *v; 321 322 if (A->factor != FACTOR_LU) SETERRQ(1,"MatSolveTrans_SeqAIJ:Not unfactored matrix"); 323 ierr = VecGetArray(bb,&b); CHKERRQ(ierr); 324 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 325 tmp = a->solve_work; 326 327 /* invert the permutations */ 328 ierr = ISInvertPermutation(isrow,&invisrow); CHKERRQ(ierr); 329 ierr = ISInvertPermutation(iscol,&inviscol); CHKERRQ(ierr); 330 331 ierr = ISGetIndices(invisrow,&r); CHKERRQ(ierr); 332 ierr = ISGetIndices(inviscol,&c); CHKERRQ(ierr); 333 334 /* copy the b into temp work space according to permutation */ 335 for ( i=0; i<n; i++ ) tmp[c[i]] = b[i]; 336 337 /* forward solve the U^T */ 338 for ( i=0; i<n; i++ ) { 339 v = aa + a->diag[i] + shift; 340 vi = aj + a->diag[i] + (!shift); 341 nz = ai[i+1] - a->diag[i] - 1; 342 tmp[i] *= *v++; 343 while (nz--) { 344 tmp[*vi++ + shift] -= (*v++)*tmp[i]; 345 } 346 } 347 348 /* backward solve the L^T */ 349 for ( i=n-1; i>=0; i-- ){ 350 v = aa + a->diag[i] - 1 + shift; 351 vi = aj + a->diag[i] - 1 + shift; 352 nz = a->diag[i] - ai[i]; 353 while (nz--) { 354 tmp[*vi-- + shift] -= (*v--)*tmp[i]; 355 } 356 } 357 358 /* copy tmp into x according to permutation */ 359 for ( i=0; i<n; i++ ) x[r[i]] = tmp[i]; 360 361 ierr = ISRestoreIndices(invisrow,&r); CHKERRQ(ierr); 362 ierr = ISRestoreIndices(inviscol,&c); CHKERRQ(ierr); 363 ierr = ISDestroy(invisrow); CHKERRQ(ierr); 364 ierr = ISDestroy(inviscol); CHKERRQ(ierr); 365 366 PLogFlops(2*a->nz-a->n); 367 return 0; 368 } 369 370 int MatSolveTransAdd_SeqAIJ(Mat A,Vec bb, Vec zz,Vec xx) 371 { 372 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 373 IS iscol = a->col, isrow = a->row, invisrow,inviscol; 374 int *r,*c, ierr, i, n = a->m, *vi, *ai = a->i, *aj = a->j; 375 int nz,shift = a->indexshift; 376 Scalar *x,*b,*tmp, *aa = a->a, *v; 377 378 if (A->factor != FACTOR_LU)SETERRQ(1,"MatSolveTransAdd_SeqAIJ:Not unfactored matrix"); 379 if (zz != xx) VecCopy(zz,xx); 380 381 ierr = VecGetArray(bb,&b); CHKERRQ(ierr); 382 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 383 tmp = a->solve_work; 384 385 /* invert the permutations */ 386 ierr = ISInvertPermutation(isrow,&invisrow); CHKERRQ(ierr); 387 ierr = ISInvertPermutation(iscol,&inviscol); CHKERRQ(ierr); 388 ierr = ISGetIndices(invisrow,&r); CHKERRQ(ierr); 389 ierr = ISGetIndices(inviscol,&c); CHKERRQ(ierr); 390 391 /* copy the b into temp work space according to permutation */ 392 for ( i=0; i<n; i++ ) tmp[c[i]] = b[i]; 393 394 /* forward solve the U^T */ 395 for ( i=0; i<n; i++ ) { 396 v = aa + a->diag[i] + shift; 397 vi = aj + a->diag[i] + (!shift); 398 nz = ai[i+1] - a->diag[i] - 1; 399 tmp[i] *= *v++; 400 while (nz--) { 401 tmp[*vi++ + shift] -= (*v++)*tmp[i]; 402 } 403 } 404 405 /* backward solve the L^T */ 406 for ( i=n-1; i>=0; i-- ){ 407 v = aa + a->diag[i] - 1 + shift; 408 vi = aj + a->diag[i] - 1 + shift; 409 nz = a->diag[i] - ai[i]; 410 while (nz--) { 411 tmp[*vi-- + shift] -= (*v--)*tmp[i]; 412 } 413 } 414 415 /* copy tmp into x according to permutation */ 416 for ( i=0; i<n; i++ ) x[r[i]] += tmp[i]; 417 418 ierr = ISRestoreIndices(invisrow,&r); CHKERRQ(ierr); 419 ierr = ISRestoreIndices(inviscol,&c); CHKERRQ(ierr); 420 ierr = ISDestroy(invisrow); CHKERRQ(ierr); 421 ierr = ISDestroy(inviscol); CHKERRQ(ierr); 422 423 PLogFlops(2*a->nz); 424 return 0; 425 } 426 /* ----------------------------------------------------------------*/ 427 int MatILUFactorSymbolic_SeqAIJ(Mat A,IS isrow,IS iscol,double f, 428 int levels,Mat *fact) 429 { 430 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data, *b; 431 IS isicol; 432 int *r,*ic, ierr, prow, n = a->m, *ai = a->i, *aj = a->j; 433 int *ainew,*ajnew, jmax,*fill, *xi, nz, *im,*ajfill,*flev; 434 int *dloc, idx, row,m,fm, nzf, nzi,len, realloc = 0; 435 int incrlev,nnz,i,shift = a->indexshift; 436 437 if (n != a->n) SETERRQ(1,"MatILUFactorSymbolic_SeqAIJ:Matrix must be square"); 438 if (!isrow) SETERRQ(1,"MatILUFactorSymbolic_SeqAIJ:Must have row permutation"); 439 if (!iscol) SETERRQ(1,"MatILUFactorSymbolic_SeqAIJ:Must have column permutation"); 440 441 ierr = ISInvertPermutation(iscol,&isicol); CHKERRQ(ierr); 442 ierr = ISGetIndices(isrow,&r); CHKERRQ(ierr); 443 ierr = ISGetIndices(isicol,&ic); CHKERRQ(ierr); 444 445 /* get new row pointers */ 446 ainew = (int *) PETSCMALLOC( (n+1)*sizeof(int) ); CHKPTRQ(ainew); 447 ainew[0] = -shift; 448 /* don't know how many column pointers are needed so estimate */ 449 jmax = (int) (f*(ai[n]+!shift)); 450 ajnew = (int *) PETSCMALLOC( (jmax)*sizeof(int) ); CHKPTRQ(ajnew); 451 /* ajfill is level of fill for each fill entry */ 452 ajfill = (int *) PETSCMALLOC( (jmax)*sizeof(int) ); CHKPTRQ(ajfill); 453 /* fill is a linked list of nonzeros in active row */ 454 fill = (int *) PETSCMALLOC( (n+1)*sizeof(int)); CHKPTRQ(fill); 455 /* im is level for each filled value */ 456 im = (int *) PETSCMALLOC( (n+1)*sizeof(int)); CHKPTRQ(im); 457 /* dloc is location of diagonal in factor */ 458 dloc = (int *) PETSCMALLOC( (n+1)*sizeof(int)); CHKPTRQ(dloc); 459 dloc[0] = 0; 460 461 for ( prow=0; prow<n; prow++ ) { 462 /* first copy previous fill into linked list */ 463 nzf = nz = ai[r[prow]+1] - ai[r[prow]]; 464 xi = aj + ai[r[prow]] + shift; 465 fill[n] = n; 466 while (nz--) { 467 fm = n; 468 idx = ic[*xi++ + shift]; 469 do { 470 m = fm; 471 fm = fill[m]; 472 } while (fm < idx); 473 fill[m] = idx; 474 fill[idx] = fm; 475 im[idx] = 0; 476 } 477 nzi = 0; 478 row = fill[n]; 479 while ( row < prow ) { 480 incrlev = im[row] + 1; 481 nz = dloc[row]; 482 xi = ajnew + ainew[row] + shift + nz; 483 flev = ajfill + ainew[row] + shift + nz + 1; 484 nnz = ainew[row+1] - ainew[row] - nz - 1; 485 if (*xi++ + shift != row) { 486 SETERRQ(1,"MatILUFactorSymbolic_SeqAIJ:zero pivot"); 487 } 488 fm = row; 489 while (nnz-- > 0) { 490 idx = *xi++ + shift; 491 if (*flev + incrlev > levels) { 492 flev++; 493 continue; 494 } 495 do { 496 m = fm; 497 fm = fill[m]; 498 } while (fm < idx); 499 if (fm != idx) { 500 im[idx] = *flev + incrlev; 501 fill[m] = idx; 502 fill[idx] = fm; 503 fm = idx; 504 nzf++; 505 } 506 else { 507 if (im[idx] > *flev + incrlev) im[idx] = *flev+incrlev; 508 } 509 flev++; 510 } 511 row = fill[row]; 512 nzi++; 513 } 514 /* copy new filled row into permanent storage */ 515 ainew[prow+1] = ainew[prow] + nzf; 516 if (ainew[prow+1] > jmax+1) { 517 /* allocate a longer ajnew */ 518 int maxadd; 519 maxadd = (int) ((f*(ai[n]+!shift)*(n-prow+5))/n); 520 if (maxadd < nzf) maxadd = (n-prow)*(nzf+1); 521 jmax += maxadd; 522 xi = (int *) PETSCMALLOC( jmax*sizeof(int) );CHKPTRQ(xi); 523 PetscMemcpy(xi,ajnew,(ainew[prow]+shift)*sizeof(int)); 524 PETSCFREE(ajnew); 525 ajnew = xi; 526 /* allocate a longer ajfill */ 527 xi = (int *) PETSCMALLOC( jmax*sizeof(int) );CHKPTRQ(xi); 528 PetscMemcpy(xi,ajfill,(ainew[prow]+shift)*sizeof(int)); 529 PETSCFREE(ajfill); 530 ajfill = xi; 531 realloc++; 532 } 533 xi = ajnew + ainew[prow] + shift; 534 flev = ajfill + ainew[prow] + shift; 535 dloc[prow] = nzi; 536 fm = fill[n]; 537 while (nzf--) { 538 *xi++ = fm - shift; 539 *flev++ = im[fm]; 540 fm = fill[fm]; 541 } 542 } 543 PETSCFREE(ajfill); 544 ierr = ISRestoreIndices(isrow,&r); CHKERRQ(ierr); 545 ierr = ISRestoreIndices(isicol,&ic); CHKERRQ(ierr); 546 ierr = ISDestroy(isicol); CHKERRQ(ierr); 547 PETSCFREE(fill); PETSCFREE(im); 548 549 PLogInfo((PetscObject)A, 550 "Info:MatILUFactorSymbolic_SeqAIJ:Realloc %d Fill ratio:given %g needed %g\n", 551 realloc,f,((double)ainew[n])/((double)ai[prow])); 552 553 /* put together the new matrix */ 554 ierr = MatCreateSeqAIJ(A->comm,n, n, 0, 0, fact); CHKERRQ(ierr); 555 b = (Mat_SeqAIJ *) (*fact)->data; 556 PETSCFREE(b->imax); 557 b->singlemalloc = 0; 558 len = (ainew[n] + shift)*sizeof(Scalar); 559 /* the next line frees the default space generated by the Create() */ 560 PETSCFREE(b->a); PETSCFREE(b->ilen); 561 b->a = (Scalar *) PETSCMALLOC( len ); CHKPTRQ(b->a); 562 b->j = ajnew; 563 b->i = ainew; 564 for ( i=0; i<n; i++ ) dloc[i] += ainew[i]; 565 b->diag = dloc; 566 b->ilen = 0; 567 b->imax = 0; 568 b->row = isrow; 569 b->col = iscol; 570 b->solve_work = (Scalar *) PETSCMALLOC( (n+1)*sizeof(Scalar)); 571 CHKPTRQ(b->solve_work); 572 /* In b structure: Free imax, ilen, old a, old j. 573 Allocate dloc, solve_work, new a, new j */ 574 PLogObjectMemory(*fact,(ainew[n]+shift-n) * (sizeof(int)+sizeof(Scalar))); 575 b->maxnz = b->nz = ainew[n] + shift; 576 (*fact)->factor = FACTOR_LU; 577 return 0; 578 } 579