#include "factorization_solve.h"
#include <vector>
#include <math.h>
#include <assert.h>
#include <lapack.h>
#include <iostream>

using namespace std;


void factorization_solve(vector<double> &A, vector<double> &b, const int32_t &n_rhs)
{
    int32_t nelem = A.size();
    int32_t N = sqrt(nelem);
    assert(N == static_cast<int32_t>(b.size())/n_rhs);

    vector<int> IPIV(N); // pivot indices
    int info_factorization;

    dgetrf_(&N, &N, A.data(), &N, IPIV.data(), &info_factorization);

    const char transp = 'N';
    int info_solve;
    dgetrs_(&transp, &N, &n_rhs, A.data(), &N, IPIV.data(), b.data(), &N, &info_solve, 1); // 1 is length of parameter 'N'

}

void print_matrix(vector<double> &A, size_t M, size_t N)
{
    for (size_t i = 0; i < M; ++i)
    {
        for (size_t j = 0; j < N; ++j)
        {
            cout << A[N*i + j] << "\t";
        }
        cout << endl;
    }
    cout << endl;
}