| 1 | |
|---|
| 2 | int KppDecomp( KPP_REAL *JVS ) |
|---|
| 3 | { |
|---|
| 4 | KPP_REAL W[KPP_NVAR]; |
|---|
| 5 | KPP_REAL a; |
|---|
| 6 | int k, kk, j, jj; |
|---|
| 7 | |
|---|
| 8 | for( k = 0; k < KPP_NVAR; k++ ) { |
|---|
| 9 | if( JVS[ LU_DIAG[k] ] == 0.0 ) return k+1; |
|---|
| 10 | for( kk = LU_CROW[k]; kk < LU_CROW[k+1]; kk++ ) |
|---|
| 11 | W[ LU_ICOL[kk] ] = JVS[kk]; |
|---|
| 12 | for( kk = LU_CROW[k]; kk < LU_DIAG[k]; kk++ ) { |
|---|
| 13 | j = LU_ICOL[kk]; |
|---|
| 14 | a = -W[j] / JVS[ LU_DIAG[j] ]; |
|---|
| 15 | W[j] = -a; |
|---|
| 16 | for( jj = LU_DIAG[j]+1; jj < LU_CROW[j+1]; jj++ ) |
|---|
| 17 | W[ LU_ICOL[jj] ] += a*JVS[jj]; |
|---|
| 18 | } |
|---|
| 19 | for( kk = LU_CROW[k]; kk < LU_CROW[k+1]; kk++ ) |
|---|
| 20 | JVS[kk] = W[ LU_ICOL[kk] ]; |
|---|
| 21 | } |
|---|
| 22 | return 0; |
|---|
| 23 | } |
|---|
| 24 | |
|---|
| 25 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 26 | Sparse LU factorization, complex |
|---|
| 27 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
|---|
| 28 | int KppDecompCmplxR( KPP_REAL JVSR[], KPP_REAL JVSI[] ) |
|---|
| 29 | { |
|---|
| 30 | KPP_REAL WR[NVAR], WI[NVAR]; |
|---|
| 31 | KPP_REAL ar, ai, den; |
|---|
| 32 | int k, kk, j, jj; |
|---|
| 33 | |
|---|
| 34 | for( k = 0; k < NVAR; k++ ) { |
|---|
| 35 | if( JVSR[ LU_DIAG[k] ] == 0.0 ) return k+1; |
|---|
| 36 | if( JVSI[ LU_DIAG[k] ] == 0.0 ) return k+1; |
|---|
| 37 | for( kk = LU_CROW[k]; kk < LU_CROW[k+1]; kk++ ) { |
|---|
| 38 | WR[ LU_ICOL[kk] ] = JVSR[kk]; |
|---|
| 39 | WI[ LU_ICOL[kk] ] = JVSI[kk]; |
|---|
| 40 | } |
|---|
| 41 | for( kk = LU_CROW[k]; kk < LU_DIAG[k]; kk++ ) { |
|---|
| 42 | j = LU_ICOL[kk]; |
|---|
| 43 | den = JVSR[LU_DIAG[j]]*JVSR[LU_DIAG[j]] + JVSI[LU_DIAG[j]]*JVSI[LU_DIAG[j]]; |
|---|
| 44 | ar = -(WR[j]*JVSR[LU_DIAG[j]] + WI[j]*JVSI[LU_DIAG[j]])/den; |
|---|
| 45 | ai = -(WI[j]*JVSR[LU_DIAG[j]] - WR[j]*JVSI[LU_DIAG[j]])/den; |
|---|
| 46 | WR[j] = -ar; |
|---|
| 47 | WI[j] = -ai; |
|---|
| 48 | for( jj = LU_DIAG[j]+1; jj < LU_CROW[j+1]; jj++ ) { |
|---|
| 49 | WR[ LU_ICOL[jj] ] = WR[ LU_ICOL[jj] ] + ar*JVSR[jj] - ai*JVSI[jj]; |
|---|
| 50 | WI[ LU_ICOL[jj] ] = WI[ LU_ICOL[jj] ] + ar*JVSI[jj] + ai*JVSR[jj]; |
|---|
| 51 | } |
|---|
| 52 | } |
|---|
| 53 | for( kk = LU_CROW[k]; kk < LU_CROW[k+1]; kk++ ) { |
|---|
| 54 | JVSR[kk] = WR[ LU_ICOL[kk] ]; |
|---|
| 55 | JVSI[kk] = WI[ LU_ICOL[kk] ]; |
|---|
| 56 | } |
|---|
| 57 | } |
|---|
| 58 | return 0; |
|---|
| 59 | } |
|---|
| 60 | |
|---|
| 61 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 62 | Complex sparse solve subroutine using indirect addressing |
|---|
| 63 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
|---|
| 64 | void KppSolveCmplxR(KPP_REAL JVSR[], KPP_REAL JVSI[], KPP_REAL XR[], KPP_REAL XI[]) |
|---|
| 65 | { |
|---|
| 66 | int i, j; |
|---|
| 67 | KPP_REAL sumr, sumi, den; |
|---|
| 68 | |
|---|
| 69 | for ( i = 0; i < NVAR; i++ ) { |
|---|
| 70 | for ( j = LU_CROW[i]; j < LU_DIAG[i]; j++ ) { |
|---|
| 71 | XR[i] = XR[i] - (JVSR[j]*XR[LU_ICOL[j]] - JVSI[j]*XI[LU_ICOL[j]]); |
|---|
| 72 | XI[i] = XI[i] - (JVSR[j]*XI[LU_ICOL[j]] + JVSI[j]*XR[LU_ICOL[j]]); |
|---|
| 73 | } |
|---|
| 74 | } |
|---|
| 75 | |
|---|
| 76 | for ( i = NVAR-1; i >= 0; i-- ) { |
|---|
| 77 | sumr = XR[i]; |
|---|
| 78 | sumi = XI[i]; |
|---|
| 79 | for ( j = LU_DIAG[i]+1; j < LU_CROW[i+1]; j++) { |
|---|
| 80 | sumr = sumr - (JVSR[j]*XR[LU_ICOL[j]] - JVSI[j]*XI[LU_ICOL[j]]); |
|---|
| 81 | sumi = sumi - (JVSR[j]*XI[LU_ICOL[j]] + JVSI[j]*XR[LU_ICOL[j]]); |
|---|
| 82 | } |
|---|
| 83 | den = JVSR[LU_DIAG[i]]*JVSR[LU_DIAG[i]] + JVSI[LU_DIAG[i]]*JVSI[LU_DIAG[i]]; |
|---|
| 84 | XR[i] = (sumr*JVSR[LU_DIAG[i]] + sumi*JVSI[LU_DIAG[i]])/den; |
|---|
| 85 | XI[i] = (sumi*JVSR[LU_DIAG[i]] - sumr*JVSI[LU_DIAG[i]])/den; |
|---|
| 86 | } |
|---|
| 87 | } |
|---|
| 88 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 89 | END FUNCTION KppSolveCmplxR |
|---|
| 90 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
|---|
