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