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