/* Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730 */

int ISAMAX ( n, x, incx )
   int n ;
   float x[] ;
   int incx ;

     /* Find index of element having maximum absolute value.*/
     /*                                                     */
     /* Jack Dongarra, LINPACK, 3/11/78.                    */
     /* Adam Fritz, C, 2/22/85.                             */

{
   float Abs() ;
   float smax ;
   int i, ix ;
   int imax ;

   imax = 0 ;
   if ( n > 0 )
   {
      if ( n > 1 )
      {
         if ( incx > 1 )
         {
                                /* incx > 1 */
            ix = 0 ;
            smax = Abs(x[0]) ;
            ix = incx ;
            for ( i = 1; i < n; i++ )
            {
               if ( Abs(x[ix]) > smax )
               {
                  imax = i ;
                  smax = Abs(x[ix]) ;
               }
               ix = ix + incx ;
            }
         }
                                /* incx = 1 */
         else
         {
            smax = Abs(x[0]) ;
            for ( i = 1; i < n; i++ )
               if ( Abs(x[i]) > smax )
               {
                  imax = i ;
                  smax = Abs(x[i]) ;
               }
         }
      }
   }
   return ( imax ) ;
}

/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/

float SASUM ( n, x, incx )
   int n ;
   float x[] ;
   int incx ;

     /* Forms sum of absolute values.                       */
     /*                                                     */
     /* Jack Dongarra, LINPACK, 3/11/78.                    */
     /* Adam Fritz, C, 2/22/85.                             */

{
   float Abs() ;
   float stemp ;
   int i, ix ;

   stemp = 0.0 ;
   if ( n > 0 ) {
      if ( incx > 1 )
      {
                                /* incx > 1 */
         ix = 0 ;
         for ( i = 0; i < n; i++ )
         {
            stemp = stemp + Abs(x[ix]) ;
            ix = ix + incx ;
         }
      }
                                /* incx = 1 */
      else
         for ( i = 0; i < n; i++ )
            stemp = stemp + Abs(x[i]) ;
   }
   return( stemp ) ;
}

/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/

SAXPY ( n, sa, x, incx, y, incy )
   int n ;
   float sa ;
   float x[] ;
   int incx ;
   float y[] ;
   int incy ;

     /* Compute constant times a vector plus a vector.      */
     /*                                                     */
     /* Jack Dongarra, LINPACK, 3/11/78.                    */
     /* Adam Fritz, C, 2/22/85.                             */

{
   int i, ix, iy ;

   if ( n > 0 )
   {
      if ( sa != 0.0 )
      {
         if ( (incx != 1) || (incy != 1) )
         {
                                /* incx != incy or incx != 1 */
            ix = 0 ;
            iy = 0 ;
            if ( incx < 0 )
               ix = (-n + 1) * incx ;
            if ( incy < 0 )
               iy = (-n + 1) * incy ;
            for ( i = 0; i < n; i++ )
            {
               y[iy] = y[iy] + sa * x[ix] ;
               ix = ix + incx ;
               iy = iy + incy ;
            }
         }
                                      /* incx, incy = 1 */
         else
            for ( i = 0; i < n; i++ )
               y[i] = y[i] + sa * x[i] ;
      }
   }
}

/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/

SCOPY ( n, x, incx, y, incy )
   int n ;
   float x[] ;
   int incx ;
   float y[] ;
   int incy ;

     /* Copies a vector, x, to a vector, y.                 */
     /*                                                     */
     /* Jack Dongarra, LINPACK, 3/11/78.                    */
     /* Adam Fritz, C, 2/22/85.                             */

{
   int i, ix, iy ;

   if ( n > 0 )
   {
      if ( (incx != 1) || (incy != 1) )
      {
                                /* incx, incy  != 1 */
         ix = 0 ;
         iy = 0 ;
         if ( incx < 0 )
            ix = (-n + 1) * incx ;
         if ( incy < 0 )
            iy = (-n + 1) * incy ;
         for ( i = 0; i < n; i++ )
         {
            y[iy] = x[ix] ;
            ix = ix + incx ;
            iy = iy + incy ;
         }
      }
                                /* incx, incx = 1 */
      else
         for ( i = 0; i < n; i++ )
            y[i] = x[i] ;
   }
}

/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/

float SDOT ( n, x, incx, y, incy )
   int n ;
   float x[] ;
   int incx ;
   float y[] ;
   int incy ;

     /* Computes dot product of two vectors.                */
     /*                                                     */
     /* Jack Dongarra, LINPACK, 3/11/78.                    */
     /* Adam Fritz, C, 2/22/85.                             */

{
   float stemp ;
   int i, ix, iy ;

   stemp = 0.0 ;
   if ( n > 0 )
   {
      if ( (incx != 1) || (incy != 1) )
      {
                                /* incx != incy or incx != 1 */
         ix = 0 ;
         iy = 0 ;
         if ( incx < 0 )
            ix = (-n + 1) * incx ;
         if ( incy < 0 )
            iy = (-n + 1) * incy ;
         for ( i = 0; i < n; i++ )
         {
            stemp = stemp + x[ix] * y[iy] ;
            ix = ix + incx ;
            iy = iy + incy ;
         }
      }
                                /* incx, incy = 1 */
      else
         for ( i = 0; i < n; i++ )
            stemp = stemp + x[i] * y[i] ;
   }
   return( stemp ) ;
}

/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/

SSCAL ( n, sa, x, incx )
   int n ;
   float sa ;
   float x[] ;
   int incx ;

     /* Computes vector scaled by a constant.               */
     /*                                                     */
     /* Jack Dongarra, LINPACK, 3/11/78.                    */
     /* Adam Fritz, C, 2/22/85.                             */

{
   int i, ix ;

   if ( n > 0 )
   {
      if ( incx != 1 )
      {
                                /* incx != 1 */
         ix = 0 ;
         for ( i = 0; i < n; i++ )
         {
            x[ix] = sa * x[ix] ;
            ix = ix + incx ;
         }
      }
                                /* incx = 1 */
      else
         for ( i = 0; i < n; i++ )
            x[i] = sa * x[i] ;
   }
}

/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/

SSWAP ( n, x, incx, y, incy )
   int n ;
   float x[] ;
   int incx ;
   float y[] ;
   int incy ;

     /* Interchanges two vectors.                           */
     /*                                                     */
     /* Jack Dongarra, LINPACK, 3/11/78.                    */
     /* Adam Fritz, C, 2/22/85.                             */

{
   float stemp ;
   int i, ix, iy ;

   if ( n > 0 )
   {
      if ( (incx != 1) || (incy != 1) )
      {
                                /* incx != incy or incx != 1 */
         ix = 0 ;
         iy = 0 ;
         if ( incx < 0 )
            ix = (-n + 1) * incx ;
         if ( incy < 0 )
            iy = (-n + 1) * incy ;
         for ( i = 0; i < n; i++ )
         {
            stemp = x[ix] ;
            x[ix] = y[iy] ;
            y[iy] = stemp ;
            ix = ix + incx ;
            iy = iy + incy ;
         }
      }
                                /* incx, incy = 1 */
      else for ( i = 0; i < n; i++ )
      {
         stemp = x[i] ;
         x[i] = y[i] ;
         y[i] = stemp ;
      }
   }
}

/* Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730 */