Ada语言矩阵求逆代码格式问题求助：文件读取与计算故障排查

我是Ada语言新手，编写的矩阵求逆代码存在问题但无法定位原因。代码目标是读取a.txt和b.txt两个文件，计算矩阵的逆并求解线性方程组。

输入文件内容 #

A.txt #

2 -2 4 -2 2 0 
1 0 4 -3 0 0
1 1 1 1 1 1 
1 0 0 -1 0 0
3 0 0 -3 0 2
1 0 0 0 -2 4

B.txt #

13
18
51
82
32
12

我的代码 #

with Ada.Text_IO, Ada.Float_Text_IO, Ada.Integer_Text_IO;
use Ada.Text_IO, Ada.Float_Text_IO, Ada.Integer_Text_IO;

procedure Matrix_Inverse is

   type Matrix is array (1 .. 6, 1 .. 6) of Float;
   type Vector is array (1 .. 6) of Float;
   
   function Get_Matrix_From_File(Filename : String) return Matrix is
      A : Matrix := (others => (others => 0.0));
      File : File_Type;
      Line : String (1 .. 20);
      Last : Natural;
   begin
      Open (File, In_File, Filename);
      for I in A'Range(1) loop
         Get_Line (File, Line, Last);
         for J in A'Range(2) loop
            A (I, J) := Float'Value (Line (1 .. Last));
            Get (File, Line);
            Get_Line (File, Line, Last);
         end loop;
         A (I, 6) := Float'Value (Line (1 .. Last));
      end loop;
      Close (File);
      return A;
   end Get_Matrix_From_File;
   
   function Get_Vector_From_File(Filename : String) return Vector is
      B : Vector := (others => 0.0);
      File : File_Type;
   begin
      Open (File, In_File, Filename);
      for I in B'Range loop
         Get (File, B (I));
      end loop;
      Close (File);
      return B;
   end Get_Vector_From_File;
   
   procedure Print_Matrix(A : Matrix) is
   begin
      for I in A'Range(1) loop
         for J in A'Range(2) loop
            Put (A (I, J), 6, 2);
            Put (" ");
         end loop;
         New_Line;
      end loop;
   end Print_Matrix;
   
   function Inverse(A : Matrix) return Matrix is
      Inv : Matrix := (others => (others => 0.0));
      Identity : Matrix := (others => (others => 0.0));
   begin
      for I in Identity'Range(1) loop
         Identity (I, I) := 1.0;
      end loop;
      
      for I in A'Range(1) loop
         if A (I, I) = 0.0 then
            for J in I + 1 .. A'Range(1)'Last loop
               if A (J, I) /= 0.0 then
                  for K in A'Range(2) loop
                     A (I, K) := A (I, K) + A (J, K);
                     Identity (I, K) := Identity (I, K) + Identity (J, K);
                  end loop;
                  exit;
               end if;
            end loop;
         end if;
    for J in A'Range(1) loop
      if J /= I then
          declare
            Factor : constant Float := A (J, I) / A (I, I);
          begin
            for K in A'Range(2) loop
                A (J, K) := A (J, K) - Factor * A (I, K);
                Identity (J, K) := Identity (J, K) - Factor * Identity (I, K);
            end loop;
          end;
      end if;
    end loop;

    for I in A'Range(1) loop
      declare
          Diagonal : constant Float := A (I, I);
      begin
          for J in A'Range(2) loop
            Identity (I, J) := Identity (I, J) / Diagonal;
          end loop;
      end;
    end loop;

      return Identity;
    end Inverse;

  function Solve(A : Matrix; B : Vector) return Vector is
    Inv_A : Matrix := Inverse(A);
    Result : Vector := (others => 0.0);
  begin
    for I in Result'Range loop
        for J in Result'Range loop
          Result (I) := Result (I) + Inv_A (I, J) * B (J);
        end loop;
    end loop;
    return Result;
  end Solve;

  A : Matrix := Get_Matrix_From_File("a.txt");
  B : Vector := Get_Vector_From_File("b.txt");

  Solution : Vector := Solve(A, B);

  begin
    Put_Line("Solution:");
    for I in Solution'Range loop
        Put(Solution(I), 6, 2);
        Put(" ");
    end loop;
    New_Line;
end Matrix_Inverse;

代码问题及修正方案 #

1. 矩阵读取逻辑完全错误 #

Get_Matrix_From_File的读取逻辑混乱，错误使用Get_Line和Get导致矩阵值读取错误。简化为直接读取数值即可：

function Get_Matrix_From_File(Filename : String) return Matrix is
   A : Matrix := (others => (others => 0.0));
   File : File_Type;
begin
   Open (File, In_File, Filename);
   for I in A'Range(1) loop
      for J in A'Range(2) loop
         Get (File, A(I, J));
      end loop;
   end loop;
   Close (File);
   return A;
end Get_Matrix_From_File;

2. 函数输入参数只读冲突 #

Ada函数的输入参数默认是只读模式，Inverse函数中直接修改参数A会导致编译错误。需创建参数副本操作：

function Inverse(Original_A : Matrix) return Matrix is
   A : Matrix := Original_A; -- 创建副本
   Identity : Matrix := (others => (others => 0.0));
begin
   -- 后续逻辑使用A而非Original_A

3. 主元处理逻辑漏洞 #

• 主元为0时，仅将下方行加到当前行而非直接交换，会导致计算异常，应改为直接互换两行：

  if abs A(I, I) < 1.0e-6 then -- 用精度阈值替代直接等于0
     for J in I + 1 .. A'Range(1)'Last loop
        if abs A(J, I) > 1.0e-6 then
           -- 交换I行和J行
           declare
              Temp : Float;
           begin
              for K in A'Range(2) loop
                 Temp := A(I, K);
                 A(I, K) := A(J, K);
                 A(J, K) := Temp;
                 Temp := Identity(I, K);
                 Identity(I, K) := Identity(J, K);
                 Identity(J, K) := Temp;
              end loop;
              exit;
           end;
        end if;
     end loop;
  end if;
  

• 高斯-约当消元步骤顺序错误，应先完成所有行的消元，再进行主元归一化，或按列先归一化再消元。

4. 浮点数精度问题 #

直接用A(I,I) = 0.0判断主元会因为精度误差误判，改用极小阈值比较：abs A(I,I) < 1.0e-6。

5. 变量初始化顺序问题 #

Solution : Vector := Solve(A, B);放在声明区虽合法，但更规范的做法是移到begin块内，避免初始化顺序潜在问题。

内容的提问来源于stack exchange，提问作者lga2