--   with MathAKD; use MathAKD;

     package body A83_016_Matrix_Solutions is
     ----------------------------------------------------------------------
     ----------------------------------------------------------------------

     function Least_Squares (
       A :  Matrix_Type;
       B :  Vector_Type )
       return Vector_Type is
     ----------------------------------------------------------------------
     begin
       return ( Inverse ( Transpose ( A ) * A ) * Transpose ( A ) ) * B;
     end Least_Squares;

--   procedure Singular_Value_Decomposition (
--     U  :     out Matrix_Type;
--     W  :     out Vector_Type;
--     V  :     out Matrix_Type;
--     A  : in      Matrix_Type ) is
--   ----------------------------------------------------------------------
--   -- Translated from "svdcmp" of _Numerical Recipes_
--   ----------------------------------------------------------------------
--     M : integer := A'last ( 1 );
--     N : integer := A'last ( 2 );
--     N_Max : constant := 100;
--     NM, L, K, JJ, ITS : integer;
--     Z, Y, X, Scale, S, H, G, F, C, A_Norm : float;
--     RV1 :  Vector_Type ( 1..N_Max );
--
--     function Sign ( A, B : float ) return float is
--     begin
--       if B >= 0.0 then
--         return abs ( A );
--       else
--         return -abs ( A );
--       end if;
--     end Sign;
--
--   begin
--     G := 0.0;
--     Scale := 0.0;
--     A_Norm := 0.0;
--     for index in A'range ( 2 ) loop
--       L := index + 1;
--       RV1 ( index ) := Scale * G;
--       G := 0.0;
--       S := 0.0;
--       Scale := 0.0;
--       if index <= M then
--         for K in index..M loop
--           Scale := Scale + abs ( A ( K, index ) );
--         end loop;
--         F := A ( index, index );
--         G := -Sign ( SqRt ( S ), F );
--         H := F * G - S;
--         A ( index, index ) := F - G;
--         if index /= N then
--           for J in A'range ( 2 ) loop
--             S := 0.0;
--             for K := index to M loop
--               S := S + A ( K, index ) * A ( K, J );
--             end loop;
--             F := S / H;
--             for K in index to M loop
--               A ( K, J ) := A ( K, J ) + F * A ( K, index );
--             end loop;
--           end loop;
--         end if;
--         W ( index ) := Scale * G;
--         G := 0.0;
--         S := 0.0;
--         Scale := 0.0;
--         if ( index <= M ) and ( index /= N ) then
--           for K in A'range ( 2 ) loop
--             Scale := Scale + abs ( A ( index, K ) );
--           end loop;
--           if Scale /= 0.0 then
--             for K := 1 to N loop
--               A ( index, K ) := A ( index, K ) / Scale;
--               S := S + A ( index, K ) * A ( index, K );
--             end loop;
--             F := A ( index, L );
--             G := -Sign ( SqRt ( S ), F );
--             H := F * G - S;
--             A ( index, L ) := F - G;
--             for K in A'range ( 2 ) loop
--               RV1 ( K ) := A ( index, K ) / H;
--             end loop;
--             if index /= M then
--               for J in A'range ( 1 ) loop
--                 S := 0.0;
--                 for K in A'range ( 2 ) loop
--                   S := S + A ( J, K ) * A ( I, K );
--                 end loop;
--                 for K in A'range ( 2 ) loop
--                   A ( J, K ) := A ( J, K ) + S * RV1 ( K );
--                 end loop;
--               end loop;
--             end if;
--             for K in A'range loop
--               A ( I, K ) := Scale * A ( I, K );
--             end loop;
--           end if;
--         end if;
--         A_Norm := Max ( A_Norm, abs ( W ( I ) ) + abs ( RV1 ( I ) ) );
--       end if;
--       for I in reverse A'range ( 2 ) loop
--         if I < N then
--           if G /= 0.0 then
--             for J in A'range ( 2 ) loop
--               V ( J, I ) := ( A ( I, J ) / A ( I, L ) ) / G;
--             end loop;
--             for J in A'range ( 2 ) loop
--               S := 0.0;
--               for K in A'range ( 2 ) loop
--                 S := S + A ( I, K ) * V ( K, J );
--               end loop;
--               for K in A'range ( 2 ) loop
--                 V ( K, J ) := V ( K, J ) + S * V ( K, I );
--               end loop;
--             end loop;
--           end if;
--           for J in A'range ( 2 ) loop
--             V ( I, J ) := 0.0;
--             V ( J, I ) := 0.0;
--           end loop;
--         end if;
--         V ( I, I ) := 1.0;
--         G := RV1 ( I );
--         L := I;
--       end loop;
--       for I in reverse A'range ( 2 ) loop
--         L := I + 1;
--         G := W ( I );
--         if I < N then
--           for J in L..N loop
--             A ( I, J ) := 0.0;
--           end loop;
--         end if;
--         if G /= 0.0 then
--           G := 1.0 / G;
--           if I /= N then
--             for J in L..N loop
--               S := 0.0;
--               for K in L..M loop
--                 S := S + A ( K, I ) * A ( K, J );
--               end loop;
--               F := ( S / A ( I, I ) ) * G;
--               for K in I..M loop
--                 A ( K, J ) := A ( K, J ) + F * A ( K, I );
--               end loop;
--             end loop;
--           end if;
--           for J in I..M loop
--             A ( J, I ) := A ( J, I ) * G;
--           end loop;
--         else
--           for J in I..M loop
--             A ( J, I ) := 0.0;
--           end loop;
--         end if;
--         A ( I, I ) := A ( I, I ) + 1.0;
--       end loop;
--
--     end loop;
--
--   end Singular_Value_Decomposition;

     ----------------------------------------------------------------------
     ----------------------------------------------------------------------
     end A83_016_Matrix_Solutions;