Procedure written on C# is in the end of article.
Program result you will find in text file "Result".
There are couple procedures like "readFile" that read matrix from text file. In files there are such matrix:
1:
1 -4 -3 8
2 -6 -4 2
3 -8 1 -4
4 2 2 -6
2:
0 1 -6 -4
3 -1 -6 -4
2 3 9 2
3 2 3 8
3:
-3 -3 4
0 -1 0
6 2 -4
4:
-3 -3 2
0 -1 0
6 2 -4
So here it is...
program Project1;
{$APPTYPE CONSOLE}
uses
SysUtils;
type TMatrix = array [1..50,1..50] of real;
TVect = array [1..50] of real;
TFile = textFile;
var A: TMatrix ; b:TVect; n, per:integer; t,f:TFile; c:byte; s:string;
procedure readFile(var A:TMatrix; var n:integer);
var i,j:integer;
begin i:=1; j:=1;
while not eof(t) do
begin j:=1;
while not eoln(t) do
begin
read(t,a[i,j]);
inc(j);
end;
readln(t);
inc(i);
end;
closeFile(t);
end;
procedure setMatrix(var A:TMatrix; var n:integer);
var i,j:integer;
begin
writeln('Input Matrix A');
for i := 1 to n do
for j := 1 to n do
begin
write ('a',i,',', j ,'=');
readln(A[i,j]);
end;
end;
procedure setVect(var b:TVect; var n:integer);
var i:integer;
begin
for i := 1 to n do
begin
write ('b',i,'=');
readln(b[i]);
end;
end;
procedure Gauss (A:TMatrix;b:TVect; n:integer);
var
i, j, k,step:integer; R,det,p:real;
x:TVect;
begin
per:=1;
for i := 1 to N - 1 do
begin
k := i;
R := Abs(A[i, i]);
for j := i + 1 to N do
if (Abs(A[j, i]) >= R) then
begin
k := j;
R := Abs(A[j, i]);
end;
if (k <> i) then
begin
R := B[k];
B[k] := B[i];
B[i] := R;
for j := i to N do
begin
R := A[k, j];
A[k, j] := A[i, j];
A[i, j] := R;
end;
end;
R := A[i, i];
if R = 0 then begin writeln('Denom <0'); readln; halt; end;
B[i] := B[i] / R;
for j := 1 to N do
begin
if R = 0 then begin writeln('Denom <0'); readln; halt; end;
A[i, j] := A[i, j] / R;
end;
for k := i + 1 to N do
begin
R := A[k, i];
B[k] := B[k] - R * B[i];
A[k, i] := 0;
for j := i + 1 to N do
A[k, j] := A[k, j] - R * A[i, j];
end;
end;
if A[N,N] = 0 then begin writeln('Denom <0');readln; halt; end;
X[N] := B[N] / A[N, N];
for i := N - 1 downto 1 do
begin
R := B[i];
for j := i + 1 to N do
R := R - A[i, j] * X[j];
X[i] := R;
per:=per+1;
end;
for i:= 1 to n do
writeln('x',i,'= ',x[i]:5:2);
end;
procedure showMatrix(A:TMatrix; n:integer);
var i,j:integer;
begin
for i := 1 to n do
begin
for j := 1 to n do
write(A[i,j]:4:2,' ');
writeln;
end;
end;
procedure Determinant;
var i,step,k:integer;
p,det:real;
begin
p:=1;
step:=1;
for i:=1 to per do
step:=step*(-1);
for k:=1 to n do
begin
p:=p * a[k,k];
det:=p * step;
end;
writeln('Determinant = ', det:5:2);
end;
procedure writeFile(x:TVect);
var i:integer;
begin
for i := 1 to n do
write(f,x[i],' ');
end;
begin
writeln('Which matrix do you want to calculate');
readln(c);
case c of
1:
begin
s:='Matrix1.txt';
n:=4;
end;
2:
begin
s:='Matrix2.txt';
n:=4
end;
3:
begin
s:='Matrix3.txt';
n:=3;
end;
4:
begin
s:='Matrix4.txt';
n:=3;
end;
end;
assignFile(t,s);
reset(t);
//setMatrix(A,n);
readFile(A,n);
showMatrix(A,n);
setVect(B,n);
AssignFile(f,'Result.txt');
Rewrite(f);
Gauss(a,b,n);
determinant;
writeFile(b);
CloseFile(f);
writeln('Results are in text file "Results"');
readln
end.
//==============================================
//Procedure written on C#
static public void Gauss (ref double[,] A, ref double[] B, int N,
ref double[] X)
{
int i, j, k;
double R;
for (i = 1; i <= N - 1; i++)
{
k = i;
R = Math.Abs(A[i, i]);
for (j = i + 1; j <= N; j++)
if (Math.Abs(A[j, i]) >= R)
{
k = j;
R = Math.Abs(A[j, i]);
}
if (k != i)
{
R = B[k];
B[k] = B[i];
B[i] = R;
for (j = i; j < N; j++)
{
R = A[k, j];
A[k, j] = A[i, j];
A[i, j] = R;
}
}
R = A[i, i];
B[i] = B[i] / R;
for (j = 1; j <= N; j++)
A[i, j] = A[i, j] / R;
for (k = i + 1; k <= N; k++)
{
R = A[k, i];
B[k] = B[k] - R * B[i];
A[k, i] = 0;
for (j = i + 1; j <= N; j++)
A[k, j] = A[k, j] - R * A[i, j];
}
}
X[N] = B[N] / A[N, N];
for (i = N - 1; i >= 1; i--)
{
R = B[i];
for (j = i + 1; j <= N; j++)
R = R - A[i, j] * X[j];
X[i] = R;
}
}
}
// program written on C++ from dreamincode.net
/******************************************************************************/
/* Perform Gauss-Jordan elimination with row-pivoting to obtain the solution to
* the system of linear equations
* A X = B
*
* Arguments:
* lhs - left-hand side of the equation, matrix A
* rhs - right-hand side of the equation, matrix B
* nrows - number of rows in the arrays lhs and rhs
* ncolsrhs- number of columns in the array rhs
*
* The function uses Gauss-Jordan elimination with pivoting. The solution X to
* the linear system winds up stored in the array rhs; create a copy to pass to
* the function if you wish to retain the original RHS array.
*
* Passing the identity matrix as the rhs argument results in the inverse of
* matrix A, if it exists.
*
* No library or header dependencies, but requires the function swaprows, which
* is included here.
*/
// swaprows - exchanges the contents of row0 and row1 in a 2d array
void swaprows(double** arr, long row0, long row1) {
double* temp;
temp=arr[row0];
arr[row0]=arr[row1];
arr[row1]=temp;
}
// gjelim
void gjelim(double** lhs, double** rhs, long nrows, long ncolsrhs) {
// augment lhs array with rhs array and store in arr2
double** arr2=new double*[nrows];
for (long row=0; row<nrows; ++row)
arr2[row]=new double[nrows+ncolsrhs];
for (long row=0; row<nrows; ++row) {
for (long col=0; col<nrows; ++col) {
arr2[row][col]=lhs[row][col];
}
for (long col=nrows; col<nrows+ncolsrhs; ++col) {
arr2[row][col]=rhs[row][col-nrows];
}
}
// perform forward elimination to get arr2 in row-echelon form
for (long dindex=0; dindex<nrows; ++dindex) {
// run along diagonal, swapping rows to move zeros in working position
// (along the diagonal) downwards
if ( (dindex==(nrows-1)) && (arr2[dindex][dindex]==0)) {
return; // no solution
} else if (arr2[dindex][dindex]==0) {
swaprows(arr2, dindex, dindex+1);
}
// divide working row by value of working position to get a 1 on the
// diagonal
if (arr2[dindex][dindex] == 0.0) {
return;
} else {
double tempval=arr2[dindex][dindex];
for (long col=0; col<nrows+ncolsrhs; ++col) {
arr2[dindex][col]/=tempval;
}
}
// eliminate value below working position by subtracting a multiple of
// the current row
for (long row=dindex+1; row<nrows; ++row) {
double wval=arr2[row][dindex];
for (long col=0; col<nrows+ncolsrhs; ++col) {
arr2[row][col]-=wval*arr2[dindex][col];
}
}
}
// backward substitution steps
for (long dindex=nrows-1; dindex>=0; --dindex) {
// eliminate value above working position by subtracting a multiple of
// the current row
for (long row=dindex-1; row>=0; --row) {
double wval=arr2[row][dindex];
for (long col=0; col<nrows+ncolsrhs; ++col) {
arr2[row][col]-=wval*arr2[dindex][col];
}
}
}
// assign result to replace rhs
for (long row=0; row<nrows; ++row) {
for (long col=0; col<ncolsrhs; ++col) {
rhs[row][col]=arr2[row][col+nrows];
}
}
for (long row=0; row<nrows; ++row)
delete[] arr2[row];
delete[] arr2;
}
No comments:
Post a Comment