//Autor: Artur Czekalski (Sator)  www.epokaY.net/artur  artur@epokaY.net  Wer. 2007-08-16d
//Liczenie układu równań metodą eliminacji Gaussa na wyrazach typu double
//z wyborem maksymalnego wyrazu do dzielenia z WIERSZY i KOLUMN.
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
//---------------------------------------------------------------------------
Function EliminacjaGaussaMaxWK(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, l, nr_W, nr_K: Integer;  w: typwyr;
 Permut: array of Integer; //permutacja do zapamiętania zamian kolumn
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]: i-ty wiersz, j-ta kolumna
 //B[Wymiar] - wektor wyrazów wolnych
 //X[Wymiar] - tablica wyników

 SetLength(Permut, Wymiar); //przydzielenie pamięci
 For i:=0 To Wymiar-1 Do Permut[i] := i; //inicjalizacja permutacji

 //===Przekształcam macierz A na trójkątną górną z jedynkami na przekątnej===
 For k:=0 To Wymiar-2 Do //macierze: wiersze (k) (bez ostatniego)
  begin
   //Znajdź największy co do wart. bezwzgl. element <> 0 w macierzy: od k-tego wiersza i k-tej kolumny
   w := 0.0; //maksimum
   For i:=k To Wymiar-1 Do //po wierszach
    For j:=k To Wymiar-1 Do //po kolumnach
     If abs(A[i*Wymiar+j]) > w Then begin w := abs(A[i*Wymiar+j]);  nr_W := i;  nr_K := j; end;
   //Musi to być element <> 0.0
   If w < MinEpsilon Then begin Permut:= NIL; RESULT := -1; exit; end; //nie ma takiego (same zera)

   //---Zamień ze sobą wiersze k-ty i nr_W oraz kolumny k-tą i nr_K
   //--Najperw wiersze
   If k <> nr_W Then //jeśli nie jest to elem. z k-tego wiersza
    begin i := k*Wymiar;  l := nr_W*Wymiar;
	   for j:=k To Wymiar-1 Do //po kolumnach; wystarczy zamieniać tylko od k-tej kolumny, bo wcześniej są wszędzie = 0.0
      begin w := A[i+j];  A[i+j] := A[l+j];  A[l+j] := w;
      end;
     w := B[k];  B[k] := B[nr_W];  B[nr_W] := w; //zamień też odpowiednie im wyrazy wolne
    end;
   //--Potem kolumny
   if k <> nr_K Then //jeśli nie jest to elem. z k-tej kolumny
    begin
     For i:=0 To Wymiar-1 Do //po wszystkich wierszach - całą kolumnę
      begin w := A[i*Wymiar+k];  A[i*Wymiar+k] := A[i*Wymiar+nr_K];  A[i*Wymiar+nr_K] := w;
      end;
     i := Permut[k];  Permut[k] := Permut[nr_K];  Permut[nr_K] := i; //zapamiętaj zmianę kolumn w permutacji
    end;

   //---Redukcja wyrazów z kolumny k-tej poniżej wiersza k-tego
   w := A[k*Wymiar+k]; //przez ten element będę dzielił elem. z wiersza
   If w <> 1.0 Then //tylko jeśli 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 będzie = 1.0
     l := k*Wymiar; //nr wiersza
     For j:=k+1 To Wymiar-1 Do A[l+j] := A[l+j] / 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 (i) poniżej. Otrzymam zera pod przekątną w kolumnie k
   For i:=k+1 To Wymiar-1 Do //po wierszach z kolumny k-tej
    begin w := A[i*Wymiar+k]; //'pierwszy' element
     If abs(w) > MinEpsilon Then //tylko te różne od 0.0 (niezredukowane)
      begin A[i*Wymiar+k] := 0.0;
	     For j:=k+1 To Wymiar-1 Do A[i*Wymiar+j] := A[i*Wymiar+j] - A[k*Wymiar+j] * w; //po kolumnach (j)
       B[i] := B[i] - B[k] * w; //oraz wyraz wolny
      end;
    end;
  end; //Koniec redukcji macierzy

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

 For k:=Wymiar-2 DownTo 0 Do //wiersze k
  begin X[Permut[k]] := B[k];
   For j:=k+1 To Wymiar-1 Do X[Permut[k]] := X[Permut[k]] - A[k*Wymiar+j] * X[Permut[j]]; //kolumny j; bez dzielenia, bo A[Wymiar*k+k]==1.0
   If abs(X[Permut[k]]) < MinEpsilon Then X[Permut[k]] := 0.0;
  end;
 Permut := NIL;
 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;
 //---
 writeln('Obliczam uklad rownan:');
 czas := GetTickCount;
 i := EliminacjaGaussaMaxWK(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
  SumaBledow := SumaBledow + abs(gOdpX[i] - gX[i]);
 writeln('Suma bezwzgledna bledow = ', SumaBledow:0:20);
 //---
 KONIEC: write('Delphi 7.0 PE: Czas=', czas);
 readln;
END.