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