source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/util/sutil.f90 @ 4620

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

Merge of branch palm4u into trunk

File size: 18.1 KB
Line 
1
2! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3SUBROUTINE KppDecomp( JVS, IER )
4! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
5!        Sparse LU factorization
6! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
7
8  USE KPP_ROOT_Parameters
9  USE KPP_ROOT_JacobianSP
10
11      INTEGER  :: IER
12      KPP_REAL :: JVS(LU_NONZERO), W(NVAR), a
13      INTEGER  :: k, kk, j, jj
14
15      a = 0. ! mz_rs_20050606
16      IER = 0
17      DO k=1,NVAR
18        ! mz_rs_20050606: don't check if real value == 0
19        ! IF ( JVS( LU_DIAG(k) ) .EQ. 0. ) THEN
20        IF ( ABS(JVS(LU_DIAG(k))) < TINY(a) ) THEN
21            IER = k
22            RETURN
23        END IF
24        DO kk = LU_CROW(k), LU_CROW(k+1)-1
25              W( LU_ICOL(kk) ) = JVS(kk)
26        END DO
27        DO kk = LU_CROW(k), LU_DIAG(k)-1
28            j = LU_ICOL(kk)
29            a = -W(j) / JVS( LU_DIAG(j) )
30            W(j) = -a
31            DO jj = LU_DIAG(j)+1, LU_CROW(j+1)-1
32               W( LU_ICOL(jj) ) = W( LU_ICOL(jj) ) + a*JVS(jj)
33            END DO
34         END DO
35         DO kk = LU_CROW(k), LU_CROW(k+1)-1
36            JVS(kk) = W( LU_ICOL(kk) )
37         END DO
38      END DO
39     
40END SUBROUTINE KppDecomp
41
42
43! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
44SUBROUTINE KppDecompCmplx( JVS, IER )
45! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
46!        Sparse LU factorization, complex
47! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
48
49  USE KPP_ROOT_Parameters
50  USE KPP_ROOT_JacobianSP
51
52      INTEGER        :: IER
53      DOUBLE COMPLEX :: JVS(LU_NONZERO), W(NVAR), a
54      KPP_REAL  :: b = 0.0
55      INTEGER        :: k, kk, j, jj
56
57      IER = 0
58      DO k=1,NVAR
59        IF ( ABS(JVS(LU_DIAG(k))) < TINY(b) ) THEN
60            IER = k
61            RETURN
62        END IF
63        DO kk = LU_CROW(k), LU_CROW(k+1)-1
64              W( LU_ICOL(kk) ) = JVS(kk)
65        END DO
66        DO kk = LU_CROW(k), LU_DIAG(k)-1
67            j = LU_ICOL(kk)
68            a = -W(j) / JVS( LU_DIAG(j) )
69            W(j) = -a
70            DO jj = LU_DIAG(j)+1, LU_CROW(j+1)-1
71               W( LU_ICOL(jj) ) = W( LU_ICOL(jj) ) + a*JVS(jj)
72            END DO
73         END DO
74         DO kk = LU_CROW(k), LU_CROW(k+1)-1
75            JVS(kk) = W( LU_ICOL(kk) )
76         END DO
77      END DO
78     
79END SUBROUTINE KppDecompCmplx
80
81
82! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
83SUBROUTINE KppDecompCmplxR( JVSR, JVSI, IER )
84! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
85!    Sparse LU factorization, complex
86!   (Real and Imaginary parts are used instead of complex data type)     
87! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
88
89  USE KPP_ROOT_Parameters
90  USE KPP_ROOT_JacobianSP
91
92      INTEGER       :: IER
93      KPP_REAL :: JVSR(LU_NONZERO), JVSI(LU_NONZERO) 
94      KPP_REAL :: WR(NVAR), WI(NVAR), ar, ai, den
95      INTEGER       :: k, kk, j, jj
96
97      IER = 0
98      ar  = 0.0
99      DO k=1,NVAR
100        IF (  ( ABS(JVSR(LU_DIAG(k))) < TINY(ar) ) .AND. &
101              ( ABS(JVSI(LU_DIAG(k))) < TINY(ar) ) )  THEN
102            IER = k
103            RETURN
104        END IF
105        DO kk = LU_CROW(k), LU_CROW(k+1)-1
106              WR( LU_ICOL(kk) ) = JVSR(kk)
107              WI( LU_ICOL(kk) ) = JVSI(kk)
108        END DO
109        DO kk = LU_CROW(k), LU_DIAG(k)-1
110            j = LU_ICOL(kk)
111            den = JVSR(LU_DIAG(j))**2 + JVSI(LU_DIAG(j))**2
112            ar = -(WR(j)*JVSR(LU_DIAG(j)) + WI(j)*JVSI(LU_DIAG(j)))/den
113            ai = -(WI(j)*JVSR(LU_DIAG(j)) - WR(j)*JVSI(LU_DIAG(j)))/den
114            WR(j) = -ar
115            WI(j) = -ai
116            DO jj = LU_DIAG(j)+1, LU_CROW(j+1)-1
117               WR( LU_ICOL(jj) ) = WR( LU_ICOL(jj) ) + ar*JVSR(jj) - ai*JVSI(jj)
118               WI( LU_ICOL(jj) ) = WI( LU_ICOL(jj) ) + ar*JVSI(jj) + ai*JVSR(jj)
119            END DO
120         END DO
121         DO kk = LU_CROW(k), LU_CROW(k+1)-1
122            JVSR(kk) = WR( LU_ICOL(kk) )
123            JVSI(kk) = WI( LU_ICOL(kk) )
124         END DO
125      END DO
126
127END SUBROUTINE KppDecompCmplxR
128
129
130! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
131SUBROUTINE KppSolveIndirect( JVS, X )
132! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
133!        Sparse solve subroutine using indirect addressing
134! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
135
136  USE KPP_ROOT_Parameters
137  USE KPP_ROOT_JacobianSP
138
139      INTEGER  :: i, j
140      KPP_REAL :: JVS(LU_NONZERO), X(NVAR), sum
141
142      DO i=1,NVAR
143         DO j = LU_CROW(i), LU_DIAG(i)-1 
144             X(i) = X(i) - JVS(j)*X(LU_ICOL(j));
145         END DO 
146      END DO
147
148      DO i=NVAR,1,-1
149        sum = X(i);
150        DO j = LU_DIAG(i)+1, LU_CROW(i+1)-1
151          sum = sum - JVS(j)*X(LU_ICOL(j));
152        END DO
153        X(i) = sum/JVS(LU_DIAG(i));
154      END DO
155     
156END SUBROUTINE KppSolveIndirect
157
158
159! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
160SUBROUTINE KppSolveTRIndirect( JVS, X )
161! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
162!        Complex sparse solve transpose subroutine using indirect addressing
163! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
164
165  USE KPP_ROOT_Parameters
166  USE KPP_ROOT_JacobianSP
167
168      INTEGER       :: i, j
169      KPP_REAL :: JVS(LU_NONZERO), X(NVAR)
170
171      DO i=1,NVAR
172        X(i) = X(i)/JVS(LU_DIAG(i))
173        ! subtract all nonzero elements in row i of JVS from X
174        DO j=LU_DIAG(i)+1,LU_CROW(i+1)-1
175          X(LU_ICOL(j)) = X(LU_ICOL(j))-JVS(j)*X(i)
176        END DO
177      END DO
178
179      DO i=NVAR, 1, -1
180        ! subtract all nonzero elements in row i of JVS from X
181        DO j=LU_CROW(i),LU_DIAG(i)-1
182          X(LU_ICOL(j)) = X(LU_ICOL(j))-JVS(j)*X(i)
183        END DO
184      END DO
185     
186END SUBROUTINE KppSolveTRIndirect
187
188
189! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
190SUBROUTINE KppSolveCmplx( JVS, X )
191! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
192!        Complex sparse solve subroutine using indirect addressing
193! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
194
195  USE KPP_ROOT_Parameters
196  USE KPP_ROOT_JacobianSP
197
198      INTEGER        :: i, j
199      DOUBLE COMPLEX :: JVS(LU_NONZERO), X(NVAR), sum
200
201      DO i=1,NVAR
202         DO j = LU_CROW(i), LU_DIAG(i)-1 
203             X(i) = X(i) - JVS(j)*X(LU_ICOL(j));
204         END DO 
205      END DO
206
207      DO i=NVAR,1,-1
208        sum = X(i);
209        DO j = LU_DIAG(i)+1, LU_CROW(i+1)-1
210          sum = sum - JVS(j)*X(LU_ICOL(j));
211        END DO
212        X(i) = sum/JVS(LU_DIAG(i));
213      END DO
214     
215END SUBROUTINE KppSolveCmplx
216
217! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
218SUBROUTINE KppSolveCmplxR( JVSR, JVSI, XR, XI )
219! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
220!   Complex sparse solve subroutine using indirect addressing
221!   (Real and Imaginary parts are used instead of complex data type)     
222! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
223
224  USE KPP_ROOT_Parameters
225  USE KPP_ROOT_JacobianSP
226
227      INTEGER       ::  i, j
228      KPP_REAL ::  JVSR(LU_NONZERO), JVSI(LU_NONZERO), XR(NVAR), XI(NVAR), sumr, sumi, den
229
230      DO i=1,NVAR
231         DO j = LU_CROW(i), LU_DIAG(i)-1 
232             XR(i) = XR(i) - (JVSR(j)*XR(LU_ICOL(j)) - JVSI(j)*XI(LU_ICOL(j)))
233             XI(i) = XI(i) - (JVSR(j)*XI(LU_ICOL(j)) + JVSI(j)*XR(LU_ICOL(j)))
234         END DO 
235      END DO
236
237      DO i=NVAR,1,-1
238        sumr = XR(i); sumi = XI(i)
239        DO j = LU_DIAG(i)+1, LU_CROW(i+1)-1
240            sumr = sumr - (JVSR(j)*XR(LU_ICOL(j)) - JVSI(j)*XI(LU_ICOL(j)))
241            sumi = sumi - (JVSR(j)*XI(LU_ICOL(j)) + JVSI(j)*XR(LU_ICOL(j)))
242        END DO
243        den   = JVSR(LU_DIAG(i))**2 + JVSI(LU_DIAG(i))**2
244        XR(i) = (sumr*JVSR(LU_DIAG(i)) + sumi*JVSI(LU_DIAG(i)))/den
245        XI(i) = (sumi*JVSR(LU_DIAG(i)) - sumr*JVSI(LU_DIAG(i)))/den
246      END DO
247     
248END SUBROUTINE KppSolveCmplxR
249
250
251! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
252SUBROUTINE KppSolveTRCmplx( JVS, X )
253! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
254!        Complex sparse solve transpose subroutine using indirect addressing
255! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
256
257  USE KPP_ROOT_Parameters
258  USE KPP_ROOT_JacobianSP
259
260      INTEGER        :: i, j
261      DOUBLE COMPLEX :: JVS(LU_NONZERO), X(NVAR)
262
263      DO i=1,NVAR
264        X(i) = X(i)/JVS(LU_DIAG(i))
265        ! subtract all nonzero elements in row i of JVS from X
266        DO j=LU_DIAG(i)+1,LU_CROW(i+1)-1
267          X(LU_ICOL(j)) = X(LU_ICOL(j))-JVS(j)*X(i)
268        END DO
269      END DO
270
271      DO i=NVAR, 1, -1
272        ! subtract all nonzero elements in row i of JVS from X
273        DO j=LU_CROW(i),LU_DIAG(i)-1
274          X(LU_ICOL(j)) = X(LU_ICOL(j))-JVS(j)*X(i)
275        END DO
276      END DO
277     
278END SUBROUTINE KppSolveTRCmplx
279
280
281! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
282SUBROUTINE KppSolveTRCmplxR( JVSR, JVSI, XR, XI )
283! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
284!   Complex sparse solve transpose subroutine using indirect addressing
285!   (Real and Imaginary parts are used instead of complex data type)     
286! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
287
288  USE KPP_ROOT_Parameters
289  USE KPP_ROOT_JacobianSP
290
291      INTEGER       ::  i, j
292      KPP_REAL ::  JVSR(LU_NONZERO), JVSI(LU_NONZERO), XR(NVAR), XI(NVAR), den
293
294      DO i=1,NVAR
295        den   = JVSR(LU_DIAG(i))**2 + JVSI(LU_DIAG(i))**2
296        XR(i) = (XR(i)*JVSR(LU_DIAG(i)) + XI(i)*JVSI(LU_DIAG(i)))/den
297        XI(i) = (XI(i)*JVSR(LU_DIAG(i)) - XR(i)*JVSI(LU_DIAG(i)))/den
298        ! subtract all nonzero elements in row i of JVS from X
299        DO j=LU_DIAG(i)+1,LU_CROW(i+1)-1
300          XR(LU_ICOL(j)) = XR(LU_ICOL(j))-(JVSR(j)*XR(i) - JVSI(j)*XI(i))
301          XI(LU_ICOL(j)) = XI(LU_ICOL(j))-(JVSI(j)*XR(i) + JVSR(j)*XI(i))
302        END DO
303      END DO
304
305      DO i=NVAR, 1, -1
306        ! subtract all nonzero elements in row i of JVS from X
307        DO j=LU_CROW(i),LU_DIAG(i)-1
308          XR(LU_ICOL(j)) = XR(LU_ICOL(j))-(JVSR(j)*XR(i) - JVSI(j)*XI(i))
309          XI(LU_ICOL(j)) = XI(LU_ICOL(j))-(JVSI(j)*XR(i) + JVSR(j)*XI(i))
310        END DO
311      END DO
312     
313END SUBROUTINE KppSolveTRCmplxR
314
315
316!
317! Next few commented subroutines perform sparse big linear algebra
318!
319!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
320!SUBROUTINE KppDecompBig( JVS, IP, IER )
321!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
322!!        Sparse LU factorization
323!!        for the Runge Kutta (3n)x(3n) linear system
324!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
325!
326!  USE KPP_ROOT_Parameters
327!  USE KPP_ROOT_JacobianSP
328!
329!      INTEGER  :: IP3(3), IER, IP(3,NVAR)
330!      KPP_REAL :: JVS(3,3,LU_NONZERO), W(3,3,NVAR), a(3,3), E(3,3)
331!      INTEGER  :: k, kk, j, jj
332!
333!      a = 0.0d0
334!      IER = 0
335!      DO k=1,NVAR
336!        DO kk = LU_CROW(k), LU_CROW(k+1)-1
337!              W( 1:3,1:3,LU_ICOL(kk) ) = JVS(1:3,1:3,kk)
338!        END DO
339!        DO kk = LU_CROW(k), LU_DIAG(k)-1
340!            j = LU_ICOL(kk)
341!            E(1:3,1:3) = JVS( 1:3,1:3,LU_DIAG(j) )
342!            ! CALL DGETRF(3,3,E,3,IP3,IER)
343!            CALL FAC3(E,IP3,IER)
344!            IF ( IER /= 0 )  RETURN
345!            ! a = W(j) / JVS( LU_DIAG(j) )
346!            a(1:3,1:3) = W( 1:3,1:3,j )
347!            ! CALL DGETRS ('N',3,3,E,3,IP3,a,3,IER)
348!            CALL SOL3('N',E,IP3,a(1,1))
349!            CALL SOL3('N',E,IP3,a(1,2))
350!            CALL SOL3('N',E,IP3,a(1,3))
351!            W(1:3,1:3,j) = a(1:3,1:3)
352!            DO jj = LU_DIAG(j)+1, LU_CROW(j+1)-1
353!               W( 1:3,1:3,LU_ICOL(jj) ) = W( 1:3,1:3,LU_ICOL(jj) ) &
354!                        - MATMUL( a(1:3,1:3) , JVS(1:3,1:3,jj) )
355!            END DO
356!         END DO
357!         DO kk = LU_CROW(k), LU_CROW(k+1)-1
358!            JVS(1:3,1:3,kk) = W( 1:3,1:3,LU_ICOL(kk) )
359!         END DO
360!      END DO
361!
362!      DO k=1,NVAR
363!         ! CALL WGEFA(JVS(1,1,LU_DIAG(k)),3,3,IP(1,k),IER)
364!         ! CALL DGETRF(3,3,JVS(1,1,LU_DIAG(k)),3,IP(1,k),IER)
365!         CALL FAC3(JVS(1,1,LU_DIAG(k)),IP(1,k),IER)
366!         IF ( IER /= 0 )  RETURN
367!      END DO
368!     
369!END SUBROUTINE KppDecompBig
370!
371!
372!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
373!SUBROUTINE KppSolveBig( JVS, IP, X )
374!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
375!!        Sparse solve subroutine using indirect addressing
376!!        for the Runge Kutta (3n)x(3n) linear system
377!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
378!
379!  USE KPP_ROOT_Parameters
380!  USE KPP_ROOT_JacobianSP
381!
382!      INTEGER  :: i, j, k, m, IP3(3), IP(3,NVAR), IER
383!      KPP_REAL :: JVS(3,3,LU_NONZERO), X(3,NVAR), sum(3)
384!
385!      DO i=1,NVAR
386!        DO j = LU_CROW(i), LU_DIAG(i)-1
387!          !X(1:3,i) = X(1:3,i) - MATMUL(JVS(1:3,1:3,j),X(1:3,LU_ICOL(j)));
388!          DO k=1,3
389!            DO m=1,3
390!              X(k,i) = X(k,i) - JVS(k,m,j)*X(m,LU_ICOL(j))
391!            END DO
392!          END DO
393!        END DO 
394!      END DO
395!
396!      DO i=NVAR,1,-1
397!        sum(1:3) = X(1:3,i);
398!        DO j = LU_DIAG(i)+1, LU_CROW(i+1)-1
399!          !sum(1:3) = sum(1:3) - MATMUL(JVS(1:3,1:3,j),X(1:3,LU_ICOL(j)));
400!          DO k=1,3
401!            DO m=1,3
402!              sum(k) = sum(k) - JVS(k,m,j)*X(m,LU_ICOL(j))
403!            END DO
404!          END DO
405!        END DO
406!        ! X(i) = sum/JVS(LU_DIAG(i));
407!        ! CALL DGETRS ('N',3,1,JVS(1:3,1:3,LU_DIAG(i)),3,IP(1,i),sum,3,0)
408!        ! CALL WGESL('N',JVS(1,1,LU_DIAG(i)),3,3,IP(1,i),sum)
409!        CALL SOL3('N',JVS(1,1,LU_DIAG(i)),IP(1,i),sum)
410!        X(1:3,i) = sum(1:3)
411!      END DO
412!     
413!END SUBROUTINE KppSolveBig
414!
415!
416!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
417!SUBROUTINE KppSolveBigTR( JVS, IP, X )
418!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
419!!        Big sparse transpose solve using indirect addressing
420!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
421!
422!  USE KPP_ROOT_Parameters
423!  USE KPP_ROOT_JacobianSP
424!
425!      INTEGER       :: i, j, k, m, IP(3,NVAR)
426!      KPP_REAL :: JVS(3,3,LU_NONZERO), X(3,NVAR)
427!
428!      DO i=1,NVAR
429!        ! X(i) = X(i)/JVS(LU_DIAG(i))
430!        CALL SOL3('T',JVS(1,1,LU_DIAG(i)),IP(1,i),X(1,i))
431!        DO j=LU_DIAG(i)+1,LU_CROW(i+1)-1
432!         !X(1:3,LU_ICOL(j)) = X(1:3,LU_ICOL(j)) &
433!          !    - MATMUL( TRANSPOSE(JVS(1:3,1:3,j)), X(1:3,i) )
434!          DO k=1,3
435!            DO m=1,3
436!              X(k,LU_ICOL(j)) = X(k,LU_ICOL(j)) - JVS(m,k,j)*X(m,i)
437!            END DO
438!          END DO
439!       END DO
440!      END DO
441!
442!      DO i=NVAR, 1, -1
443!        DO j=LU_CROW(i),LU_DIAG(i)-1
444!         !X(1:3,LU_ICOL(j)) = X(1:3,LU_ICOL(j)) &
445!          !   - MATMUL( TRANSPOSE(JVS(1:3,1:3,j)), X(1:3,i) )
446!          DO k=1,3
447!            DO m=1,3
448!              X(k,LU_ICOL(j)) = X(k,LU_ICOL(j)) - JVS(m,k,j)*X(m,i)
449!            END DO
450!          END DO
451!       END DO
452!      END DO
453!     
454!END SUBROUTINE KppSolveBigTR
455!
456!
457!
458!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
459!SUBROUTINE FAC3(A,IPVT,INFO)
460!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
461!!     FAC3 FACTORS THE MATRIX A (3,3) BY
462!!           GAUSS ELIMINATION WITH PARTIAL PIVOTING
463!!     LINPACK - LIKE
464!!
465!!     Remove comments to perform pivoting
466!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
467!!
468!      KPP_REAL :: A(3,3)
469!      INTEGER       :: IPVT(3),INFO
470!!      INTEGER       :: L
471!!      KPP_REAL :: t, dmax, da, TMP(3)
472!      KPP_REAL, PARAMETER :: ZERO = 0.0, ONE = 1.0
473!
474!      info = 0
475!!      t = TINY(da)
476!!     
477!!      da = ABS(A(1,1)); L = 1
478!!      IF ( ABS(A(2,1))>da ) THEN
479!!        da = ABS(A(2,1)); L = 2
480!!        IF ( ABS(A(3,1))>da ) THEN
481!!          L = 3
482!!        END IF 
483!!      END IF 
484!!      IPVT(1)  = L
485!!      IF (L /=1 ) THEN
486!!         TMP(1:3) = A(L,1:3)
487!!         A(L,1:3) = A(1,1:3)
488!!         A(1,1:3) = TMP(1:3)
489!!      END IF
490!!      IF (ABS(A(1,1)) < t) THEN
491!!         info = 1
492!!         return
493!!      END IF   
494!!
495!      A(2,1) = A(2,1)/A(1,1)
496!      A(2,2) = A(2,2) - A(2,1)*A(1,2)
497!      A(2,3) = A(2,3) - A(2,1)*A(1,3)
498!      A(3,1) = A(3,1)/A(1,1)
499!      A(3,2) = A(3,2) - A(3,1)*A(1,2)
500!      A(3,3) = A(3,3) - A(3,1)*A(1,3)
501!     
502!!      IPVT(2)  = 2
503!!      IF (ABS(A(3,2))>ABS(A(2,2))) THEN
504!!         IPVT(2)  = 3
505!!         TMP(2:3) = A(3,2:3)
506!!         A(3,2:3) = A(2,2:3)
507!!         A(2,2:3) = TMP(2:3)
508!!      END IF
509!!      IF (ABS(A(2,2)) < t) THEN
510!!         info = 1
511!!         return
512!!      END IF   
513!!     
514!      A(3,2)   = A(3,2)/A(2,2)
515!      A(3,3)   = A(3,3) - A(3,2)*A(2,3)
516!      IPVT(3)  = 3
517!     
518!END SUBROUTINE FAC3
519!
520!
521!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
522!SUBROUTINE SOL3(Trans,A,IPVT,b)
523!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
524!!     SOL3 solves the system 3x3
525!!     A * x = b  or  trans(a) * x = b
526!!     using the factors computed by WGEFA.
527!!
528!!     Trans      = 'N'   to solve  A*x = b ,
529!!                = 'T'   to solve  transpose(A)*x = b
530!!     LINPACK - LIKE
531!!
532!!     Remove comments to use pivoting
533!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
534!
535!      CHARACTER     :: Trans
536!      KPP_REAL :: a(3,3),b(3)
537!      INTEGER       :: IPVT(3)
538!!      INTEGER       :: L
539!!      KPP_REAL :: TMP
540!     
541!      SELECT CASE (Trans)
542!
543!      CASE ('n','N')  !  Solve  A * x = b
544!
545!!     Solve  L*y = b
546!!         L = IPVT(1)
547!!         IF (L /= 1) THEN
548!!            TMP = B(1); B(1) = B(L); B(L) = TMP
549!!         END IF
550!         b(2) = b(2)-A(2,1)*b(1)
551!         b(3) = b(3)-A(3,1)*b(1)
552!         
553!!         L = IPVT(2)
554!!         IF (L /= 2) THEN
555!!            TMP = B(2); B(2) = B(L); B(L) = TMP
556!!         END IF
557!         b(3) = b(3)-A(3,2)*b(2)
558!
559!!     Solve  U*x = y
560!         b(3) = b(3)/A(3,3)
561!         b(2) = (b(2)-A(2,3)*b(3))/A(2,2)
562!         b(1) = (b(1)-A(1,3)*b(3)-A(1,2)*b(2))/A(1,1)
563!     
564!     
565!      CASE ('t','T')  !  Solve transpose(A) * x = b
566!
567!!      Solve transpose(U)*y = b
568!         b(1) = b(1)/A(1,1)
569!         b(2) = (b(2)-A(1,2)*b(1))/A(2,2)
570!         b(3) = (b(3)-A(1,3)*b(1)-A(2,3)*b(2))/A(3,3)
571!
572!!      Solve transpose(L)*x = y
573!         b(2) = b(2)-A(3,2)*b(3)
574!!         L = ipvt(2)
575!!         IF (L /= 2) THEN
576!!            TMP = B(2); B(2) = B(L); B(L) = TMP
577!!         END IF
578!         b(1) = b(1)-A(3,1)*b(3)-A(2,1)*b(2)
579!!         L = ipvt(1)
580!!         IF (L /= 1) THEN
581!!            TMP = B(1); B(1) = B(L); B(L) = TMP
582!!         END IF
583!   
584!      END SELECT
585!
586!END SUBROUTINE SOL3
587
Note: See TracBrowser for help on using the repository browser.