task_sem3_3.1

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include "eigen.h"
#include "linsys.h"

#define de 1e-12

double val(int i, int n) {
    return pow((2.0 * (1.0 - cos((M_PI * (2.0 * i - 1.0)) / (2.0 * (double)n + 1.0)))), -1.0);
}

double f(double x, double *matrix, int n) {
    // fill matrix
    double *matrixp;
    matrixp = matrix;
    for (int i = 0; i < n; i++, matrix += 80)
        for (int j = 0; j < n; j++)
            if (i == j)
                matrixp[j] = (double)n - (double)i - x; // diagonal elements
            else
                matrixp[j] = (double)n - (double)j; // others

    // decompose && calculate
    if (decomp(matrix, matrix, n, 80) == FAILURE)
        return 0.0; // we've found solution
    else
        return determ(matrix, n, 80); // otherwise
}

double sign(double x) {
    return (x / fabs(x));
}

double solvex(double a, double b, double *matrix, int n) {
    double c, x;

    do {
        c = a + (b - a) / 2;

        double t = f(c, matrix, n); // get the determinant

        if ((t == 0.0) || (fabs(t) <= de) || ((b - a) <= de)) { // some checks
            x = c;
            break;
        }

        if (sign(f(a, matrix, n)) * sign(t) < 0) // divide interval
            b = c;
        else
            a = c;
    } while (1);

    return x;
}

int main(int argc, char *argv[]) {
    clock_t start;
    double points[80]; // own values
    double jacobip[80]; // Jacobi own values
    double *a, *ap; // the source matrix
    int n; // dimensions

    a = (double *)malloc(sizeof(double) * 80 * 80);

    for (int e = 0; e <= 3; e++) {
        n = 10 * (int)pow(2.0, (double)e);
        printf("n = %d\n**********\n", n);

        // out method
        start = clock();
        for (int i = 0; i < n; i++)
            points[i] = val(i + 1, n);

        double min1 = points[0];
        for (int i = 0; i < n; i++)
            if (points[i] < min1)
                min1 = points[i];

        double min2 = min1;
        for (int i = 0; i < n; i++)
            if ((points[i] < min2) && (points[i] != min1))
                min2 = points[i];

        double delta = fabs(min2 - min1) * 0.95;

        for (int i = 0; i < n; i++)
            /*printf("#%d\texact value = %.16lf\tour value = %.16lf\n", i + 1, points[i], solvex(points[i] - delta, points[i] + delta, a, n));*/
            solvex(points[i] - delta, points[i] + delta, a, n);

        double elapsed = (double)(clock() - start) / CLOCKS_PER_SEC;
        printf("Our method = %lf", elapsed);

        // Jacobi method
        ap = a;
        for (int i = 0; i < n; i++, ap += 80)
            for (int j = 0; j < n; j++)
                ap[j] = (double)n - (double)i;

        start = clock();
        jacobi(a, n, 80, jacobip, NULL, 0);
        elapsed = (double)(clock() - start) / CLOCKS_PER_SEC;
        printf("\tJacobi method = %lf\n", elapsed);

        for (int i = 0; i < n; i++)
            printf("#%d\tjacobi value = %.16lf\tour value = %.16lf\texact value = %.16lf\n", i + 1, jacobip[i], solvex(points[i] - delta, points[i] + delta, a, n), points[i]);

        putchar('\n');
    }

    free(a);

    return 0;
}