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