with ConsAK; use ConsAK;
     with ErroAK; use ErroAK;
     with FileAK; use FileAK; -- for Load and Save
     with InteAK; use InteAK;
     with FloaAK; use FloaAK;
     with MathAK; use MathAK;
     with Text_IO; use Text_IO;

     package body A83_015_Matrix is
     ---------------------------------------------------------------------------
     ---------------------------------------------------------------------------

     function "*" (
       M : Matrix_Type;
       F : float )
       return Matrix_Type is
     ---------------------------------------------------------------------------
       MR : Matrix_Type ( M'range ( 1 ), M'range ( 2 ) ) := M;
     begin
       for Row in MR'range ( 1 ) loop
	 for Col in MR'range ( 2 ) loop
	   MR ( Row, Col ) := MR ( Row, Col ) * F;
	 end loop;
       end loop;
       return MR;
     end "*";

     function "*" (
       M1 : Matrix_Type;
       M2 : Matrix_Type )
       return Matrix_Type is
     ---------------------------------------------------------------------------
       R : Matrix_Type ( M1'range ( 1 ), M2'range ( 2 ) )
	 := ( others => ( others => 0.0 ) );
     begin
--     Put_Line ( M1 );
--     Pause;
--     Put_Line ( M2 );
--     Pause;
       if ( M1'last ( 2 ) - M1'first ( 2 ) + 1 )
	 /= ( M2'last ( 1 ) - M2'first ( 1 ) + 1 ) then
	   Put_Line ( "# of columns of 1st matrix ("
	     & integer'image ( M1'last ( 2 ) - M1'first ( 2 ) + 1 )
	     & ") must equal # of rows of 2nd matrix ("
	     & integer'image ( M2'last ( 1 ) - M2'first ( 1 ) + 1 )
	     & ")." );
	   Pause;
	   raise constraint_error;
       end if;
       for Row in R'range ( 1 ) loop
	 for Col in R'range ( 2 ) loop
	   R ( Row, Col ) := 0.0;
	   for R_Col in M1'range ( 2 ) loop
	     R ( Row, Col ) := R ( Row, Col )
	       + M1 ( Row, R_Col ) * M2 ( R_Col, Col );
	   end loop;
	 end loop;
       end loop;
       return R;
     end "*";

     function "*" (
       M : Matrix_Type;
       V : Vector_Type )
       return Vector_Type is
     ----------------------------------------------------------------------
     begin
       return Vector_Convert ( M * Matrix_Convert ( V ) );
     end "*";

     function "/" (
       M : Matrix_Type;
       F : float )
       return Matrix_Type is
     ---------------------------------------------------------------------------
       MR : Matrix_Type ( M'range ( 1 ), M'range ( 2 ) ) := M;
     begin
       for Row in MR'range ( 1 ) loop
	 for Col in MR'range ( 2 ) loop
	   MR ( Row, Col ) := MR ( Row, Col ) / F;
	 end loop;
       end loop;
       return MR;
     end "/";

     function "-" (
       M : Matrix_Type )
       return Matrix_Type is
     ----------------------------------------------------------------------
       Temp : Matrix_Type ( M'range ( 1 ), M'range ( 2 ) );
     begin
       for Row in M'range ( 1 ) loop
	 for Col in M'range ( 2 ) loop
	   Temp ( Row, Col ) := - M ( Row, Col );
	 end loop;
       end loop;
       return Temp;
     end "-";

     function "<=" (
       M1, M2 : Matrix_Type )
       return boolean is
     ---------------------------------------------------------------------------
     begin
       for Row in M1'range ( 1 ) loop
	 for Col in M1'range ( 2 ) loop
	   if M1 ( Row, Col ) > M2 ( Row, Col ) then
	     return false;
	   end if;
	 end loop;
       end loop;
       return true;
     end "<=";

     function Adjoint (
       M : Matrix_Type )
       return Matrix_Type is
     ---------------------------------------------------------------------------
       R : Matrix_Type ( M'range ( 1 ), M'range ( 2 ) );
     begin
       for Row in R'range ( 1 ) loop
	 for Col in R'range ( 2 ) loop
	   R ( Row, Col ) := Cofactor ( M, Row, Col );
	 end loop;
       end loop;
       return Transpose ( R );
     end Adjoint;

     function Ask (
       Row_Default : in positive := 2;
       Col_Default : in positive := 2 )
       return Matrix_Type is
     ----------------------------------------------------------------------
       Rows, Cols : positive;
       Good : boolean;
     begin
       loop
	 Good := true;
	 Rows := InteAK.Ask_Pos ( "How many rows are in your matrix?    ",
	   Row_Default, 1 );
	 Cols := InteAK.Ask_Pos ( "How many columns are in your matrix? ",
	   Col_Default, 1 );
	 declare
	   Matrix : Matrix_Type ( 1..Rows, 1..Cols );
	 begin
	   Put_Line (
	     "Please enter all values in floating-point format." );
	   Put_Line ( "For example, instead of entering 0, enter 0.0." );
	   for Row in 1..Rows loop
	     Put_Line ( "Please enter row" & positive'image ( Row )
	       & " of your matrix (example:  1.0 -3.4 1.0e-10)." );
	     for Col in 1..Cols loop
	       begin
		 Get ( Matrix ( Row, Col ) );
	       exception
		 when others =>
     --          Had to be changed when I added Load and Save.
     --          when Data_Error =>
		   Put_Line (
		     "Please enter all values in floating-point format." );
		   Put_Line (
		     "For example, instead of entering 0, enter 0.0." );
		   Pause;
		   Put_Line ( "" );
		   Good := false;
	       end;
	       exit when not Good;
	     end loop;
	     exit when not Good;
	   end loop;
	   if Good then
	     New_Line;
	     Put_Line ( Matrix );
	     Flush;
	     Ask ( Good, "Is this the matrix you meant? ", true );
	     if Good then
	       return Matrix;
	     end if;
	   end if;
	 end;
       end loop;
       exception
	 when others =>
	   Notify ( "function MatrAK.Ask return Matrix_Type" );
	   raise;
     end Ask;

     procedure Ask (
       M      : in out Matrix_Type;
       Prompt : in     string ) is
     ----------------------------------------------------------------------
       Good : boolean;
     begin
       loop
	 Good := true;
	 Put_Line ( Prompt );
	 for Row in M'range ( 1 ) loop
	   Put_Line ( "Please enter row" & integer'image ( Row ) & "." );
	   for Col in M'range ( 2 ) loop
	     begin
	       Get ( M ( Row, Col ) );
	     exception
	       when others =>
--             when Data_Error =>
		 Put_Line (
		   "Please enter all values in floating-point format." );
		 Put_Line (
		   "For example, instead of entering 0, enter 0.0." );
		 Pause;
		 Put_Line ( "" );
		 Good := false;
	     end;
	     exit when not Good;
	   end loop;
	   Flush;
	   exit when not Good;
	 end loop;
	 exit when Good;
       end loop;
     end Ask;

     function Cofactor (
       M        : Matrix_Type;
       Row, Col : integer )
       return float is
     ---------------------------------------------------------------------------
     begin
       return ( ( - 1.0 ) ** ( Row + Col ) ) * Minor ( M, Row, Col );
     exception
       when others =>
	 Notify ( "MatrAK.Cofactor" );
	 raise;
     end Cofactor;

     function Column (
       M   : Matrix_Type;
       Col : integer )
       return Vector_Type is
     ----------------------------------------------------------------------
       V : Vector_Type ( M'range ( 1 ) );
     begin
       for Row in V'range loop
	 V ( Row ) := M ( Row, Col );
       end loop;
       return V;
     end Column;

     function Concatenate (
       M1, M2 : Matrix_Type;
       Col    : boolean := true )
       return Matrix_Type is
     ---------------------------------------------------------------------------
     -- If input Col is true, then concatenates two matrices to make one bigger
     -- matrix with same number of rows but with the number of columns equal to
     -- sum of columns of the two inputs.
     --
     -- M1 := ( ( 1, 2 ),    M2 := ( ( 5, 6 ),
     --         ( 3, 4 ) );          ( 7, 8 ) );
     --
     -- M3 := Concatenate ( M1, M2, true ) = ( ( 1, 2, 5, 6 ),
     --                                        ( 3, 4, 7, 8 ) );
     --
     -- If input Col is false, then concatenates two matrices to make one bigger
     -- matrix with same number of columns but with the number of rows equal
     -- to sum of rows of the two inputs.
     --
     -- M1 := ( ( 1, 2 ),    M2 := ( ( 5, 6 ),
     --         ( 3, 4 ) );          ( 7, 8 ) );
     --
     -- M3 := Concatenate ( M1, M2, false ) = ( ( 1, 2 ),
     --                                         ( 3, 4 ),
     --                                         ( 5, 6 ),
     --                                         ( 7, 8 ) );
     ---------------------------------------------------------------------------
     begin
       if Col then
	 declare
	   M3 : Matrix_Type ( M1'range ( 1 ),
	     M1'first ( 2 )..( M1'last ( 2 ) + M2'last ( 2 ) ) );
	 begin
	   -- Rows of M1 and M2 must be equal.
	   if M1'last ( 1 ) /= M2'last ( 1 ) then
	     raise constraint_error;
	   end if;
	   for Row in M1'range ( 1 ) loop
	     for Col in M1'range ( 2 ) loop
	       M3 ( Row, Col ) := M1 ( Row, Col );
	     end loop;
	     for Col in ( M1'last ( 2 ) + 1 )..M3'last ( 2 ) loop
	       M3 ( Row, Col ) := M2 ( Row, Col - M1'last ( 2 ) );
	     end loop;
	   end loop;
	   return M3;
	 end;
       else
	 declare
	   M3 : Matrix_Type (
	     M1'first ( 1 )..( M1'last ( 1 ) + M2'last ( 1 ) ),
	     M1'range ( 2 ) );
	 begin
	   -- Columns of M1 and M2 must be equal.
	   if M1'last ( 2 ) /= M2'last ( 2 ) then
	     raise constraint_error;
	   end if;
	   for Col in M3'range ( 2 ) loop
	     for Row in M1'range ( 1 ) loop
	       M3 ( Row, Col ) := M1 ( Row, Col );
	     end loop;
	     for Row in ( M1'last ( 1 ) + 1 )..M3'last ( 1 ) loop
	       M3 ( Row, Col ) := M2 ( Row - M1'last ( 1 ), Col );
	     end loop;
	   end loop;
	   return M3;
	 end;
       end if;
     end Concatenate;

     function Concatenate (
       V1, V2 : Vector_Type )
       return Matrix_Type is
     ----------------------------------------------------------------------
     -- Example:  V1 = ( 1, 2, 3 ), V2 = ( 4, 5, 6 )
     --   returns ( ( 1, 4 ),
     --             ( 2, 5 ),
     --             ( 3, 6 ) )
     ----------------------------------------------------------------------
     begin
       return Concatenate ( Matrix_Convert ( V1 ), Matrix_Convert ( V2 ) );
     end Concatenate;

     procedure Demo is
     ---------------------------------------------------------------------------
       Default :          natural :=  1;
       Min     : constant natural :=  0;
       Max     : constant natural := 20;
       Option  :          natural;
     begin
       loop
	 Put_Line ( Copyright );
	 Put_Line ( Description );
	 New_Line;
	 Put_Line ( " 0 = Quit" );
	 Put_Line ( " 1 = ""*"" (multiplies a matrix by a float)" );
	 Put_Line ( " 2 = ""*"" (multiplies matrices)" );
	 Put_Line ( " 3 = ""/"" (divides matrix by float)" );
	 Put_Line ( " 4 = ""<=""" );
	 Put_Line ( " 5 = Adjoint" );
	 Put_Line ( " 6 = Ask" );
	 Put_Line ( " 7 = Cofactor" );
	 Put_Line ( " 8 = Concatenate (matrix)" );
	 Put_Line ( " 9 = Determinant" );
	 Put_Line ( "10 = Fill (matrix)" );
	 Put_Line ( "11 = Identity" );
	 Put_Line ( "12 = Inverse" );
	 Put_Line ( "13 = Is_Square" );
	 Put_Line ( "14 = Load" );
	 Put_Line ( "15 = Minor" );
	 Put_Line ( "16 = Put_Line" );
	 Put_Line ( "17 = Save" );
	 Put_Line ( "18 = Slice" );
	 Put_Line ( "19 = Sub_Matrix" );
	 Put_Line ( "20 = Transpose" );
	 New_Line;
	 Option := Ask_Nat ( "Option ", Default, Min, Max );
	 Default := Option + 1;
	 if Default > Max then
	   Default := Min;
	 end if;
	 case Option is
	   when  0 => exit;
	   when  1 => Put_Line ( A83_015_Matrix.Ask * Ask (
	     "Please enter a floating-point value: " ) );
	   when  2 =>
	     Put_Line ( "Warning:  please enter the second matrix first!" );
	     -- Don't know why this is necessary, but it is.
	     -- Possibly compiler dependent bug.
	     Pause;
	     Put_Line ( A83_015_Matrix.Ask * A83_015_Matrix.Ask );
	   when  3 => Put_Line ( A83_015_Matrix.Ask / Ask (
	     "Please enter a non-zero floating-point value: " ) );
	   when  4 => Put_Line ( A83_015_Matrix.Ask <= A83_015_Matrix.Ask );
	   when  5 => Put_Line ( Adjoint ( Ask ) );
	   when  6 => Put_Line ( A83_015_Matrix.Ask );
--         when  7 => cofactor
--         when  8 => concatenate
	   when  9 => Put_Line ( Determinant ( Ask ) );
	   when 10 => Put_Line ( Fill ( Ask ) ); -- need more
	   when 11 => Put_Line ( Identity ( Ask ) );
	   when 12 => Put_Line ( Inverse ( Ask ) );
	   when 13 => Put_Line ( Is_Square ( Ask ) );
	   when 14 =>
	     declare
	       Rows, Cols : positive := 10;
	       File_Name : string ( 1..12 ) := "Matrix.Mat  ";
	     begin
	       Rows := Ask_Pos ( "Rows    ", Rows );
	       Cols := Ask_Pos ( "Columns ", Cols );
	       declare
		 M : Matrix_Type ( 1..Rows, 1..Cols );
	       begin
		 Load ( M, File_Name );
		 Put_Line ( M );
	       end;
	     end;
     --    when 15 => Put_Line ( Minor ( Ask ) );
	   when 16 => Put_Line ( A83_015_Matrix.Ask );
	   when 17 =>
	     declare
	       File_Name : string ( 1..12 ) := "Matrix.Mat  ";
	     begin
	       Ask ( File_Name, "File Name ", File_Name );
	       Save ( Ask, File_Name );
	     end;
	   when 18 => Put_Line ( Slice ( Ask, Ask ) );
     --    when 19 => Put_Line ( Sub_Matrix ( Ask ) );
	   when 20 => Put_Line ( Transpose ( Ask ) );
	   when others =>
	     Put_Line (
	       "The demonstration of that option is currently not available." );
	 end case;
	 Pause;
	 New_Line;
       end loop;
     end Demo;

     function Determinant (
       M : Matrix_Type )
       return float is
     ---------------------------------------------------------------------------
       Sum : float := 0.0;
     begin
       if not Is_Square ( M ) then
	 raise constraint_error;
       end if;
       if M'length = 2 then
	 Sum := M ( M'first ( 1 ), M'first ( 2 ) )
	      * M ( M'last  ( 1 ), M'last  ( 2 ) )
	      - M ( M'last  ( 1 ), M'first ( 2 ) )
	      * M ( M'first ( 1 ), M'last  ( 2 ) );
       else
	 for Row in M'range loop
	   Sum := Sum + M ( Row, 1 ) * Cofactor ( M, Row, 1 );
	 end loop;
       end if;
       return Sum;
     end Determinant;

     function Fill (
       M : Matrix_Type;
       F : float := 0.0 )
       return Matrix_Type is
     ---------------------------------------------------------------------------
     -- Returns the matrix filled with F.
     ---------------------------------------------------------------------------
       R : Matrix_Type ( M'range ( 1 ), M'range ( 2 ) )
	 := ( others => ( others => F ) );
     begin
       return R;
     end Fill;

     function Identity (
       M : Matrix_Type )
       return Matrix_Type is
     ---------------------------------------------------------------------------
     -- Must be square.
     ---------------------------------------------------------------------------
       B : Matrix_Type ( M'range ( 1 ), M'range ( 2 ) )
	 := ( others => ( others => 0.0 ) );
     begin
       if not Is_Square ( M ) then
	 raise constraint_error;
       end if;
       for index in B'range loop
	 B ( index, index ) := 1.0;
       end loop;
       return B;
     end Identity;

     function Inverse (
       M : Matrix_Type )
       return Matrix_Type is
     ---------------------------------------------------------------------------
     -- M * Inverse ( M ) = I
     -- Inverse is not the transpose.
     ---------------------------------------------------------------------------
     begin
       return ( Adjoint ( M ) / Determinant ( M ) );
     end Inverse;

     function Is_Square (
       M : Matrix_Type )
       return boolean is
     ---------------------------------------------------------------------------
     begin
       if ( M'first ( 1 ) /= M'first ( 2 ) )
	 or ( M'last ( 1 ) /= M'last ( 2 ) ) then
	   return false;
       else
	 return true;
       end if;
     end Is_Square;

     procedure Load (
       M         :    out Matrix_Type;
       File_Name : in     string ) is
     ----------------------------------------------------------------------
       File : File_Type;
     begin
       M := ( others => ( others => 0.0 ) );
       Open ( File, In_File, File_Name );
       for Row in M'range ( 1 ) loop
	 for Col in M'range ( 2 ) loop
	   exit when End_Of_File ( File );
	   Get ( File, M ( Row, Col ) );
	 end loop;
       end loop;
       Close ( File );
     end Load;

     function Matrix_Convert (
       V : Vector_Type )
       return Matrix_Type is
     ----------------------------------------------------------------------
       Temp : Matrix_Type ( V'range, 1..1 );
     begin
       for Row in V'range loop
	 Temp ( Row, 1 ) := V ( Row );
       end loop;
       return Temp;
     end Matrix_Convert;

     function Minor (
       A        : Matrix_Type;
       Row, Col : integer )
       return float is
     ---------------------------------------------------------------------------
     -- "The determinants of the square submatrices of any matrix A, obtained
     -- by striking out certain rows or columns, or both, are called the
     -- determinants or minors of A."
     -- Eshbach, "Handbook of Engineering Fundamentals", 3rd Ed., 1975, p230.
     ---------------------------------------------------------------------------
     begin
       return Determinant ( Sub_Matrix ( A, Row, Col ) );
     end;

     procedure Put_Line (
       M : in Matrix_Type ) is
     ---------------------------------------------------------------------------
     begin
       for Row in M'range ( 1 ) loop
	 for Col in M'range ( 2 ) loop
	   Put ( M ( Row, Col ), 3, 2 );
	 end loop;
	 New_Line;
       end loop;
     end Put_Line;

     function Row (
       M   : Matrix_Type;
       Row : integer )
       return Vector_Type is
     ----------------------------------------------------------------------
       V : Vector_Type ( M'range ( 2 ) );
     begin
       for Col in V'range loop
	 V ( Col ) := M ( Row, Col );
       end loop;
       return V;
     end Row;

     procedure Save (
       M         : in     Matrix_Type;
       File_Name : in     string ) is
     ----------------------------------------------------------------------
       File : File_Type;
     begin
       Create ( File, Out_File, File_Name );
       for Row in M'range ( 1 ) loop
	 for Col in M'range ( 2 ) loop
	   Put ( File, M ( Row, Col ) );
	   Put ( File, ' ' );
	 end loop;
	 Put_Line ( File, "" );
       end loop;
       Close ( File );
     end Save;

     function Slice (
       Tiny, Huge : Matrix_Type )
       return Matrix_Type is
     ---------------------------------------------------------------------------
     -- Returns a matrix of the ranges and dimensions of Tiny from Huge.
     ---------------------------------------------------------------------------
       R : Matrix_Type ( Tiny'range ( 1 ), Tiny'range ( 2 ) );
     begin
       for Row in R'range ( 1 ) loop
	 for Col in R'range ( 2 ) loop
	   R ( Row, Col ) := Huge ( Row, Col );
	 end loop;
       end loop;
       return R;
     end Slice;

     function Sub_Matrix (
       A        : Matrix_Type;
       Row, Col : integer )
       return Matrix_Type is
     ---------------------------------------------------------------------------
     -- Returns the matrix generated by striking the row and column that
     -- contains the element ( Row, Col ).
     ---------------------------------------------------------------------------
       S : Matrix_Type ( A'first ( 1 )..( A'last ( 1 ) - 1 ),
			 A'first ( 2 )..( A'last ( 2 ) - 1 ) );
     begin
       for R in S'range ( 1 ) loop
	 for C in S'range ( 2 ) loop
	   S ( R , C ) := A ( R + boolean'pos ( R >= Row ),
	     C + boolean'pos ( C >= Col ) );
	 end loop;
       end loop;
       return S;
     end Sub_Matrix;

     function Transpose (
       M : Matrix_Type )
       return Matrix_Type is
     ---------------------------------------------------------------------------
       B : Matrix_Type ( M'range ( 2 ), M'range ( 1 ) );
     begin
       for Row in B'range ( 1 ) loop
	 for Col in B'range ( 2 ) loop
	   B ( Row, Col ) := M ( Col, Row );
	 end loop;
       end loop;
       return B;
     end Transpose;

     function Vector_Convert (
       M : Matrix_Type )
       return Vector_Type is
     ----------------------------------------------------------------------
       Temp : Vector_Type ( M'range ( 1 ) );
     begin
       if M'first ( 2 ) /= M'last ( 2 ) then
	 raise Constraint_Error;
       end if;
       for Row in M'range ( 1 ) loop
	 Temp ( Row ) := M ( Row, M'first ( 2 ) );
       end loop;
       return Temp;
     end Vector_Convert;

     ---------------------------------------------------------------------------
     ---------------------------------------------------------------------------
     end A83_015_Matrix;