function isamax ( n, x, incx )
!
!*******************************************************************************
!
!! ISAMAX finds the index of the vector element of maximum absolute value.
!
!
!  Modified:
!
!    08 April 1999
!
!  Reference:
!
!    Lawson, Hanson, Kincaid, Krogh,
!    Basic Linear Algebra Subprograms for Fortran Usage,
!    Algorithm 539, 
!    ACM Transactions on Mathematical Software, 
!    Volume 5, Number 3, September 1979, pages 308-323.
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real X(*), the vector to be examined.
!
!    Input, integer INCX, the increment between successive entries of SX.
!
!    Output, integer ISAMAX, the index of the element of SX of maximum
!    absolute value.
!
  implicit none
!
  integer i
  integer incx
  integer isamax
  integer ix
  integer n
  real samax
  real x(*)
!
  if ( n <= 0 ) then

    isamax = 0

  else if ( n == 1 ) then

    isamax = 1

  else if ( incx == 1 ) then

    isamax = 1
    samax = abs ( x(1) )

    do i = 2, n

      if ( abs ( x(i) ) > samax ) then
        isamax = i
        samax = abs ( x(i) )
      end if

    end do

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( - n + 1 ) * incx + 1
    end if

    isamax = 1
    samax = abs ( x(ix) )

    ix = ix + incx

    do i = 2, n
      if ( abs ( x(ix) ) > samax ) then
        isamax = i
        samax = abs ( x(ix) )
      end if
      ix = ix + incx
    end do

  end if

  return
end
function isamin ( n, x, incx )
!
!*******************************************************************************
!
!! ISAMIN finds the index of the vector element having minimum absolute value.
!
!
!  Modified:
!
!    08 April 1999
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real X(*), the vector to be examined.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Output, integer ISAMIN, the index of the element of SX of minimum
!    absolute value.
!
  implicit none
!
  integer i
  integer incx
  integer isamin
  integer ix
  integer n
  real samin
  real x(*)
!
  if ( n <= 0 ) then

    isamin = 0

  else if ( n == 1 ) then

    isamin = 1

  else if ( incx == 1 ) then

    isamin = 1
    samin = abs ( x(1) )

    do i = 2, n

      if ( abs ( x(i) ) < samin ) then
        isamin = i
        samin = abs ( x(i) )
      end if

    end do

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( - n + 1 ) * incx + 1
    end if

    isamin = 1
    samin = abs ( x(ix) )

    ix = ix + incx

    do i = 2, n
      if ( abs ( x(ix) ) < samin ) then
        isamin = i
        samin = abs ( x(ix) )
      end if
      ix = ix + incx
    end do

  end if

  return
end
function ismax ( n, x, incx )
!
!*******************************************************************************
!
!! ISMAX finds the index of the vector element having maximum value.
!
!
!  Modified:
!
!    08 April 1999
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real X(*), the vector to be examined.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Output, integer ISMAX, the index of the element of X of maximum
!    value.
!
  implicit none
!
  integer i
  integer incx
  integer ismax
  integer ix
  integer n
  real smax
  real x(*)
!
  if ( n <= 0 ) then

    ismax = 0

  else if ( n == 1 ) then

    ismax = 1

  else if ( incx == 1 ) then

    ismax = 1
    smax = x(1)

    do i = 2, n
      if ( x(i) > smax ) then
        ismax = i
        smax = x(i)
      end if
    end do

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( - n + 1 ) * incx + 1
    end if

    ismax = 1
    smax = x(ix)
    ix = ix + incx

    do i = 2, n
      if ( x(ix) > smax ) then
        ismax = i
        smax = x(ix)
      end if
      ix = ix + incx
    end do

  end if

  return
end
function ismin ( n, x, incx )
!
!*******************************************************************************
!
!! ISMIN finds the index of the vector element having minimum value.
!
!
!  Modified:
!
!    09 April 1999
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real X(*), the vector to be examined.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Output, integer ISMIN, the index of the element of X of minimum
!    value.
!
  implicit none
!
  integer i
  integer incx
  integer ismin
  integer ix
  integer n
  real smin
  real x(*)
!
  if ( n <= 0 ) then

    ismin = 0

  else if ( n == 1 ) then

    ismin = 1

  else if ( incx == 1 ) then

    ismin = 1
    smin = x(1)

    do i = 2, n
      if ( x(i) < smin ) then
        ismin = i
        smin = x(i)
      end if
    end do

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( - n + 1 ) * incx + 1
    end if

    ismin = 1
    smin = x(ix)
    ix = ix + incx

    do i = 2, n
      if ( x(ix) < smin ) then
        ismin = i
        smin = x(ix)
      end if
      ix = ix + incx
    end do

  end if

  return
end
function samax ( n, x, incx )
!
!*******************************************************************************
!
!! SAMAX returns the maximum absolute value of the entries in a vector.
!
!
!  Modified:
!
!    08 April 1999
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real X(*), the vector to be examined.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Output, real SAMAX, the maximum absolute value of an element of X.
!
  implicit none
!
  integer i
  integer incx
  integer ix
  integer n
  real samax
  real x(*)
!
  if ( n <= 0 ) then

    samax = 0.0E+00

  else if ( n == 1 ) then

    samax = abs ( x(1) )

  else if ( incx == 1 ) then

    samax = abs ( x(1) )

    do i = 2, n
      if ( abs ( x(i) ) > samax ) then
        samax = abs ( x(i) )
      end if
    end do

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( - n + 1 ) * incx + 1
    end if

    samax = abs ( x(ix) )
    ix = ix + incx

    do i = 2, n
      if ( abs ( x(ix) ) > samax ) then
        samax = abs ( x(ix) )
      end if
      ix = ix + incx
    end do

  end if

  return
end
function samin ( n, x, incx )
!
!*******************************************************************************
!
!! SAMIN returns the minimum absolute value of entries in a vector.
!
!
!  Modified:
!
!    08 April 1999
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real X(*), the vector to be examined.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Output, real SAMIN, the minimum absolute value of an element of X.
!
  implicit none
!
  integer i
  integer incx
  integer ix
  integer n
  real samin
  real x(*)
!
  if ( n <= 0 ) then

    samin = 0.0E+00

  else if ( n == 1 ) then

    samin = abs ( x(1) )

  else if ( incx == 1 ) then

    samin = abs ( x(1) )

    do i = 2, n
      if ( abs ( x(i) ) < samin ) then
        samin = abs ( x(i) )
      end if
    end do

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( - n + 1 ) * incx + 1
    end if

    samin = abs ( x(ix) )
    ix = ix + incx

    do i = 2, n
      if ( abs ( x(ix) ) < samin ) then
        samin = abs ( x(ix) )
      end if
      ix = ix + incx
    end do

  end if

  return
end
function sasum ( n, x, incx )
!
!*******************************************************************************
!
!! SASUM sums the absolute values of the entries of a vector.
!
!
!  Modified:
!
!    15 February 2001
!
!  Reference:
!
!    Lawson, Hanson, Kincaid, Krogh,
!    Basic Linear Algebra Subprograms for Fortran Usage,
!    Algorithm 539, 
!    ACM Transactions on Mathematical Software, 
!    Volume 5, Number 3, September 1979, pages 308-323.
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real X(*), the vector to be examined.
!
!    Input, integer INCX, the increment between successive entries of X.
!    INCX must not be negative.
!
!    Output, real SASUM, the sum of the absolute values of X.
!
  implicit none
!
  integer incx
  integer n
  real sasum
  real x(*)
!
  sasum = sum ( abs ( x(1:1+(n-1)*incx:incx) ) )

  return
end
subroutine saxpy ( n, sa, x, incx, y, incy )
!
!*******************************************************************************
!
!! SAXPY adds a constant times one vector to another.
!
!
!  Modified:
!
!    08 April 1999
!
!  Reference:
!
!    Lawson, Hanson, Kincaid, Krogh,
!    Basic Linear Algebra Subprograms for Fortran Usage,
!    Algorithm 539, 
!    ACM Transactions on Mathematical Software, 
!    Volume 5, Number 3, September 1979, pages 308-323.
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real SA, the multiplier.
!
!    Input, real X(*), the vector to be scaled and added to Y.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Input/output, real Y(*), the vector to which a multiple of X is to 
!    be added.
!
!    Input, integer INCY, the increment between successive entries of Y.
!
  implicit none
!
  integer i
  integer incx
  integer incy
  integer ix
  integer iy
  integer n
  real sa
  real x(*)
  real y(*)
!
  if ( n <= 0 ) then

  else if ( sa == 0.0E+00 ) then

  else if ( incx == 1 .and. incy == 1 ) then

    y(1:n) = y(1:n) + sa * x(1:n)

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( - n + 1 ) * incx + 1
    end if

    if ( incy >= 0 ) then
      iy = 1
    else
      iy = ( - n + 1 ) * incy + 1
    end if

    do i = 1, n
      y(iy) = y(iy) + sa * x(ix)
      ix = ix + incx
      iy = iy + incy
    end do

  end if

  return
end
subroutine saxpyi ( nz, sa, x, indx, y )
!
!*******************************************************************************
!
!! SAXPYI adds a scalar times a sparse vector to a vector.
!
!
!  Modified:
!
!    08 April 1999
!
!  Parameters:
!
!    Input, integer NZ, the number of entries in the sparse vector.
!
!    Input, real SA, the multiplier.
!
!    Input, integer X(NZ), the sparse vector to be scaled and added to Y.
!
!    Input, integer INDX(NZ), the full storage indices of the sparse vector.
!
!    Input/output, real Y(*), the full storage vector.  On output,
!    a scaled copy of X has been added to Y.
!
  implicit none
!
  integer nz
!
  integer i
  integer indx(nz)
  real sa
  real x(nz)
  real y(*)
!
  if ( sa /= 0.0E+00 ) then

    do i = 1, nz
      y(indx(i))  = y(indx(i)) + sa * x(i)
    end do

  end if

  return
end
subroutine saxpyx ( n, sa, x, incx, y, incy )
!
!*******************************************************************************
!
!! SAXPYX replaces a vector by a constant multiple plus another vector.
!
!
!  Modified:
!
!    08 April 1999
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vectors.
!
!    Input, real SA, the multiplier.
!
!    Input/output, real X(*), the vector to be scaled and increased by Y.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Input, real Y(*), the vector to be added to X.
!
!    Input, integer INCY, the increment between successive elements of Y.
!
  implicit none
!
  integer i
  integer incx
  integer incy
  integer ix
  integer iy
  integer n
  real sa
  real x(*)
  real y(*)
!
  if ( n <= 0 ) then

  else if ( incx == 1 .and. incy == 1 ) then

    x(1:n) = y(1:n) + sa * x(1:n)

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( - n + 1 ) * incx + 1
    end if

    if ( incy >= 0 ) then
      iy = 1
    else
      iy = ( - n + 1 ) * incy + 1
    end if

    do i = 1, n
      x(ix) = y(iy) + sa * x(ix)
      ix = ix + incx
      iy = iy + incy
    end do

  end if

  return
end
subroutine scopy ( n, x, incx, y, incy )
!
!*******************************************************************************
!
!! SCOPY copies one real vector into another.
!
!
!  Modified:
!
!    08 April 1999
!
!  Reference:
!
!    Lawson, Hanson, Kincaid, Krogh,
!    Basic Linear Algebra Subprograms for Fortran Usage,
!    Algorithm 539, 
!    ACM Transactions on Mathematical Software, 
!    Volume 5, Number 3, September 1979, pages 308-323.
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real X(*), the vector to be copied into Y.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Output, real Y(*), the copy of X.
!
!    Input, integer INCY, the increment between successive elements of Y.
!
  implicit none
!
  integer i
  integer incx
  integer incy
  integer ix
  integer iy
  integer n
  real x(*)
  real y(*)
!
  if ( n <= 0 ) then

  else if ( incx == 1 .and. incy == 1 ) then

    y(1:n) = x(1:n)

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( - n + 1 ) * incx + 1
    end if

    if ( incy >= 0 ) then
      iy = 1
    else
      iy = ( - n + 1 ) * incy + 1
    end if

    do i = 1, n
      y(iy) = x(ix)
      ix = ix + incx
      iy = iy + incy
    end do

  end if

  return
end
function sdot ( n, x, incx, y, incy )
!
!*******************************************************************************
!
!! SDOT forms the dot product of two vectors.
!
!
!  Modified:
!
!    02 June 2000
!
!  Reference:
!
!    Lawson, Hanson, Kincaid, Krogh,
!    Basic Linear Algebra Subprograms for Fortran Usage,
!    Algorithm 539, 
!    ACM Transactions on Mathematical Software, 
!    Volume 5, Number 3, September 1979, pages 308-323.
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vectors.
!
!    Input, real X(*), one of the vectors to be multiplied.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Input, real Y(*), one of the vectors to be multiplied.
!
!    Input, integer INCY, the increment between successive elements of Y.
!
!    Output, real SDOT, the dot product of X and Y.
!
  implicit none
!
  integer i
  integer incx
  integer incy
  integer ix
  integer iy
  integer n
  real sdot
  real stemp
  real x(*)
  real y(*)
!
  if ( n <= 0 ) then

    sdot = 0.0E+00

  else if ( incx == 1 .and. incy == 1 ) then

    sdot = dot_product ( x(1:n), y(1:n) )

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( - n + 1 ) * incx + 1
    end if

    if ( incy >= 0 ) then
      iy = 1
    else
      iy = ( - n + 1 ) * incy + 1
    end if

    stemp = 0.0E+00
    do i = 1, n
      stemp = stemp + x(ix) * y(iy)
      ix = ix + incx
      iy = iy + incy
    end do

    sdot = stemp

  end if

  return
end
function sdoti ( nz, x, indx, y )
!
!*******************************************************************************
!
!! SDOTI forms the dot product of a sparse vector and a full vector.
!
!
!  Modified:
!
!    08 April 1999
!
!  Parameters:
!
!    Input, integer NZ, the number of entries in the sparse vector.
!
!    Input, real X(NZ), the sparse vector.
!
!    Input, integer INDX(NZ), the full storage indices of the sparse vector.
!
!    Input, real Y(*), the full storage vector.
!
!    Output, real SDOTI, the dot product of X and Y.
!
  implicit none
!
  integer nz
!
  integer i
  integer indx(nz)
  real sdoti
  real x(nz)
  real y(*)
!
  sdoti = 0.0E+00

  do i = 1, nz
    sdoti = sdoti + x(i) * y(indx(i))
  end do

  return
end
function sdsdot ( n, x, incx, y, incy )
!
!*******************************************************************************
!
!! SDSDOT forms the dot product of two vectors using higher precision.
!
!
!  Modified:
!
!    02 June 2000
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vectors.
!
!    Input, real X(*), one of the vectors to be multiplied.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Input, real Y(*), one of the vectors to be multiplied.
!
!    Input, integer INCY, the increment between successive elements of Y.
!
!    Output, real SDSDOT, the dot product of X and Y.
!
  implicit none
!
  double precision dsdot
  integer i
  integer incx
  integer incy
  integer ix
  integer iy
  integer n
  real sdsdot
  real x(*)
  real y(*)
!
  if ( n <= 0 ) then

    dsdot = 0.0D+00

  else if ( incx == 1 .and. incy == 1 ) then

    dsdot = dot_product ( dble ( x(1:n) ), dble ( y(1:n) ) )

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( - n + 1 ) * incx + 1
    end if

    if ( incy >= 0 ) then
      iy = 1
    else
      iy = ( - n + 1 ) * incy + 1
    end if

    dsdot = 0.0D+00
    do i = 1, n
      dsdot = dsdot + dble ( x(ix) ) * dble ( y(iy) )
      ix = ix + incx
      iy = iy + incy
    end do

  end if

  sdsdot = sngl ( dsdot )

  return
end
subroutine sgthr ( nz, y, x, indx )
!
!*******************************************************************************
!
!! SGTHR gathers specified elements of a full vector into a sparse vector.
!
!
!  Modified:
!
!    08 April 1999
!
!  Parameters:
!
!    Input, integer NZ, the number of entries in the sparse vector.
!
!    Input, real Y(*), the full storage vector.
!
!    Output, real X(NZ), the sparse vector gathered from Y.
!
!    Input, integer INDX(NZ), the full storage indices of the sparse vector.
!
  implicit none
!
  integer nz
!
  integer i
  integer indx(nz)
  real x(nz)
  real y(*)
!
  do i = 1, nz
    x(i) = y(indx(i))
  end do

  return
end
subroutine sgthrz ( nz, y, x, indx )
!
!*******************************************************************************
!
!! SGTHRZ gathers elements of a full vector into a sparse one, and zeroes the originals.
!
!
!  Modified:
!
!    08 April 1999
!
!  Parameters:
!
!    Input, integer NZ, the number of entries in the sparse vector.
!
!    Input/output, real Y(*), the full storage vector.  On output,
!    the indexed entries of Y have been zeroed out.
!
!    Output, real X(NZ), the sparse vector gathered from Y.
!
!    Input, integer INDX(NZ), the full storage indices of the sparse vector.
!
  implicit none
!
  integer nz
!
  integer i
  integer indx(nz)
  real x(nz)
  real y(*)
!
  do i = 1, nz
    x(i) = y(indx(i))
    y(indx(i)) = 0.0E+00
  end do

  return
end
subroutine sinit ( n, sa, x, incx )
!
!*******************************************************************************
!
!! SINIT initializes a real vector to a constant.
!
!
!  Modified:
!
!    08 April 1999
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real SA, the constant value to be used to initialize X.
!
!    Output, real X(*), the vector to be initialized.
!
!    Input, integer INCX, the increment between successive entries of X.
!
  implicit none
!
  integer i
  integer incx
  integer ix
  integer n
  real sa
  real x(*)
!
  if ( n <= 0 ) then

  else if ( incx == 1 ) then

    x(1:n) = sa

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( - n + 1 ) * incx + 1
    end if

    do i = 1, n
      x(ix) = sa
      ix = ix + incx
    end do

  end if

  return
end
function smach ( job )
!
!*******************************************************************************
!
!! SMACH computes real machine parameters.
!
!
!  Discussion:
!
!    Assuming the computer has
!
!      B = base of arithmetic
!      T = number of base  b  digits
!      L = smallest possible exponent
!      U = largest possible exponent
!
!    then
!
!      EPS = b**(1-T)
!      TEENY = 100.0 * B**(-L+T)
!      BIG = 0.01 * B**(U-T)
!
!  Modified:
!
!    08 April 1999
!
!  Parameters:
!
!    Input, integer JOB, specifies the quantity to be computed.
!    1, EPS
!    2, TINY
!    3, BIG
!
!    Output, real SMACH, the requested quantity.
!
  implicit none
!
  real big
  real eps
  integer job
  real s
  real smach
  real teeny
!
  eps = epsilon ( eps )

  s = 1.0E+00

  do
    teeny = s
    s = s / 16.0E+00
    if ( s * 100.0E+00 == 0.0E+00 ) then
      exit
    end if
  end do

  teeny = ( teeny / eps ) * 100.0E+00
  big = 1.0E+00 / teeny

  if ( job == 1 ) then
    smach = eps
  else if ( job == 2 ) then
    smach = teeny
  else if ( job == 3 ) then
    smach = big
  end if

  return
end
function smax ( n, x, incx )
!
!*******************************************************************************
!
!! SMAX returns the maximum value of all entries in a vector.
!
!
!  Modified:
!
!    08 April 1999
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real X(*), the vector to be examined.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Output, real SMAX, the maximum entry of X.
!
  implicit none
!
  integer i
  integer incx
  integer ix
  integer n
  real smax
  real x(*)
!
  if ( n <= 0 ) then

    smax = 0.0E+00

  else if ( n == 1 ) then

    smax = x(1)

  else if ( incx == 1 ) then

    smax = x(1)
    do i = 2, n
      if ( x(i) > smax ) then
        smax = x(i)
      end if
    end do

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( - n + 1 ) * incx + 1
    end if

    smax = x(ix)
    ix = ix + incx
    do i = 2, n
      if ( x(ix) > smax ) then
        smax = x(ix)
      end if
      ix = ix + incx
    end do

  end if

  return
end
function smin ( n, x, incx )
!
!*******************************************************************************
!
!! SMIN returns the minimum value of all entries in a vector.
!
!
!  Modified:
!
!    08 April 1999
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real X(*), the vector to be examined.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Output, real SMIN, the value of the minimum entry of X.
!
  implicit none
!
  integer i
  integer incx
  integer ix
  integer n
  real smin
  real x(*)
!
  if ( n <= 0 ) then

    smin = 0.0E+00

  else if ( n == 1 ) then

    smin = x(1)

  else if ( incx == 1 ) then

    smin = x(1)
    do i = 2, n
      if ( x(i) < smin ) then
        smin = x(i)
      end if
    end do

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( - n + 1 ) * incx + 1
    end if

    smin = x(ix)
    ix = ix + incx
    do i = 2, n
      if ( x(ix) < smin ) then
        smin = x(ix)
      end if
      ix = ix + incx
    end do

  end if

  return
end
function snrm2 ( n, x, incx )
!
!*******************************************************************************
!
!! SNRM2 computes the Euclidean norm of a vector.
!
!
!  Discussion:
!
!    The original SNRM2 algorithm is accurate but written in a bizarre,
!    unreadable and obsolete format.  This version goes for clarity.
!
!  Modified:
!
!    01 June 2000
!
!  Reference:
!
!    Lawson, Hanson, Kincaid, Krogh,
!    Basic Linear Algebra Subprograms for Fortran Usage,
!    Algorithm 539, 
!    ACM Transactions on Mathematical Software, 
!    Volume 5, Number 3, September 1979, pages 308-323.
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real X(*), the vector whose norm is to be computed.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Output, real SNRM2, the Euclidean norm of X.
!
  implicit none
!
  integer i
  integer incx
  integer ix
  integer n
  real samax
  real snrm2
  real stemp
  real x(*)
  real xmax
!
  if ( n <= 0 ) then

    snrm2 = 0.0E+00

  else

    xmax = samax ( n, x, incx )

    if ( xmax == 0.0E+00 ) then

      snrm2 = 0.0E+00

    else

      if ( incx >= 0 ) then
        ix = 1
      else
        ix = ( - n + 1 ) * incx + 1
      end if

      stemp = 0.0E+00
      do i = 1, n
        stemp = stemp + ( x(ix) / xmax )**2
        ix = ix + incx
      end do

      snrm2 = xmax * sqrt ( stemp )

    end if

  end if

  return
end
subroutine spaxpy ( nz, sa, x, y, indx )
!
!*******************************************************************************
!
!! SPAXPY adds a scalar times a sparse vector to a full one.
!
!
!  Modified:
!
!    08 April 1999
!
!  Parameters:
!
!    Input, integer NZ, the number of entries in the sparse vector.
!
!    Input, real SA, the multiplier.
!
!    Input, real X(NZ), the sparse vector to be scaled and added to Y.
!
!    Input/output, real Y(*), the full storage vector to which a scaled
!    copy of X is added.
!
!    Input, integer INDX(NZ), the full storage indices of the sparse vector.
!
  implicit none
!
  integer nz
!
  integer i
  integer indx(nz)
  real sa
  real x(nz)
  real y(*)
!
  do i = 1, nz
    y(indx(i)) = sa * x(i) + y(indx(i))
  end do

  return
end
function spdot ( nz, y, indx, x )
!
!*******************************************************************************
!
!! SPDOT computes the product of a sparse vector and a full one.
!
!
!  Modified:
!
!    08 April 1999
!
!  Parameters:
!
!    Input, integer NZ, the number of entries in the sparse vector.
!
!    Input, real Y(*), the full storage vector.
!
!    Input, integer INDX(NZ), the full storage indices of the sparse vector.
!
!    Input, real X(NZ), the sparse vector.
!
  implicit none
!
  integer nz
!
  integer i
  integer indx(nz)
  real spdot
  real sum
  real x(nz)
  real y(*)
!
  sum = 0.0E+00
  do i = 1, nz
    sum = sum + y(indx(i)) * x(i)
  end do

  spdot = sum

  return
end
subroutine srot ( n, x, incx, y, incy, c, s )
!
!*******************************************************************************
!
!! SROT applies a plane rotation.
!
!
!  Modified:
!
!    08 April 1999
!
!  Reference:
!
!    Lawson, Hanson, Kincaid, Krogh,
!    Basic Linear Algebra Subprograms for Fortran Usage,
!    Algorithm 539, 
!    ACM Transactions on Mathematical Software, 
!    Volume 5, Number 3, September 1979, pages 308-323.
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vectors.
!
!    Input/output, real X(*), one of the vectors to be rotated.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Input/output, real Y(*), one of the vectors to be rotated.
!
!    Input, integer INCY, the increment between successive elements of Y.
!
!    Input, real C, S, parameters (presumably the cosine and sine of
!    some angle) that define a plane rotation.
!
  implicit none
!
  real c
  integer i
  integer incx
  integer incy
  integer ix
  integer iy
  integer n
  real s
  real stemp
  real x(*)
  real y(*)
!
  if ( n <= 0 ) then

  else if ( incx == 1 .and. incy == 1 ) then

    do i = 1, n
      stemp = c * x(i) + s * y(i)
      y(i) = c * y(i) - s * x(i)
      x(i) = stemp
    end do

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( - n + 1 ) * incx + 1
    end if

    if ( incy >= 0 ) then
      iy = 1
    else
      iy = ( - n + 1 ) * incy + 1
    end if

    do i = 1, n
      stemp = c * x(ix) + s * y(iy)
      y(iy) = c * y(iy) - s * x(ix)
      x(ix) = stemp
      ix = ix + incx
      iy = iy + incy
    end do

  end if

  return
end
subroutine srotg ( sa, sb, c, s )
!
!*******************************************************************************
!
!! SROTG constructs a Givens plane rotation.
!
!
!  Modified:
!
!    08 April 1999
!
!  Reference:
!
!    Lawson, Hanson, Kincaid, Krogh,
!    Basic Linear Algebra Subprograms for Fortran Usage,
!    Algorithm 539, 
!    ACM Transactions on Mathematical Software, 
!    Volume 5, Number 3, September 1979, pages 308-323.
!
!  Parameters:
!
!    Input/output, real SA, SB, ...
!
!    Output, real C, S, ...
!
  implicit none
!
  real c
  real r
  real roe
  real s
  real sa
  real sb
  real scale
  real z
!
  if ( abs ( sa ) > abs ( sb ) ) then
    roe = sa
  else
    roe = sb
  end if

  scale = abs ( sa ) + abs ( sb )

  if ( scale == 0.0E+00 ) then
    c = 1.0E+00
    s = 0.0E+00
    r = 0.0E+00
  else
    r = scale * sqrt ( ( sa / scale )**2 + ( sb / scale )**2 )
    r = sign ( 1.0E+00, roe ) * r
    c = sa / r
    s = sb / r
  end if

  if ( abs ( c ) > 0.0E+00 .and. abs ( c ) <= s ) then
    z = 1.0E+00 / c
  else
    z = s
  end if

  sa = r
  sb = z

  return
end
subroutine sroti ( nz, x, indx, y, c, s )
!
!*******************************************************************************
!
!! SROTI applies a Givens rotation to a sparse vector.
!
!
!  Modified:
!
!    01 June 2000
!
!  Parameters:
!
!    Input, integer NZ, the number of entries in the sparse vector.
!
!    Input/output, real X(NZ), the sparse vector to be rotated.
!
!    Input, integer INDX(NZ), the full storage indices of the sparse vector.
!
!    Input/output, real Y(*), the full storage vector to be rotated.
!
!    Input, real C, S, the parameters that define the rotation.
!
  implicit none
!
  integer nz
!
  real c
  integer i
  integer indx(nz)
  real s
  real temp
  real x(nz)
  real y(*)
!
  if ( c /= 1.0E+00 .or. s /= 0.0E+00 ) then

    do i = 1, nz
      temp       = - s * x(i)  +  c * y(indx(i))
      x(i)       =   c * x(i)  +  s * y(indx(i))
      y(indx(i)) = temp
    end do

  end if

  return
end
subroutine srotm ( n, x, incx, y, incy, sparam )
!
!*******************************************************************************
!
!! SROTM applies a modified Givens plane rotation.
!
!
!  Discussion:
!
!    The modified Givens rotation matrix H is applied to the 2 by N matrix
!
!      ( x(1) ... x(n) )
!      ( y(1) ... y(n) )
!
!  Modified:
!
!    28 December 2000
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vectors.
!
!    Input/output, real X(*), one of the vectors to be rotated.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Input/output, real Y(*), one of the vectors to be rotated.
!
!    Input, integer INCY, the increment between successive elements of Y.
!
!    Input, real SPARAM(5), defines the rotation matrix H.
!    H takes one of the following forms, determined by SPARAM(1):
!
!    sparam(1) = -2.0, then
!
!       h11 = 1.0        h12 = 0.0
!       h21 = 0.0        h22 = 1.0
!
!    sparam(1) = -1.0, then
!
!       h11 = sparam(2)  h12 = sparam(4)
!       h21 = sparam(3)  h22 = sparam(5)
!
!    sparam(1) = 0.0, then
!
!       h11 = 1.0        h12 = sparam(4)
!       h21 = sparam(3)  h22 = 1.0
!
!    sparam(1) = 1.0, then
!
!       h11 = sparam(2)  h12 = 1.0
!       h21 = -1.0       h22 = sparam(5)
!
  implicit none
!
  integer i
  integer incx
  integer incy
  integer ix
  integer iy
  integer n
  real sh11
  real sh12
  real sh21
  real sh22
  real sparam(5)
  real x(*)
  real y(*)
  real w
  real z
!
  if ( n <= 0 ) then
!
!  Identity matrix.
!
  else if ( sparam(1) == -2.0E+00 ) then
!
!  S(2) S(4)
!  S(3) S(5)
!
  else if ( sparam(1) == -1.0E+00 ) then

    sh11 = sparam(2)
    sh12 = sparam(4)
    sh21 = sparam(3)
    sh22 = sparam(5)

    if ( incx == 1 .and. incy == 1 ) then

      do i = 1, n
        w = x(i)
        z = y(i)
        x(i) = w * sh11 + z * sh12
        y(i) = w * sh21 + z * sh22
      end do

    else

      if ( incx >= 0 ) then
        ix = 1
      else
        ix = 1 + ( 1 - n ) * incx
      end if

      if ( incy >= 0 ) then
        iy = 1
      else
        iy = 1 + ( 1 - n ) * incy
      end if

      do i = 1, n
        w = x(ix)
        z = y(iy)
        x(ix) = w*sh11 + z * sh12
        y(iy) = w*sh21 + z * sh22
        ix = ix + incx
        iy = iy + incy
      end do

    end if
!
!  1.0  S(4)
!  S(3) 1.0
!
  else if ( sparam(1) ==  0.0E+00 ) then

    sh12 = sparam(4)
    sh21 = sparam(3)

    if ( incx == 1 .and. incy == 1 ) then

      do i = 1, n
        w = x(i)
        z = y(i)
        x(i) = w        + z * sh12
        y(i) = w * sh21 + z
      end do

    else

      if ( incx >= 0 ) then
        ix = 1
      else
        ix = 1 + ( 1 - n ) * incx
      end if

      if ( incy >= 0 ) then
        iy = 1
      else
        iy = 1 + ( 1 - n ) * incy
      end if

      do i = 1, n
        w = x(ix)
        z = y(iy)
        x(ix) = w        + z * sh12
        y(iy) = w * sh21 + z
        ix = ix + incx
        iy = iy + incy
      end do

    end if
!
!  S(2) 1.0
!  -1.0 S(5)
!
  else if ( sparam(1) == +1.0E+00 ) then

    sh11 = sparam(2)
    sh22 = sparam(5)

    if ( incx == 1 .and. incy == 1 ) then

      do i = 1, n
        w = x(i)
        z = y(i)
        x(i) =  w * sh11 +        z
        y(i) = -w        + sh22 * z
      end do

    else

      if ( incx >= 0 ) then
        ix = 1
      else
        ix = 1 + ( 1 - n ) * incx
      end if

      if ( incy >= 0 ) then
        iy = 1
      else
        iy = 1 + ( 1 - n ) * incy
      end if

      do i = 1, n
        w = x(ix)
        z = y(iy)
        x(ix) =  w * sh11 +        z
        y(iy) = -w        + sh22 * z
        ix = ix + incx
        iy = iy + incy
      end do

    end if

  else

    write ( *, * ) ' '
    write ( *, * ) 'SROTM - Fatal error!'
    write ( *, * ) '  Illegal input value of SPARAM(1) = ', sparam(1)
    stop

  end if

  return
end
subroutine srotmg ( sd1, sd2, sx1, sy1, sparam )
!
!*******************************************************************************
!
!! SROTMG constructs a modified Givens plane rotation.
!
!
!  Discussion:
!
!    SROTMG constructs a modified Givens rotation H and updates the
!    scale factors SD1 and SD2 which zero SY1.  
!
!    on input,  
!
!      sw1  =  sqrt(sd1)*sx1
!      sz1  =  sqrt(sd2)*sy1.
!
!    on output, 
!
!      ( c  s ) (sw1) = (c*sw1 + s*sz1) = (sqrt(sd1)*sx1)
!      (-s  c ) (sz1) = (      0      ) = (      0       )
!
!    where C and S define a Givens rotation.
!
!  Modified:
!
!    08 April 1999
!
!  Parameters:
!
!    Input/output, real SD1, a scale factor.
!    On input, SD1 contains the first scale factor.  On
!    output, SD1 contains the updated scale factor.
!
!    Input/output, real SD2, a scale factor. 
!    on input, sd2 contains the second scale factor.  on
!    output, sd2 contains the updated scale factor.
!
!    Input/output, real SX1.  On input, sx1 conatins the first component 
!    of the vector to be rotated.  On output sx1 contains the rotated value
!    of the first component.
!
!    Input, real SY1, the second component of the vector to be rotated.
!    Since this component is zeroed by the rotation, it is
!    left unchanged in storage.
!
!    Input/output, real SPARAM(5), defines the rotation matrix H.  
!
!    Depending on the value of SPARAM(1), H takes the form:
!
!     sparam(1) = -2.0
!       sparam(2) = unchanged   sparam(4) = unchanged
!       sparam(3) = unchanged   sparam(5) = unchanged
!
!     sparam(1) = -1.0
!       sparam(2) = h11         sparam(4) = h12
!       sparam(3) = h21         sparam(5) = h22
!
!     sparam(1) = 0.0
!       sparam(2) = unchanged   sparam(4) = h12
!       sparam(3) = h21         sparam(5) = unchanged
!
!     sparam(1) = 1.0
!       sparam(2) = h11         sparam(4) = unchanged
!       sparam(3) = unchanged   sparam(5) = h22
!
  implicit none
!
  real, save :: gam = 4096.0E+00
  real, save :: gamsq = 1.678E+07
  integer igo
  real, save :: rgamsq = 5.96E-08
  real sd1
  real sd2
  real sflag
  real sh11
  real sh12
  real sh21
  real sh22
  real sp1
  real sp2
  real sparam(5)
  real sq1
  real sq2
  real stemp
  real su
  real sx1
  real sy1
!
  if ( sd1 < 0.0E+00 ) then
    go to 60
  end if
!
!  go zero-h-d-and-sx1..
!
   10 continue
!
!  case-sd1-nonnegative
!
  sp2 = sd2 * sy1

  if ( sp2 == 0.0E+00 ) then
    sflag = -2.0E+00
    sparam(1) = sflag
    return
  end if
!
!  regular-case..
!
   20 continue
  sp1 = sd1*sx1
  sq2 = sp2*sy1
  sq1 = sp1*sx1

  if ( .not. abs(sq1) > abs(sq2) ) go to 40
  sh21 = -sy1/sx1
  sh12 = sp2/sp1

  su = 1.0E+00 - sh12*sh21

  if ( su <= 0.0E+00 ) then
    go to 60
  end if
!
!  go zero-h-d-and-sx1..
!
30 continue

  sflag = 0.0E+00
  sd1 = sd1/su
  sd2 = sd2/su
  sx1 = sx1*su
!
!  go scale-check..
!
  go to 100

40 continue

  if ( .not. sq2 < 0.0E+00 ) go to 50
!
!  go zero-h-d-and-sx1..
!
  go to 60

50 continue

  sflag = 1.0E+00
  sh11 = sp1/sp2
  sh22 = sx1/sy1
  su = 1.0E+00 + sh11*sh22
  stemp = sd2/su
  sd2 = sd1/su
  sd1 = stemp
  sx1 = sy1*su
!
!  go scale-check
!
  go to 100
!
!  procedure..zero-h-d-and-sx1..
!
60 continue

  sflag = -1.0E+00
  sh11 = 0.0E+00
  sh12 = 0.0E+00
  sh21 = 0.0E+00
  sh22 = 0.0E+00
  sd1 = 0.0E+00
  sd2 = 0.0E+00
  sx1 = 0.0E+00
  go to 230
!
!  procedure..fix-h..
!
70 continue

  if (.not.sflag >= 0.0E+00 ) go to 90

  if (.not.sflag == 0.0E+00 ) go to 80
  sh11 = 1.0E+00
  sh22 = 1.0E+00
  sflag = -1.0E+00
  go to 90

80 continue

  sh21 = -1.0E+00
  sh12 = 1.0E+00
  sflag = -1.0E+00

90 continue

  go to igo, (130, 160, 190, 220)
!
!  procedure..scale-check
!
100 continue

110 continue

  if (.not.sd1 <= rgamsq) go to 140
  120 continue
  if (sd1 == 0.0E+00 ) go to 170
  assign 130 to igo
!
!  fix-h..
!
  go to 70

130 continue
  sd1 = sd1*gam**2
  sx1 = sx1/gam
  sh11 = sh11/gam
  sh12 = sh12/gam
  go to 110
  140 continue
  150 continue
  if (.not.sd1 >= gamsq) go to 170
  assign 160 to igo
!
!  fix-h..
!
  go to 70
  160 continue
  sd1 = sd1/gam**2
  sx1 = sx1*gam
  sh11 = sh11*gam
  sh12 = sh12*gam
  go to 150
  170 continue
  180 continue
  if (.not.abs(sd2) <= rgamsq) go to 200
  if (sd2 == 0.0E+00 ) go to 230
  assign 190 to igo
!
!  fix-h..
!
  go to 70
  190 continue
  sd2 = sd2 * gam**2
  sh21 = sh21 / gam
  sh22 = sh22 / gam
  go to 180
  200 continue
  210 continue
  if ( .not. abs ( sd2 ) >= gamsq ) go to 230
  assign 220 to igo
!
!  fix-h..
!
  go to 70

  220 continue
  sd2 = sd2 / gam**2
  sh21 = sh21 * gam
  sh22 = sh22 * gam
  go to 210

  230 continue

  if ( sflag == 0.0E+00 ) then
    sparam(3) = sh21
    sparam(4) = sh12
  else if ( sflag > 0.0E+00 ) then
    sparam(2) = sh11
    sparam(5) = sh22
  else if ( sflag < 0.0E+00 ) then
    sparam(2) = sh11
    sparam(3) = sh21
    sparam(4) = sh12
    sparam(5) = sh22
  end if

  sparam(1) = sflag

  return
end
subroutine sscal ( n, sa, x, incx )
!
!*******************************************************************************
!
!! SSCAL scales a vector by a constant.
!
!
!  Modified:
!
!    08 April 1999
!
!  Reference:
!
!    Lawson, Hanson, Kincaid, Krogh,
!    Basic Linear Algebra Subprograms for Fortran Usage,
!    Algorithm 539, 
!    ACM Transactions on Mathematical Software, 
!    Volume 5, Number 3, September 1979, pages 308-323.
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real SA, the multiplier.
!
!    Input/output, real X(*), the vector to be scaled.
!
!    Input, integer INCX, the increment between successive entries of X.
!
  implicit none
!
  integer i
  integer incx
  integer ix
  integer m
  integer n
  real sa
  real x(*)
!
  if ( n <= 0 ) then

  else if ( incx == 1 ) then

    m = mod ( n, 5 )

    x(1:m) = sa * x(1:m)

    do i = m+1, n, 5
      x(i)   = sa * x(i)
      x(i+1) = sa * x(i+1)
      x(i+2) = sa * x(i+2)
      x(i+3) = sa * x(i+3)
      x(i+4) = sa * x(i+4)
    end do

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( - n + 1 ) * incx + 1
    end if

    do i = 1, n
      x(ix) = sa * x(ix)
      ix = ix + incx
    end do

  end if

  return
end
subroutine ssctr ( nz, x, indx, y )
!
!*******************************************************************************
!
!! SSCTR scatters the values of a sparse vector into a full one.
!
!
!  Modified:
!
!    08 April 1999
!
!  Parameters:
!
!    Input, integer NZ, the number of entries in the sparse vector.
!
!    Input, real X(NZ), the sparse vector to be scattered into Y.
!
!    Input, integer INDX(NZ), the full storage indices of the sparse vector.
!
!    Output, real Y(*), a full storage vector into which the values of
!    X have been copied.
!
  implicit none
!
  integer nz
!
  integer i
  integer indx(nz)
  real x(nz)
  real y(*)
!
  do i = 1, nz
    y(indx(i)) = x(i)
  end do

  return
end
function ssum ( n, x, incx )
!
!*******************************************************************************
!
!! SSUM sums the entries of a vector.
!
!
!  Modified:
!
!    02 June 2000
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vector.
!
!    Input, real X(*), the vector to be summed.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Output, real SSUM, the sum of the entries of X.
!
  implicit none
!
  integer i
  integer incx
  integer ix
  integer m
  integer n
  real ssum
  real stemp
  real x(*)
!
  if ( n <= 0 ) then

    stemp = 0.0E+00

  else if ( incx == 1 ) then

    stemp = sum ( x(1:n) )

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( - n + 1 ) * incx + 1
    end if

    stemp = 0.0E+00
    do i = 1, n
      stemp = stemp + x(ix)
      ix = ix + incx
    end do

  end if

  ssum = stemp

  return
end
subroutine sswap ( n, x, incx, y, incy )
!
!*******************************************************************************
!
!! SSWAP interchanges two vectors.
!
!
!  Modified:
!
!    08 April 1999
!
!  Reference:
!
!    Lawson, Hanson, Kincaid, Krogh,
!    Basic Linear Algebra Subprograms for Fortran Usage,
!    Algorithm 539, 
!    ACM Transactions on Mathematical Software, 
!    Volume 5, Number 3, September 1979, pages 308-323.
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vectors.
!
!    Input/output, real X(*), one of the vectors to swap.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Input/output, real Y(*), one of the vectors to swap.
!
!    Input, integer INCY, the increment between successive elements of Y.
!
  implicit none
!
  integer i
  integer incx
  integer incy
  integer ix
  integer iy
  integer m
  integer n
  real temp
  real x(*)
  real y(*)
!
  if ( n <= 0 ) then

  else if ( incx == 1 .and. incy == 1 ) then

    m = mod ( n, 3 )

    do i = 1, m
      temp = x(i)
      x(i) = y(i)
      y(i) = temp
    end do

    do i = m+1, n, 3

      temp = x(i)
      x(i) = y(i)
      y(i) = temp

      temp = x(i+1)
      x(i+1) = y(i+1)
      y(i+1) = temp

      temp = x(i+2)
      x(i+2) = y(i+2)
      y(i+2) = temp

    end do

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( - n + 1 ) * incx + 1
    end if

    if ( incy >= 0 ) then
      iy = 1
    else
      iy = ( - n + 1 ) * incy + 1
    end if

    do i = 1, n
      temp = x(ix)
      x(ix) = y(iy)
      y(iy) = temp
      ix = ix + incx
      iy = iy + incy
    end do

  end if

  return
end
subroutine svcal ( n, sa, x, incx, y, incy )
!
!*******************************************************************************
!
!! SVCAL sets a vector to a multiple of another vector.
!
!
!  Modified:
!
!    08 April 1999
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vectors.
!
!    Input, real SA, the multiplier.
!
!    Input, real X(*), the array to be scaled and copied into Y.
!
!    Input, integer INCX, the increment between successive entries of X.
!
!    Output, real Y(*), the scaled copy of X.
!
!    Input, integer INCY, the increment between successive elements of Y.
!
  implicit none
!
  integer i
  integer incx
  integer incy
  integer ix
  integer iy
  integer n
  real sa
  real x(*)
  real y(*)
!
  if ( n <= 0 ) then

  else if ( incx == 1 .and. incy == 1 ) then

    y(1:n) = sa * x(1:n)

  else

    if ( incx >= 0 ) then
      ix = 1
    else
      ix = ( - n + 1 ) * incx + 1
    end if

    if ( incy >= 0 ) then
      iy = 1
    else
      iy = ( - n + 1 ) * incy + 1
    end if

    do i = 1, n
      y(iy) = sa * x(ix)
      ix = ix + incx
      iy = iy + incy
    end do

  end if

  return
end
function svsame ( n, x, y )
!
!*******************************************************************************
!
!! SVSAME determines if two real vectors are equal.
!
!
!  Modified:
!
!    10 November 2000
!
!  Parameters:
!
!    Input, integer N, the number of entries in the vectors.
!
!    Input, real X(N), Y(N), the two vectors to compare.
!
!    Output, logical SVSAME, is TRUE if the vectors are identical,
!    and FALSE if at least one pair of entries differ.
!
  implicit none
!
  integer n
!
  logical svsame
  real x(n)
  real y(n)
!
  svsame = all ( x(1:n) == y(1:n) )

  return
end