Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/util/blas.c @ 4002

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

Merge of branch palm4u into trunk
File size: 6.9 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 
9--------------------------------------------------------------*/
10
11#define ZERO  (KPP_REAL)0.0
12#define ONE   (KPP_REAL)1.0
13#define HALF  (KPP_REAL)0.5
14#define TWO   (KPP_REAL)2.0
15#define MOD(A,B) (int)((A)%(B))
16
17/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
18void WCOPY(int N, KPP_REAL X[], int incX, KPP_REAL Y[], int incY)
19/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
20    copies a vector, x, to a vector, y:  y <- x
21    only for incX=incY=1
22    after BLAS
23    replace this by the function from the optimized BLAS implementation:
24        CALL  SCOPY(N,X,1,Y,1)   or   CALL  DCOPY(N,X,1,Y,1)
25~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
26{
27      int i, M;
28      if (N <= 0) return;
29
30      M = MOD(N,8);
31      if( M != 0 ) {
32        for ( i = 0; i < M; i++ )
33          Y[i] = X[i];
34        if( N < 8 ) return;
35      } /* end if */   
36      for ( i = M; i<N; i+=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 for */
46
47} /* end function WCOPY */
48
49
50/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
51void WAXPY(int N, KPP_REAL Alpha, 
52         KPP_REAL X[], int incX, 
53         KPP_REAL Y[], int incY )
54/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
55    constant times a vector plus a vector: y <- y + Alpha*x
56    only for incX=incY=1
57    after BLAS
58    replace this by the function from the optimized BLAS implementation:
59        CALL SAXPY(N,Alpha,X,1,Y,1) or  CALL DAXPY(N,Alpha,X,1,Y,1)
60~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
61{
62      int i, M;
63
64      if (Alpha == ZERO) return;
65      if (N  <=  0) return;
66
67      M = MOD(N,4);
68      if( M != 0 ) {
69        for ( i = 0; i < M; i++ )
70          Y[i] = Y[i] + Alpha*X[i];
71        if ( N < 4 ) return;
72      } /* end if */
73     
74      for ( i = M; i < N; i += 4 ) {
75        Y[i] = Y[i] + Alpha*X[i];
76        Y[i + 1] = Y[i + 1] + Alpha*X[i + 1];
77        Y[i + 2] = Y[i + 2] + Alpha*X[i + 2];
78        Y[i + 3] = Y[i + 3] + Alpha*X[i + 3];
79      } /* end for */
80
81} /* end function  WAXPY */
82
83
84
85/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
86void WSCAL(int N, KPP_REAL Alpha, KPP_REAL X[], int incX)
87/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
88    constant times a vector: x(1:N) <- Alpha*x(1:N)
89    only for incX=incY=1
90    after BLAS
91    replace this by the function from the optimized BLAS implementation:
92        CALL SSCAL(N,Alpha,X,1) or  CALL DSCAL(N,Alpha,X,1)
93~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
94{
95      int i, M;
96     
97      if (Alpha == ONE) return;
98      if (N  <=  0) return;
99
100      M = MOD(N,5);
101      if( M  !=  0 ) {
102        if (Alpha == (-ONE))
103          for ( i = 0; i < M; i++ )  X[i] = -X[i];
104        else {
105          if (Alpha == ZERO)
106            for ( i = 0; i < M; i++ ) X[i] = ZERO;
107          else
108            for ( i = 0; i < M; i++ ) X[i] = Alpha*X[i];
109        } /* end else */
110        if( N < 5 ) return;
111      } /* end if */
112     
113      if (Alpha == (-ONE))
114        for ( i = M; i<N; i+=5 ) {
115          X[i]     = -X[i];
116          X[i + 1] = -X[i + 1];
117          X[i + 2] = -X[i + 2];
118          X[i + 3] = -X[i + 3];
119          X[i + 4] = -X[i + 4];
120        } /* end for */
121      else {
122        if (Alpha == ZERO)
123          for ( i = M; i < N; i += 5 ) {
124            X[i]     = ZERO;
125            X[i + 1] = ZERO;
126            X[i + 2] = ZERO;
127            X[i + 3] = ZERO;
128            X[i + 4] = ZERO;
129          } /* end for */
130        else
131          for ( i = M; i < N; i += 5 ) {
132            X[i]     = Alpha*X[i];
133            X[i + 1] = Alpha*X[i + 1];
134            X[i + 2] = Alpha*X[i + 2];
135            X[i + 3] = Alpha*X[i + 3];
136            X[i + 4] = Alpha*X[i + 4];
137           } /* end for */
138      }  /* else  */
139     
140} /* end function WSCAL */
141   
142     
143/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
144KPP_REAL WLAMCH_ADD( KPP_REAL  A, KPP_REAL  B )
145{
146      return (A + B);
147} /* end function  WLAMCH_ADD */
148
149/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
150KPP_REAL WLAMCH( char 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{
158      int i;
159      KPP_REAL Suma;
160      static KPP_REAL Eps;
161      static char First = 1;
162     
163      if (First) {
164        First = 0;
165        Eps = pow(HALF,16);
166        for ( i = 17; i <= 80; i++ ) {
167          Eps = Eps*HALF;
168          Suma = WLAMCH_ADD(ONE,Eps);
169          if (Suma <= ONE) break;
170        } /* end for */
171        if (i==80) {
172           printf("\nERROR IN WLAMCH. Very small EPS = %g\n",Eps);
173           return (double)2.2e-16;
174        }
175        Eps *= TWO; i--;     
176      } /* end if First */
177
178      return Eps;
179
180} /* end function WLAMCH */
181
182
183/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
184        copies zeros into the vector y:  y <- 0 after BLAS
185  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
186void Set2Zero(int N, KPP_REAL Y[])
187{
188   int i,M;
189       
190   if (N <= 0) return;
191
192   M = MOD(N,8);
193   if  (M != 0)
194   {
195        for (i = 0; i < M; i++)
196        {
197           Y[i] = ZERO;
198        }
199        if (N < 8) 
200           return;
201   } /* end if */
202   for (i = M; i < N; i += 8)
203   {
204        Y[i] = ZERO;
205        Y[i + 1] = ZERO;
206        Y[i + 2] = ZERO;
207        Y[i + 3] = ZERO;
208        Y[i + 4] = ZERO;
209        Y[i + 5] = ZERO;
210        Y[i + 6] = ZERO;
211        Y[i + 7] = ZERO;
212   } /* end for */
213} /* end function Set2Zero */
214
215
216/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
217        adds two vectors: z <- x + y     BLAS - like
218  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
219void WADD(int N, KPP_REAL X[], KPP_REAL Y[], KPP_REAL Z[])
220{
221        int i, M;
222
223        if ( N <= 0 ) return;
224       
225        M = MOD(N,5);
226        if ( M != 0 )
227        {
228           for(i = 0; i < M; i++)
229           {
230                Z[i] = X[i] + Y[i];
231           }
232           if ( N < 5 ) return;
233        } /* end if */
234        for (i = M; i < N; i += 5)
235        {
236           Z[i] = X[i] + Y[i];
237           Z[i + 1] = X[i + 1] + Y[i + 1];
238           Z[i + 2] = X[i + 2] + Y[i + 2];
239           Z[i + 3] = X[i + 3] + Y[i + 3];
240           Z[i + 4] = X[i + 4] + Y[i + 4];
241        } /* end for */
242} /* end function WADD */
Note: See TracBrowser for help on using the repository browser.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |