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