// (c) 2003. Dinko Korunic 'kreator'
//
// This program is free software; you can redistribute it and/or modify it
// under the terms of the GNU General Public License as published by the
// Free Software Foundation; either version 2 of the License, or (at your
// option) any later version.

#include <stl.h>
#include <iostream>
#include <fstream>
#include <string>

// makroi za skipanje razlicitih delimitera..
#define SkipDelimiterNL(ifs) SkipChar(ifs, " \n")
#define SkipDelimiters(ifs) SkipChar(ifs, " \t")

// E ili preciznost
#define E 0.0001

using namespace std;

typedef vector<double> Double1D;  // jednodimenzionalni STL vektor
typedef vector<Double1D> Double2D;  // dvodimenzionalni STL vektor

// osnovna klasa
class Matrica
{
  protected:
    int redaka, stupaca; // podaci o podacima :-)
    Double2D podaci; // 2D STL vektor sa elementima
  public:
    Matrica(); // osnovni konstruktor
    Matrica(const int, const int); // konstruktor sa dimenzijama
    Matrica(const Matrica &); // konstruktor sa kopiranjem
    double Dohvati(int, int); // dohvati element
    double Postavi(int, int, const double); // postavi element
    Matrica& operator= (const Matrica &); // nadogr. pridruzivanje
    Matrica operator+ (const Matrica &); // nadogr. zbrajanje
    Matrica operator- (const Matrica &); // nadogr. oduzimanje
    Matrica operator* (const Matrica &); // nadogr. mnozenje
    Matrica operator* (double); // nadogr. mnozenje skalarom
    Matrica& operator+= (const Matrica &); // nadogr. +=
    Matrica& operator-= (const Matrica &); // nadogr. -=
    Matrica& Transponiraj(); // nadogradjeno transponiranje
    Matrica LUdekompozicija(); // LU
    Matrica LUPdekompozicija(Matrica &); // LUP
    Matrica& SupstUnaprijed(Matrica &); // unaprijed
    Matrica& SupstUnazad(Matrica &); // unazad
    friend ostream& operator<< (ostream &, Matrica &); // nadogr. <<
    friend istream& operator>> (istream &, Matrica &); // nadogr. >>
};

class LinearniSustav : public Matrica 
{
  private:
    Matrica A, b;
  public:
    Matrica SolveLU();
    Matrica SolveLUP();
    LinearniSustav(Matrica &, Matrica &);
};

// prazni konstruktor - treba nam samo za retke/stupce
Matrica::Matrica()
{
  redaka = stupaca = 0;
}

// konstruktor za matricu iz matrice
//      ulaz: referenca na izvornu matricu
//      izlaz: nista
Matrica::Matrica(const Matrica &izvor)
{
  *this = izvor;
}

// konstruktor za stvaranje nove matrice proizvoljnih dimenzija
//      ulaz: retci, stupci
//      izlaz: nista
Matrica::Matrica(const int _redaka, const int _stupaca) :
  redaka(_redaka), stupaca(_stupaca)
{
  Double2D tmp(redaka, stupaca);
  podaci = tmp;
}

// dohvacanje proizvoljnog elementa iz matrice
//      ulaz: stupac, redak
//      izlaz: dohvaceni element
double Matrica::Dohvati(int stupac, int redak)
{
  if (!(--stupac < stupaca))
    throw "Dohvati: stupac izvan matrice";
  if (!(--redak < redaka))
    throw "Dohvati: redak izvan matrice";

  return podaci[redak][stupac];
}

// postavljanje proizvoljnog elementa u matrici
//      ulaz: stupac, redak, podatak
//      izlaz: postavljeni element
double Matrica::Postavi(int stupac, int redak, const double podatak)
{
  if (!(--stupac < stupaca))
    throw "Postavi: stupac izvan matrice";
  if (!(--redak < redaka))
    throw "Postavi: redak izvan matrice";

  return (podaci[redak][stupac] = podatak);
}

// operator pridruzivanja
//      ulaz: referenca na matricu
//      izlaz: referenca na vlastiti objekt
Matrica& Matrica::operator= (const Matrica &izvor)
{
  stupaca = izvor.stupaca;
  redaka = izvor.redaka; // izjednacimo stupce/retke
  podaci = izvor.podaci; // i onda same podatke (nadogr. =)

  return *this;
}

// operator zbrajanja
//      ulaz: referenca na matricu
//      izlaz: matrica (novi objekt)
Matrica Matrica::operator+ (const Matrica &izvor)
{
  if ((izvor.stupaca != stupaca)
      || (izvor.redaka != redaka))
    throw "operator+: matrice nemaju iste dimenzije";

  Matrica rezultat(*this); // koristimo nadogradjeno instanciranje
 
  int i, j;
  
  for (i = 0; i < redaka; ++i)
    for (j = 0; j < stupaca; ++j)
      rezultat.podaci[i][j] += izvor.podaci[i][j]; // i za sve elem.
  
  return rezultat;
}

// operator oduzimanja
//      ulaz: referenca na matricu
//      izlaz: matrica (novi objekt)
Matrica Matrica::operator- (const Matrica &izvor)
{
  if ((izvor.stupaca != stupaca)
      || (izvor.redaka != redaka))
    throw "operator+: matrice nemaju iste dimenzije";

  Matrica rezultat(*this);
 
  int i, j;
  
  for (i = 0; i < redaka; ++i)
    for (j = 0; j < stupaca; ++j)
      rezultat.podaci[i][j] -= izvor.podaci[i][j]; // za sve..
  
  return rezultat;
}

// operator mnozenja skalarom
//      ulaz: skalar
//      izlaz: matrica (novi objekt)
Matrica Matrica::operator* (const double izvor)
{
  Matrica rezultat(*this);
 
  int i, j;
  
  for (i = 0; i < redaka; ++i)
    for (j = 0; j < stupaca; ++j)
      rezultat.podaci[i][j] *= izvor; // za sve elemente
  
  return rezultat;
}

// operator mnozenja matricom
//      ulaz: referenca na matricu
//      izlaz: matrica (novi objekt)
Matrica Matrica::operator* (const Matrica &izvor)
{
  if (stupaca != izvor.redaka)
    throw "operator*: matrice nisu ulancane";

  Matrica rezultat(redaka, izvor.stupaca);

  int i, j, k;
  
  for (i = 0; i < redaka; ++i)
    for (j = 0; j < izvor.stupaca; ++j)
      for (k = 0; k < stupaca; ++k)
        rezultat.podaci[i][j] += podaci[i][k] * izvor.podaci[k][j];
        // prema formuli A[i][j] = Sum a[i][j]b[k][j], k=1..n

  return rezultat;
}

// operator zbrajanja i izjednacenja
//      ulaz: referenca na matricu
//      izlaz: referenca na matricu
Matrica& Matrica::operator+= (const Matrica &izvor)
{
  return (*this = *this + izvor); // najjednostavnije
}

// operator oduzimanja i izjednacenja
//      ulaz: referenca na matricu
//      izlaz: referenca na matricu
Matrica& Matrica::operator-= (const Matrica &izvor)
{
  return (*this = *this - izvor); // najjednostavnije
}

// transponiranje matrice
//      ulaz: nista
//      izlaz: referenca na matricu
Matrica& Matrica::Transponiraj()
{
  Matrica rezultat(stupaca, redaka);

  int i, j;
  for (i = 0; i < redaka; ++i)
    for (j = 0; j < stupaca; ++j)
      rezultat.podaci[j][i] = podaci[i][j];
      // napravimo novi objekt tocnih dimenzija
      // i napunimo ga podacima
  
  return (*this = rezultat);
}

// supstitucija unaprijed
//  ulaz: matrica
//  izlaz: matrica
Matrica& Matrica::SupstUnaprijed(Matrica &b)
{
  int i, j, n = b.redaka;

  for (i = 0; i < n - 1; ++i)
    for (j = i + 1; j < n; ++j)
      b.podaci[j][0] -= podaci[j][i] * b.podaci[i][0];

  return b;
}

// supstitucija unazad
//  ulaz: matrica
//  izlaz: matrica
Matrica& Matrica::SupstUnazad(Matrica &b)
{
  int i, j, n = b.redaka;

  for (i = n - 1; i >= 0; --i)
  {
    b.podaci[i][0] /= podaci[i][i];
    for (j = 0; j < i; ++j)
      b.podaci[j][0] -= podaci[j][i] * b.podaci[i][0];
  }

  return b;
}

// LU dekompozicija matrica
//  ulaz: nista
//  izlaz: nova LU matrica
Matrica Matrica::LUdekompozicija()
{
  if (redaka != stupaca)
    throw "Matrica A nije kvadratna";

  Matrica A(*this);

  int i, j, k, n = redaka;
  for (i = 0; i < n - 1; ++i)
    for (j = i + 1; j < n; ++j)
    {
      if (fabs(A.podaci[i][i]) < E)
        throw "LU dekompozicija: dijeljenje sa elementom bliskim 0";
      A.podaci[j][i] /= A.podaci[i][i];
      for (k = i + 1; k < n; ++k)
        A.podaci[j][k] -= A.podaci[j][i] * A.podaci[i][k];
    }

  return A;
}

// LUP dekompozicija matrica
//  ulaz: matrica b
//  izlaz: nova LU matrica
//  NB: matrica[redak][stupac]
Matrica Matrica::LUPdekompozicija(Matrica &b)
{
  if (redaka != stupaca)
    throw "Matrica A nije kvadratna";
  if (b.stupaca != 1)
    throw "Matrica B nije [n, 1]";

  Matrica A(*this);

  int i, j, k, pivot_x = 0, n = redaka;
  double pivot, tmp_pivot;

  for (i = 0; i < n - 1; ++i)
  {
    pivot = 0.0;
    for (k = i; k < n; ++k)
    {
      tmp_pivot = fabs(A.podaci[k][i]);
      if (pivot < tmp_pivot)
      {
        pivot = tmp_pivot;
        pivot_x = k;
      }
    }
    if (pivot < E)
      throw "LUP dekompozicija: singularna matrica";

    for (k = 0; k < n; ++k)
      swap(A.podaci[i][k], A.podaci[pivot_x][k]);
    swap(b.podaci[i][0], b.podaci[pivot_x][0]);

    for (j = i + 1; j < n; ++j)
    {
      A.podaci[j][i] /= A.podaci[i][i];
      for (k = i + 1; k < n; ++k)
        A.podaci[j][k] -= A.podaci[j][i] * A.podaci[i][k];
    }
  }

  return A;
}

// rjesavanje linearnog sustava LU metodom
//  ulaz: nista
//  izlaz: matrica rjesenja
Matrica LinearniSustav::SolveLU()
{
  Matrica rez, b_orig(b);

  rez = A.LUdekompozicija();
  cout << "LU:" << endl << rez << endl;
  rez.SupstUnaprijed(b);
  rez = rez.SupstUnazad(b);
  b = b_orig;

  return rez;
}

// rjesavanje linearnog sustava LUP metodom
//  ulaz: nista
//  izlaz: matrica rjesenja
Matrica LinearniSustav::SolveLUP()
{
  Matrica rez, b_orig(b);

  rez = A.LUPdekompozicija(b);
  //cout << "LUP:" << endl << rez << endl;
  //cout << "b': " << endl << b << endl;
  rez.SupstUnaprijed(b);
  //cout << "y: " << endl << b << endl;
  rez = rez.SupstUnazad(b);
  //cout << "x: " << endl << b << endl;
  b = b_orig;

  return rez;
}

LinearniSustav::LinearniSustav(Matrica &A_, Matrica &b_) : A(A_), b(b_)
{
  cout << "Matrica A: " << endl << A << endl;
  cout << "Matrica b: " << endl << b << endl;
}

// preskacemo zeljene znakove, ma koliko ih god bilo
//      ulaz: input stream, znak za ignoriranje
//      izlaz: nista
void SkipChar(istream &ifs, char *igme)
{
  for (; index(igme, ifs.peek()) != NULL; ifs.ignore());
}

// operator >>
//      ulaz: input stream, referenca na matricu
//      izlaz: referenca na input stream
istream& operator>> (istream &ifs, Matrica &MOb)
{
  while (ifs.peek() != EOF) // dok mozes citaj
  {
    double tmp;

    MOb.stupaca = 0; // najjednostavnije resetirati svaki put
    Double1D tmpRedak; // privremeni redak

    while (ifs && ifs.peek() != '\n')
    {
      ifs >> tmp; // ucitaj jedan znak
      tmpRedak.push_back(tmp); // i jedan po jedan u vektor
      MOb.stupaca++;
      SkipDelimiters(ifs); // preskoci " " ili "\t"
    }
    MOb.podaci.push_back(tmpRedak); // redak u kompletnu tablicu
    SkipDelimiterNL(ifs); // preskoci "\n" i sl.
    MOb.redaka++;
  }

  return ifs;
}

// operator <<
//      ulaz: output stream, referenca na matricu
//      izlaz: referenca na output stream
ostream& operator<< (ostream &ofs, Matrica &MOb)
{
  Double2D::iterator it;

  for (it = MOb.podaci.begin(); it != MOb.podaci.end(); ++it)
  {
    ostream_iterator<double> out(ofs, "\t");
    copy(it->begin(), it->end(), out); // iteriraj sve elemente
    ofs << endl;
  }

  return ofs;
}

// glavna funkcija
int main()
{
  try
  {
    ifstream prva("lab2-prva");
    ifstream druga("lab2-druga");

    if (prva && druga)
    {
      Matrica A;
      Matrica b;

      prva >> A;
      druga >> b;

      LinearniSustav Prvi(A, b);
      Matrica Rjesenje;
      
      Rjesenje = Prvi.SolveLU();
      cout << "Rjesenje LU: " << endl << Rjesenje << endl;

      Rjesenje = Prvi.SolveLUP();
      cout << "Rjesenje LUP: " << endl << Rjesenje << endl;
    }
  }

  catch (const char *msg)
  {
    cout << "Greska: " << msg << endl;
    return EXIT_FAILURE;
  }

  return EXIT_SUCCESS;
}