//Autor: Artur Czekalski (Sator)  www.epokaY.net/artur  artur@epokaY.net  Wer. 2007-08-15d
//Liczenie ukladu rownan metoda eliminacji Gaussa na wyrazach typu double
//z wyborem maksymalnego wyrazu do dzielenia z WIERSZY i KOLUMN.
import std.c.windows.windows; //GetTickCount
import std.stdio; //writefln
import std.math; ///fabs
//---------------------------------------------------------------------------
typedef double typwyr;
const double MinEpsilon = 1.0E-9; //dla spr czy == 0, bo bledy zaokraglen !!!WAZNE!!!
const int MAXWymiar = 2000; //2000 dla 995
double gM[]; //macierz wspolczynnikow
double gB[MAXWymiar]; //wektor wyrazow wolnych
double gX[MAXWymiar]; //dla wyniku
double gOdpX[MAXWymiar]; //Odpowiedni wynik
//---------------------------------------------------------------------------
int EliminacjaGaussaMaxWK(const int Wymiar, typwyr A[], typwyr B[], typwyr X[]) //zmienia A[] !
{//Rozwiazuje uklad rownan liniowych metoda 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 wyrazow wolnych
 //X[Wymiar] - tablica wynikow
 int i, j, k, l, nr_W, nr_K;  typwyr w;
 int Permut[] = new int[Wymiar]; //permutacja do zapamietania zamian kolumn
 for (i=0; i<Wymiar; ++i) Permut[i] = i; //inicjalizacja permutacji

 //===Przeksztalcam macierz A na trojkatna gorna z jedynkami na przekatnej===
 for (k=0; k < Wymiar-1; ++k) //macierze: wiersze (k) (bez ostatniego)
  {//Znajdz najwiekszy co do wart. bezwzgl. element <> 0 w macierzy: od k-tego wiersza i k-tej kolumny 
   w = cast(typwyr)0.0; //maksimum
   for (i=k; i<Wymiar; ++i) //po wierszach
    for (j=k; j<Wymiar; ++j) //po kolumnach
     if (fabs(A[i*Wymiar+j]) > w) {w = cast(typwyr)fabs(A[i*Wymiar+j]);  nr_W = i;  nr_K = j;}
   //Musi to byc element <> 0.0
   if (w < MinEpsilon) {delete Permut;  return -1;} //nie ma takiego (same zera)

   //---Zamien ze soba wiersze k-ty i nr_W oraz kolumny k-ta i nr_K
   //--Najperw wiersze
   if (k != nr_W) //jesli nie jest to elem. z k-tego wiersza
    {i = k*Wymiar;  l = nr_W*Wymiar;
	   for (j=k; j<Wymiar; ++j) //po kolumnach; wystarczy zamieniac tylko od k-tej kolumny, bo wczesniej sa wszedzie = 0.0
      {w = A[i+j];  A[i+j] = A[l+j];  A[l+j] = w;
      }
     w = B[k];  B[k] = B[nr_W];  B[nr_W] = w; //zamien tez odpowiednie im wyrazy wolne
    }
   //--Potem kolumny
   if (k != nr_K) //jeali nie jest to elem. z k-tej kolumny
    {for (i=0; i<Wymiar; ++i) //po wszystkich wierszach - cala kolumne
      {w = A[i*Wymiar+k];  A[i*Wymiar+k] = A[i*Wymiar+nr_K];  A[i*Wymiar+nr_K] = w;
      }
     i = Permut[k];  Permut[k] = Permut[nr_K];  Permut[nr_K] = i; //zapamietaj zmiane kolumn w permutacji
    }

   //---Redukcja wyrazow z kolumny k-tej ponizej wiersza k-tego
   w = A[k*Wymiar+k]; //przez ten element bede dzielil elem. z wiersza
   if (w != 1.0) //tylko jesli rozny od 1.0;  if (fabs(w-1.0) > MinEpsilon)
    //Dziele przez 'pierwszy' tzn. k-ty wyraz wszystkie elementy tego wiersza
    {A[k*Wymiar+k] = cast(typwyr)1.0; //juz nie dziele; wiadomo bedzie = 1.0
     l = k*Wymiar; //nr wiersza
     for (j=k+1; j<Wymiar; ++j) A[l+j] /= w; //wystarczy te od kolumny k-tej, bo wczesniejsze sa juz = 0.0
     B[k] /= w; //oraz wyraz wolny
    }
   //--Redukuje wszystkie wiersze (i) ponizej. Otrzymam zera pod przekatna w kolumnie k
   for (i=k+1; i<Wymiar; ++i) //po wierszach z kolumny k-tej
    {w = A[i*Wymiar+k]; //'pierwszy' element
     if (fabs(w) > MinEpsilon) //tylko te rozne od 0.0 (niezredukowane)
      {A[i*Wymiar+k] = cast(typwyr)0.0;
	   for (j=k+1; j<Wymiar; ++j) A[i*Wymiar+j] -= A[k*Wymiar+j] * w; //po kolumnach (j)
       B[i] -= B[k] * w; //oraz wyraz wolny
      }
    }
  } //Koniec redukcji macierzy

 //===Moge juz latwo wyliczyc X[]===
 if (fabs(A[Wymiar*Wymiar-1]) < MinEpsilon) {delete Permut;  return -2;} //ten wiersz jest linowa kombinacja innego wiersza
 X[Permut[Wymiar-1]] = B[Wymiar-1] / A[Wymiar*Wymiar-1];
 if (fabs(X[Permut[Wymiar-1]]) < MinEpsilon) X[Permut[Wymiar-1]] = cast(typwyr)0.0; //niedokladnosc reprezentacji liczb rzeczywistych

 for (k=Wymiar-2; k>=0; --k) //wiersze k
  {X[Permut[k]] = B[k];
   for (j=k+1; j<Wymiar; ++j) X[Permut[k]] -= A[k*Wymiar+j] * X[Permut[j]]; //kolumny j; bez dzielenia, bo A[Wymiar*k+k]==1.0
   if (fabs(X[Permut[k]]) < MinEpsilon) X[Permut[k]] = cast(typwyr)0.0;
  }
 delete Permut;
 return 0; //wsz. OK.
}
//===========================================================================
int main(char[][] args)
{gM = new double[MAXWymiar*MAXWymiar];
 int i, j, czas;
 //Ustalenie macierzy wspolczynnikow
 for (i=0; i<MAXWymiar*MAXWymiar; ++i) gM[i] = 1+i-(i%3)-(i%13)+(i%5)+i/3- i*i/995-((i*5)%11); //995

 //Ustalenie wektora poprawnego wyniku (dla spr.)
 for (i=0; i<MAXWymiar; ++i) gOdpX[i] = (i+1)*10;

 //Ustalenie wektora wyrazow wolnych zgodnie z gM[] i gOdpX[]
 for (j=0; j<MAXWymiar; ++j) //wiersze
 {gB[j] = 0.0;
  for (i=0; i<MAXWymiar; ++i) //kolumny
   gB[j] += gM[j*MAXWymiar+i] * gOdpX[i];
 }
 //---
 printf("Obliczam uklad rownan:\n");
 czas = GetTickCount();
 i = EliminacjaGaussaMaxWK(MAXWymiar, cast(typwyr[])gM, cast(typwyr[])gB, cast(typwyr[])gX); //zmienia gM
 czas = GetTickCount() - czas;
 //---
 double SumaBledow = 0.0;
 printf("i=%d\n", i);
 if (i < 0) {printf("Nie mozna obliczyc tego ukladu.\n");  goto KONIEC;}
 //---Wynik
 for (i=0; i<MAXWymiar; ++i)
  SumaBledow += fabs(gOdpX[i]-gX[i]);
 printf("Suma bezwzgledna bledow = %.20f\n", SumaBledow);
 //---
 KONIEC: writef("DigitalMars D v2.0: Czas=%d", czas);
 delete gM;
 getchar();  return 0;
}