#include <iostream>
#include <sstream>
#include <vector>
#include <cmath>
#include <lapacke.h>
using namespace std;

int n = 10;
double threshold = 1.0/sqrt(2);
double lambda(double x)
{
    if(x < threshold)
    {
        return 1;
    }
    return 10;
}

double fillOffDiagonal(unsigned int i)
{
    
    double x_l = i/double(n);
    double x_u = (i+1)/double(n);

    if(x_l < threshold && x_u > threshold)
    {
        return -n*n*(threshold-x_l+10*(x_u-threshold));
    }
    return -n*lambda(x_l);
}

double fillDiagonal(int i)
{
    double x_l = (i-1)/double(n);
    double x_u = (i+1)/double(n);
    if(x_l < threshold && x_u > threshold)
    {
        return n*n*(threshold-x_l+10*(x_u-threshold));
    }
    return 2*n*lambda(x_l);
}

int main()
{   
    //Ignored the first Row of K cause we now that u(0)=0
    
    int rhs = 1;
    unsigned int nodes = n+1;


    //Lapacke overwrites upper/lower diagonal so store twice 
    vector<double> upperdiagonal(n,0.0);
    vector<double> lowerdiagonal(n,0.0);
    vector<double> diagonal(nodes,1.0);
    vector<double> F(nodes,0.0);

    for(int i=0; i < n; i++)
    {
        upperdiagonal[i] = fillOffDiagonal(i);
        lowerdiagonal[i] = fillOffDiagonal(i);
        diagonal[i+1] = fillDiagonal(i+1);
    }

    upperdiagonal[0]=0.0;
    lowerdiagonal[n-1]=0.0;
    diagonal[n] = 1.0;
    F[n] = 1.0;
    

    //Tridiagonal
    LAPACKE_dgtsv(LAPACK_COL_MAJOR,nodes,rhs,lowerdiagonal.data(),diagonal.data(),upperdiagonal.data(),F.data(),nodes);
    

    cout << "Solution u_h at nodes x_0..x_" << n << ":\n";
    
    for (unsigned int i = 0; i < nodes; ++i) {

        double xi = double(i) /double(n);
        cout << "u_h(" << xi << ") = " << F[i] << "\n";
        // F will be overwritten with sol
    }
    return 0;

    
}