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