ChM_03

#include <stdio.h>  //(+0.65 -0.93 +0.45 | -0.72)
#include <stdlib.h> //(+1.15 +0.43 -0.72 | +1.24)
#include <math.h>   //(+0.56 -0.18 +1.03 | +2.15)
#include <conio.h>

double ** CreateMassA(const int SLAR, double **massA);
double * CreateMassB(const int SLAR, double *massB);
void CreateMass1(const int SLAR, double **mass);
double ** CreateMass2(const int SLAR, double **mass);
void swap (double *a, double *b);
void temp (const int SLAR, double **massA, double **tempA, const double *massB, double * tempB);
void PrintMass(const int SLAR, double **massA, const double *massB);
void PrintLU (const int SLAR, double **massL, double **massU);
void PrintX(const int SLAR, const double *massX);
void Gauss(const int SLAR, double **massA, double *massB);
void LU(const int SLAR, double **massA, double *massB);

int main()
{
    double print;
    int SLAR;
    do
    {
        printf ("Enter the number of equations: ");
        scanf_s ("%lf", &print);

    } while (print < 1);
    SLAR = print;
    double **massA = CreateMassA(SLAR, massA);
    double *massB = CreateMassB(SLAR, massB);
    double **tempA = CreateMass2(SLAR, tempA);
    double *tempB;
    CreateMass1 (SLAR, &tempB);
    temp (SLAR, massA, tempA, massB, tempB);
    Gauss(SLAR, massA, massB);
    temp (SLAR, tempA, massA, tempB, massB);
    LU(SLAR, massA, massB);
    for (int i = 0; i < SLAR; ++i)
        free (massA[i]);
    free (massA);
    for (int i = 0; i < SLAR; ++i)
        free (tempA[i]);
    free (tempA);
    tempA = NULL;
    free (massB);
    massB = NULL;
    free (tempB);
    tempB = NULL;
    system ("pause");
    return 0;
}

double ** CreateMassA(const int SLAR, double **massA)
{
    massA = (double **) malloc (SLAR*sizeof (double *));
    if (!massA)
    {
        printf ("Not enough memory for massA!\n");
        system ("pause");
        exit (1);
    }
    for (int i = 0; i < SLAR; ++i)
    {
        massA[i] = (double *)malloc (SLAR*sizeof (double));
        if (!massA[i])
        {
            printf ("Not enough memory for massA[%d]!\n", i+1);
            system ("pause");
            exit (1);
        }
        for (int j = 0; j < SLAR; ++j)
        {
            printf ("Enter the a[%i][%i]: ", i+1, j+1);
            scanf_s ("%lf", &massA[i][j]);
        }
    }
    return massA;
}
double * CreateMassB(const int SLAR, double *massB)
{
    massB = (double *) malloc (SLAR*sizeof (double));
    if(!massB)
    {
        printf ("Not enough memory for massB!\n");
        system ("pause");
        exit (1);
    }
    for (int i = 0; i < SLAR; ++i)
    {
        printf ("Enter the b[%d]: ", i + 1);
        scanf ("%lf", &massB[i]);
    }
    return massB;
}

double ** CreateMass2(const int SLAR, double **mass)
{
    mass = (double **) malloc (SLAR*sizeof (double *));
    if (!mass)
    {
        printf ("Not enough memory for mass[][]!\n");
        system ("pause");
        exit (1);
    }
    for (int i = 0; i < SLAR; ++i)
    {
        mass[i] = (double *)malloc (SLAR*sizeof (double));
        if (!mass[i])
        {
            printf ("Not enough memory for mass[%d]!\n", i+1);
            system ("pause");
            exit (1);
        }
    }
    return mass;
}
void CreateMass1(const int SLAR, double **mass)
{
    *mass = (double *) malloc (SLAR*sizeof (double));
    if(!mass)
    {
        printf ("Not enough memory for mass[]!\n");
        system ("pause");
        exit (1);
    }
}

void swap (double *a, double *b)
{
    double temp = *a;
    *a = *b;
    *b = temp;
}

void temp (const int SLAR, double **massA, double **tempA, const double *massB, double * tempB)
{
    for (int i = 0; i < SLAR; ++i)
    {
        for (int j = 0; j < SLAR; ++j)
            tempA[i][j] = massA[i][j];
        tempB[i] = massB[i];
    }
}

void PrintMass(const int SLAR, double **massA, const double *massB)
{
    int k = (SLAR+1)*8+5;
    for (int i = 0; i < k; ++i)
        putchar ('=');
    putchar('\n');
    for (int i = 0; i < SLAR; ++i)
    {
        putchar ('(');
        for (int j = 0; j < SLAR; ++j)
            printf (" %+7.2lf", massA[i][j]);
        printf (" | %+7.2lf )\n", massB[i]);
    }
}

void PrintLU (const int SLAR, double **massL, double **massU)
{
    int k = SLAR*8+6;
    int m = SLAR/2;
    for (int i = 0; i < k; ++i)
        putchar ('=');
    putchar('\n');
    for (int i = 0; i < SLAR; ++i)
    {
        if (i == m) printf ("L = (");
        else printf ("    (");
        for (int j = 0; j < SLAR; ++j)
            printf (" %+7.2lf", massL[i][j]);
        putchar(')'); putchar('\n');
    }
    for (int i = 0; i < k; ++i)
        putchar ('-');
    putchar('\n');
    for (int i = 0; i < SLAR; ++i)
    {
        if (i == m) printf ("U = (");
        else printf ("    (");
        for (int j = 0; j < SLAR; ++j)
            printf (" %+7.2lf", massU[i][j]);
        putchar(')'); putchar('\n');
    }
}

void PrintX(const int SLAR, const double *massX)
{
    putchar ('\n');
    for (int i = 0; i < SLAR; ++i)
        printf ("x[%d] = %+.2lf\n", i+1, massX[i]);
    putchar ('\n');
}

void Gauss(const int SLAR, double **massA, double *massB)
{
    double *massX;
    CreateMass1(SLAR, &massX);
    printf ("Press 1 if you want to see step output of method of Gauss\n");
    char ch = _getch();
    if (ch == '1') PrintMass(SLAR, massA, massB);
    for (int k = 0; k < SLAR - 1; ++k)  //k - номер рівняння, над яким працюємо
    {
        int check = 0;
        double MAX = fabs(massA[k][k]);
        for (int i = k + 1; i < SLAR; ++i)  //і - номер рівняння, що іде після k-го рівняння
            if (MAX < fabs(massA[i][k]))    //визначаємо найбільше за модулем число у стовпці
            {
                check = 1;
                for (int a = 0; a < SLAR; ++a) //а - номер елемента в рівнянні
                    swap (&massA[k][a], &massA[i][a]);
                swap (&massB[k], &massB[i]);
            }
        if (MAX == 0)
        {
            for (int j = k; j < SLAR-1; ++j)
            {
                int check = 0;
                for (int i = 0; i < SLAR - 1; ++i)
                {
                    if (massA[j][i] != 0)
                    {
                        check = 1; continue;
                    }
                }
                if (massB[j] != 0)
                {
                    printf ("\nERROR!  System is indefinite\n", k+1);
                    system ("pause");
                    exit(1);
                }
            }
            printf ("\nERROR!  System is compatible\n", k+1);
            system ("pause");
            exit(1);
        }
        if ((ch == '1') && (check == 1)) PrintMass(SLAR, massA, massB);
        for (int j = k + 1; j < SLAR; ++j) //j - номер рівняння, над яким працюємо
        {                                  //занулюємо елементи під головною діагоналлю
            double a = massA[j][k] / massA[k][k]; // a - число, на яке ділиться рівняння
            for (int n = k; n < SLAR; ++n) //n - номер елемента в рівнянні
            {
                massA[j][n] = massA[j][n]-a*massA[k][n];
            }
            massB[j] = massB[j] - a*massB[k];
        }
        if (ch == '1') PrintMass(SLAR, massA, massB);
    }

    if (massA[SLAR-1][SLAR-1] == 0)
    {
        if (massB[SLAR-1] != 0) printf ("\nERROR! System is indefinite!\n");
        else printf ("\nERROR! System is compatible!\n");
        system ("pause");
        exit (1);
    }

    for (int i = SLAR-1; i >= 0; --i) //визначаємо корені рівняння
    {
        double x = 0;
        for (int j = SLAR-1; j > i; --j)
            x = x + massA[i][j]*massX[j];
        massX[i] = (massB[i] - x) / massA[i][i];
    }
    PrintX(SLAR, massX);
    free (massX);
    massX = NULL;
}

void LU(const int SLAR, double **massA, double *massB)
{
    double *massX;
    double **massL = CreateMass2(SLAR, massL);
    double **massU = CreateMass2(SLAR, massU);
    CreateMass1(SLAR, &massX);
    for (int i = 0; i < SLAR; ++i)
        for (int j = 0; j < SLAR; ++j)
        {
            if (i == j) {massU[i][j] = 1; massL[i][j] = 0;}
            else {massL[i][j] = 0; massU[i][j] = 0;}
        }
    massL[0][0] = massA[0][0];
    for (int i = 1; i < SLAR; ++i)
    {
        massU[0][i] = massA[0][i] / massL[0][0];
        massL[i][0] = massA[i][0];  //+2;1;1|1
    }                               //+4;1;0|-2
    for (int i = 1; i < SLAR; ++i)  //-2;2;1|7
    {
        massL[i][0] = massA[i][0];
        for (int k = 1; k < SLAR; ++k)
        {
            if (k <= i)
            {
                for (int m = 0; m < k; ++m)
                    massL[i][k] += (massL[i][m]*massU[m][k]);
                massL[i][k] = massA[i][k] - massL[i][k];
            }
            else
            {
                for (int m = 0; m < i; ++m)
                    massU[i][k] += (massL[i][m]*massU[m][k]);
            massU[i][k] = (massA[i][k] - massU[i][k]) / massL[i][i];
            }
        }
    }
    printf ("Press 1 if you want to see step output of LU-method\n");
    char ch = _getch();
    if (ch == '1') PrintLU (SLAR, massL, massU);
    for (int i = 0; i < SLAR; ++i) //визначаємо massY
    {
        double x = 0;
        for (int j = 0; j < i; ++j)
            x += massL[i][j]*massX[j];
        massX[i] = (massB[i] - x) / massL[i][i];
    }
    for (int i = SLAR-1; i >= 0; --i) //визначаємо корені рівняння
    {
        double x = 0;
        for (int j = SLAR-1; j > i; --j)
            x += massU[i][j]*massX[j];
        massX[i] = massX[i] - x;
    }
    PrintX(SLAR, massX);
    for (int i = 0; i < SLAR; ++i)
        free (massL[i]);
    free (massL);
    massL = NULL;
    for (int i = 0; i < SLAR; ++i)
        free (massU[i]);
    free (massU);
    massU = NULL;
    free (massX);
    massX = NULL;
}