|
DirectSolverN Library "JINRLIB" F099 
(N=2,3,4,5,6)
Author: A.P.Sapojnikov
Language: Fortran
GENERAL TYPE LINEAR SYSTEM SOLVING WITH USAGE
OF A DIRECT METHOD (BY KRAMER'S RULE)
The function DirectSolverN (N=2,3,4,5,6) solves a general linear
system A*X = B using a direct method, i.e. by Kramer's rule:
X(i) = DX(i) / D .
It is well-known that for the solving of the system of the N order
by methods of Gauss-Gordan type, O (N ** 3) operations are required.
At the same time, the direct use of Kramer's rule requires to
calculate N+1 determinant, i.e. (N+1)! operations.
However, at small N both estimations are quite commensurable, and
at N < 4 simply (N+1)! < N ** 3.
Moreover, it was experimentally found that down to the 6-th order,
the system is solved by a direct method much faster (at N=6 -
twice faster).
Therefore, we present the family of the functions DirectSolverN (A, B, X)
(N=2,3,4,5,6) for the solving of linear system A*X=B,
where A - a system matrix NxN; B - an array of length N;
X - an array of length N, where the solution will be written.
Each function returns the determinant of the matrix or 0
for the singular matrix.
The feature of these programs is that they were generated
automatically with usage of a specially constructed recursive
Pascal-written procedure DeterGen.
Indeed, who wants to "write away" manually the determinant of the
6-th order! And it is so easy to make a mistake here. The code
of the mentioned procedure is published below.
Notes:
------
Because of some compiler restrictions the maximal value
of system dimension N for Microsoft Fortran 5.00 is equal to 5.
Example:
--------
Implicit Real*8 (a-h,o-z)
dimension a(4,4),b(4),x(4)
do i=1,4
b(i)=0.0d0
do j=1,4
a(i,j)=1.0d0/(i+j) ! Hilbert matrix
b(i)=b(i)+a(i,j) ! solution = all 1
enddo
enddo
d=DirectSolver4(a,b,x)
write(*,*) ' Determinant for Hilbert matrix = ',d
Result:
---------
Determinant for Hilbert matrix = 2.362055933475061E-009
Procedure DeterGen:
-------------------
Procedure DeterGen(N:integer; Lines:TStrings);
Type Matrix=array[1..9,1..9] of string[6];
var i,j,k,N0:integer;
a:Matrix; Indent:string;
Function Dig(d:integer):char;
begin Result:=Chr(d+Ord('0')); end;
Procedure Determ(n:integer; a:Matrix; Lines:TStrings);
var s,s1,sg:string; i,j,k:integer;
b:Matrix;
begin sg:=''; s:='D'+Dig(n);
if (n=2) and (N0=2) then begin Lines.Add(Indent+'D2 = '+a[1][1]+'*'+a[2][2]+'-'+a[1][2]+'*'+a[2][1]); Exit; end;
if n=3 then
begin s1:=Indent+' '; s1[6]:='&';
Lines.Add(Indent+'D3 = '+a[1][1]+'*('+a[2][2]+'*'+a[3][3]+'-'+a[3][2]+'*'+a[2][3]+')-'+a[1][2]+'*('+a[2][1]+'*'+a[3][3]+'-');
Lines.Add(s1+a[2][3]+'*'+a[3][1]+')+'+a[1][3]+'*('+a[2][1]+'*'+a[3][2]+'-'+a[2][2]+'*'+a[3][1]+')');
Exit;
end;
for j:=1 to n do
begin
for i:=1 to n-1 do for k:=1 to n do b[i][k]:=a[i+1][k];
for i:=1 to n-1 do for k:=j to n-1 do b[i][k]:=b[i][k+1];
Determ(n-1,b,Lines);
s1:=Indent+s+' = '; if j<>1 then s1:=s1+s+sg;
Lines.Add(s1+a[1][j]+' * D'+Dig(n-1));
if sg=' - ' then sg:=' + ' else sg:=' - ';
end;
end;
begin N0:=N; Indent:=' ';
Lines.Clear; if n>9 then n:=9;
Lines.Add(' Function DirectSolver'+Dig(n)+'(A,B,X) ! for Linear System A*X=B');
Lines.Add('C *** Generated by "DeterGen" procedure ***');
Lines.Add(' Implicit Real*8 (A-H,O-Z)');
Lines.Add(' Dimension A('+Dig(n)+','+Dig(n)+') ! system Matrix');
Lines.Add(' Dimension B('+Dig(n)+') ! system Right Part');
Lines.Add(' Dimension X('+Dig(n)+') ! system Solution');
for i:=1 to n do for j:=1 to n do a[i][j]:='A('+Dig(i)+','+Dig(j)+')';
Determ(n,a,Lines);
Lines.Add(' D = D'+Dig(n));
Lines.Add(' if(D.ne.0) then');
for k:=1 to n do
begin Indent:=' ';
for i:=1 to n do a[i][k]:='B('+Dig(i)+')';
Determ(n,a,Lines);
Lines.Add(Indent+'X('+Dig(k)+') = D'+Dig(n)+' / D');
for i:=1 to n do a[i][k]:='A('+Dig(i)+','+Dig(k)+')';
end;
Lines.Add(' endif');
Lines.Add(' DirectSolver'+Dig(n)+'=D');
Lines.Add(' return');
Lines.Add(' End');
end;
|
