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