//Autor: Artur Czekalski (Sator)  www.epokaY.net/artur  artur@epokaY.net  Wer. 2007-08-04
//Liczenie ukladu rownan - wersja z wyborem max. wyrazu do dzielenia
import std.c.windows.windows; //GetTickCount
import std.stdio; //writefln
import std.math; ///fabs
//---------------------------------------------------------------------------
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 EliminacjaGaussaMax(int Wymiar, double A[], double B[], double 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]
 //B[Wymiar] - wektor wyrazow wolnych
 //X[Wymiar] - tablica wynikow
 int i, j, k, nr_max;  double w, max;

 //---Przeksztalcam macierz A na trojkatna gorna z jedynkami na przekatnej---
 for (k=0; k < Wymiar-1; ++k) //wiersze (k) (bez ostatniego)
  {
   //Znajdz najwiekszy element != 0 w kolumnie k (bo mniejsze beda bledy przy dzieleniu przez mala liczbe)
   i = k;  max = fabs(A[i*Wymiar+k]);  nr_max = i;
   while (++i < Wymiar)
    if (fabs(A[i*Wymiar+k]) > max) {max = fabs(A[i*Wymiar+k]);  nr_max = i;}
   //Musi to by wiersz z NIEZEROWYM elementem w kolumnie k (bo nim bede dzielil)
   if (max < MinEpsilon) return -1; //nie ma takiego -kolumna (k) jest linowa kombinacja innej kolumny (ma same zera)

   if (k != nr_max) //Jesli nie jest to elem. z k-tego wiersza, to zamien wiersz k-ty z nr_max
    {
     for (j=k; j<Wymiar; ++j) //kolumny; wystarczy zamieniac tylko od k-tej kolumny
      {w = A[k*Wymiar+j];  A[k*Wymiar+j] = A[nr_max*Wymiar+j];  A[nr_max*Wymiar+j] = w;
      }
     w = B[k];  B[k] = B[nr_max];  B[nr_max] = w; //zamien tez odpowiednie im wyrazy wolne
    }
   //---------------
   w = A[k*Wymiar+k]; //przez ten element bede dzielil
   if (w != 1.0) //tylko gdy 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] = 1.0; //juz nie dziele; wiadomo ma byc =1.0
     for (i=k+1; i<Wymiar; ++i) A[k*Wymiar+i] /= w; //wystarczy te od kolumny k-tej, bo wczesniejsze sa juz =0.0
     B[k] /= w; //oraz wyraz wolny
    }
   //Redukuje wszystkie wiersze (j) ponizej. Otrzymam zera pod przekatna w kolumnie k
   for (j=k+1; j<Wymiar; ++j)
    {w = A[j*Wymiar+k];
     if (fabs(w) > MinEpsilon) //tylko te rozne od 0.0 (niezredukowane)
      {for (i=k; i<Wymiar; ++i) A[j*Wymiar+i] -= A[k*Wymiar+i] * w; //kolumny (i)
       B[j] -= B[k] * w; //oraz wyraz wolny
      }
    }
  }
 //---Moge juz latwo wyliczyc X[]---
 if (fabs(A[Wymiar*Wymiar-1]) < MinEpsilon) return -2; //ten wiersz jest linowa kombinacja innego wiersza
 X[Wymiar-1] = B[Wymiar-1] / A[Wymiar*Wymiar-1];
 if (fabs(X[Wymiar-1]) < MinEpsilon) X[Wymiar-1] = 0.0; //niedokladnosc reprezentacji liczb rzeczywistych

 for (k=Wymiar-2; k>=0; --k) //wiersze k
  {X[k] = B[k];
   for (j=k+1; j<Wymiar; ++j) X[k] -= A[k*Wymiar+j] * X[j]; //kolumny j; bez dzielenia, bo A[Wymiar*k+k]==1.0
   if (fabs(X[k]) < MinEpsilon) X[k] = 0.0;
  }
 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 = EliminacjaGaussaMax(MAXWymiar, gM, gB, 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;
}