Matrix-exponentiation operator: Difference between revisions
→{{header|Ada}}
m (→{{header|J}}) |
|||
Line 226:
begin
for K in RL'Range (1) loop
Sum := Sum + X (
end loop;
R (I, J) := Sum;
Line 237:
begin
for I in A'Range (1) loop
for J in A'Range (
Put (A (I, J));
end loop;
Line 251:
end Test_Matrix;</syntaxhighlight>
This solution is not tested, because the available version of GNAT GPL Ada compiler (20070405-41) does not provide an implementation of the standard library.
''(Another person is talking here:)'' I have made small corrections and tested this in 2023, and it did not work as I expected. However, I have questions about the mathematical libraries (to put it mildly). I tried both GCC 12 and GCC 13. What might be needed here is one's own eigensystem routine.
On the other hand, I did get a version working to raise a real matrix to a natural number power, thus demonstrating the correctness of the approach:
<syntaxhighlight lang="ada">
with Ada.Text_IO; use Ada.Text_IO;
with Ada.Float_Text_IO; use Ada.Float_Text_IO;
with Ada.Numerics.Real_Arrays; use Ada.Numerics.Real_Arrays;
procedure Test_Matrix is
procedure Put (A : Real_Matrix) is
begin
for I in A'Range (1) loop
for J in A'Range (2) loop
Put (" ");
Put (A (I, J));
end loop;
New_Line;
end loop;
end Put;
function "**" (A : Real_Matrix; Power : Integer) return Real_Matrix is
L : Real_Vector (A'Range (1));
X : Real_Matrix (A'Range (1), A'Range (2));
R : Real_Matrix (A'Range (1), A'Range (2));
RL : Real_Vector (A'Range (1));
begin
Eigensystem (A, L, X);
for I in L'Range loop
RL (I) := L (I) ** Power;
end loop;
for I in R'Range (1) loop
for J in R'Range (2) loop
declare
Sum : Float := 0.0;
begin
for K in RL'Range loop
Sum := Sum + X (I, K) * RL (K) * X (J, K);
end loop;
R (I, J) := Sum;
end;
end loop;
end loop;
return R;
end "**";
M : Real_Matrix (1..2, 1..2) := ((3.0, 2.0), (2.0, 1.0));
begin
Put_Line ("M ="); Put (M);
Put_Line ("M**0 ="); Put (M**0);
Put_Line ("M**1 ="); Put (M**1);
Put_Line ("M**2 ="); Put (M**2);
Put_Line ("M**3 ="); Put (M**3);
Put_Line ("M**50 ="); Put (M**50);
end Test_Matrix;
</syntaxhighlight>
{{out}}
<pre>
M =
3.00000E+00 2.00000E+00
2.00000E+00 1.00000E+00
M**0 =
1.00000E+00 0.00000E+00
0.00000E+00 1.00000E+00
M**1 =
3.00000E+00 2.00000E+00
2.00000E+00 1.00000E+00
M**2 =
1.30000E+01 8.00000E+00
8.00000E+00 5.00000E+00
M**3 =
5.50000E+01 3.40000E+01
3.40000E+01 2.10000E+01
M**50 =
1.61305E+31 9.96919E+30
9.96919E+30 6.16130E+30
</pre>
The method works because there are many products of the eigenvector matrices that cancel out.
=={{header|ALGOL 68}}==
| |||