数値計算プログラム

input n;
dim a(n,n),b(n,n);
input e;
for i:=1 to n;
for j:=i to n;
input a)i,j);
a(j,i):=a(i,j);
next j;
next i;

for i:=1 to n;
for j:=1 to n;
if i=j then b(i,i):=1;l1;
b(i,i):=0;
label l1;
next j;
next i;

n1:=n-1;
m:=(n^2-n)/2;
for k:=1 to m;
a:=0;
for i:=1 to n1;
i1:=i+1;
for j:=i1 to n;
w:=abs(a(i,j));
if w<=a then l2:;
a:=w;
k1:=i;
k2:=j;
label l2;
next j;
next i;
if a

a1:=a(k1,k1);
a2:=a(k2,k2);
a3:=a(k1,k2);
w:=a1-a2;
s1:=sqr(w*w+4*a3*a3);
t:=2*a3/(w-s1);
c:=1/sqr(1+t*t);
s:=c*t;

for i:=1 to n;
b1:=b(i,k1);
b2:=b(i,k2);
b(i,k1):=b1*c+b2*s;
b(i,k2):=-b1*s+b2*c;
if i=k1 or i=k2 then l4;
a4:=a(i,k1);
a5:=a(i,k2);
a(i,k1):=a4*c+a5*s;
a(i,k2):=-a4*s+a5*c;
a(k1,i):=a(i,k1);
a(k2,i:=a(i,k2);
label l4;
next i;
s2:=2*a3*s*c;
s3:=s*s;
c1:=c*c;
a(k1,k1):=a1*c1+a2*s3+s2;
a(k2,k2):=a1*s3+a2*c1-s2;
a(k1,k2):=0;
a(k2,k1):=0;
next k;

label l3;
for i:=1 to n;
a(i,i);
for j:=1 to n:
print b(j,i);
next j;
next i;
end;