#include "factorization_solve_tests.h"
#include "factorization_solve.h"
#include "benchmarks.h"
#include <chrono>
#include <iostream>

using namespace std::chrono; 

void CheckCorrectness()
{
    cout.precision(4);

    size_t N = 5;
    size_t n_rhs = 5;
    vector<double> A(N*N);
    vector<double> b(N*n_rhs);

    // initialize A
    for (size_t i = 0; i < N; ++i)
    {
        for (size_t j = 0; j < N; ++j)
            if (i == j)
                A[N*i + j] = 4.0;
            else
                A[N*i + j] = 1.0/((i - j)*(i - j));
    }
    const vector<double> A_copy = A;
        
    // initialize b as identity matrix
    for (size_t j = 0; j < N; ++j)
    {
        for (size_t l = 0; l < n_rhs; ++l)
            if (l == j)
                b[N*l + j] = 1.0;
            else
                b[N*l + j] = 0.0;
    }
    
    cout << "A = " << endl;
    print_matrix(A, N, N);
    cout << "b = " << endl;
    print_matrix(b, N, n_rhs);

    // solve system
    factorization_solve(A, b, n_rhs);

    vector<double> check_matrix_identity = MatMat_cBLAS(A_copy, b, N);
    cout << "A*A^{-1} = " << endl;
    print_matrix(check_matrix_identity, N, n_rhs);
}


void CheckDuration(const size_t &N)
{
    for (size_t n_rhs : {1, 2, 4, 8, 16, 32})
    {
        auto time_start = system_clock::now();

        vector<double> A(N*N);
        vector<double> b(N*n_rhs);

        // initialize A
        for (size_t i = 0; i < N; ++i)
        {
            for (size_t j = 0; j < N; ++j)
                if (i == j)
                    A[N*i + j] = 4.0;
                else
                    A[N*i + j] = 1.0/((i - j)*(i - j));
        }

        // initialize b
        for (size_t j = 0; j < N; ++j)
            for (size_t l = 0; l < n_rhs; ++l)
                b[N*l + j] = 1.0;


        // solve system
        factorization_solve(A, b, n_rhs);



        auto time_end = system_clock::now();
        auto duration = duration_cast<microseconds>(time_end - time_start);        // duration in microseconds
        double t_diff = static_cast<double>(duration.count()) / 1e6;               // overall duration in seconds

        cout << "n_rhs = " << n_rhs << "\tTime per rhs: " << t_diff/n_rhs << endl;
    }
}