with MathAKD; use MathAKD;
     with MatrAK ; use MatrAK ;

     package body DOS_010_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 (
       W  :    out Vector_Type;
       V  :    out Matrix_Type;
       A  : in out 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 Max ( A, B : float ) return float is
       begin
	 if A > B then
	   return A;
	 else
	   return B;
	 end if;
       end 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 I in A'range ( 2 ) loop
	 L := I + 1;
	 RV1 ( I ) := Scale * G;
	 G := 0.0;
	 S := 0.0;
	 Scale := 0.0;
	 if I <= M then
	   for K in I..M loop
	     Scale := Scale + abs ( A ( K, I ) );
	   end loop;
	   F := A ( I, I );
	   G := -Sign ( SqRt ( S ), F );
	   H := F * G - S;
	   A ( I, I ) := F - G;
	   if I /= N then
	     for J in A'range ( 2 ) loop
	       S := 0.0;
	       for K in I..M loop
		 S := S + A ( K, I ) * A ( K, J );
	       end loop;
	       F := S / H;
	       for K in I..M loop
		 A ( K, J ) := A ( K, J ) + F * A ( K, I );
	       end loop;
	     end loop;
	   end if;
	   W ( I ) := Scale * G;
	   G := 0.0;
	   S := 0.0;
	   Scale := 0.0;
	   if ( I <= M ) and ( I /= N ) then
	     for K in A'range ( 2 ) loop
	       Scale := Scale + abs ( A ( I, K ) );
	     end loop;
	     if Scale /= 0.0 then
	       for K in 1..N loop
		 A ( I, K ) := A ( I, K ) / Scale;
		 S := S + A ( I, K ) * A ( I, K );
	       end loop;
	       F := A ( I, L );
	       G := -Sign ( SqRt ( S ), F );
	       H := F * G - S;
	       A ( I, L ) := F - G;
	       for K in A'range ( 2 ) loop
		 RV1 ( K ) := A ( I, K ) / H;
	       end loop;
	       if I /= 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 DOS_010_Matrix_Solutions;