//Autor: Artur Czekalski (Sator)  www.epokaY.net/artur  artur@epokaY.net  Wer. 2007-08-15d
//Liczenie układu równań metodą eliminacji Gaussa na wyrazach typu double
//z wyborem maksymalnego wyrazu do dzielenia z WIERSZY i KOLUMN.
#include <stdio.h> //printf
#include <windows.h> //GetTickCount()
#include <math.h> //fabs
//---------------------------------------------------------------------------
typedef double typwyr;
const typwyr MinEpsilon = 1.0e-9; //dla spr czy == 0, bo błędy zaokrągleń !!!WAŻNE!!!
const int MAXWymiar = 2000; //2000 dla 995
typwyr gM[MAXWymiar*MAXWymiar]; //macierz współczynników
typwyr gB[MAXWymiar]; //wektor wyrazów wolnych
typwyr gX[MAXWymiar]; //dla wyniku
typwyr gOdpX[MAXWymiar]; //Odpowiedni wynik
//---------------------------------------------------------------------------
inline int EliminacjaGaussaMaxWK(const int Wymiar, typwyr *A, typwyr *B, typwyr *X) //zmienia A[] !
{//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
 int i, j, k, l, nr_W, nr_K;  typwyr w;
 int* Permut = new int[Wymiar]; //permutacja do zapamiętania zamian kolumn
 for (i=0; i<Wymiar; ++i) Permut[i] = i; //inicjalizacja permutacji

 //===Przekształcam macierz A na trójkątną górną z jedynkami na przekątnej===
 for (k=0; k < Wymiar-1; ++k) //macierze: wiersze (k) (bez ostatniego)
  {//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; i<Wymiar; ++i) //po wierszach
    for (j=k; j<Wymiar; ++j) //po kolumnach
     if (fabs(A[i*Wymiar+j]) > w) {w = fabs(A[i*Wymiar+j]);  nr_W = i;  nr_K = j;}
   //Musi to być element <> 0.0
   if (w < MinEpsilon) {delete[] Permut;  return -1;} //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) //jeśli nie jest to elem. z k-tego wiersza
    {i = k*Wymiar;  l = nr_W*Wymiar;
	 for (j=k; j<Wymiar; ++j) //po kolumnach; wystarczy zamieniać tylko od k-tej kolumny, bo wcześniej są wszędzie = 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; //zamień też odpowiednie im wyrazy wolne
    }
   //--Potem kolumny
   if (k != nr_K) //jeśli nie jest to elem. z k-tej kolumny
    {for (i=0; i<Wymiar; ++i) //po wszystkich wierszach - całą kolumnę
      {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; //zapamiętaj zmianę kolumn w permutacji
    }

   //---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) //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
    {A[k*Wymiar+k] = 1.0; //już nie dzielę; wiadomo będzie = 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 wcześniejsze są już = 0.0
     B[k] /= w; //oraz wyraz wolny
    }
   //--Redukuję wszystkie wiersze (i) poniżej. Otrzymam zera pod przekątną 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 różne od 0.0 (niezredukowane)
      {A[i*Wymiar+k] = 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

 //===Mogę już łatwo wyliczyć X[]===
 if (fabs(A[Wymiar*Wymiar-1]) < MinEpsilon) {delete[] Permut;  return -2;} //ten wiersz jest linową kombinacją innego wiersza
 X[Permut[Wymiar-1]] = B[Wymiar-1] / A[Wymiar*Wymiar-1];
 if (fabs(X[Permut[Wymiar-1]]) < MinEpsilon) X[Permut[Wymiar-1]] = 0.0; //niedokładność 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]] = 0.0;
  }
 delete[] Permut;
 return 0; //wsz. OK.
}
//===========================================================================
int main(int, char* [])
{int i, j, czas;
 //Ustalenie macierzy współczynników
 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 wyrazów 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, gM, gB, gX); //zmienia gM
 czas = GetTickCount() - czas;
 //---
 typwyr 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: printf("MinGW Dev. Studio C++ 2.05 (GCC 3.4.2): Czas=%d", czas);
 getchar();  return 0;
}