source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/util/blas.f @ 4127

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

Merge of branch palm4u into trunk

File size: 9.2 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     
22      INTEGER i,incX,incY,M,MP1,N
23      KPP_REAL X(N),Y(N)
24
25      IF (N.LE.0) RETURN
26
27      M = MOD(N,8)
28      IF( M .NE. 0 ) THEN
29        DO i = 1,M
30          Y(i) = X(i)
31        END DO
32        IF( N .LT. 8 ) RETURN
33      END IF   
34      MP1 = M+1
35      DO i = MP1,N,8
36        Y(i) = X(i)
37        Y(i + 1) = X(i + 1)
38        Y(i + 2) = X(i + 2)
39        Y(i + 3) = X(i + 3)
40        Y(i + 4) = X(i + 4)
41        Y(i + 5) = X(i + 5)
42        Y(i + 6) = X(i + 6)
43        Y(i + 7) = X(i + 7)
44      END DO
45      RETURN
46      END
47
48
49!--------------------------------------------------------------
50      SUBROUTINE WAXPY(N,Alpha,X,incX,Y,incY)
51!--------------------------------------------------------------
52!     constant times a vector plus a vector: y <- y + Alpha*x
53!     only for incX=incY=1
54!     after BLAS
55!     replace this by the function from the optimized BLAS implementation:
56!         CALL SAXPY(N,Alpha,X,1,Y,1) or  CALL DAXPY(N,Alpha,X,1,Y,1)
57!--------------------------------------------------------------
58
59      INTEGER i,incX,incY,M,MP1,N
60      KPP_REAL X(N),Y(N),Alpha
61      KPP_REAL ZERO
62      PARAMETER( ZERO = 0.0d0 )
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      RETURN
82      END
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 ZERO, ONE
99      PARAMETER( ZERO = 0.0d0 )
100      PARAMETER( ONE  = 1.0d0 )
101
102      IF (Alpha .EQ. ONE) RETURN
103      IF (N .LE. 0) RETURN
104
105      M = MOD(N,5)
106      IF( M .NE. 0 ) THEN
107        IF (Alpha .EQ. (-ONE)) THEN
108          DO i = 1,M
109            X(i) = -X(i)
110          END DO
111        ELSEIF (Alpha .EQ. ZERO) THEN
112          DO i = 1,M
113            X(i) = ZERO
114          END DO
115        ELSE
116          DO i = 1,M
117            X(i) = Alpha*X(i)
118          END DO
119        END IF
120        IF( N .LT. 5 ) RETURN
121      END IF
122      MP1 = M + 1
123      IF (Alpha .EQ. (-ONE)) THEN
124        DO i = MP1,N,5
125          X(i)     = -X(i)
126          X(i + 1) = -X(i + 1)
127          X(i + 2) = -X(i + 2)
128          X(i + 3) = -X(i + 3)
129          X(i + 4) = -X(i + 4)
130        END DO
131      ELSEIF (Alpha .EQ. ZERO) THEN
132        DO i = MP1,N,5
133          X(i)     = ZERO
134          X(i + 1) = ZERO
135          X(i + 2) = ZERO
136          X(i + 3) = ZERO
137          X(i + 4) = ZERO
138        END DO
139      ELSE
140        DO i = MP1,N,5
141          X(i)     = Alpha*X(i)
142          X(i + 1) = Alpha*X(i + 1)
143          X(i + 2) = Alpha*X(i + 2)
144          X(i + 3) = Alpha*X(i + 3)
145          X(i + 4) = Alpha*X(i + 4)
146        END DO
147      END IF
148      RETURN
149      END
150
151!--------------------------------------------------------------
152      KPP_REAL FUNCTION WLAMCH( C )
153!--------------------------------------------------------------
154!     returns epsilon machine
155!     after LAPACK
156!     replace this by the function from the optimized LAPACK implementation:
157!          CALL SLAMCH('E') or CALL DLAMCH('E')
158!--------------------------------------------------------------
159
160      CHARACTER C
161      INTEGER   i
162      KPP_REAL  ONE, HALF, Eps, Suma
163      PARAMETER (ONE  = 1.0d0)
164      PARAMETER (HALF = 0.5d0)
165      LOGICAL   First
166      SAVE     First, Eps
167      DATA     First /.TRUE./
168     
169      IF (First) THEN
170        First = .FALSE.
171        Eps = HALF**(16)
172        DO i = 17, 80
173          Eps = Eps*HALF
174          CALL WLAMCH_ADD(ONE,Eps,Suma)
175          IF (Suma.LE.ONE) GOTO 10
176        END DO
177        PRINT*,'ERROR IN WLAMCH. EPS < ',Eps
178        RETURN
17910      Eps = Eps*2
180        i = i-1     
181      END IF
182
183      WLAMCH = Eps
184
185      RETURN
186      END
187     
188      SUBROUTINE WLAMCH_ADD( A, B, Suma )
189      KPP_REAL A, B, Suma
190      Suma = A + B
191      RETURN
192      END
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
204      IF (N.LE.0) RETURN
205
206      M = MOD(N,8)
207      IF( M .NE. 0 ) THEN
208        DO i = 1,M
209          Y(i) = 0.0d0
210        END DO
211        IF( N .LT. 8 ) RETURN
212      END IF   
213      MP1 = M+1
214      DO i = MP1,N,8
215        Y(i)     = 0.0d0
216        Y(i + 1) = 0.0d0
217        Y(i + 2) = 0.0d0
218        Y(i + 3) = 0.0d0
219        Y(i + 4) = 0.0d0
220        Y(i + 5) = 0.0d0
221        Y(i + 6) = 0.0d0
222        Y(i + 7) = 0.0d0
223      END DO
224
225      END SUBROUTINE SET2ZERO
226
227
228!--------------------------------------------------------------
229      KPP_REAL FUNCTION WDOT (N, DX, incX, DY, incY) 
230!--------------------------------------------------------------
231!     dot produce: wdot = x(1:N)*y(1:N)
232!     only for incX=incY=1
233!     after BLAS
234!     replace this by the function from the optimized BLAS implementation:
235!         CALL SDOT(N,X,1,Y,1) or  CALL DDOT(N,X,1,Y,1)
236!--------------------------------------------------------------
237!      USE messy_mecca_kpp_Precision
238!--------------------------------------------------------------
239      IMPLICIT NONE
240      INTEGER  N, incX, incY
241      KPP_REAL DX(N), DY(N) 
242
243      INTEGER  i, IX, IY, M, MP1, NS
244                                 
245      WDOT = 0.0D0 
246      IF (N .LE. 0) RETURN
247      IF (incX .EQ. incY) IF (incX-1) 5,20,60 
248!                                                                       
249!     Code for unequal or nonpositive increments.                       
250!                                                                       
251    5 IX = 1 
252      IY = 1 
253      IF (incX .LT. 0) IX = (-N+1)*incX + 1 
254      IF (incY .LT. 0) IY = (-N+1)*incY + 1 
255      DO i = 1,N 
256        WDOT = WDOT + DX(IX)*DY(IY) 
257        IX = IX + incX 
258        IY = IY + incY 
259      END DO
260      RETURN 
261!                                                                       
262!     Code for both increments equal to 1.                             
263!                                                                       
264!     Clean-up loop so remaining vector length is a multiple of 5.     
265!                                                                       
266   20 M = MOD(N,5) 
267      IF (M .EQ. 0) GO TO 40 
268      DO i = 1,M 
269         WDOT = WDOT + DX(i)*DY(i) 
270      END DO
271      IF (N .LT. 5) RETURN
272   40 MP1 = M + 1 
273      DO i = MP1,N,5 
274          WDOT = WDOT + DX(i)*DY(i) 
275     &            + 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     
297      INTEGER  i, M, MP1, N
298      KPP_REAL  X(N),Y(N),Z(N)
299
300      IF (N.LE.0) RETURN
301
302      M = MOD(N,5)
303      IF( M .NE. 0 ) THEN
304         DO i = 1,M
305            Z(i) = X(i) + Y(i)
306         END DO
307         IF( N .LT. 5 ) RETURN
308      END IF   
309      MP1 = M+1
310      DO i = MP1,N,5
311         Z(i)     = X(i)     + Y(i)
312         Z(i + 1) = X(i + 1) + Y(i + 1)
313         Z(i + 2) = X(i + 2) + Y(i + 2)
314         Z(i + 3) = X(i + 3) + Y(i + 3)
315         Z(i + 4) = X(i + 4) + Y(i + 4)
316      END DO
317
318      END SUBROUTINE WADD
319!--------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.