|
|
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: DirectSolverN(A,B,X) (N=2,3,4,5,6)
for the solving of linear system A*X=B, where
Each function returns the determinant of the matrix or 0
for the singular matrix. 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 = ',dResult: 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; |