{Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730 }

function ISAMAX (      n : integer ;
                  var sx : real ;
                    incx : integer ) : integer ;
     {                                                     }
     { Find index of element having maximum absolute value.}
     {                                                     }
     { Jack Dongarra, LINPACK, 3/11/78.                    }
     { Adam Fritz, TURBO Pascal, 2/22/85.                  }
     {                                                     }

var
     smax           : real ;
     i, ix          : integer ;
     x              : RowPointer ;

begin
   ISAMAX := 0 ;
   if n > 0 then begin
      ISAMAX := 1 ;
      if n > 1 then begin
         x := Ptr(Addr(sx)) ;
         if incx > 1 then begin
                                { incx > 1 }
            ix := 1 ;
            smax := Abs(x^.s[1]) ;
            ix := ix + incx ;
            for i := 2 to n do begin
               if Abs(x^.s[ix]) > smax then begin
                  ISAMAX := i ;
                  smax := Abs(x^.s[ix])
               end ;
               ix := ix + incx
            end
         end
                                { incx = 1 }
         else begin
            smax := Abs(x^.s[1]) ;
            for i := 2 to n do
               if Abs(x^.s[i]) > smax then begin
                  ISAMAX := i ;
                  smax := Abs(x^.s[i])
               end
         end
      end
   end
end ;

{~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~}

function SASUM (      n : integer ;
                 var sx : real ;
                   incx : integer ) : real ;
     {                                                     }
     { Forms sum of absolute values.                       }
     {                                                     }
     { Jack Dongarra, LINPACK, 3/11/78.                    }
     { Adam Fritz, TURBO Pascal, 2/22/85.                  }
     {                                                     }

var
     stemp          : real ;
     i, ix          : integer ;
     x              : RowPointer ;

begin
   stemp := 0.0 ;
   if n > 0 then begin
      x := Ptr(Addr(sx)) ;
      if incx > 1 then begin
                                { incx > 1 }
         ix := 1 ;
         for i := 1 to n do begin
            stemp := stemp + Abs(x^.s[ix]) ;
            ix := ix + incx
         end
      end
                                { incx = 1 }
      else
         for i := 1 to n do
            stemp := stemp + Abs(x^.s[i])
   end ;

   SASUM := stemp

end ;

{~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~}

procedure SAXPY (      n : integer ;
                      sa : real ;
                  var sx : real  ;
                    incx : integer ;
                  var sy : real  ;
                    incy : integer ) ;
     {                                                     }
     { Compute constant times a vector plus a vector.      }
     {                                                     }
     { Jack Dongarra, LINPACK, 3/11/78.                    }
     { Adam Fritz, TURBO Pascal, 2/22/85.                  }
     {                                                     }

var
     i, ix, iy      : integer ;
     x, y           : RowPointer ;

begin
   if n > 0 then begin
      if sa <> 0.0 then begin
         x := Ptr(Addr(sx)) ;
         y := Ptr(Addr(sy)) ;
         if (incx <> 1) or (incy <> 1) then begin
                                { incx <> incy or incx <> 1 }
            ix := 1 ;
            iy := 1 ;
            if incx < 0 then
               ix := (-n + 1) * incx + 1 ;
            if incy < 0 then
               iy := (-n + 1) * incy + 1 ;
            for i := 1 to n do begin
               y^.s[iy] := y^.s[iy] + sa * x^.s[ix] ;
               ix := ix + incx ;
               iy := iy + incy
            end
         end
                                      { incx, incy = 1 }
         else
            for i := 1 to n do
               y^.s[i] := y^.s[i] + sa * x^.s[i]
      end
   end
end ;

{~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~}

procedure SCOPY (      n : integer ;
                  var sx : real ;
                    incx : integer ;
                  var sy : real ;
                    incy : integer ) ;
     {                                                     }
     { Copies a vector, x, to a vector, y.                 }
     {                                                     }
     { Jack Dongarra, LINPACK, 3/11/78.                    }
     { Adam Fritz, TURBO Pascal, 2/22/85.                  }
     {                                                     }

var
     i, ix, iy      : integer ;
     x, y           : RowPointer ;

begin
   if n > 0 then begin
      x := Ptr(Addr(sx)) ;
      y := Ptr(Addr(sy)) ;
      if (incx <> 1) or (incy <> 1) then begin
                                { incx, incy  <> 1 }
         ix := 1 ;
         iy := 1 ;
         if incx < 0 then
            ix := (-n + 1) * incx + 1 ;
         if incy < 0 then
            iy := (-n + 1) * incy + 1 ;
         for i := 1 to n do begin
            y^.s[iy] := x^.s[ix] ;
            ix := ix + incx ;
            iy := iy + incy
         end
      end
                                { incx, incx = 1 }
      else
         for i := 1 to n do
            y^.s[i] := x^.s[i]
   end
end ;

{~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~}

function SDOT (      n : integer ;
                var sx : real ;
                  incx : integer ;
                var sy : real ;
                  incy : integer ) : real ;
     {                                                     }
     { Computes dot product of two vectors.                }
     {                                                     }
     { Jack Dongarra, LINPACK, 3/11/78.                    }
     { Adam Fritz, TURBO Pascal, 2/22/85.                  }
     {                                                     }

var
     stemp          : real ;
     i, ix, iy      : integer ;
     x, y           : RowPointer ;

begin
   stemp := 0.0 ;
   if n > 0 then begin
      x := Ptr(Addr(sx)) ;
      y := Ptr(Addr(sy)) ;
      if (incx <> 1) or (incy <> 1) then begin
                                { incx <> incy or incx <> 1 }
         ix := 1 ;
         iy := 1 ;
         if incx < 0 then
            ix := (-n + 1) * incx + 1 ;
         if incy < 0 then
            iy := (-n + 1) * incy + 1 ;
         for i := 1 to n do begin
            stemp := stemp + x^.s[ix] * y^.s[iy] ;
            ix := ix + incx ;
            iy := iy + incy
         end
      end
                                { incx, incy = 1 }
      else
         for i := 1 to n do
            stemp := stemp + x^.s[i] * y^.s[i]
   end ;
   SDOT := stemp
end ;

{~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~}

procedure SSCAL (      n : integer ;
                      sa : real ;
                  var sx : real ;
                    incx : integer ) ;
     {                                                     }
     { Computes vector scaled by a constant.               }
     {                                                     }
     { Jack Dongarra, LINPACK, 3/11/78.                    }
     { Adam Fritz, TURBO Pascal, 2/22/85.                  }
     {                                                     }

var
     i, ix          : integer ;
     x              : RowPointer ;

begin
   if n > 0 then begin
      x := Ptr(Addr(sx)) ;
      if incx <> 1 then begin
                                { incx <> 1 }
         ix := 1 ;
         for i := 1 to n do begin
            x^.s[ix] := sa * x^.s[ix] ;
            ix := ix + incx
         end
      end
                                { incx = 1 }
      else
         for i := 1 to n do
            x^.s[i] := sa * x^.s[i]
   end
end ;

{~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~}

procedure SSWAP (      n : integer ;
                  var sx : real ;
                    incx : integer ;
                  var sy : real ;
                    incy : integer ) ;
     {                                                     }
     { Interchanges two vectors.                           }
     {                                                     }
     { Jack Dongarra, LINPACK, 3/11/78.                    }
     { Adam Fritz, TURBO Pascal, 2/22/85.                  }
     {                                                     }
var
     stemp          : real ;
     i, ix, iy      : integer ;
     x, y           : RowPointer ;

begin
   if n > 0 then begin
      x := Ptr(Addr(sx)) ;
      y := Ptr(Addr(sy)) ;
      if (incx <> 1) or (incy <> 1) then begin
                                { incx <> incy or incx <> 1 }
         ix := 1 ;
         iy := 1 ;
         if incx < 0 then
            ix := (-n + 1) * incx + 1 ;
         if incy < 0 then
            iy := (-n + 1) * incy + 1 ;
         for i := 1 to n do begin
            stemp := x^.s[ix] ;
            x^.s[ix] := y^.s[iy] ;
            y^.s[iy] := stemp ;
            ix := ix + incx ;
            iy := iy + incy
         end
      end
                                { incx, incy = 1 }
      else for i := 1 to n do begin
         stemp := x^.s[i] ;
         x^.s[i] := y^.s[i] ;
         y^.s[i] := stemp
      end
   end
end ;

{Copyright (C) 1985 Adam Fritz, 133 Main St., Afton, NY 13730 }