fpc/packages/numlib/tests/eigbs3te.pas
2016-03-05 19:16:54 +00:00

81 lines
2.0 KiB
ObjectPascal

program eigbs3te;
uses
typ,
iom,
eig,
omv;
const
n1 = -100;
n2 = 100;
i1 = -10;
i2 = 10;
rwx = i2 - i1 + 1;
var
ex, nex, nel, ind, p, q, r, s, n, i, j, b, term: ArbInt;
a: array[n1..n2] of ArbFloat;
lam: array[i1..i2] of ArbFloat;
x, mat, e: array[i1..i2, i1..i2] of ArbFloat;
begin
Write(' program results eigbs3te');
case sizeof(ArbFloat) of
4: writeln('(single)');
6: writeln('(real)');
8: writeln('(double)');
end;
Read(nex);
writeln;
writeln('number of examples', nex: 2);
writeln;
for ex := 1 to nex do
begin
writeln('example number', ex: 2);
writeln;
Read(p, q, r, s, n, b);
nel := n * (b + 1) - (b * (b + 1)) div 2;
iomrev(input, a[p], nel);
eigbs3(a[p], n, b, lam[q], x[r, s], rwx, term);
writeln(' A = ');
iomwrv(output, a[p], nel, numdig);
writeln;
writeln('term=', term: 2);
if term = 1 then
begin
writeln;
writeln('lambda=');
iomwrv(output, lam[q], n, numdig);
writeln;
writeln('X=');
iomwrm(output, x[r, s], n, n, rwx, numdig);
ind := p;
for i := 1 to n do
for j := 1 to i do
if j < i - b then
mat[i + r - 1, j + s - 1] := 0
else
begin
mat[i + r - 1, j + s - 1] := a[ind];
ind := ind + 1;
end;
for i := 1 to n do
for j := i + 1 to n do
mat[i + r - 1, j + s - 1] := mat[j + r - 1, i + s - 1];
writeln;
writeln(' matrix A =');
iomwrm(output, mat[r, s], n, n, rwx, numdig);
writeln;
writeln('Ax-lambda.x = ');
omvmmm(mat[r, s], n, n, rwx, x[r, s], n, rwx, e[r, s], rwx);
for j := 1 to n do
for i := 1 to n do
e[i + r - 1, j + s - 1] := e[i + r - 1, j + s - 1] - lam[q + j - 1] * x[i + r - 1, j + s - 1];
iomwrm(output, e[r, s], n, n, rwx, numdig);
end;
writeln;
writeln('-------------------------------------------');
end;
Close(input);
Close(output);
end.