task_sem3_3.2

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

double *a; // main matrix
double *at; // matrix backup
double *v; // vectors matrix

int cmpitems_p(const void *p1, const void *p2) {
    if (*(double *)p1 < *(double *)p2)
        return -1;
    else if (*(double *)p1 > *(double *)p2)
        return 1;

    return 0;
}

double val(int i, int n) {
    double s;

    s = 4 * pow(sin((i + 1.0) * M_PI / (2.0 * n + 2.0)), 2.0);
    return (pow(s, 3.0) - 5.0 * pow(s, 2.0) + 8.0 * s);
}

void getA(int n) {
    double *c = (double *)malloc(sizeof(double) * n * n);
    double *t = (double *)malloc(sizeof(double) * n * n);
    double *t1 = (double *)malloc(sizeof(double) * n * n);

    a = (double *)malloc(sizeof(double) * n * n);
    at = (double *)malloc(sizeof(double) * n * n);
    v = (double *)malloc(sizeof(double) * n * n);

    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++) {
            if (i == j)
                *(c + i * n + j) = 2;
            else if ((j == (i + 1)) || (j == (i - 1)))
                *(c + i * n + j) = 1;
            else
                *(c + i * n + j) = 0;

            *(t + i * n + j) = 0;
            *(t1 + i * n + j) = 0;
        }

    for (int i = 0; i < n; i++) // A = 8C
        for (int j = 0; j < n; j++)
            *(a + i * n + j) = 8 * *(c + i * n + j);

    for (int i = 0; i < n; i++) // T = C^2
        for (int j = 0; j < n; j++)
            for (int x = 0; x < n; x++)
                *(t + i * n + j) += *(c + i * n + x) * *(c + x * n + j);

    for (int i = 0; i < n; i++) // A = 8C - 5C^2
        for (int j = 0; j < n; j++)
            *(a + i * n + j) -= 5 * *(t + i * n + j);

    for (int i = 0; i < n; i++) // T_1 = C^3 = TC^2
        for (int j = 0; j < n; j++)
            for (int x = 0; x < n; x++)
                *(t1 + i * n + j) += *(t + i * n + x) * *(c + x * n + j);

    for (int i = 0; i < n; i++) // A = 8C - 5C^2 + C^3
        for (int j = 0; j < n; j++)
            *(a + i * n + j) += *(t1 + i * n + j);

    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            *(at + i * n + j) = *(a + i * n + j);

    free(c);
    free(t);
    free(t1);
}

int main(int argc, char *argv[]) {
    int n;
    double *vals;
    double *j_vals;
    int *order;
    double *t;
    double *j_vecs;
    double *vecs;

    printf("n = "); scanf("%d", &n);
    getA(n);

    vals = (double *)malloc(sizeof(double) * n);
    j_vals = (double *)malloc(sizeof(double) * n);
    j_vecs = (double *)malloc(sizeof(double) * n * n);
    vecs = (double *)malloc(sizeof(double) * n * n);
    order = (int *)malloc(sizeof(int) * n);
    t = (double *)malloc(sizeof(double) * n);

    for (int i = 0; i < n; i++)
        *(vals + i) = val(i, n);

    clock_t start;
    start = clock();

    jacobi(a, n, n, j_vals, v, 1);

    double elapsed = (double)(clock() - start) / CLOCKS_PER_SEC;

    qsort(vals, n, sizeof(double), cmpitems_p);

    for (int i = 0; i < n; i++)
        *(t + i) = *(j_vals + i);

    qsort(j_vals, n, sizeof(double), cmpitems_p);

    printf("Exact values\t\tJacobi values\n");
    for (int i = 0; i < n; i++)
        printf("%.16lf\t%.16lf\n", *(vals + i), *(j_vals + i));

    putchar('\n');

    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            if (*(t + j) == *(j_vals + i))
                *(order + i) = j;

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

    for (int i = 0; i < n; i++) // Av
        for (int j = 0; j < n; j++)
            for (int x = 0; x < n; x++)
                *(j_vecs + i * n + j) += *(at + j * n + x) * *(v + i * n + x);

    for (int i = 0; i < n; i++) // lambda v
        for (int j = 0; j < n; j++)
            *(vecs + i * n + j) = *(vals + i) * *(v + *(order + i) * n + j);

    printf("A * v\n");

    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            printf("%.16lf%c", *(j_vecs + i * n + j), ((j) == ((n) - 1)) ? '\n' : ' ');

    putchar('\n');

    printf("lambda * v\n");
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            printf("%.16lf%c", *(vecs + i * n + j), ((j) == ((n) - 1)) ? '\n' : ' ');

    putchar('\n');
    printf("Time elapsed = %.16lf\n", elapsed);

    free(a);
    free(at);
    free(v);
    free(vals);
    free(j_vals);
    free(order);
    free(t);
    free(j_vecs);
    free(vecs);

    return 0;
}