home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS - Coast to Coast / simteldosarchivecoasttocoast2.iso / biology / gsrc208a.zip / DAXPY.C < prev    next >
Text File  |  1992-05-01  |  2KB  |  94 lines
  1. /*
  2. From tam@dragonfly.wri.com Wed Apr 24 15:48:31 1991
  3. Return-Path: <tam>
  4. Date: Wed, 24 Apr 91 17:48:43 CDT
  5. From: tam@dragonfly.wri.com
  6. To: whitbeck@sanjuan.wrc.unr.edu
  7. */
  8.  
  9. void daxpy( n, da, dx, incx, dy, incy )
  10.  
  11. double da, *dx, *dy;
  12. int n, incx, incy;
  13.  
  14. /*
  15.    Purpose : To compute
  16.  
  17.    dy = da * dx + dy
  18.  
  19.  
  20.    --- Input ---
  21.  
  22.    n    : number of elements in input vector(s)
  23.    da   : double scalar multiplier
  24.    dx   : double vector with n+1 elements, dx[0] is not used
  25.    incx : storage spacing between elements of dx
  26.    dy   : double vector with n+1 elements, dy[0] is not used
  27.    incy : storage spacing between elements of dy
  28.  
  29.  
  30.    --- Output ---
  31.  
  32.    dy = da * dx + dy, unchanged if n <= 0
  33.  
  34.  
  35.    For i = 0 to n-1, replace dy[ly+i*incy] with
  36.    da*dx[lx+i*incx] + dy[ly+i*incy], where lx = 1
  37.    if  incx >= 0, else lx = (-incx)*(n-1)+1 and ly is
  38.    defined in a similar way using incy.
  39.  
  40. */
  41.  
  42. {
  43.    int ix, iy, i, m;
  44.  
  45.    if ( n < 0 || da == 0. )
  46.       return;
  47.  
  48. /* Code for nonequal or nonpositive increments.  */
  49.  
  50.    if ( incx != incy || incx < 1 ) {
  51.       ix = 1;
  52.       iy = 1;
  53.       if ( incx < 0 )
  54.          ix = ( -n + 1 ) * incx + 1;
  55.       if ( incy < 0 )
  56.          iy = ( -n + 1 ) * incy + 1;
  57.       for ( i = 1 ; i <= n ; i++ ) {
  58.          dy[iy] = dy[iy] + da * dx[ix];
  59.          ix = ix + incx;
  60.          iy = iy + incy;
  61.       }
  62.       return;
  63.    }
  64.  
  65. /* Code for both increments equal to 1.   */
  66.  
  67. /* Clean-up loop so remaining vector length is a multiple of 4.  */
  68.  
  69.    if ( incx == 1 ) {
  70.       m = n % 4;
  71.       if ( m != 0 ) {
  72.          for ( i = 1 ; i <= m ; i++ )
  73.             dy[i] = dy[i] + da * dx[i];
  74.          if ( n < 4 )
  75.             return;
  76.       }
  77.       for ( i = m + 1 ; i <= n ; i = i + 4 ) {
  78.          dy[i] = dy[i] + da * dx[i];
  79.          dy[i+1] = dy[i+1] + da * dx[i+1];
  80.          dy[i+2] = dy[i+2] + da * dx[i+2];
  81.          dy[i+3] = dy[i+3] + da * dx[i+3];
  82.       }
  83.       return;
  84.    }
  85.  
  86. /* Code for equal, positive, nonunit increments.   */
  87.  
  88.    for ( i = 1 ; i <= n * incx ; i = i + incx )
  89.       dy[i] = da * dx[i] + dy[i];
  90.    return;
  91.  
  92. }
  93.  
  94.