Math.lab2

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define MAXS 3

int read_matrix(char *filename, float matrix[MAXS][MAXS+1]);
void print_matrix(float matrix[MAXS][MAXS+1], int size);
void diagonal_dominate(float matrix[MAXS][MAXS+1], int size);
void cnt(float matrix[MAXS][MAXS+1], int size);

int main(int argc, char **argv)
{
    char *filename = "input.txt";
    float matrix[MAXS][MAXS+1];
    int size;

    printf("\tJacob's method\n");
    size = read_matrix(filename, matrix);
    printf("\n\tOriginal matrix:\n");
    print_matrix(matrix, size);
    diagonal_dominate(matrix, size);
    printf("\tAfter making diagonal elements dominate:\n");
    print_matrix(matrix, size);

    cnt(matrix, size);
    //print_matrix(matrix, size);

    return 0;
}

void diagonal_dominate(float matrix[MAXS][MAXS+1], int size)
{
    int j;

    for (j = 0; j < size+1; j++) {
        matrix[0][j] += matrix[2][j];
        matrix[2][j] -= matrix[0][j];
        matrix[0][j] *= 2;
        matrix[1][j] -= matrix[0][j];

    }
}

void cnt(float matrix[MAXS][MAXS+1], int size)
{
    int i, j, k;
    float factor, pure[MAXS][MAXS+1], B[MAXS], C[MAXS], X[MAXS], norm_B = 0.0, norm_C = 0.0;

    for (j = 0; j < size; j++) {
        B[j] = 0.0;
        C[j] = matrix[j][size];
    }

    for (i = 0; i < size; i++)
        for (j = 0; j < size+1; j++)
            pure[i][j] = matrix[i][j];

    for (i = 0; i < size; i++) {
        factor = matrix[i][i];
        for (j = 0; j < size+1; j++) {
            matrix[i][j] /= factor;
            if ((i != j) && (j != size)) {
                B[i] += fabsf(matrix[i][j]);
                if (norm_B < B[i])
                    norm_B = B[i];
            }
            if ((j == size) && (norm_C < matrix[i][j]))
                norm_C = matrix[i][j];
        }
    }
    k = abs(log10(0.001) + log10(0.12) - log10(2.16) - log10(0.88) - 1);


    printf("\tAfter dividing elements by main diagonal's factor:\n");
    print_matrix(matrix, size);
    printf("\t||B|| = %.2f\n\n\t||C|| = %.2f\n\n\tepsilon = 0.001\n\n", norm_B, norm_C);
    printf("\tk > %d\n\n", k);

    for (j = 0; j < size; j++)
        X[j] = C[j];
    printf("\tk =  0    x1 = %8.4f    x2 = %8.4f    x3 = %8.4f\n", X[0], X[1], X[2]);

    for (i = 1; i <= k+3; i++) {
        if (i == k+2)
            printf("\t---------------------------------------------------------\n");
        X[0] = -1/pure[0][0] * ( (matrix[0][1] + matrix[0][2]) * X[0] + C[0] );
        X[1] = -1/pure[1][1] * ( (matrix[1][0] + matrix[1][2]) * X[1] + C[1] );
        X[2] = -1/pure[2][2] * ( (matrix[2][0] + matrix[2][1]) * X[2] + C[2] );
        printf("\tk = %2d    x1 = %8.4f    x2 = %8.4f    x3 = %8.4f\n", i, X[0], X[1], X[2]);
    }

}

int read_matrix(char *filename, float matrix[MAXS][MAXS+1])
{
    FILE *ifp;
    int i, j, size = 0;

    ifp = fopen(filename, "r");
    fscanf(ifp, "%d", &size);
    for (i = 0; i < size; i++)
        for (j = 0; j < size+1; j++)
            fscanf(ifp, "%f", &matrix[i][j]);
    fclose(ifp);

    return size;
}

void print_matrix(float matrix[MAXS][MAXS+1], int size)
{
    int i, j;

    putchar('\n');
    for (i = 0; i < size; i++) {
        for (j = 0; j < size+1; j++)
            printf("  %10.2f",matrix[i][j]);
        putchar('\n');
    }
    putchar('\n');
}