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