#define mat_type double

const Eps = 0.1;

int stepen(int a, int b) {
int r = 1;
for (int l = 1; l <= b; l++)
r *= a;

return r;
}

class Vektor {
private:
int n; //razmer vektora
mat_type *data;

public:
Vektor(int in);
Vektor(int in, float a[]);
Vektor(Vektor *a);
~Vektor();

int getSize();
mat_type get(int x);
float norma();

void put(int x, mat_type d);
void print();
};

Vektor::Vektor(int in) {
n = in;

if ((data = (mat_type *) malloc(n * sizeof(mat_type))) == NULL) {
cout << "Not enough memory to allocate buffer\n";
exit(1); /* terminate program if out of memory */
}

mat_type *temp = data;

for(int i = 0; i < n; i++) {
*temp = 0;
temp++;
}
}

Vektor::Vektor(int in, float a[]) {
n = in;

if ((data = (mat_type *) malloc(n * sizeof(mat_type))) == NULL) {
cout << "Not enough memory to allocate buffer\n";
exit(1); /* terminate program if out of memory */
}

mat_type *temp = data;

for(int i = 0; i < n; i++) {
*temp = a[i];
temp++;
}
}

Vektor::Vektor(Vektor *a) {
n = a->getSize();

if ((data = (mat_type *) malloc(n * sizeof(mat_type))) == NULL) {
cout << "Not enough memory to allocate buffer\n";
exit(1); /* terminate program if out of memory */
}

for(int i = 1; i <= n; i++) {
*data = a->get(i);
data++;
}

data -= n;
}

Vektor::~Vektor() {
free(data);
}

int Vektor::getSize() {
return n;
}

mat_type Vektor::get(int x) {
return *(data + x - 1);
}

float Vektor::norma() {
float sum = 0;
for (int i = 1; i <= n; i++) {
// cout << "\n*** : " << get(i);
sum += get(i) * get(i);
}

return sqrt(sum);
}

void Vektor::put(int x, mat_type d) {
if ((x < 1) || (x > n)) exit(1);

*(data + x - 1) = d;
}

void Vektor::print() {
cout << "(";
for (int i = 1; i <= n; i++) {
cout << *data;
if (i != n) cout << ", ";
data++;
}
cout << ")";
data -= n;
}

class Matrix {
private:
int n; //razmer matrici
mat_type *data;

public:
Matrix(int in); //input matrix size for get mem
Matrix(Matrix *a);
~Matrix();

void put(int x, int y, mat_type el);
int getSize();
void trans();

mat_type get(int x, int y);
mat_type minor(int x, int y); //absolyutni sdvig
mat_type norma();
Vektor *getRow(int i);
Vektor *getCol(int i);

void swapCol(int y, int y1);
void swapRow(int x, int x1);

float det();
int getMax(int y); //Index max elem v stroke y
Matrix *copy();
Matrix *obrat();

void print();
};

Matrix::Matrix(int in) {
n = in;

if ((data = (mat_type *) malloc(n * n * sizeof(mat_type))) == NULL) {
cout << "Not enough memory to allocate buffer\n";
exit(1); /* terminate program if out of memory */
}

mat_type *temp = data;

for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++) {
*temp = 0;
temp++;
}

}

Matrix::Matrix(Matrix *a) {
n = a->getSize();

if ((data = (mat_type *) malloc(n * n * sizeof(mat_type))) == NULL) {
cout << "Not enough memory to allocate buffer\n";
exit(1); /* terminate program if out of memory */
}

for (int j = 1; j<=n; j++)
for (int i = 1; i<=n; i++) {
*data = a->get(i, j);
data++;
}

data -= n * n;
}

Matrix::~Matrix() {
free(data);
}

void Matrix::put(int x, int y, mat_type el) {
*(data + (y - 1) * n + (x - 1)) = el;
}

mat_type Matrix::get(int x, int y) {
return *(data + (y - 1) * n + (x - 1));
}

float Matrix::det() {
if (n == 1) return get(1, 1);

mat_type d = 0;
int y;

Matrix new_mat(n - 1);

for(int l = 1; l <= n; l++) {

for(int j = 1; j <= n-1; j++)
for(int i = 1; i <= n-1; i++)
if (i < l) new_mat.put(i, j, get(i, j+1));
else new_mat.put(i, j, get(i+1, j+1));

d += stepen(-1, l + 1) * get(l, 1) * new_mat.det();
}

return d;
}

int Matrix::getMax(int y) {
int max = 1;
mat_type *temp = data + (y - 1) * n;

for(int j = 1; j <= n; j++) {
if (get(max, y) < *temp) max = j;
temp++;
}

return max;
}

mat_type Matrix::minor(int x, int y) {

int sx, sy;
Matrix temp(n - 1);

sx = 0;
sy = 0;
for(int j = 1; j <= n-1; j++) {
if (j >= y) sy = 1;
sx = 0;
for(int i = 1; i <= n-1; i++) {
if (i >= x) sx = 1;
temp.put(i, j, get(i + sx, j + sy));
}
}

return temp.det();
}

void Matrix::swapCol(int x, int x1) {
mat_type buf;
mat_type *temp = data + (x - 1);
mat_type *temp1 = data + (x1 - 1);

for(int j = 0; j < n; j++) {
buf = *temp;
*temp = *temp1;
*temp1 = buf;
temp += n;
temp1 += n;
}
}

void Matrix::swapRow(int y, int y1) {
mat_type buf;
mat_type *temp = data + (y - 1) * n;
mat_type *temp1 = data + (y1 - 1) * n;

for(int i = 0; i < n; i++) {
buf = *temp;
*temp = *temp1;
*temp1 = buf;
temp++;
temp1++;
}
}

int Matrix::getSize() {
return n;
}

void Matrix::trans() {
mat_type d;

for (int i = 2; i <= n; i++)
for (int j = i-1; j >= 1; j--)
if (i != j) {
d = get(i, j);
put(i, j, get(j, i));
put(j, i, d);
}

}

Matrix *Matrix::obrat() {
Matrix *S = new Matrix(n);
float d = det();

for (int i = 1; i <= n; i++)
for (int j = 1; j <= n; j++)
S->put(i, j, minor(i, j) / d * stepen(-1, i+j));

S->trans();
return S;
}

void Matrix::print() {

mat_type *temp = data;

for(int i = 0; i < n; i++) {
for(int j = 0; j < n; j++) {
cout << *temp << "\t";
temp++;
}
cout << "\n";
}
}

mat_type Matrix::norma() {
mat_type r = 0;

for (int i = 1; i <= n; i++)
for (int j = 1; j <= n; j++)
r += get(i, j) * get(i, j);

return sqrt(r);
}

Vektor *Matrix::getRow(int i) {
Vektor *temp = new Vektor(n);

for (int j = 1; j <= n; j++)
temp->put(j, get(i, j));

return temp;
}

Vektor *Matrix::getCol(int i) {
Vektor *temp = new Vektor(n);

for (int j = 1; j <= n; j++)
temp->put(j, get(j, i));

return temp;
}

Vektor *minVV(Vektor *a, Vektor *x) {
if (a->getSize() != x->getSize()) {
cout << "Sorry, size diffrent.\n";
exit(1);
}

int n = a->getSize();

Vektor *b = new Vektor(n);

for (int j = 1; j <= n; j++) {
b->put(j, a->get(j) - x->get(j));
}

return b;
}

Vektor *mulMV(Matrix *a, Vektor *x) {
if (a->getSize() != x->getSize()) {
cout << "Sorry, size diffrent.\n";
exit(1);
}

int n = a->getSize();
mat_type summ;

Vektor *b = new Vektor(n);

for (int j = 1; j <= n; j++) {
summ = 0;
for (int i = 1; i <= n; i++) {
summ += x->get(i) * a->get(i, j);
}
b->put(j, summ);
}

return b;
}

Matrix *mulMM(Matrix *a, Matrix *b) {
if (a->getSize() != b->getSize()) {
cout << "Sorry, size diffrent.\n";
exit(1);
}

int n = a->getSize();

Matrix *c = new Matrix(n);
mat_type summ;

for (int j = 1; j <= n; j++)
for (int i = 1; i <= n; i++) {
summ = 0;

for (int k = 1; k <= n; k++)
summ += a->get(k, j) * b->get(i, k);

c->put(i, j, summ);
}

return c;
}

Vektor *nevyzka(Matrix *a, Vektor *x, Vektor *b) {
Vektor *c = mulMV(a, x);

for (int i = 1; i <= c->getSize(); i++)
c->put(i, c->get(i) - b->get(i));

return c;
}

mat_type mulVV(Vektor *a, Vektor *b) {
mat_type sum = 0;

for (int i = 1; i <= b->getSize(); i++)
sum += a->get(i) * b->get(i);

return sum;
}

Vektor *SimpleX(Matrix *a, Vektor *b, Vektor *x0, int **step) { //x1 = A * x0 + b

**step += 1;

float sum = 0;

Vektor *x1 = mulMV(a, x0);

for (int i = 1; i <= b->getSize(); i++) {
x1->put(i, x1->get(i) + b->get(i));
sum += (x1->get(i) - x0->get(i)) * (x1->get(i) - x0->get(i));
}

if (**step != 1) delete x0;

if (sqrt(sum) > 0.000001) return SimpleX(a, b, x1, step);
return x1;
}

Vektor *Gauss(Matrix *a, Vektor *b, Vektor *x0, int **step) {

**step += 1;

float pogr;

Vektor *x1 = new Vektor(x0);
Vektor *av;

for (int i = 1; i <= x0->getSize(); i++) {
av = a->getCol(i);
x1->put(i, mulVV(av, x1) + b->get(i));
delete av;
}

av = minVV(x0, x1);
pogr = av->norma();

delete av;
if (**step != 1) delete x0;

if (pogr > 0.000001) return Gauss(a, b, x1, step);
return x1;
}

Vektor *Spusk(Matrix *a, Vektor *b, Vektor *x0, int **step) {
**step += 1;

float pogr, m;
Vektor *x1;
Matrix *aT = new Matrix(a); aT->trans();

Vektor *v = nevyzka(a, x0, b);
Vektor *v1 = mulMV(aT, v);

Vektor *v2 = mulMV(a, v1);

m = mulVV(v, v2) / mulVV(v2, v2);

for (int i = 1; i <= v1->getSize(); i++)
v1->put(i, v1->get(i) * m);

x1 = minVV(x0, v1);

delete v;
delete v1;
delete v2;

v1 = minVV(x0, x1);
pogr = v1->norma();

delete v1;
delete aT;

delete x0;

if (pogr > 0.000001) return Spusk(a, b, x1, step);
return x1;
}

int CheckSimpleX(Matrix *a) { //Proverka shodimosti metoda
int i, j;
float sumI, sumJ, indexI = 0, indexJ = 0;

for (i = 1; i <= a->getSize(); i++) {
sumI = 0;
sumJ = 0;
for (j = 1; j <= a->getSize(); j++) {
sumI += a->get(i, j);
sumJ += a->get(j, i);
}
if (sumI < 1) indexI++;
if (sumJ < 1) indexJ++;
}
if ((indexI == a->getSize()) || (indexJ == a->getSize())) return 1;

return 0;
}

void laba4(Matrix *a, Vektor *b) {
int *step;
Matrix *a1 = new Matrix(a);
Vektor *b1 = new Vektor(b);
Vektor *x, *r;

cout << "\tMatrix A:\n"; a->print();
cout << "Privedenie systemi A * x = b -> x = B * x + g: \n";

for (int j = 1; j <= a->getSize(); j++) {
for (int i = 1; i <= a->getSize(); i++)
a1->put(i, j, (a->get(i, j) / a->get(j, j)) * (-1));
a1->put(j, j, 0);
b1->put(j, b->get(j) / a->get(j, j));
}

Vektor *x0 = new Vektor(b1); //nachalnoe priblejenie

cout << "\tMatrix B:\n"; a1->print();
cout << "\n\tNorma B: " << a1->norma();
cout << "\tVektor g: "; b1->print();

cout << "\n\n\t Metod prostoi itteracii:\n";
*step = 0;
if (CheckSimpleX(a1) == 1) {
x = SimpleX(a1, b1, x0, &step);
r = nevyzka(a, x, b);
cout << " Reshenie: \t"; x->print();
cout << "\nChislo itteracii: \t" << *step << "\n";
cout << " Vektor nevyzki: \t"; r->print();
delete(x);
delete(r);
} else
cout << "Metod rashoditsya.\n";

cout << "\n\n\t Metod Gausa-Zeidelya:\n";
*step = 0;
if (CheckSimpleX(a1) == 1) {
x = Gauss(a1, b1, x0, &step);
r = nevyzka(a, x, b);
cout << " Reshenie: \t"; x->print();
cout << "\nChislo itteracii: \t" << *step << "\n";
cout << " Vektor nevyzki: \t"; r->print();
delete(x);
delete(r);
} else
cout << "Metod rashoditsya.\n";

cout << "\n\n\t Metod Skoreishego spuska:\n";
*step = 0;
x = Spusk(a, b, new Vektor(b), &step);
r = nevyzka(a, x, b);
cout << " Reshenie: \t"; x->print();
cout << "\nChislo itteracii: \t" << *step << "\n";
cout << " Vektor nevyzki: \t"; r->print();
delete(x);
delete(r);

delete(a1);
delete(b1);
delete(x0);

}

int main() {
clrscr();

unsigned long StartMem = coreleft();

int n = 4;

float v[4] = {3, 1, 3, 1};

Vektor *b = new Vektor(n, v);

Matrix *a3 = new Matrix(n);

/* a1->put(1, 1, 10); a1->put(2, 1, 1); a1->put(3, 1, -1);
a1->put(1, 2, 1); a1->put(2, 2, 10); a1->put(3, 2, -1);
a1->put(1, 3, -1); a1->put(2, 3, 1); a1->put(3, 3, 10);*/

/* a1->put(1, 1, 2); a1->put(2, 1, 1);
a1->put(1, 2, 1); a1->put(2, 2, -2);*/

a3->put(1, 1, 5.88); a3->put(2, 1, 0.04); a3->put(3, 1, 1.28); a3->put(4, 1, 0.82);
a3->put(1, 2, 0.04); a3->put(2, 2, 7.87); a3->put(3, 2, 0.12); a3->put(4, 2, 1.32);
a3->put(1, 3, 1.28); a3->put(2, 3, 0.12); a3->put(3, 3, 3.19); a3->put(4, 3, 0.82);
a3->put(1, 4, 0.82); a3->put(2, 4, 1.32); a3->put(3, 4, 0.82); a3->put(4, 4, 6.22);

laba4(a3, b);

delete a3;
delete b;

cout << "\n\nLost memory: " << StartMem - coreleft();

return 0;
};

/* a1->put(1, 1, 81); a1->put(2, 1, -45); a1->put(3, 1, 45);
a1->put(1, 2, -45); a1->put(2, 2, 50); a1->put(3, 2, -15);
a1->put(1, 3, 45); a1->put(2, 3, -15); a1->put(3, 3, 38);

a2->put(1, 1, 1.06); a2->put(2, 1, 2.14); a2->put(3, 1, 0.08); a2->put(4, 1, 1.17);
a2->put(1, 2, 3.13); a2->put(2, 2, 4.17); a2->put(3, 2, 0.10); a2->put(4, 2, 9.12);
a2->put(1, 3, 5.04); a2->put(2, 3, 2.20); a2->put(3, 3, 6.11); a2->put(4, 3, 7.19);
a2->put(1, 4, 1.14); a2->put(2, 4, 3.19); a2->put(3, 4, 2.22); a2->put(4, 4, 0.07);

a3->put(1, 1, 5.88); a3->put(2, 1, 0.04); a3->put(3, 1, 1.28); a3->put(4, 1, 0.82);
a3->put(1, 2, 0.04); a3->put(2, 2, 7.87); a3->put(3, 2, 0.12); a3->put(4, 2, 1.32);
a3->put(1, 3, 1.28); a3->put(2, 3, 0.12); a3->put(3, 3, 3.19); a3->put(4, 3, 0.82);
a3->put(1, 4, 0.82); a3->put(2, 4, 1.32); a3->put(3, 4, 0.82); a3->put(4, 4, 6.22);*/

�����