source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/util/sutil.c @ 4272

Last change on this file since 4272 was 2696, checked in by kanani, 7 years ago

Merge of branch palm4u into trunk

File size: 2.9 KB
Line 
1
2int KppDecomp( KPP_REAL *JVS )
3{
4KPP_REAL W[KPP_NVAR];
5KPP_REAL a;
6int 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  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
28int 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  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
64void 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  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
Note: See TracBrowser for help on using the repository browser.