TUPO LEGKO

#define _CRT_SECURE_NO_WARNINGS
#include <malloc.h>
#include <string.h>
#include <stdio.h>
#include <math.h>
#include <time.h>

#define EPS 0.000000001

struct matrix
{
        size_t   w;
        size_t   h;
        double** elems;
};

struct vector
{
        size_t  size;
        double* elems;
};

struct matrix* create_matrix(size_t w, size_t h)
{
        struct matrix* m = malloc(sizeof(*m));
        m->elems = malloc(sizeof(double*) * h);
        for (size_t i = 0; i < h; i++)
                m->elems[i] = malloc(sizeof(double) * w);
        m->w = w;
        m->h = h;
        return m;
}

void delete_matrix(struct matrix* m)
{
        for (size_t i = 0; i < m->h; i++)
                free(m->elems[i]);
        free(m);
}

struct vector* create_vector(size_t size)
{
        struct vector* v = malloc(sizeof(*v));
        v->elems = malloc(sizeof(double) * size);
        v->size = size;
        return v;
}

void delete_vector(struct vector* v)
{
        free(v->elems);
        free(v);
}

void print_matrix(const struct matrix* m)
{
        for (size_t i = 0; i < m->h; i++)
        {
                for (size_t j = 0; j < m->w; j++)
                        printf("%.3lf ", m->elems[i][j]);
                printf("\n");
        }
        printf("\n");
}

void print_vector(const struct vector* v)
{
        for (size_t i = 0; i < v->size; i++)
                printf("%.3lf ", v->elems[i]);
        printf("\n\n");
}

struct vector* mult(const struct matrix* m, const struct vector* v)
{
        struct vector* r = create_vector(v->size);
        for (size_t i = 0; i < m->h; i++)
        {
                r->elems[i] = 0.0;
                for (size_t j = 0; j < m->w; j++)
                        r->elems[i] += m->elems[i][j] * v->elems[j];
        }
        return r;
}

void divide(struct vector* v, double val)
{
        for (size_t i = 0; i < v->size; i++)
                v->elems[i] /= val;
}

struct matrix* sub(const struct matrix* m, double val)
{
        struct matrix* r = create_matrix(m->w, m->h);
        for (size_t i = 0; i < m->h; i++)
                for (size_t j = 0; j < m->w; j++)
                        r->elems[i][j] = m->elems[i][j] - val;
        return r;
}

double norm(const struct vector* v)
{
        double sum = 0.0;
        for (size_t i = 0; i < v->size; i++)
                sum += v->elems[i] * v->elems[i];
        return sqrt(sum);
}

struct vector* ones(int size)
{
        struct vector* v = create_vector(size);
        for (int i = 0; i < size; i++)
                v->elems[i] = 1.0;
        return v;
}

struct matrix* hilb(int size)
{
        struct matrix* m = create_matrix(size, size);
        for (int i = 0; i < size; i++)
                for (int j = 0; j < size; j++)
                        m->elems[i][j] = 1.0 / (i + j + 1);
        return m;
}

struct vector* hilbv(int size)
{
        struct vector* v = create_vector(size);
        for (int i = 0; i < size; i++)
        {
                v->elems[i] = 0;
                for (int j = 0; j < size; j++)
                        v->elems[i] += 1.0 / (i + j + 1);
        }
        return v;
}

struct vector* solve_linear_equation(const struct matrix* m, const struct vector* v)
{
        if (m->h != v->size || m->w != v->size)
                return NULL;
        struct matrix* aug = create_matrix(m->w + 1, m->h);
        size_t h = aug->h;
        size_t w = aug->w;
        double** elems = aug->elems;
        for (size_t i = 0; i < h; i++)
        {
                memcpy(aug->elems[i], m->elems[i], sizeof(double) * m->w);
                elems[i][w - 1] = v->elems[i];
        }
        for (size_t i = 0; i < h - 1; i++)
        {
                size_t max = i;
                for (size_t j = i; j < h; j++)
                        if (fabs(elems[j][i]) > fabs(elems[max][i]))
                                max = j;
                double* t = elems[i];
                elems[i] = elems[max];
                elems[max] = t;
                for (size_t j = i + 1; j < h; j++)
                {
                        double mult = -elems[j][i] / elems[i][i];
                        for (size_t k = 0; k < w; k++)
                                aug->elems[j][k] += elems[i][k] * mult;
                }
        }
        struct vector* res = create_vector(h);
        res->elems[h - 1] = elems[h - 1][w - 1] / elems[h - 1][w - 2];
        for (size_t i = h - 2; i != -1; i--)
        {
                double sum = 0.0;
                for (size_t j = i + 1; j < w - 1; j++)
                        sum += aug->elems[i][j] * res->elems[j];
                res->elems[i] = (elems[i][w - 1] - sum) / aug->elems[i][i];
        }
        delete_matrix(aug);
        return res;
}

double min_eigen_value(const struct matrix* m)
{
        struct vector* cur_guess = ones(m->w);
        for (size_t i = 0; i < 50; i++)
        {
                struct vector* prev = cur_guess;
                cur_guess = solve_linear_equation(m, cur_guess);
                delete_vector(prev);
                divide(cur_guess, norm(cur_guess));
        }
        struct vector* prev = cur_guess;
        cur_guess = mult(m, cur_guess);
        delete_vector(prev);
        return norm(cur_guess);
}

int main()
{
        struct matrix* m = create_matrix(2, 2);


        m->elems[0][0] = 1;
        m->elems[0][1] = 2;
        m->elems[1][0] = 3;
        m->elems[1][1] = 4;

        print_matrix(m);
        double min = min_eigen_value(m);
        printf("min: %.5lf", min);
        getchar();
        return 0;
}