//Autor: Artur Czekalski (Sator)  www.epokaY.net/artur  artur@epokaY.net  Wer. 2007-08-10
//Liczenie układu równań - wersja z wyborem min. wyrazu do dzielenia
program Project1;
{$APPTYPE CONSOLE}
uses SysUtils, windows;
type typwyr = type Double; 
//---------------------------------------------------------------------------
const MinEpsilon: typwyr = 1.0e-9; //dla spr czy = 0, bo błędy zaokrągleń !!!WAŻNE!!!
const MAXWymiar: Integer = 2000;
var
 gM: array [0..2000*2000-1] of typwyr; //macierz współczynników
 gB: array [0..2000-1] of typwyr; //wektor wyrazów wolnych
 gX: array [0..2000-1] of typwyr; //dla wyniku
 gOdpX: array [0..2000-1] of typwyr; //Odpowiedni wynik
//---------------------------------------------------------------------------
Procedure WypiszUklad(Wymiar: Integer; var M: array of typwyr; B: array of typwyr);
var i, j: Integer;
begin
 For j:=0 To Wymiar-1 Do //wiersze
  begin
   for i:=0 To Wymiar-1 Do //kolumny
    write(M[j*Wymiar+i]:3:0, ' ');
   writeln(': ', B[j]:3:0);
  end;
end;
//---------------------------------------------------------------------------
Function EliminacjaGaussaMax(const Wymiar: Integer; var A: array of typwyr;
          var B: array of typwyr; var X: array of typwyr): Integer; //zmienia A[] !
var
 i, j, k, nr_max: Integer;   w, max: typwyr;
begin
 //Rozwiązuje układ równań liniowych metodą eliminacji Gaussa -wsz. przypadki
 //A[Wymiar*Wymiar] - macierz kwadratowa stopnia Wymiar; A[i,j] = A[i*Wymiar+j]
 //B[Wymiar] - wektor wyrazów wolnych
 //X[Wymiar] - tablica wyników
 //---Przekształcam macierz A na trójkątną górną z jedynkami na przekątnej---
 For k:=0 To Wymiar-2 Do //wiersze (k) (bez ostatniego)
  begin
   //Znajdź największy element <> 0 w kolumnie k (bo mniejsze będą błędy przy dzieleniu przez małą liczbę)
   i := k;  max := abs(A[i*Wymiar+k]);  nr_max := i;
   Inc(i);
   while i < Wymiar Do
    begin
     If abs(A[i*Wymiar+k]) > max Then
      begin max := abs(A[i*Wymiar+k]);  nr_max := i;
      end;
     Inc(i);
    end;
   //Musi to być wiersz z NIEZEROWYM elementem w kolumnie k (bo nim będę dzielił)
   //writeln('Max= ', max:0:20, ' nr_max= ', nr_max);
   If max < MinEpsilon Then begin RESULT := -1; exit; end; //nie ma takiego -kolumna (k) jest linową kombinacją innej kolumny (ma same zera)

   If k <> nr_max Then //Jeśli nie jest to elem. z k-tego wiersza, to zamień wiersz k-ty z nr_max
    begin
     for j:=k To Wymiar-1 Do //kolumny; wystarczy zamieniać tylko od k-tej kolumny
      begin w := A[k*Wymiar+j];  A[k*Wymiar+j] := A[nr_max*Wymiar+j];  A[nr_max*Wymiar+j] := w;
      end;
     w := B[k];  B[k] := B[nr_max];  B[nr_max] := w; //zamień też odpowiednie im wyrazy wolne
    end;
   //---------------
   w := A[k*Wymiar+k]; //przez ten element będę dzielił
   If w <> 1.0 Then //tylko gdy różny od 1.0;  if (fabs(w-1.0) > MinEpsilon)
    //Dzielę przez 'pierwszy' tzn. k-ty wyraz wszystkie elementy tego wiersza
    begin A[k*Wymiar+k] := 1.0; //już nie dzielę; wiadomo ma być =1.0
     For i:=k+1 To Wymiar-1 Do A[k*Wymiar+i] := A[k*Wymiar+i] / w; //wystarczy te od kolumny k-tej, bo wcześniejsze są już =0.0
     B[k]:= B[k] / w; //oraz wyraz wolny
    end;
   //Redukuję wszystkie wiersze (j) poniżej. Otrzymam zera pod przekątną w kolumnie k
   For j:=k+1 To Wymiar-1 Do
    begin w := A[j*Wymiar+k];
     If abs(w) > MinEpsilon Then //tylko te różne od 0.0 (niezredukowane)
      begin For i:=k To Wymiar-1 Do A[j*Wymiar+i] := A[j*Wymiar+i] - A[k*Wymiar+i] * w; //kolumny (i)
       B[j] := B[j] - B[k] * w; //oraz wyraz wolny
      end;
    end;
  end;

 //---Mogę już łatwo wyliczyć X[]---
 If abs(A[Wymiar*Wymiar-1]) < MinEpsilon Then begin RESULT := -2; exit; end; //ten wiersz jest linową kombinacją innego wiersza
 X[Wymiar-1] := B[Wymiar-1] / A[Wymiar*Wymiar-1];
 If abs(X[Wymiar-1]) < MinEpsilon Then X[Wymiar-1] := 0.0; //niedokładność reprezentacji liczb rzeczywistych

 For k:=Wymiar-2 DownTo 0 Do //wiersze k
  begin X[k] := B[k];
   For j:=k+1 To Wymiar-1 Do X[k] := X[k] - A[k*Wymiar+j] * X[j]; //kolumny j; bez dzielenia, bo A[Wymiar*k+k]==1.0
   If abs(X[k]) < MinEpsilon Then X[k] := 0.0;
  end;
 RESULT := 0; //wsz. OK.
end;
//===========================================================================
VAR czas: Cardinal;  i, j: Integer;  SumaBledow: typwyr;
 label KONIEC;
BEGIN
 //Ustalenie macierzy współczynników
 For i:=0 To MAXWymiar*MAXWymiar-1 Do gM[i] := 1+i-(i mod 3)-(i mod 13)+(i mod 5)+i div 3- i*i div 995-((i*5) mod 11);
 //Ustalenie wektora poprawnego wyniku (dla spr.)
 For i:=0 To MAXWymiar-1 Do gOdpX[i] := (i+1)*10;

 //Ustalenie wektora wyrazów wolnych zgodnie z gM[] i gOdpX[]
 For j:=0 To MAXWymiar-1 Do //wiersze
  begin gB[j] := 0.0;
   for i:=0 To MAXWymiar-1 Do //kolumny
    gB[j] := gB[j] + gM[j*MAXWymiar+i] * gOdpX[i];
  end;
 //---
 //for i:=0 To MAXWymiar-1 Do writeln('gOdpX[', i, ']= ', gOdpX[i]:23:19);
 //WypiszUklad(MAXWymiar, gM, gB);
 //---
 writeln('Obliczam uklad rownan:');
 czas := GetTickCount;
 i := EliminacjaGaussaMax(MAXWymiar, gM, gB, gX); //zmienia gM
 czas := GetTickCount - czas;
 //---
 SumaBledow := 0.0;
 writeln('i=', i);
 If i < 0 Then begin writeln('Nie mozna obliczyc tego ukladu.');  goto KONIEC; end;
 //---Wynik
 For i:=0 To MAXWymiar-1 Do
  begin
   //writeln('x[', i, ']= ',  gX[i]:23:19, ' Blad:', (gOdpX[i]-gX[i]):30:27);
   SumaBledow := SumaBledow + abs(gOdpX[i] - gX[i]);
  end;
 writeln('Suma bezwzgledna bledow = ', SumaBledow:0:20);
 //---
 KONIEC: write('Delphi 7.0 PE: Czas=', czas);
 readln;
END.