chem/KPP/kpp/kpp-2.1/util/blas.c

   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 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  18 void 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 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  51 void 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 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  86 void 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 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
 144 KPP_REAL WLAMCH_ADD( KPP_REAL  A, KPP_REAL  B )
 145 {
 146       return (A + B);
 147 } /* end function  WLAMCH_ADD */
 148
 149 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
 150 KPP_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 Sum;
 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           Sum = WLAMCH_ADD(ONE,Eps);
 169           if (Sum <= 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 */