source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/util/blas.f90 @ 3655

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

Merge of branch palm4u into trunk

File size: 12.6 KB
Line 
1!--------------------------------------------------------------
2!
3! BLAS/LAPACK-like subroutines used by the integration algorithms
4! It is recommended to replace them by calls to the optimized
5!      BLAS/LAPACK library for your machine
6!
7!  (C) Adrian Sandu, Aug. 2004
8!      Virginia Polytechnic Institute and State University
9!--------------------------------------------------------------
10
11
12!--------------------------------------------------------------
13      SUBROUTINE WCOPY(N,X,incX,Y,incY)
14!--------------------------------------------------------------
15!     copies a vector, x, to a vector, y:  y <- x
16!     only for incX=incY=1
17!     after BLAS
18!     replace this by the function from the optimized BLAS implementation:
19!         CALL  SCOPY(N,X,1,Y,1)   or   CALL  DCOPY(N,X,1,Y,1)
20!--------------------------------------------------------------
21!     USE KPP_ROOT_Precision
22     
23      INTEGER  :: i,incX,incY,M,MP1,N
24      KPP_REAL :: X(N),Y(N)
25
26      IF (N.LE.0) RETURN
27
28      M = MOD(N,8)
29      IF( M .NE. 0 ) THEN
30        DO i = 1,M
31          Y(i) = X(i)
32        END DO
33        IF( N .LT. 8 ) RETURN
34      END IF   
35      MP1 = M+1
36      DO i = MP1,N,8
37        Y(i) = X(i)
38        Y(i + 1) = X(i + 1)
39        Y(i + 2) = X(i + 2)
40        Y(i + 3) = X(i + 3)
41        Y(i + 4) = X(i + 4)
42        Y(i + 5) = X(i + 5)
43        Y(i + 6) = X(i + 6)
44        Y(i + 7) = X(i + 7)
45      END DO
46
47      END SUBROUTINE WCOPY
48
49
50!--------------------------------------------------------------
51      SUBROUTINE WAXPY(N,Alpha,X,incX,Y,incY)
52!--------------------------------------------------------------
53!     constant times a vector plus a vector: y <- y + Alpha*x
54!     only for incX=incY=1
55!     after BLAS
56!     replace this by the function from the optimized BLAS implementation:
57!         CALL SAXPY(N,Alpha,X,1,Y,1) or  CALL DAXPY(N,Alpha,X,1,Y,1)
58!--------------------------------------------------------------
59
60      INTEGER  :: i,incX,incY,M,MP1,N
61      KPP_REAL :: X(N),Y(N),Alpha
62      KPP_REAL, PARAMETER :: ZERO = 0.0_dp
63
64      IF (Alpha .EQ. ZERO) RETURN
65      IF (N .LE. 0) RETURN
66
67      M = MOD(N,4)
68      IF( M .NE. 0 ) THEN
69        DO i = 1,M
70          Y(i) = Y(i) + Alpha*X(i)
71        END DO
72        IF( N .LT. 4 ) RETURN
73      END IF
74      MP1 = M + 1
75      DO i = MP1,N,4
76        Y(i) = Y(i) + Alpha*X(i)
77        Y(i + 1) = Y(i + 1) + Alpha*X(i + 1)
78        Y(i + 2) = Y(i + 2) + Alpha*X(i + 2)
79        Y(i + 3) = Y(i + 3) + Alpha*X(i + 3)
80      END DO
81     
82      END SUBROUTINE WAXPY
83
84
85
86!--------------------------------------------------------------
87      SUBROUTINE WSCAL(N,Alpha,X,incX)
88!--------------------------------------------------------------
89!     constant times a vector: x(1:N) <- Alpha*x(1:N)
90!     only for incX=incY=1
91!     after BLAS
92!     replace this by the function from the optimized BLAS implementation:
93!         CALL SSCAL(N,Alpha,X,1) or  CALL DSCAL(N,Alpha,X,1)
94!--------------------------------------------------------------
95
96      INTEGER  :: i,incX,M,MP1,N
97      KPP_REAL  :: X(N),Alpha
98      KPP_REAL, PARAMETER  :: ZERO=0.0_dp, ONE=1.0_dp
99
100      IF (Alpha .EQ. ONE) RETURN
101      IF (N .LE. 0) RETURN
102
103      M = MOD(N,5)
104      IF( M .NE. 0 ) THEN
105        IF (Alpha .EQ. (-ONE)) THEN
106          DO i = 1,M
107            X(i) = -X(i)
108          END DO
109        ELSEIF (Alpha .EQ. ZERO) THEN
110          DO i = 1,M
111            X(i) = ZERO
112          END DO
113        ELSE
114          DO i = 1,M
115            X(i) = Alpha*X(i)
116          END DO
117        END IF
118        IF( N .LT. 5 ) RETURN
119      END IF
120      MP1 = M + 1
121      IF (Alpha .EQ. (-ONE)) THEN
122        DO i = MP1,N,5
123          X(i)     = -X(i)
124          X(i + 1) = -X(i + 1)
125          X(i + 2) = -X(i + 2)
126          X(i + 3) = -X(i + 3)
127          X(i + 4) = -X(i + 4)
128        END DO
129      ELSEIF (Alpha .EQ. ZERO) THEN
130        DO i = MP1,N,5
131          X(i)     = ZERO
132          X(i + 1) = ZERO
133          X(i + 2) = ZERO
134          X(i + 3) = ZERO
135          X(i + 4) = ZERO
136        END DO
137      ELSE
138        DO i = MP1,N,5
139          X(i)     = Alpha*X(i)
140          X(i + 1) = Alpha*X(i + 1)
141          X(i + 2) = Alpha*X(i + 2)
142          X(i + 3) = Alpha*X(i + 3)
143          X(i + 4) = Alpha*X(i + 4)
144        END DO
145      END IF
146
147      END SUBROUTINE WSCAL
148
149!--------------------------------------------------------------
150      KPP_REAL FUNCTION WLAMCH( C )
151!--------------------------------------------------------------
152!     returns epsilon machine
153!     after LAPACK
154!     replace this by the function from the optimized LAPACK implementation:
155!          CALL SLAMCH('E') or CALL DLAMCH('E')
156!--------------------------------------------------------------
157!      USE KPP_ROOT_Precision
158
159      CHARACTER ::  C
160      INTEGER    :: i
161      KPP_REAL, SAVE  ::  Eps
162      KPP_REAL  ::  Suma
163      KPP_REAL, PARAMETER  ::  ONE=1.0_dp, HALF=0.5_dp
164      LOGICAL, SAVE   ::  First=.TRUE.
165     
166      IF (First) THEN
167        First = .FALSE.
168        Eps = HALF**(16)
169        DO i = 17, 80
170          Eps = Eps*HALF
171          CALL WLAMCH_ADD(ONE,Eps,Suma)
172          IF (Suma.LE.ONE) GOTO 10
173        END DO
174        PRINT*,'ERROR IN WLAMCH. EPS < ',Eps
175        RETURN
17610      Eps = Eps*2
177        i = i-1     
178      END IF
179
180      WLAMCH = Eps
181
182      END FUNCTION WLAMCH
183     
184      SUBROUTINE WLAMCH_ADD( A, B, Suma )
185!      USE KPP_ROOT_Precision
186     
187      KPP_REAL A, B, Suma
188      Suma = A + B
189
190      END SUBROUTINE WLAMCH_ADD
191!--------------------------------------------------------------
192
193
194!--------------------------------------------------------------
195      SUBROUTINE SET2ZERO(N,Y)
196!--------------------------------------------------------------
197!     copies zeros into the vector y:  y <- 0
198!     after BLAS
199!--------------------------------------------------------------
200     
201      INTEGER ::  i,M,MP1,N
202      KPP_REAL ::  Y(N)
203      KPP_REAL, PARAMETER :: ZERO = 0.0d0
204
205      IF (N.LE.0) RETURN
206
207      M = MOD(N,8)
208      IF( M .NE. 0 ) THEN
209        DO i = 1,M
210          Y(i) = ZERO
211        END DO
212        IF( N .LT. 8 ) RETURN
213      END IF   
214      MP1 = M+1
215      DO i = MP1,N,8
216        Y(i)     = ZERO
217        Y(i + 1) = ZERO
218        Y(i + 2) = ZERO
219        Y(i + 3) = ZERO
220        Y(i + 4) = ZERO
221        Y(i + 5) = ZERO
222        Y(i + 6) = ZERO
223        Y(i + 7) = ZERO
224      END DO
225
226      END SUBROUTINE SET2ZERO
227
228
229!--------------------------------------------------------------
230      KPP_REAL FUNCTION WDOT (N, DX, incX, DY, incY) 
231!--------------------------------------------------------------
232!     dot produce: wdot = x(1:N)*y(1:N)
233!     only for incX=incY=1
234!     after BLAS
235!     replace this by the function from the optimized BLAS implementation:
236!         CALL SDOT(N,X,1,Y,1) or  CALL DDOT(N,X,1,Y,1)
237!--------------------------------------------------------------
238!      USE messy_mecca_kpp_Precision
239!--------------------------------------------------------------
240      IMPLICIT NONE
241      INTEGER :: N, incX, incY
242      KPP_REAL :: DX(N), DY(N) 
243
244      INTEGER :: i, IX, IY, M, MP1, NS
245                                 
246      WDOT = 0.0D0 
247      IF (N .LE. 0) RETURN
248      IF (incX .EQ. incY) IF (incX-1) 5,20,60 
249!                                                                       
250!     Code for unequal or nonpositive increments.                       
251!                                                                       
252    5 IX = 1 
253      IY = 1 
254      IF (incX .LT. 0) IX = (-N+1)*incX + 1 
255      IF (incY .LT. 0) IY = (-N+1)*incY + 1 
256      DO i = 1,N 
257        WDOT = WDOT + DX(IX)*DY(IY) 
258        IX = IX + incX 
259        IY = IY + incY 
260      END DO
261      RETURN 
262!                                                                       
263!     Code for both increments equal to 1.                             
264!                                                                       
265!     Clean-up loop so remaining vector length is a multiple of 5.     
266!                                                                       
267   20 M = MOD(N,5) 
268      IF (M .EQ. 0) GO TO 40 
269      DO i = 1,M 
270         WDOT = WDOT + DX(i)*DY(i) 
271      END DO
272      IF (N .LT. 5) RETURN
273   40 MP1 = M + 1 
274      DO i = MP1,N,5 
275          WDOT = WDOT + DX(i)*DY(i) + DX(i+1)*DY(i+1) + DX(i+2)*DY(i+2) +  &
276                   DX(i+3)*DY(i+3) + DX(i+4)*DY(i+4)                   
277      END DO
278      RETURN 
279!                                                                       
280!     Code for equal, positive, non-unit increments.                   
281!                                                                       
282   60 NS = N*incX 
283      DO i = 1,NS,incX 
284        WDOT = WDOT + DX(i)*DY(i) 
285      END DO
286
287      END FUNCTION WDOT                                         
288
289
290!--------------------------------------------------------------
291      SUBROUTINE WADD(N,X,Y,Z)
292!--------------------------------------------------------------
293!     adds two vectors: z <- x + y
294!     BLAS - like
295!--------------------------------------------------------------
296!     USE KPP_ROOT_Precision
297     
298      INTEGER :: i, M, MP1, N
299      KPP_REAL :: X(N),Y(N),Z(N)
300
301      IF (N.LE.0) RETURN
302
303      M = MOD(N,5)
304      IF( M /= 0 ) THEN
305         DO i = 1,M
306            Z(i) = X(i) + Y(i)
307         END DO
308         IF( N < 5 ) RETURN
309      END IF   
310      MP1 = M+1
311      DO i = MP1,N,5
312         Z(i)     = X(i)     + Y(i)
313         Z(i + 1) = X(i + 1) + Y(i + 1)
314         Z(i + 2) = X(i + 2) + Y(i + 2)
315         Z(i + 3) = X(i + 3) + Y(i + 3)
316         Z(i + 4) = X(i + 4) + Y(i + 4)
317      END DO
318
319      END SUBROUTINE WADD
320     
321     
322     
323!--------------------------------------------------------------
324      SUBROUTINE WGEFA(N,A,Ipvt,info)
325!--------------------------------------------------------------
326!     WGEFA FACTORS THE MATRIX A (N,N) BY
327!           GAUSS ELIMINATION WITH PARTIAL PIVOTING
328!     LINPACK - LIKE
329!--------------------------------------------------------------
330!
331      INTEGER       :: N,Ipvt(N),info
332      KPP_REAL :: A(N,N)
333      KPP_REAL :: t, dmax, da
334      INTEGER       :: j,k,l
335      KPP_REAL, PARAMETER :: ZERO = 0.0, ONE = 1.0
336
337      info = 0
338
339size: IF (n > 1) THEN
340     
341col:  DO k = 1, n-1
342
343!        find l = pivot index
344!        l = idamax(n-k+1,A(k,k),1) + k - 1
345         l = k; dmax = abs(A(k,k))
346         DO j = k+1,n
347            da = ABS(A(j,k))
348            IF (da > dmax) THEN
349              l = j; dmax = da
350            END IF
351         END DO
352         Ipvt(k) = l
353
354!        zero pivot implies this column already triangularized
355         IF (ABS(A(l,k)) < TINY(ZERO)) THEN
356            info = k
357            return
358         ELSE   
359            IF (l /= k) THEN
360               t = A(l,k); A(l,k) = A(k,k); A(k,k) = t
361            END IF
362            t = -ONE/A(k,k)
363            CALL WSCAL(n-k,t,A(k+1,k),1)
364            DO j = k+1, n
365               t = A(l,j)
366               IF (l /= k) THEN
367                  A(l,j) = A(k,j); A(k,j) = t
368               END IF
369               CALL WAXPY(n-k,t,A(k+1,k),1,A(k+1,j),1)
370            END DO         
371         END IF
372         
373       END DO col
374       
375      END IF size
376     
377      Ipvt(N) = N
378      IF (ABS(A(N,N)) == ZERO) info = N
379     
380      END SUBROUTINE WGEFA
381
382
383!--------------------------------------------------------------
384      SUBROUTINE WGESL(Trans,N,A,Ipvt,b)
385!--------------------------------------------------------------
386!     WGESL solves the system
387!     a * x = b  or  trans(a) * x = b
388!     using the factors computed by WGEFA.
389!
390!     Trans      = 'N'   to solve  A*x = b ,
391!                = 'T'   to solve  transpose(A)*x = b
392!     LINPACK - LIKE
393!--------------------------------------------------------------
394
395      INTEGER       :: N,Ipvt(N)
396      CHARACTER     :: trans
397      KPP_REAL :: A(N,N),b(N)
398      KPP_REAL :: t
399      INTEGER       :: k,kb,l
400
401     
402      SELECT CASE (Trans)
403
404      CASE ('n','N')  !  Solve  A * x = b
405
406!        first solve  L*y = b
407         IF (n >= 2) THEN
408          DO k = 1, n-1
409            l = Ipvt(k)
410            t = b(l)
411            IF (l /= k) THEN
412               b(l) = b(k)
413               b(k) = t
414            END IF
415            CALL WAXPY(n-k,t,a(k+1,k),1,b(k+1),1)
416          END DO
417         END IF
418!        now solve  U*x = y
419         DO kb = 1, n
420            k = n + 1 - kb
421            b(k) = b(k)/a(k,k)
422            t = -b(k)
423            CALL WAXPY(k-1,t,a(1,k),1,b(1),1)
424         END DO
425     
426      CASE ('t','T')  !  Solve transpose(A) * x = b
427
428!        first solve  trans(U)*y = b
429         DO k = 1, n
430            t = WDOT(k-1,a(1,k),1,b(1),1)
431            b(k) = (b(k) - t)/a(k,k)
432         END DO
433!        now solve trans(L)*x = y
434         IF (n >= 2) THEN
435         DO kb = 1, n-1
436            k = n - kb
437            b(k) = b(k) + WDOT(n-k,a(k+1,k),1,b(k+1),1)
438            l = Ipvt(k)
439            IF (l /= k) THEN
440               t = b(l); b(l) = b(k); b(k) = t
441            END IF
442         END DO
443         END IF
444   
445      END SELECT
446
447      END SUBROUTINE WGESL
Note: See TracBrowser for help on using the repository browser.