//Autor: Artur Czekalski (Sator)  www.epokaY.net/artur  artur@epokaY.net  Wer. 2007-08-04
//Liczenie układu równań - wersja z wyborem max. wyrazu do dzielenia
#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
//---------------------------------------------------------------------------
int EliminacjaGaussaMax(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]
 //B[Wymiar] - wektor wyrazów wolnych
 //X[Wymiar] - tablica wyników
 int i, j, k, nr_max;  typwyr w, max;

 //---Przekształcam macierz A na trójkątną górną z jedynkami na przekątnej---
 for (k=0; k < Wymiar-1; ++k) //wiersze (k) (bez ostatniego)
  {
   //Znajdź największy element != 0 w kolumnie k (bo mniejsze będą błędy przy dzieleniu przez małą liczbę)
   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 będę dzielił)
   if (max < MinEpsilon) return -1; //nie ma takiego -kolumna (k) jest linową kombinacją innej kolumny (ma same zera)

   if (k != nr_max) //Jeśli nie jest to elem. z k-tego wiersza, to zamień wiersz k-ty z nr_max
    {
     for (j=k; j<Wymiar; ++j) //kolumny; wystarczy zamieniać 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; //zamień też odpowiednie im wyrazy wolne
    }
   //---------------
   w = A[k*Wymiar+k]; //przez ten element będę dzielił
   if (w != 1.0) //tylko gdy 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 ma być =1.0
     for (i=k+1; i<Wymiar; ++i) A[k*Wymiar+i] /= w; //wystarczy te od kolumny k-tej, bo wcześniejsze są już =0.0
     B[k] /= w; //oraz wyraz wolny
    }
   //Redukuję wszystkie wiersze (j) poniżej. Otrzymam zera pod przekątną w kolumnie k
   for (j=k+1; j<Wymiar; ++j)
    {w = A[j*Wymiar+k];
     if (fabs(w) > MinEpsilon) //tylko te różne 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
      }
    }
  }
 //---Mogę już łatwo wyliczyć X[]---
 if (fabs(A[Wymiar*Wymiar-1]) < MinEpsilon) return -2; //ten wiersz jest linową kombinacją innego wiersza
 X[Wymiar-1] = B[Wymiar-1] / A[Wymiar*Wymiar-1];
 if (fabs(X[Wymiar-1]) < MinEpsilon) X[Wymiar-1] = 0.0; //niedokładność 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(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 = EliminacjaGaussaMax(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;
}