JOINT   INSTITUTE   FOR   NUCLEAR   RESEARCH
lit
  Program   Library  JINRLIB rus

DirectSolverN (N=2,3,4,5,6) - general type linear system solving with usage
of a direct method (by Kramer's rule)

F099

Author: A.P.Sapojnikov Language: Fortran

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)! < N3.
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;


home up e-mail