#include "matr_bvv.hpp"

// konstructor

matrix::matrix(int n, int m, int s) {
    row = n;
    col = m;
    int j, i;
    sparse = s;
    y = new double*[row];
    for (i = 0; i < n; i++) y[i] = new double[col];
    for (i = 0; i < row; i++)
        for (j = 0; j < col; j++)
            y[i][j] = 0.0;
}

// ���������� �����������

matrix::matrix(const matrix &a) {
    row = a.row;
    col = a.col;
    int j, i;
    sparse = a.sparse;
    y = new double*[row];
    for (i = 0; i < row; i++) y[i] = new double[col];
    for (i = 0; i < row; i++)
        for (j = 0; j < col; j++)
            y[i][j] = a.y[i][j];
}

// ����������

matrix::~matrix() {
    int i;
    for (i = 0; i < row; i++)delete y[i];
    delete [row]y;
    //[row] in this position is need. Destructor is called for every element of array.
}

//�������� ������������

matrix matrix::operator = (const matrix &a){
    int j, i;
    for (i = 0; i < row; i++) delete y[i];
    delete [row]y;
    row = a.row;
    col = a.col;
    y = new double*[row];
    for (i = 0; i < row; i++) y[i] = new double[col];
    for (i = 1; i <= row; i++)
        for (j = 1; j <= col; j++) {
            (*this)(i, j) = a(i, j);
        }
    return *this;
}

// �������� ������

matrix matrix::operator + (const matrix & a) const {
    if ((row != a.row) || (col != a.col)) {
        //sound_tone(1500,1331,1331);
        //        sound(2000);
        //        delay(500);
        //        nosound();

        printf("matrix dimensions must agree for operator+");
        exit(1);
    }
    int j, i;
    matrix t(row, col);
    for (i = 0; i < row; i++)
        for (j = 0; j < col; j++)
            t.y[i][j] = y[i][j] + a.y[i][j];
    return t;
}

// ��������� ������

matrix matrix::operator - (const matrix & a) const {
    if ((row != a.row) || (col != a.col)) {
        //sound_tone(1500,1331,1331);
        //        sound(2000);
        //        delay(500);
        //        nosound();

        printf("matrix dimensions must agree for operator+");
        exit(1);
    }
    matrix t(row, col);
    int j, i;
    for (i = 0; i < row; i++)
        for (j = 0; j < col; j++)
            t.y[i][j] = y[i][j] - a.y[i][j];
    return t;
}

// ������

void matrix::print() {
    int j, i;
    for (i = 0; i < row; i++) {
        printf("\n");
        for (j = 0; j < col; j++) printf(" %5.2f", y[i][j]);
    }
    printf("\n");
}

// ������������ ������

matrix matrix::operator*(const matrix &a) const {
    if (col != a.row) {
        //                sound_tone(1500,1331,1331);
        //        sound(2000);
        //        delay(500);
        //        nosound();

        printf("matrix dimensions must agree for operator*");
        exit(1);
    }
    matrix t(row, a.col);
    int j, k, i;
    for (i = 1; i <= t.row; i++) // ������ �� ������� �������������� �������
        for (j = 1; j <= t.col; j++) // ������ �� �������� �������������� �������
            // ��������� ������� i,j �������������� �������
            for (k = 1; k <= col; k++) t(i, j) = t(i, j)+(*this)(i, k) * a(k, j);
    return (t);
}

// ������� ������ �������� ���������

matrix matrix::operator | (const matrix &b) const {
    if (col != b.row) {
        //                sound_tone(1500,1331,1331);
        //        sound(2000);
        //        delay(500);
        //        nosound();

        printf("matrix dimensions must agree for operator|\n");
        exit(1);
    }
    if (b.col > 1) {
        //                sound_tone(1500,1331,1331);
        //        sound(2000);
        //        delay(500);
        //        nosound();
        printf("number columns for second argument must be 1 for operator|\n");
        exit(1);
    }
    matrix a1(row, col), c(row, col), x(row, 1), b1(row, 1);
    int i, j, h, n1, k;
    double v, r;
    a1 = *this;
    b1 = b;
    c = +a1;
    n1 = row - 1;
    //recalculation b1
    for (i = 1; i <= n1; i++)
        for (h = i + 1; h <= row; h++)
            b1(h, 1) = b1(h, 1) - c(h, i) * b1(i, 1);
    //solve triangular system
    int n = row, i1;
    x(n, 1) = b1(n, 1) / a1(n, n);
    i = n1;
l110:
    v = 0.0;
    i1 = i + 1;
    for (h = i1; h <= n; h++)v += a1(i, h) * x(h, 1);
    x(i, 1) = (b1(i, 1) - v) / a1(i, i);
    i--;
    if (i >= 1)goto l110;
    return x;
}
//**********************member function for triangulation********************

matrix matrix::operator + () {
    matrix c(row, col);
    int i, j, h, n1, k;
    double v, r;
    n1 = row - 1;
    //triangulation
    for (i = 1; i <= n1; i++) {
        h = i + 1;
        if ((*this)(i, i) != 0)goto l50;
        //find nonzero element
        for (k = h; k <= row; k++)if ((*this)(k, i) != 0)goto l30;
        //nonzero element not found
        //                sound_tone(1500,1331,1331);
        //        sound(2000);
        //        delay(500);
        //        nosound();
        printf("nonzero element not found in Gauss elemination");
        exit(1);

        //replace rows
l30:
        for (j = i; j <= row; j++) {
            v = (*this)(i, j);
            (*this)(i, j) = (*this)(k, j);
            (*this)(k, j) = v;
        }
        //erase a(h,i)
l50:
        if ((*this)(h, i) == 0.0)goto l80;
        r = -(*this)(h, i) / (*this)(i, i);
        c(h, i) = -r;
        for (j = i; j <= row; j++)(*this)(h, j) = (*this)(h, j) + r * (*this)(i, j);
l80:
        h++;
        if (h <= row)goto l50;
    }
    //end triangulation
    return c;
}
//************************determinant**********************************

double matrix::det() {
    matrix a1(row, col);
    double d;
    int i, j, h, n1, k, num = 0;
    double v, r;
    a1 = *this;
    n1 = row - 1;
    //triangulation
    for (i = 1; i <= n1; i++) {
        h = i + 1;
        if (a1(i, i) != 0)goto l50;
        //find nonzero element
        for (k = h; k <= row; k++)if (a1(k, i) != 0)goto l30;
        //nonzero element not found
        //                sound_tone(1500,1331,1331);
        //        sound(2000);
        //        delay(500);
        //        nosound();
        printf("nonzero element not found in Gauss elemination");
        exit(1);

        //replace rows
l30:
        for (j = i; j <= row; j++) {
            v = a1(i, j);
            a1(i, j) = a1(k, j);
            a1(k, j) = v;
        }
        num++;
        //erase a(h,i)
l50:
        if (a1(h, i) == 0.0)goto l80;
        r = -a1(h, i) / a1(i, i);
        for (j = i; j <= row; j++)a1(h, j) = a1(h, j) + r * a1(i, j);
l80:
        h++;
        if (h <= row)goto l50;
    }
    //end triangulation
    d = 1;
    for (i = 1; i <= row; i++)d *= a1(i, i);
    d *= pow(-1, num);
    return d;
}

//*********************** norm ***************************************

double matrix::norm() const {
    double a = 0;
    int i;
    for (i = 1; i <= row; i++)a = a + (*this)(i, 1)*(*this)(i, 1);
    return sqrt(a);
}

matrix matrix::column(int n) {
    if (n > col) {
        printf("number column exceed matrix dimension\n");
        exit(1);
    }
    matrix t(row, 1);
    int i;
    for (i = 1; i <= row; i++)t(i, 1) = (*this)(i, n);
    return t;
}

matrix operator*(double c, matrix &a) {
    matrix t(a.row, a.col);
    int j, i;
    for (i = 1; i <= a.row; i++)
        for (j = 1; j <= a.col; j++)
            t(i, j) = c * a(i, j);
    return t;
}

matrix operator*(matrix &a, double c) {
    matrix t(a.row, a.col);
    int j, i;
    for (i = 1; i <= a.row; i++)
        for (j = 1; j <= a.col; j++)
            t(i, j) = c * a(i, j);
    return t;
}

double & matrix::operator()(int n, int m) const {
    if ((n > row) || (m > col) || (n < 1) || (m < 1)) {
        //        sound_tone(1500,1331,1331);
        //        sound(2000);
        //        delay(500);
        //        nosound();
        printf("index exceeds matrix dimension");
        exit(1);
    }
    return y[n - 1][m - 1];
}

double & matrix::operator()(int n) const {
    if ((n > row * col) || (n < 1)) {
        //        sound_tone(1500,1331,1331);
        //        sound(2000);
        //        delay(500);
        //        nosound();
        printf("index exceeds matrix length");
        exit(1);
    }
    int n1, m1;
    m1 = n / row + 1;
    n1 = n % row;
    return y[n1 - 1][m1 - 1];
}

int eigJ(matrix &a, double eps, double *d, matrix &v) {
    int n = a.row, i, j, k, p, q;
    matrix u(n, n), ut(n, n), a1(n, n);
    double max, s, fi, t;
    for (i = 1; i <= n; i++) v(i, i) = u(i, i) = ut(i, i) = 1.0;
    a1 = a;
    t = 1;
    int n1 = n * (n - 1) / 2;
    while (t > eps) {
        max = fabs(a1(1, 2));
        p = 1;
        q = 2;
        for (i = 1; i <= n; i++)
            for (j = 1; j <= n; j++)if (i != j)
                    if (fabs(a1(i, j)) > max) {
                        p = i;
                        q = j;
                        max = fabs(a1(i, j));
                    }
        fi = 0.5 * atan(2 * a1(p, q) / (a1(p, p) - a1(q, q)));
        u(p, p) = u(q, q) = ut(p, p) = ut(q, q) = cos(fi);
        u(p, q) = ut(q, p) = -sin(fi);
        u(q, p) = ut(p, q) = sin(fi);
        a1 = ut * a1*u;
        v = v*u;
        t = 0.0;
        for (i = 1; i <= n; i++)
            for (j = i + 1; j <= n; j++) t += fabs(a1(i, j));
        t /= n1;
        u(p, p) = u(q, q) = ut(p, p) = ut(q, q) = 1.0;
        u(p, q) = ut(q, p) = 0.0;
        u(q, p) = ut(p, q) = 0.0;
    }
    for (i = 0; i < n; i++)d[i] = a1(i + 1, i + 1);
    return 1;
}