Home

Context Navigation

source: palm/trunk/UTIL/chemistry/gasphase_preproc/tmp_kpp4palm/chem_gasphase_mod_LinearAlgebra.f90 @ 2696

Last change on this file since 2696 was 2696, checked in by kanani, 6 years ago
Merge of branch palm4u into trunk
File size: 34.3 KB

Line
1	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2	!
3	! Linear Algebra Data and Routines File
4	!
5	! Generated by KPP-2.2.3 symbolic chemistry Kinetics PreProcessor
6	! (http://www.cs.vt.edu/~asandu/Software/KPP)
7	! KPP is distributed under GPL, the general public licence
8	! (http://www.gnu.org/copyleft/gpl.html)
9	! (C) 1995-1997, V. Damian & A. Sandu, CGRER, Univ. Iowa
10	! (C) 1997-2005, A. Sandu, Michigan Tech, Virginia Tech
11	! With important contributions from:
12	! M. Damian, Villanova University, USA
13	! R. Sander, Max-Planck Institute for Chemistry, Mainz, Germany
14	!
15	! File : chem_gasphase_mod_LinearAlgebra.f90
16	! Time : Fri Dec 1 18:10:53 2017
17	! Working directory : /data/kanani/branches/palm4u/GASPHASE_PREPROC/tmp_kpp4palm
18	! Equation file : chem_gasphase_mod.kpp
19	! Output root filename : chem_gasphase_mod
20	!
21	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
22
23
24
25	MODULE chem_gasphase_mod_LinearAlgebra
26
27	USE chem_gasphase_mod_Parameters
28	USE chem_gasphase_mod_JacobianSP
29
30	IMPLICIT NONE
31
32	CONTAINS
33
34
35	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
36	!
37	! SPARSE_UTIL - SPARSE utility functions
38	! Arguments :
39	!
40	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
41
42
43	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
44	SUBROUTINE KppDecomp( JVS, IER )
45	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
46	! Sparse LU factorization
47	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
48
49	USE chem_gasphase_mod_Parameters
50	USE chem_gasphase_mod_JacobianSP
51
52	INTEGER :: IER
53	REAL(kind=dp) :: JVS(LU_NONZERO), W(NVAR), a
54	INTEGER :: k, kk, j, jj
55
56	a = 0. ! mz_rs_20050606
57	IER = 0
58	DO k=1,NVAR
59	! mz_rs_20050606: don't check if real value == 0
60	! IF ( JVS( LU_DIAG(k) ) .EQ. 0. ) THEN
61	IF ( ABS(JVS(LU_DIAG(k))) < TINY(a) ) THEN
62	IER = k
63	RETURN
64	END IF
65	DO kk = LU_CROW(k), LU_CROW(k+1)-1
66	W( LU_ICOL(kk) ) = JVS(kk)
67	END DO
68	DO kk = LU_CROW(k), LU_DIAG(k)-1
69	j = LU_ICOL(kk)
70	a = -W(j) / JVS( LU_DIAG(j) )
71	W(j) = -a
72	DO jj = LU_DIAG(j)+1, LU_CROW(j+1)-1
73	W( LU_ICOL(jj) ) = W( LU_ICOL(jj) ) + a*JVS(jj)
74	END DO
75	END DO
76	DO kk = LU_CROW(k), LU_CROW(k+1)-1
77	JVS(kk) = W( LU_ICOL(kk) )
78	END DO
79	END DO
80
81	END SUBROUTINE KppDecomp
82
83
84	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
85	SUBROUTINE KppDecompCmplx( JVS, IER )
86	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
87	! Sparse LU factorization, complex
88	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
89
90	USE chem_gasphase_mod_Parameters
91	USE chem_gasphase_mod_JacobianSP
92
93	INTEGER :: IER
94	DOUBLE COMPLEX :: JVS(LU_NONZERO), W(NVAR), a
95	REAL(kind=dp) :: b = 0.0
96	INTEGER :: k, kk, j, jj
97
98	IER = 0
99	DO k=1,NVAR
100	IF ( ABS(JVS(LU_DIAG(k))) < TINY(b) ) THEN
101	IER = k
102	RETURN
103	END IF
104	DO kk = LU_CROW(k), LU_CROW(k+1)-1
105	W( LU_ICOL(kk) ) = JVS(kk)
106	END DO
107	DO kk = LU_CROW(k), LU_DIAG(k)-1
108	j = LU_ICOL(kk)
109	a = -W(j) / JVS( LU_DIAG(j) )
110	W(j) = -a
111	DO jj = LU_DIAG(j)+1, LU_CROW(j+1)-1
112	W( LU_ICOL(jj) ) = W( LU_ICOL(jj) ) + a*JVS(jj)
113	END DO
114	END DO
115	DO kk = LU_CROW(k), LU_CROW(k+1)-1
116	JVS(kk) = W( LU_ICOL(kk) )
117	END DO
118	END DO
119
120	END SUBROUTINE KppDecompCmplx
121
122
123	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
124	SUBROUTINE KppDecompCmplxR( JVSR, JVSI, IER )
125	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
126	! Sparse LU factorization, complex
127	! (Real and Imaginary parts are used instead of complex data type)
128	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
129
130	USE chem_gasphase_mod_Parameters
131	USE chem_gasphase_mod_JacobianSP
132
133	INTEGER :: IER
134	REAL(kind=dp) :: JVSR(LU_NONZERO), JVSI(LU_NONZERO)
135	REAL(kind=dp) :: WR(NVAR), WI(NVAR), ar, ai, den
136	INTEGER :: k, kk, j, jj
137
138	IER = 0
139	ar = 0.0
140	DO k=1,NVAR
141	IF ( ( ABS(JVSR(LU_DIAG(k))) < TINY(ar) ) .AND. &
142	( ABS(JVSI(LU_DIAG(k))) < TINY(ar) ) ) THEN
143	IER = k
144	RETURN
145	END IF
146	DO kk = LU_CROW(k), LU_CROW(k+1)-1
147	WR( LU_ICOL(kk) ) = JVSR(kk)
148	WI( LU_ICOL(kk) ) = JVSI(kk)
149	END DO
150	DO kk = LU_CROW(k), LU_DIAG(k)-1
151	j = LU_ICOL(kk)
152	den = JVSR(LU_DIAG(j))2 + JVSI(LU_DIAG(j))2
153	ar = -(WR(j)JVSR(LU_DIAG(j)) + WI(j)JVSI(LU_DIAG(j)))/den
154	ai = -(WI(j)JVSR(LU_DIAG(j)) - WR(j)JVSI(LU_DIAG(j)))/den
155	WR(j) = -ar
156	WI(j) = -ai
157	DO jj = LU_DIAG(j)+1, LU_CROW(j+1)-1
158	WR( LU_ICOL(jj) ) = WR( LU_ICOL(jj) ) + arJVSR(jj) - aiJVSI(jj)
159	WI( LU_ICOL(jj) ) = WI( LU_ICOL(jj) ) + arJVSI(jj) + aiJVSR(jj)
160	END DO
161	END DO
162	DO kk = LU_CROW(k), LU_CROW(k+1)-1
163	JVSR(kk) = WR( LU_ICOL(kk) )
164	JVSI(kk) = WI( LU_ICOL(kk) )
165	END DO
166	END DO
167
168	END SUBROUTINE KppDecompCmplxR
169
170
171	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
172	SUBROUTINE KppSolveIndirect( JVS, X )
173	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
174	! Sparse solve subroutine using indirect addressing
175	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
176
177	USE chem_gasphase_mod_Parameters
178	USE chem_gasphase_mod_JacobianSP
179
180	INTEGER :: i, j
181	REAL(kind=dp) :: JVS(LU_NONZERO), X(NVAR), sum
182
183	DO i=1,NVAR
184	DO j = LU_CROW(i), LU_DIAG(i)-1
185	X(i) = X(i) - JVS(j)*X(LU_ICOL(j));
186	END DO
187	END DO
188
189	DO i=NVAR,1,-1
190	sum = X(i);
191	DO j = LU_DIAG(i)+1, LU_CROW(i+1)-1
192	sum = sum - JVS(j)*X(LU_ICOL(j));
193	END DO
194	X(i) = sum/JVS(LU_DIAG(i));
195	END DO
196
197	END SUBROUTINE KppSolveIndirect
198
199
200	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
201	SUBROUTINE KppSolveTRIndirect( JVS, X )
202	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
203	! Complex sparse solve transpose subroutine using indirect addressing
204	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
205
206	USE chem_gasphase_mod_Parameters
207	USE chem_gasphase_mod_JacobianSP
208
209	INTEGER :: i, j
210	REAL(kind=dp) :: JVS(LU_NONZERO), X(NVAR)
211
212	DO i=1,NVAR
213	X(i) = X(i)/JVS(LU_DIAG(i))
214	! subtract all nonzero elements in row i of JVS from X
215	DO j=LU_DIAG(i)+1,LU_CROW(i+1)-1
216	X(LU_ICOL(j)) = X(LU_ICOL(j))-JVS(j)*X(i)
217	END DO
218	END DO
219
220	DO i=NVAR, 1, -1
221	! subtract all nonzero elements in row i of JVS from X
222	DO j=LU_CROW(i),LU_DIAG(i)-1
223	X(LU_ICOL(j)) = X(LU_ICOL(j))-JVS(j)*X(i)
224	END DO
225	END DO
226
227	END SUBROUTINE KppSolveTRIndirect
228
229
230	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
231	SUBROUTINE KppSolveCmplx( JVS, X )
232	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
233	! Complex sparse solve subroutine using indirect addressing
234	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
235
236	USE chem_gasphase_mod_Parameters
237	USE chem_gasphase_mod_JacobianSP
238
239	INTEGER :: i, j
240	DOUBLE COMPLEX :: JVS(LU_NONZERO), X(NVAR), sum
241
242	DO i=1,NVAR
243	DO j = LU_CROW(i), LU_DIAG(i)-1
244	X(i) = X(i) - JVS(j)*X(LU_ICOL(j));
245	END DO
246	END DO
247
248	DO i=NVAR,1,-1
249	sum = X(i);
250	DO j = LU_DIAG(i)+1, LU_CROW(i+1)-1
251	sum = sum - JVS(j)*X(LU_ICOL(j));
252	END DO
253	X(i) = sum/JVS(LU_DIAG(i));
254	END DO
255
256	END SUBROUTINE KppSolveCmplx
257
258	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
259	SUBROUTINE KppSolveCmplxR( JVSR, JVSI, XR, XI )
260	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
261	! Complex sparse solve subroutine using indirect addressing
262	! (Real and Imaginary parts are used instead of complex data type)
263	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
264
265	USE chem_gasphase_mod_Parameters
266	USE chem_gasphase_mod_JacobianSP
267
268	INTEGER :: i, j
269	REAL(kind=dp) :: JVSR(LU_NONZERO), JVSI(LU_NONZERO), XR(NVAR), XI(NVAR), sumr, sumi, den
270
271	DO i=1,NVAR
272	DO j = LU_CROW(i), LU_DIAG(i)-1
273	XR(i) = XR(i) - (JVSR(j)XR(LU_ICOL(j)) - JVSI(j)XI(LU_ICOL(j)))
274	XI(i) = XI(i) - (JVSR(j)XI(LU_ICOL(j)) + JVSI(j)XR(LU_ICOL(j)))
275	END DO
276	END DO
277
278	DO i=NVAR,1,-1
279	sumr = XR(i); sumi = XI(i)
280	DO j = LU_DIAG(i)+1, LU_CROW(i+1)-1
281	sumr = sumr - (JVSR(j)XR(LU_ICOL(j)) - JVSI(j)XI(LU_ICOL(j)))
282	sumi = sumi - (JVSR(j)XI(LU_ICOL(j)) + JVSI(j)XR(LU_ICOL(j)))
283	END DO
284	den = JVSR(LU_DIAG(i))2 + JVSI(LU_DIAG(i))2
285	XR(i) = (sumrJVSR(LU_DIAG(i)) + sumiJVSI(LU_DIAG(i)))/den
286	XI(i) = (sumiJVSR(LU_DIAG(i)) - sumrJVSI(LU_DIAG(i)))/den
287	END DO
288
289	END SUBROUTINE KppSolveCmplxR
290
291
292	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
293	SUBROUTINE KppSolveTRCmplx( JVS, X )
294	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
295	! Complex sparse solve transpose subroutine using indirect addressing
296	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
297
298	USE chem_gasphase_mod_Parameters
299	USE chem_gasphase_mod_JacobianSP
300
301	INTEGER :: i, j
302	DOUBLE COMPLEX :: JVS(LU_NONZERO), X(NVAR)
303
304	DO i=1,NVAR
305	X(i) = X(i)/JVS(LU_DIAG(i))
306	! subtract all nonzero elements in row i of JVS from X
307	DO j=LU_DIAG(i)+1,LU_CROW(i+1)-1
308	X(LU_ICOL(j)) = X(LU_ICOL(j))-JVS(j)*X(i)
309	END DO
310	END DO
311
312	DO i=NVAR, 1, -1
313	! subtract all nonzero elements in row i of JVS from X
314	DO j=LU_CROW(i),LU_DIAG(i)-1
315	X(LU_ICOL(j)) = X(LU_ICOL(j))-JVS(j)*X(i)
316	END DO
317	END DO
318
319	END SUBROUTINE KppSolveTRCmplx
320
321
322	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
323	SUBROUTINE KppSolveTRCmplxR( JVSR, JVSI, XR, XI )
324	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
325	! Complex sparse solve transpose subroutine using indirect addressing
326	! (Real and Imaginary parts are used instead of complex data type)
327	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
328
329	USE chem_gasphase_mod_Parameters
330	USE chem_gasphase_mod_JacobianSP
331
332	INTEGER :: i, j
333	REAL(kind=dp) :: JVSR(LU_NONZERO), JVSI(LU_NONZERO), XR(NVAR), XI(NVAR), den
334
335	DO i=1,NVAR
336	den = JVSR(LU_DIAG(i))2 + JVSI(LU_DIAG(i))2
337	XR(i) = (XR(i)JVSR(LU_DIAG(i)) + XI(i)JVSI(LU_DIAG(i)))/den
338	XI(i) = (XI(i)JVSR(LU_DIAG(i)) - XR(i)JVSI(LU_DIAG(i)))/den
339	! subtract all nonzero elements in row i of JVS from X
340	DO j=LU_DIAG(i)+1,LU_CROW(i+1)-1
341	XR(LU_ICOL(j)) = XR(LU_ICOL(j))-(JVSR(j)XR(i) - JVSI(j)XI(i))
342	XI(LU_ICOL(j)) = XI(LU_ICOL(j))-(JVSI(j)XR(i) + JVSR(j)XI(i))
343	END DO
344	END DO
345
346	DO i=NVAR, 1, -1
347	! subtract all nonzero elements in row i of JVS from X
348	DO j=LU_CROW(i),LU_DIAG(i)-1
349	XR(LU_ICOL(j)) = XR(LU_ICOL(j))-(JVSR(j)XR(i) - JVSI(j)XI(i))
350	XI(LU_ICOL(j)) = XI(LU_ICOL(j))-(JVSI(j)XR(i) + JVSR(j)XI(i))
351	END DO
352	END DO
353
354	END SUBROUTINE KppSolveTRCmplxR
355
356
357	!
358	! Next few commented subroutines perform sparse big linear algebra
359	!
360	!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
361	!SUBROUTINE KppDecompBig( JVS, IP, IER )
362	!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
363	!! Sparse LU factorization
364	!! for the Runge Kutta (3n)x(3n) linear system
365	!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
366	!
367	! USE chem_gasphase_mod_Parameters
368	! USE chem_gasphase_mod_JacobianSP
369	!
370	! INTEGER :: IP3(3), IER, IP(3,NVAR)
371	! REAL(kind=dp) :: JVS(3,3,LU_NONZERO), W(3,3,NVAR), a(3,3), E(3,3)
372	! INTEGER :: k, kk, j, jj
373	!
374	! a = 0.0d0
375	! IER = 0
376	! DO k=1,NVAR
377	! DO kk = LU_CROW(k), LU_CROW(k+1)-1
378	! W( 1:3,1:3,LU_ICOL(kk) ) = JVS(1:3,1:3,kk)
379	! END DO
380	! DO kk = LU_CROW(k), LU_DIAG(k)-1
381	! j = LU_ICOL(kk)
382	! E(1:3,1:3) = JVS( 1:3,1:3,LU_DIAG(j) )
383	! ! CALL DGETRF(3,3,E,3,IP3,IER)
384	! CALL FAC3(E,IP3,IER)
385	! IF ( IER /= 0 ) RETURN
386	! ! a = W(j) / JVS( LU_DIAG(j) )
387	! a(1:3,1:3) = W( 1:3,1:3,j )
388	! ! CALL DGETRS ('N',3,3,E,3,IP3,a,3,IER)
389	! CALL SOL3('N',E,IP3,a(1,1))
390	! CALL SOL3('N',E,IP3,a(1,2))
391	! CALL SOL3('N',E,IP3,a(1,3))
392	! W(1:3,1:3,j) = a(1:3,1:3)
393	! DO jj = LU_DIAG(j)+1, LU_CROW(j+1)-1
394	! W( 1:3,1:3,LU_ICOL(jj) ) = W( 1:3,1:3,LU_ICOL(jj) ) &
395	! - MATMUL( a(1:3,1:3) , JVS(1:3,1:3,jj) )
396	! END DO
397	! END DO
398	! DO kk = LU_CROW(k), LU_CROW(k+1)-1
399	! JVS(1:3,1:3,kk) = W( 1:3,1:3,LU_ICOL(kk) )
400	! END DO
401	! END DO
402	!
403	! DO k=1,NVAR
404	! ! CALL WGEFA(JVS(1,1,LU_DIAG(k)),3,3,IP(1,k),IER)
405	! ! CALL DGETRF(3,3,JVS(1,1,LU_DIAG(k)),3,IP(1,k),IER)
406	! CALL FAC3(JVS(1,1,LU_DIAG(k)),IP(1,k),IER)
407	! IF ( IER /= 0 ) RETURN
408	! END DO
409	!
410	!END SUBROUTINE KppDecompBig
411	!
412	!
413	!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
414	!SUBROUTINE KppSolveBig( JVS, IP, X )
415	!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
416	!! Sparse solve subroutine using indirect addressing
417	!! for the Runge Kutta (3n)x(3n) linear system
418	!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
419	!
420	! USE chem_gasphase_mod_Parameters
421	! USE chem_gasphase_mod_JacobianSP
422	!
423	! INTEGER :: i, j, k, m, IP3(3), IP(3,NVAR), IER
424	! REAL(kind=dp) :: JVS(3,3,LU_NONZERO), X(3,NVAR), sum(3)
425	!
426	! DO i=1,NVAR
427	! DO j = LU_CROW(i), LU_DIAG(i)-1
428	! !X(1:3,i) = X(1:3,i) - MATMUL(JVS(1:3,1:3,j),X(1:3,LU_ICOL(j)));
429	! DO k=1,3
430	! DO m=1,3
431	! X(k,i) = X(k,i) - JVS(k,m,j)*X(m,LU_ICOL(j))
432	! END DO
433	! END DO
434	! END DO
435	! END DO
436	!
437	! DO i=NVAR,1,-1
438	! sum(1:3) = X(1:3,i);
439	! DO j = LU_DIAG(i)+1, LU_CROW(i+1)-1
440	! !sum(1:3) = sum(1:3) - MATMUL(JVS(1:3,1:3,j),X(1:3,LU_ICOL(j)));
441	! DO k=1,3
442	! DO m=1,3
443	! sum(k) = sum(k) - JVS(k,m,j)*X(m,LU_ICOL(j))
444	! END DO
445	! END DO
446	! END DO
447	! ! X(i) = sum/JVS(LU_DIAG(i));
448	! ! CALL DGETRS ('N',3,1,JVS(1:3,1:3,LU_DIAG(i)),3,IP(1,i),sum,3,0)
449	! ! CALL WGESL('N',JVS(1,1,LU_DIAG(i)),3,3,IP(1,i),sum)
450	! CALL SOL3('N',JVS(1,1,LU_DIAG(i)),IP(1,i),sum)
451	! X(1:3,i) = sum(1:3)
452	! END DO
453	!
454	!END SUBROUTINE KppSolveBig
455	!
456	!
457	!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
458	!SUBROUTINE KppSolveBigTR( JVS, IP, X )
459	!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
460	!! Big sparse transpose solve using indirect addressing
461	!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
462	!
463	! USE chem_gasphase_mod_Parameters
464	! USE chem_gasphase_mod_JacobianSP
465	!
466	! INTEGER :: i, j, k, m, IP(3,NVAR)
467	! REAL(kind=dp) :: JVS(3,3,LU_NONZERO), X(3,NVAR)
468	!
469	! DO i=1,NVAR
470	! ! X(i) = X(i)/JVS(LU_DIAG(i))
471	! CALL SOL3('T',JVS(1,1,LU_DIAG(i)),IP(1,i),X(1,i))
472	! DO j=LU_DIAG(i)+1,LU_CROW(i+1)-1
473	! !X(1:3,LU_ICOL(j)) = X(1:3,LU_ICOL(j)) &
474	! ! - MATMUL( TRANSPOSE(JVS(1:3,1:3,j)), X(1:3,i) )
475	! DO k=1,3
476	! DO m=1,3
477	! X(k,LU_ICOL(j)) = X(k,LU_ICOL(j)) - JVS(m,k,j)*X(m,i)
478	! END DO
479	! END DO
480	! END DO
481	! END DO
482	!
483	! DO i=NVAR, 1, -1
484	! DO j=LU_CROW(i),LU_DIAG(i)-1
485	! !X(1:3,LU_ICOL(j)) = X(1:3,LU_ICOL(j)) &
486	! ! - MATMUL( TRANSPOSE(JVS(1:3,1:3,j)), X(1:3,i) )
487	! DO k=1,3
488	! DO m=1,3
489	! X(k,LU_ICOL(j)) = X(k,LU_ICOL(j)) - JVS(m,k,j)*X(m,i)
490	! END DO
491	! END DO
492	! END DO
493	! END DO
494	!
495	!END SUBROUTINE KppSolveBigTR
496	!
497	!
498	!
499	!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
500	!SUBROUTINE FAC3(A,IPVT,INFO)
501	!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
502	!! FAC3 FACTORS THE MATRIX A (3,3) BY
503	!! GAUSS ELIMINATION WITH PARTIAL PIVOTING
504	!! LINPACK - LIKE
505	!!
506	!! Remove comments to perform pivoting
507	!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
508	!!
509	! REAL(kind=dp) :: A(3,3)
510	! INTEGER :: IPVT(3),INFO
511	!! INTEGER :: L
512	!! REAL(kind=dp) :: t, dmax, da, TMP(3)
513	! REAL(kind=dp), PARAMETER :: ZERO = 0.0, ONE = 1.0
514	!
515	! info = 0
516	!! t = TINY(da)
517	!!
518	!! da = ABS(A(1,1)); L = 1
519	!! IF ( ABS(A(2,1))>da ) THEN
520	!! da = ABS(A(2,1)); L = 2
521	!! IF ( ABS(A(3,1))>da ) THEN
522	!! L = 3
523	!! END IF
524	!! END IF
525	!! IPVT(1) = L
526	!! IF (L /=1 ) THEN
527	!! TMP(1:3) = A(L,1:3)
528	!! A(L,1:3) = A(1,1:3)
529	!! A(1,1:3) = TMP(1:3)
530	!! END IF
531	!! IF (ABS(A(1,1)) < t) THEN
532	!! info = 1
533	!! return
534	!! END IF
535	!!
536	! A(2,1) = A(2,1)/A(1,1)
537	! A(2,2) = A(2,2) - A(2,1)*A(1,2)
538	! A(2,3) = A(2,3) - A(2,1)*A(1,3)
539	! A(3,1) = A(3,1)/A(1,1)
540	! A(3,2) = A(3,2) - A(3,1)*A(1,2)
541	! A(3,3) = A(3,3) - A(3,1)*A(1,3)
542	!
543	!! IPVT(2) = 2
544	!! IF (ABS(A(3,2))>ABS(A(2,2))) THEN
545	!! IPVT(2) = 3
546	!! TMP(2:3) = A(3,2:3)
547	!! A(3,2:3) = A(2,2:3)
548	!! A(2,2:3) = TMP(2:3)
549	!! END IF
550	!! IF (ABS(A(2,2)) < t) THEN
551	!! info = 1
552	!! return
553	!! END IF
554	!!
555	! A(3,2) = A(3,2)/A(2,2)
556	! A(3,3) = A(3,3) - A(3,2)*A(2,3)
557	! IPVT(3) = 3
558	!
559	!END SUBROUTINE FAC3
560	!
561	!
562	!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
563	!SUBROUTINE SOL3(Trans,A,IPVT,b)
564	!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
565	!! SOL3 solves the system 3x3
566	!! A * x = b or trans(a) * x = b
567	!! using the factors computed by WGEFA.
568	!!
569	!! Trans = 'N' to solve A*x = b ,
570	!! = 'T' to solve transpose(A)*x = b
571	!! LINPACK - LIKE
572	!!
573	!! Remove comments to use pivoting
574	!! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
575	!
576	! CHARACTER :: Trans
577	! REAL(kind=dp) :: a(3,3),b(3)
578	! INTEGER :: IPVT(3)
579	!! INTEGER :: L
580	!! REAL(kind=dp) :: TMP
581	!
582	! SELECT CASE (Trans)
583	!
584	! CASE ('n','N') ! Solve A * x = b
585	!
586	!! Solve L*y = b
587	!! L = IPVT(1)
588	!! IF (L /= 1) THEN
589	!! TMP = B(1); B(1) = B(L); B(L) = TMP
590	!! END IF
591	! b(2) = b(2)-A(2,1)*b(1)
592	! b(3) = b(3)-A(3,1)*b(1)
593	!
594	!! L = IPVT(2)
595	!! IF (L /= 2) THEN
596	!! TMP = B(2); B(2) = B(L); B(L) = TMP
597	!! END IF
598	! b(3) = b(3)-A(3,2)*b(2)
599	!
600	!! Solve U*x = y
601	! b(3) = b(3)/A(3,3)
602	! b(2) = (b(2)-A(2,3)*b(3))/A(2,2)
603	! b(1) = (b(1)-A(1,3)b(3)-A(1,2)b(2))/A(1,1)
604	!
605	!
606	! CASE ('t','T') ! Solve transpose(A) * x = b
607	!
608	!! Solve transpose(U)*y = b
609	! b(1) = b(1)/A(1,1)
610	! b(2) = (b(2)-A(1,2)*b(1))/A(2,2)
611	! b(3) = (b(3)-A(1,3)b(1)-A(2,3)b(2))/A(3,3)
612	!
613	!! Solve transpose(L)*x = y
614	! b(2) = b(2)-A(3,2)*b(3)
615	!! L = ipvt(2)
616	!! IF (L /= 2) THEN
617	!! TMP = B(2); B(2) = B(L); B(L) = TMP
618	!! END IF
619	! b(1) = b(1)-A(3,1)b(3)-A(2,1)b(2)
620	!! L = ipvt(1)
621	!! IF (L /= 1) THEN
622	!! TMP = B(1); B(1) = B(L); B(L) = TMP
623	!! END IF
624	!
625	! END SELECT
626	!
627	!END SUBROUTINE SOL3
628
629	! End of SPARSE_UTIL function
630	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
631
632
633	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
634	!
635	! KppSolve - sparse back substitution
636	! Arguments :
637	! JVS - sparse Jacobian of variables
638	! X - Vector for variables
639	!
640	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
641
642	SUBROUTINE KppSolve ( JVS, X )
643
644	! JVS - sparse Jacobian of variables
645	REAL(kind=dp) :: JVS(LU_NONZERO)
646	! X - Vector for variables
647	REAL(kind=dp) :: X(NVAR)
648
649	X(3) = X(3)/JVS(3)
650	X(2) = X(2)/JVS(2)
651	X(1) = X(1)/JVS(1)
652
653	END SUBROUTINE KppSolve
654
655	! End of KppSolve function
656	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
657
658
659	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
660	!
661	! KppSolveTR - sparse, transposed back substitution
662	! Arguments :
663	! JVS - sparse Jacobian of variables
664	! X - Vector for variables
665	! XX - Vector for output variables
666	!
667	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
668
669	SUBROUTINE KppSolveTR ( JVS, X, XX )
670
671	! JVS - sparse Jacobian of variables
672	REAL(kind=dp) :: JVS(LU_NONZERO)
673	! X - Vector for variables
674	REAL(kind=dp) :: X(NVAR)
675	! XX - Vector for output variables
676	REAL(kind=dp) :: XX(NVAR)
677
678	XX(1) = X(1)/JVS(1)
679	XX(2) = X(2)/JVS(2)
680	XX(3) = X(3)/JVS(3)
681	XX(3) = XX(3)
682	XX(2) = XX(2)
683	XX(1) = XX(1)
684
685	END SUBROUTINE KppSolveTR
686
687	! End of KppSolveTR function
688	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
689
690
691	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
692	!
693	! BLAS_UTIL - BLAS-LIKE utility functions
694	! Arguments :
695	!
696	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
697
698	!--------------------------------------------------------------
699	!
700	! BLAS/LAPACK-like subroutines used by the integration algorithms
701	! It is recommended to replace them by calls to the optimized
702	! BLAS/LAPACK library for your machine
703	!
704	! (C) Adrian Sandu, Aug. 2004
705	! Virginia Polytechnic Institute and State University
706	!--------------------------------------------------------------
707
708
709	!--------------------------------------------------------------
710	SUBROUTINE WCOPY(N,X,incX,Y,incY)
711	!--------------------------------------------------------------
712	! copies a vector, x, to a vector, y: y <- x
713	! only for incX=incY=1
714	! after BLAS
715	! replace this by the function from the optimized BLAS implementation:
716	! CALL SCOPY(N,X,1,Y,1) or CALL DCOPY(N,X,1,Y,1)
717	!--------------------------------------------------------------
718	! USE chem_gasphase_mod_Precision
719
720	INTEGER :: i,incX,incY,M,MP1,N
721	REAL(kind=dp) :: X(N),Y(N)
722
723	IF (N.LE.0) RETURN
724
725	M = MOD(N,8)
726	IF( M .NE. 0 ) THEN
727	DO i = 1,M
728	Y(i) = X(i)
729	END DO
730	IF( N .LT. 8 ) RETURN
731	END IF
732	MP1 = M+1
733	DO i = MP1,N,8
734	Y(i) = X(i)
735	Y(i + 1) = X(i + 1)
736	Y(i + 2) = X(i + 2)
737	Y(i + 3) = X(i + 3)
738	Y(i + 4) = X(i + 4)
739	Y(i + 5) = X(i + 5)
740	Y(i + 6) = X(i + 6)
741	Y(i + 7) = X(i + 7)
742	END DO
743
744	END SUBROUTINE WCOPY
745
746
747	!--------------------------------------------------------------
748	SUBROUTINE WAXPY(N,Alpha,X,incX,Y,incY)
749	!--------------------------------------------------------------
750	! constant times a vector plus a vector: y <- y + Alpha*x
751	! only for incX=incY=1
752	! after BLAS
753	! replace this by the function from the optimized BLAS implementation:
754	! CALL SAXPY(N,Alpha,X,1,Y,1) or CALL DAXPY(N,Alpha,X,1,Y,1)
755	!--------------------------------------------------------------
756
757	INTEGER :: i,incX,incY,M,MP1,N
758	REAL(kind=dp) :: X(N),Y(N),Alpha
759	REAL(kind=dp), PARAMETER :: ZERO = 0.0_dp
760
761	IF (Alpha .EQ. ZERO) RETURN
762	IF (N .LE. 0) RETURN
763
764	M = MOD(N,4)
765	IF( M .NE. 0 ) THEN
766	DO i = 1,M
767	Y(i) = Y(i) + Alpha*X(i)
768	END DO
769	IF( N .LT. 4 ) RETURN
770	END IF
771	MP1 = M + 1
772	DO i = MP1,N,4
773	Y(i) = Y(i) + Alpha*X(i)
774	Y(i + 1) = Y(i + 1) + Alpha*X(i + 1)
775	Y(i + 2) = Y(i + 2) + Alpha*X(i + 2)
776	Y(i + 3) = Y(i + 3) + Alpha*X(i + 3)
777	END DO
778
779	END SUBROUTINE WAXPY
780
781
782
783	!--------------------------------------------------------------
784	SUBROUTINE WSCAL(N,Alpha,X,incX)
785	!--------------------------------------------------------------
786	! constant times a vector: x(1:N) <- Alpha*x(1:N)
787	! only for incX=incY=1
788	! after BLAS
789	! replace this by the function from the optimized BLAS implementation:
790	! CALL SSCAL(N,Alpha,X,1) or CALL DSCAL(N,Alpha,X,1)
791	!--------------------------------------------------------------
792
793	INTEGER :: i,incX,M,MP1,N
794	REAL(kind=dp) :: X(N),Alpha
795	REAL(kind=dp), PARAMETER :: ZERO=0.0_dp, ONE=1.0_dp
796
797	IF (Alpha .EQ. ONE) RETURN
798	IF (N .LE. 0) RETURN
799
800	M = MOD(N,5)
801	IF( M .NE. 0 ) THEN
802	IF (Alpha .EQ. (-ONE)) THEN
803	DO i = 1,M
804	X(i) = -X(i)
805	END DO
806	ELSEIF (Alpha .EQ. ZERO) THEN
807	DO i = 1,M
808	X(i) = ZERO
809	END DO
810	ELSE
811	DO i = 1,M
812	X(i) = Alpha*X(i)
813	END DO
814	END IF
815	IF( N .LT. 5 ) RETURN
816	END IF
817	MP1 = M + 1
818	IF (Alpha .EQ. (-ONE)) THEN
819	DO i = MP1,N,5
820	X(i) = -X(i)
821	X(i + 1) = -X(i + 1)
822	X(i + 2) = -X(i + 2)
823	X(i + 3) = -X(i + 3)
824	X(i + 4) = -X(i + 4)
825	END DO
826	ELSEIF (Alpha .EQ. ZERO) THEN
827	DO i = MP1,N,5
828	X(i) = ZERO
829	X(i + 1) = ZERO
830	X(i + 2) = ZERO
831	X(i + 3) = ZERO
832	X(i + 4) = ZERO
833	END DO
834	ELSE
835	DO i = MP1,N,5
836	X(i) = Alpha*X(i)
837	X(i + 1) = Alpha*X(i + 1)
838	X(i + 2) = Alpha*X(i + 2)
839	X(i + 3) = Alpha*X(i + 3)
840	X(i + 4) = Alpha*X(i + 4)
841	END DO
842	END IF
843
844	END SUBROUTINE WSCAL
845
846	!--------------------------------------------------------------
847	REAL(kind=dp) FUNCTION WLAMCH( C )
848	!--------------------------------------------------------------
849	! returns epsilon machine
850	! after LAPACK
851	! replace this by the function from the optimized LAPACK implementation:
852	! CALL SLAMCH('E') or CALL DLAMCH('E')
853	!--------------------------------------------------------------
854	! USE chem_gasphase_mod_Precision
855
856	CHARACTER :: C
857	INTEGER :: i
858	REAL(kind=dp), SAVE :: Eps
859	REAL(kind=dp) :: Suma
860	REAL(kind=dp), PARAMETER :: ONE=1.0_dp, HALF=0.5_dp
861	LOGICAL, SAVE :: First=.TRUE.
862
863	IF (First) THEN
864	First = .FALSE.
865	Eps = HALF**(16)
866	DO i = 17, 80
867	Eps = Eps*HALF
868	CALL WLAMCH_ADD(ONE,Eps,Suma)
869	IF (Suma.LE.ONE) GOTO 10
870	END DO
871	PRINT*,'ERROR IN WLAMCH. EPS < ',Eps
872	RETURN
873	10 Eps = Eps*2
874	i = i-1
875	END IF
876
877	WLAMCH = Eps
878
879	END FUNCTION WLAMCH
880
881	SUBROUTINE WLAMCH_ADD( A, B, Suma )
882	! USE chem_gasphase_mod_Precision
883
884	REAL(kind=dp) A, B, Suma
885	Suma = A + B
886
887	END SUBROUTINE WLAMCH_ADD
888	!--------------------------------------------------------------
889
890
891	!--------------------------------------------------------------
892	SUBROUTINE SET2ZERO(N,Y)
893	!--------------------------------------------------------------
894	! copies zeros into the vector y: y <- 0
895	! after BLAS
896	!--------------------------------------------------------------
897
898	INTEGER :: i,M,MP1,N
899	REAL(kind=dp) :: Y(N)
900	REAL(kind=dp), PARAMETER :: ZERO = 0.0d0
901
902	IF (N.LE.0) RETURN
903
904	M = MOD(N,8)
905	IF( M .NE. 0 ) THEN
906	DO i = 1,M
907	Y(i) = ZERO
908	END DO
909	IF( N .LT. 8 ) RETURN
910	END IF
911	MP1 = M+1
912	DO i = MP1,N,8
913	Y(i) = ZERO
914	Y(i + 1) = ZERO
915	Y(i + 2) = ZERO
916	Y(i + 3) = ZERO
917	Y(i + 4) = ZERO
918	Y(i + 5) = ZERO
919	Y(i + 6) = ZERO
920	Y(i + 7) = ZERO
921	END DO
922
923	END SUBROUTINE SET2ZERO
924
925
926	!--------------------------------------------------------------
927	REAL(kind=dp) FUNCTION WDOT (N, DX, incX, DY, incY)
928	!--------------------------------------------------------------
929	! dot produce: wdot = x(1:N)*y(1:N)
930	! only for incX=incY=1
931	! after BLAS
932	! replace this by the function from the optimized BLAS implementation:
933	! CALL SDOT(N,X,1,Y,1) or CALL DDOT(N,X,1,Y,1)
934	!--------------------------------------------------------------
935	! USE messy_mecca_kpp_Precision
936	!--------------------------------------------------------------
937	IMPLICIT NONE
938	INTEGER :: N, incX, incY
939	REAL(kind=dp) :: DX(N), DY(N)
940
941	INTEGER :: i, IX, IY, M, MP1, NS
942
943	WDOT = 0.0D0
944	IF (N .LE. 0) RETURN
945	IF (incX .EQ. incY) IF (incX-1) 5,20,60
946	!
947	! Code for unequal or nonpositive increments.
948	!
949	5 IX = 1
950	IY = 1
951	IF (incX .LT. 0) IX = (-N+1)*incX + 1
952	IF (incY .LT. 0) IY = (-N+1)*incY + 1
953	DO i = 1,N
954	WDOT = WDOT + DX(IX)*DY(IY)
955	IX = IX + incX
956	IY = IY + incY
957	END DO
958	RETURN
959	!
960	! Code for both increments equal to 1.
961	!
962	! Clean-up loop so remaining vector length is a multiple of 5.
963	!
964	20 M = MOD(N,5)
965	IF (M .EQ. 0) GO TO 40
966	DO i = 1,M
967	WDOT = WDOT + DX(i)*DY(i)
968	END DO
969	IF (N .LT. 5) RETURN
970	40 MP1 = M + 1
971	DO i = MP1,N,5
972	WDOT = WDOT + DX(i)DY(i) + DX(i+1)DY(i+1) + DX(i+2)*DY(i+2) + &
973	DX(i+3)DY(i+3) + DX(i+4)DY(i+4)
974	END DO
975	RETURN
976	!
977	! Code for equal, positive, non-unit increments.
978	!
979	60 NS = N*incX
980	DO i = 1,NS,incX
981	WDOT = WDOT + DX(i)*DY(i)
982	END DO
983
984	END FUNCTION WDOT
985
986
987	!--------------------------------------------------------------
988	SUBROUTINE WADD(N,X,Y,Z)
989	!--------------------------------------------------------------
990	! adds two vectors: z <- x + y
991	! BLAS - like
992	!--------------------------------------------------------------
993	! USE chem_gasphase_mod_Precision
994
995	INTEGER :: i, M, MP1, N
996	REAL(kind=dp) :: X(N),Y(N),Z(N)
997
998	IF (N.LE.0) RETURN
999
1000	M = MOD(N,5)
1001	IF( M /= 0 ) THEN
1002	DO i = 1,M
1003	Z(i) = X(i) + Y(i)
1004	END DO
1005	IF( N < 5 ) RETURN
1006	END IF
1007	MP1 = M+1
1008	DO i = MP1,N,5
1009	Z(i) = X(i) + Y(i)
1010	Z(i + 1) = X(i + 1) + Y(i + 1)
1011	Z(i + 2) = X(i + 2) + Y(i + 2)
1012	Z(i + 3) = X(i + 3) + Y(i + 3)
1013	Z(i + 4) = X(i + 4) + Y(i + 4)
1014	END DO
1015
1016	END SUBROUTINE WADD
1017
1018
1019
1020	!--------------------------------------------------------------
1021	SUBROUTINE WGEFA(N,A,Ipvt,info)
1022	!--------------------------------------------------------------
1023	! WGEFA FACTORS THE MATRIX A (N,N) BY
1024	! GAUSS ELIMINATION WITH PARTIAL PIVOTING
1025	! LINPACK - LIKE
1026	!--------------------------------------------------------------
1027	!
1028	INTEGER :: N,Ipvt(N),info
1029	REAL(kind=dp) :: A(N,N)
1030	REAL(kind=dp) :: t, dmax, da
1031	INTEGER :: j,k,l
1032	REAL(kind=dp), PARAMETER :: ZERO = 0.0, ONE = 1.0
1033
1034	info = 0
1035
1036	size: IF (n > 1) THEN
1037
1038	col: DO k = 1, n-1
1039
1040	! find l = pivot index
1041	! l = idamax(n-k+1,A(k,k),1) + k - 1
1042	l = k; dmax = abs(A(k,k))
1043	DO j = k+1,n
1044	da = ABS(A(j,k))
1045	IF (da > dmax) THEN
1046	l = j; dmax = da
1047	END IF
1048	END DO
1049	Ipvt(k) = l
1050
1051	! zero pivot implies this column already triangularized
1052	IF (ABS(A(l,k)) < TINY(ZERO)) THEN
1053	info = k
1054	return
1055	ELSE
1056	IF (l /= k) THEN
1057	t = A(l,k); A(l,k) = A(k,k); A(k,k) = t
1058	END IF
1059	t = -ONE/A(k,k)
1060	CALL WSCAL(n-k,t,A(k+1,k),1)
1061	DO j = k+1, n
1062	t = A(l,j)
1063	IF (l /= k) THEN
1064	A(l,j) = A(k,j); A(k,j) = t
1065	END IF
1066	CALL WAXPY(n-k,t,A(k+1,k),1,A(k+1,j),1)
1067	END DO
1068	END IF
1069
1070	END DO col
1071
1072	END IF size
1073
1074	Ipvt(N) = N
1075	IF (ABS(A(N,N)) == ZERO) info = N
1076
1077	END SUBROUTINE WGEFA
1078
1079
1080	!--------------------------------------------------------------
1081	SUBROUTINE WGESL(Trans,N,A,Ipvt,b)
1082	!--------------------------------------------------------------
1083	! WGESL solves the system
1084	! a * x = b or trans(a) * x = b
1085	! using the factors computed by WGEFA.
1086	!
1087	! Trans = 'N' to solve A*x = b ,
1088	! = 'T' to solve transpose(A)*x = b
1089	! LINPACK - LIKE
1090	!--------------------------------------------------------------
1091
1092	INTEGER :: N,Ipvt(N)
1093	CHARACTER :: trans
1094	REAL(kind=dp) :: A(N,N),b(N)
1095	REAL(kind=dp) :: t
1096	INTEGER :: k,kb,l
1097
1098
1099	SELECT CASE (Trans)
1100
1101	CASE ('n','N') ! Solve A * x = b
1102
1103	! first solve L*y = b
1104	IF (n >= 2) THEN
1105	DO k = 1, n-1
1106	l = Ipvt(k)
1107	t = b(l)
1108	IF (l /= k) THEN
1109	b(l) = b(k)
1110	b(k) = t
1111	END IF
1112	CALL WAXPY(n-k,t,a(k+1,k),1,b(k+1),1)
1113	END DO
1114	END IF
1115	! now solve U*x = y
1116	DO kb = 1, n
1117	k = n + 1 - kb
1118	b(k) = b(k)/a(k,k)
1119	t = -b(k)
1120	CALL WAXPY(k-1,t,a(1,k),1,b(1),1)
1121	END DO
1122
1123	CASE ('t','T') ! Solve transpose(A) * x = b
1124
1125	! first solve trans(U)*y = b
1126	DO k = 1, n
1127	t = WDOT(k-1,a(1,k),1,b(1),1)
1128	b(k) = (b(k) - t)/a(k,k)
1129	END DO
1130	! now solve trans(L)*x = y
1131	IF (n >= 2) THEN
1132	DO kb = 1, n-1
1133	k = n - kb
1134	b(k) = b(k) + WDOT(n-k,a(k+1,k),1,b(k+1),1)
1135	l = Ipvt(k)
1136	IF (l /= k) THEN
1137	t = b(l); b(l) = b(k); b(k) = t
1138	END IF
1139	END DO
1140	END IF
1141
1142	END SELECT
1143
1144	END SUBROUTINE WGESL
1145	! End of BLAS_UTIL function
1146	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1147
1148
1149
1150	END MODULE chem_gasphase_mod_LinearAlgebra
1151

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |