interpolation1d

done single-threded 1d interpolations; linear, rational, cubic spline, polynomial and barycentric. Still not done test functions yet. Still missing multi-core options.
2025-09-18 20:18:27 +02:00
parent 92437e5ef1
commit 3a53b6ebf7
29 changed files with 982 additions and 15 deletions
@@ -0,0 +1,40 @@
+#pragma once
+
+#include "utils/vector.h"
+#include "modules/grid1d.h"
+
+
+namespace fvm {
+
+	template <typename T>
+
+	struct Field1D{
+
+	    const Grid1D* grid = nullptr;   // not owning
+	    utils::Vector<T> u; 			// size = grid->N
+		                          
+		Field1D() = default;
+
+		explicit Field1D(const Grid1D& g, double init = 0.0) : grid(&g), u(g.N){
+			
+		}
+
+		void generate_vertices(){
+			vertices.resize(N_vertices);
+			vertices[0] = centers[0] - ((centers[1] - centers[0])*0.5);	
+			vertices[N_vertices-1] = centers[N_centers-1] + ((centers[N_centers-1] - centers[N_centers-2])*0.5);
+			for (uint64_t i = 1; i < N_centers; ++i){
+				vertices[i] = (centers[i-1] + centers[i])*0.5;
+			}
+		}
+
+		T dx(uint64_t i) const { check(i); return vertices(i+1) - vertices(i); }
+		T center(uint64_t i) const { check(i); return centers(i); }
+private:
+	    void check(uint64_t i) const {
+	        if (i >= N_centers) throw std::runtime_error("Grid1D: cell index out of range");
+	    }
+	};
+
+
+}
@@ -0,0 +1,58 @@
+#pragma once
+
+#include "modules/grid1d.h"
+
+
+namespace fd1d {
+
+	// -----------------------------------------------------------------------------
+	// Second derivative (u_xx) at interior cell i, central difference
+	// Works on NON-uniform grids
+	// On uniform: (u[i-1] - 2 u[i] + u[i+1]) / dx^2
+	// -----------------------------------------------------------------------------
+	template <typename T>
+	void inplace_Build_CentralDerivative_Matrix(const fvm::Grid1D<T>& g, utils::Matrix<T>& A, utils::Vector<T>& b, const utils::Vector<T>& s, const T& c){
+
+		for (uint64_t i = 1; i < g.center_idx; ++i){
+
+			A(i,i-1) = -(c/(g.centers[i] - g.centers[i-1]));
+			A(i,i) = -((c/(g.centers[i+1] - g.centers[i])) + (c/(g.centers[i] - g.centers[i-1])));
+			A(i,i+1) = -(c/(g.centers[i+1] - g.centers[i]));
+			b[i] = -s[i]*(g.vertices[i+1] - g.vertices[i]);
+		}
+    }
+	template <typename T>
+	utils::Matrix<T> Build_CentralDerivative_Matrix(const fvm::Grid1D<T>& g, utils::Vector<T>& b, const utils::Vector<T>& s, const T& c){
+		utils::Matrix<T> A(g.center_idx+1, g.center_idx+1, T{0});
+		inplace_Build_CentralDerivative_Matrix(g, A, b, s, c);
+		return A;
+    }
+
+
+	template <typename T>
+	void inplace_BC_Dirichlet(const fvm::Grid1D<T>& g, utils::Matrix<T>& A, utils::Vector<T>& b, const utils::Vector<T>& s, const T& c){
+
+		A(0,0) = -((c/(g.centers[1] - g.centers[0])) + (c/(g.centers[0] - g.vertices[0])));
+		A(0,1) = c/(g.centers[1] - g.centers[0]);
+
+		A(g.center_idx, g.center_idx-1) = c/(g.centers[g.center_idx]-g.centers[g.center_idx-1]);
+		A(g.center_idx, g.center_idx) = -((c/(g.vertices[g.vertices_idx] - g.centers[g.center_idx])) + (c/(g.centers[g.center_idx] - g.centers[g.center_idx-1])));
+    }
+
+
+	template <typename T>
+	void inplace_BC_Neumann(const fvm::Grid1D<T>& g, utils::Matrix<T>& A, const T& c){
+
+		A(0,0) = -c/(g.centers[1]-g.centers[0]);
+		A(0,1) = c/(g.centers[1]-g.centers[0]);
+
+		A(g.center_idx, g.center_idx-1) = c/(g.centers[g.center_idx]-g.centers[g.center_idx-1]);
+		A(g.center_idx, g.center_idx) = -c/(g.centers[g.center_idx]-g.centers[g.center_idx-1]);
+    }
+
+
+
+
+
+
+}
@@ -0,0 +1,43 @@
+#pragma once
+
+#include "utils/vector.h"
+
+
+namespace fvm {
+
+	template <typename T>
+
+	struct Grid1D{
+		uint64_t center_idx;				// max cell index
+		uint64_t vertices_idx;				// max vertice index
+		utils::Vector<T> centers;			// size N (unknowns at cell centers)
+		utils::Vector<T> vertices;			// size N+1 (face positions)
+
+		Grid1D() = default;
+
+		explicit Grid1D(const utils::Vector<T>& midpoints){
+			centers = midpoints;
+			center_idx = centers.size()-1;
+			vertices_idx = centers.size();		
+		}
+
+		void generate_vertices(T left, T right){
+			vertices.resize(vertices_idx+1);
+			vertices[0] = left;	
+			vertices[vertices_idx] = right;
+			for (uint64_t i = 1; i < center_idx+1; ++i){
+				vertices[i] = (centers[i-1] + centers[i])*0.5;
+			}
+
+		}
+
+		T dx(uint64_t i) const { check(i); return vertices[i+1] - vertices[i]; }
+		T center(uint64_t i) const { check(i); return centers[i]; }
+
+	    void check(uint64_t i) const {
+	        if (i > center_idx) throw std::runtime_error("Grid1D: cell index out of range");
+	    }
+	};
+
+
+}
@@ -0,0 +1,23 @@
+#pragma once
+
+
+#include "./utils/vector.h"
+#include "./utils/matrix.h"
+
+
+namespace numerics{
+	
+	template <typename T>
+	T abs(const T a){
+
+		if(a < 0){
+			return -a;
+		}else{
+			return a;
+		}
+	}
+
+
+
+} // namespace numerics
+
@@ -0,0 +1,88 @@
+#pragma once
+
+#include "./numerics/interpolation1d_base.h"
+
+#include "./utils/vector.h"
+#include "./numerics/min.h"
+#include "./numerics/max.h"
+
+
+namespace numerics{
+	
+	template <typename T>
+	struct interp_barycentric : Base_interp<T> {
+		using Base   = Base_interp<T>;
+
+	  	// bring base data members into scope (or use this->xx / this->yy below)
+	  	using Base::xx;
+	  	using Base::yy;
+	  	using Base::n;
+
+	  	utils::Vector<T> w;
+	  	int64_t d;
+
+		interp_barycentric(const utils::Vector<T> &xv, const utils::Vector<T> &yv, uint64_t dd)
+		: Base_interp<T>(xv, &yv[0], xv.size()), w(n,T{0}), d(dd) {
+			// Constructor arguments are x and y vectors of length n, and order d of desired approximation.
+			if (n <= d){
+				throw std::invalid_argument("d too large for number of points in interp_barycentric");
+			}
+			for (int64_t k = 0; k < n; ++k){
+				int64_t imin = numerics::max(k-d, static_cast<int64_t>(0));
+				int64_t imax;
+				if (k >= n - d) {
+				    imax = n - d - 1;
+				} else {
+				    imax = k;
+				}
+				T temp;
+				if ( (imin & 1) != 0 ) {   // odd?
+				    temp = T{-1};
+				} else {                   // even
+				    temp = T{1};
+				}
+				T sum = T{0};
+
+				for (int64_t i = imin; i <= imax; ++i){
+					int64_t jmax = numerics::min(i+d, n-1);
+					T term = T{1};
+					for (int64_t j = i; j <= jmax; ++j){
+						if (j == k){
+							continue;
+						}
+						term *= (xx[k] - xx[j]);
+					}
+					term = temp/term;
+					temp = -temp;
+					sum += term;
+				}
+				w[k] = sum;
+			}
+		}
+
+		T rawinterp(int64_t jl, T x) override{
+
+			T num{T{0}}, den{T{0}};
+			
+			for (int64_t i = 0; i < n; ++i){
+				T h = x - xx[i];
+				if (h == T{0}){
+					return yy[i];
+				}else{
+					T temp = w[i]/h;
+					num += temp*yy[i];
+					den += temp;
+				}
+			}
+			return num/den;
+		}
+
+
+		T interp(T x) {
+			return rawinterp(1, x);
+		}
+
+	};
+
+} // namespace numerics
+
@@ -0,0 +1,86 @@
+#pragma once
+
+#include "./numerics/interpolation1d_base.h"
+
+//#include "./numerics/abs.h"
+#include "./utils/vector.h"
+
+
+namespace numerics{
+	
+	template <typename T>
+	struct interp_cubic_spline : Base_interp<T> {
+		using Base   = Base_interp<T>;
+	  	// bring base data members into scope (or use this->xx / this->yy below)
+	  	using Base::xx;
+	  	using Base::yy;
+	  	//using Base::mm;
+
+	  	utils::Vector<T> y2;
+
+
+		interp_cubic_spline(utils::Vector<T> &xv, utils::Vector<T> &yv, T yp1=T{1.e99}, T ypn=T{1.e99})
+		: Base_interp<T>(xv, &yv[0], 2), y2(xv.size(),T{0}) {
+			sety2(&xv[0], &yv[0], yp1, ypn);
+		}
+
+		interp_cubic_spline(utils::Vector<T> &xv, const T *yv, T yp1=T{1.e99}, T ypn=T{1.e99})
+		: Base_interp<T>(xv, yv, 2), y2(xv.size(),T{0}) {
+			sety2(&xv[0], yv, yp1, ypn);
+		}
+
+
+		void sety2(const T *xv, const T *yv, T yp1, T ypn){
+
+			T p, qn, sig, un;
+			uint64_t n = y2.size();
+			utils::Vector<T> u(n-1, T{0});
+
+			if (yp1 > static_cast<T>(0.99e99)){						// The lower boundary condition is set either to be “natural”
+				y2[0] = u[0] = T{0};
+			}else{													// or else to have a speciﬁed ﬁrst derivative.
+				y2[0] = T{-0.5};
+				u[0] = (3.0/(xv[1]-xv[0]))*(((yv[1]-yv[0])/(xv[1]-xv[0]))-yp1);
+			}
+			for (uint64_t i = 1; i < n-1; ++i){						// This is the decomposition loop of the tridiagonal algorithm
+				sig = (xv[i]-xv[i-1])/(xv[i+1]-xv[i-1]);
+				p = sig*y2[i-1]+T{2};
+				y2[i] = (sig - T{1})/p;								// y2 and u are used for temporary storage of the decomposed factors.
+				u[i]=((yv[i+1]-yv[i])/(xv[i+1]-xv[i])) - ((yv[i]-yv[i-1])/(xv[i]-xv[i-1]));
+				u[i]=((T{6}*u[i]/(xv[i+1]-xv[i-1])) - sig*u[i-1])/p;
+			}
+			if (ypn > static_cast<T>(0.99e99)){						// The upper boundary condition is set either to be “natural”
+				qn = un = T{0};
+			}else{													// or else to have a speciﬁed ﬁrst derivative.
+				qn = T{0.5};
+				un = (T{3}/(xv[n-1]-xv[n-2]))*(ypn-((yv[n-1]-yv[n-2])/(xv[n-1]-xv[n-2])));
+			}
+			y2[n-1] = (un-(qn*u[n-2]))/((qn*y2[n-2])+T{1});
+			for (int64_t k = n-2; k >= 0; --k){
+				y2[k] = y2[k] * y2[k+1]+u[k];
+			}
+		}
+
+
+		T rawinterp(int64_t jl, T x) override{
+
+			int64_t klo=jl, khi=jl+1;
+			T y, h, b, a;
+
+			h = xx[khi] - xx[klo];
+			if (h == T{0}){											// The xa’s must be distinct.
+				throw std::invalid_argument("interp_cubic_spline: Bad input to routine splint");
+			}
+
+			a = (xx[khi] - x)/h;									// Cubic spline polynomial is now evaluated.
+			b = (x - xx[klo])/h;
+			y = a*yy[klo] + b*yy[khi] + ( ((a*a*a) - a)*y2[klo] + ((b*b*b) - b)*y2[khi] ) * (h*h) / T{6};
+
+
+			return y;
+		}
+
+	};
+
+} // namespace numerics
+
@@ -0,0 +1,29 @@
+#pragma once
+
+#include "./numerics/interpolation1d_base.h"
+
+
+namespace numerics{
+	
+	template <typename T>
+	struct interp_linear : Base_interp<T> {
+		using Base   = Base_interp<T>;
+	  	// bring base data members into scope (or use this->xx / this->yy below)
+	  	using Base::xx;
+	  	using Base::yy;
+
+
+		interp_linear(const utils::Vector<T> &xv, const utils::Vector<T> &yv): Base_interp<T>(xv, &yv[0], 2){}
+
+		T rawinterp(int64_t j, T x) override{
+			if (xx[j]==xx[j+1]){
+				return yy[j]; 			// Table is defective, but we can recover.
+			}else {
+				return (yy[j] + ((x-xx[j])/(xx[j+1]-xx[j]))*(yy[j+1]-yy[j]));
+			}
+		}
+
+	};
+
+} // namespace numerics
+
@@ -0,0 +1,75 @@
+#pragma once
+
+#include "./numerics/interpolation1d_base.h"
+
+#include "./numerics/abs.h"
+#include "./utils/vector.h"
+
+
+namespace numerics{
+	
+	template <typename T>
+	struct interp_polynomial : Base_interp<T> {
+		using Base   = Base_interp<T>;
+	  	// bring base data members into scope (or use this->xx / this->yy below)
+	  	using Base::xx;
+	  	using Base::yy;
+	  	using Base::mm;
+
+	  	T dy;
+
+
+		interp_polynomial(const utils::Vector<T> &xv, const utils::Vector<T> &yv, uint64_t m)
+		: Base_interp<T>(xv, &yv[0], m), dy(T{0}){}
+
+		T rawinterp(int64_t jl, T x) override{
+
+			int64_t ns=0;
+			T y, den, dif, dift, ho, hp, w;
+			const T *xa = &xx[jl], *ya = &yy[jl];
+			utils::Vector<T> c(mm,0), d(mm,0);
+			dif = numerics::abs(x-xa[0]);
+
+			for (int64_t i = 0; i < mm; ++i){					// Here we ﬁnd the index ns of the closest table entry,
+				dift = numerics::abs(x-xa[i]);
+				if (dift < dif){
+					ns = i;
+					dif=dift;
+				}
+				c[i]=ya[i];									// and initialize the tableau of c’s and d’s.
+				d[i]=ya[i];
+			}
+			y = ya[ns];										// This is the initial approximation to y.
+			ns -= 1;
+
+			for (int64_t m = 1; m < mm; ++m){				// For each column of the tableau,
+				for (int64_t i = 0; i < mm-m; ++i){		// we loop over the current c’s and d’s and update them.
+					ho = xa[i]-x;
+					hp = xa[i+m]-x;
+					w = c[i+1]-d[i];
+					den = ho-hp;
+					if (den == T{0.0}){
+						throw std::invalid_argument("interp_polynomial error"); // This error can occur only if two input xa’s are (to within roundoff identical.
+					}
+					den = w/den;							// Here the c’s and d’s are updated.
+					d[i] = hp*den;
+					c[i] = ho*den;
+				}
+				bool take_left = 2 * (ns + 1) < (mm - m);
+
+				if (take_left) {
+				    dy = c[ns + 1];
+				    y += dy;
+				} else {
+				    dy = d[ns];   
+				    y += dy;
+				    ns -= 1;
+				}
+			}
+			return y;
+		}
+
+	};
+
+} // namespace numerics
+
@@ -0,0 +1,79 @@
+#pragma once
+
+#include "./numerics/interpolation1d_base.h"
+
+#include "./utils/vector.h"
+#include "./numerics/abs.h"
+
+namespace numerics{
+	
+	template <typename T>
+	struct interp_rational : Base_interp<T> {
+		using Base   = Base_interp<T>;
+
+	  	// bring base data members into scope (or use this->xx / this->yy below)
+	  	using Base::xx;
+	  	using Base::yy;
+	  	using Base::mm;
+
+	  	T dy;
+
+
+
+		interp_rational(const utils::Vector<T> &xv, const utils::Vector<T> &yv, uint64_t m)
+		: Base_interp<T>(xv, &yv[0], m), dy(T{0}){}
+
+		T rawinterp(int64_t jl, T x) override{
+
+			const T TINY = T{1.0e-99};
+			int64_t ns=0;
+			T y, w, t, hh, h, dd;
+			const T *xa = &xx[jl], *ya = &yy[jl];
+			utils::Vector<T> c(mm, T{0}), d(mm, T{0});
+
+			hh = numerics::abs(x - xa[0]);
+
+			for (int64_t i = 0; i < mm; ++i){
+				h = numerics::abs(x-xa[i]);
+				if (h == T{0}){
+					dy = T{0};
+					return ya[i];
+				}else if (h < hh){
+					ns = i;
+					hh = h;
+				}
+				c[i] = ya[i];
+				d[i] = ya[i] + TINY;					// The TINY part is needed to prevent a rare zero-over-zero condition.
+			}
+			y = ya[ns];
+			ns -= 1;
+			for (int64_t m = 1; m < mm; ++m){
+				for (int64_t i = 0; i < mm-m; ++i){
+					w = c[i+1] - d[i];
+					h = xa[i+m] - x;					// h will never be zero, since this was tested in the initializing loop.
+					t = (xa[i] - x)*d[i]/h;
+					dd = t - c[i+1];
+					if (dd == T{0}){					// This error condition indicates that the interpolating function has a pole at the requested value of x.
+						throw std::invalid_argument("Error in routine interp_rational"); // 
+					}
+					dd = w/dd;
+					d[i] = c[i+1]*dd;
+					c[i] = t*dd;
+				}
+				const bool take_left = (2 * (ns + 1) < (mm - m));
+
+				if (take_left) {
+					dy = c[ns + 1];
+				} else {
+					dy = d[ns];
+					ns -= 1;
+				}
+				y += dy;
+			}
+			return y;
+		}
+
+	};
+
+} // namespace numerics
+
@@ -0,0 +1,153 @@
+#pragma once
+
+#include "./numerics/min.h"
+#include "./numerics/max.h"
+#include "./numerics/abs.h"
+
+#include "./utils/vector.h"
+
+
+namespace numerics{
+	
+	template <typename T>
+	struct Base_interp{
+
+		int64_t n, mm;
+		int64_t jsav, dj;
+		bool cor;
+		const T *xx, *yy;
+
+
+		Base_interp(const utils::Vector<T>& x, const T *y, uint64_t m)
+			:n(x.size()), mm(m), jsav(0), cor(false), xx(&x[0]), yy(y){
+			//dj = numerics::min(static_cast<int64_t>(1), static_cast<int64_t>(std::pow(static_cast<T>(n), 0.25))); // from NR
+			dj = numerics::max(static_cast<int64_t>(1), static_cast<int64_t>(std::pow(static_cast<T>(n), 0.25))); // from chatbot
+
+			if (mm < 2 || n < mm) throw std::invalid_argument("Base_interp: invalid mm or n");
+		  	if (!xx || !yy)       throw std::invalid_argument("Base_interp: null data pointers");
+		  	if (n < 2)            throw std::invalid_argument("Base_interp: need at least 2 points");
+
+			bool asc = false;
+			if (xx[0] < xx[1]){
+				asc = true;
+			}
+			for (int64_t i = 1; i < n; ++i){
+				if (!(xx[i] > xx[i-1]) && asc) {
+					throw std::invalid_argument("x must be strictly increasing");
+				} else if (!(xx[i] < xx[i-1]) && !asc){
+        			throw std::invalid_argument("x must be strictly decreasing");
+				}
+			}
+		}
+
+		T interp(T x){
+			int64_t jlo;
+			if (cor){
+				std::cout << "Uses hunt()" << std::endl;
+				jlo = hunt(x);
+			}
+			else{
+				std::cout << "Uses locate()" << std::endl;
+				jlo = locate(x);
+			}
+			return rawinterp(jlo,x);
+		}
+
+		// Derived classes provide this as the actual interpolation method.
+		T virtual rawinterp(int64_t jlo, T x) = 0;
+
+
+		int64_t locate(const T x){
+			int64_t ju, jl;
+			int64_t jm;
+
+			if (n < 2 || mm < 2 || mm > n){
+				throw std::runtime_error("Interpolate: locate size error");	
+			}
+
+			bool ascnd = (xx[n-1] >= xx[0]);			// True if ascending order of table, false otherwise.
+			jl = 0;										// Initialize lower
+			ju = n-1;									// and upper limits.
+			while (ju - jl > 1) {						// If we are not yet done,
+				jm = (ju+jl) >> 1;						// compute a midpoint,
+				if ((x >= xx[jm]) == ascnd){
+					jl=jm;								// and replace either the lower limit
+				}else{
+					ju=jm;								// or the upper limit, as appropriate.
+				}
+			}											// Repeat until the test condition is satisﬁed.
+
+			if (std::abs(jl - jsav) > dj){				// Decide whether to use hunt or locate next time.
+				cor = false;
+			}else{
+				cor = true;
+			}
+			jsav = jl;
+			return numerics::max(static_cast<int64_t>(0), numerics::min(n-mm, jl-((mm-2)>>1)));
+		}
+
+		int64_t hunt(const T x){
+			int64_t jl=jsav, jm, ju, inc=1;
+
+			if (n < 2 || mm < 2 || mm > n){
+				throw std::runtime_error("Interpolate: hunt size error");	
+			}
+			bool ascnd=(xx[n-1] >= xx[0]);				// True if ascending order of table, false otherwise.
+			if (jl < 0 || jl > n-1) {					// Input guess not useful. Go immediately to bisection.
+				jl=0;
+				ju=n-1;
+			}else{
+				if ((x >= xx[jl]) == ascnd){			// Hunt up:
+					for (;;){
+						ju = jl + inc;
+						if (ju >= n-1){
+							ju = n-1;
+							break;						// Off end of table.
+						}else if((x < xx[ju]) == ascnd){
+							break;						// Found bracket.
+						}else{							// Not done, so double the increment and try again.
+							jl = ju;
+							inc += inc;
+						}
+					}	
+				}else{									// Hunt down:
+					ju = jl;
+					for (;;){
+						jl = jl - inc;
+						if (jl <= 0){					//Off end of table.
+							jl = 0;	
+							break;
+						}else if((x >= xx[jl]) == ascnd){
+							break;						// Found bracket.
+						}
+						else{							// Not done, so double the increment and try again.
+							ju = jl;
+							inc += inc;
+						}
+					}
+				}
+			}
+
+
+			while(ju-jl > 1){							// Hunt is done, so begin the ﬁnal bisection phase:
+				jm = (ju+jl) >> 1;
+				if ((x >= xx[jm]) == ascnd){
+					jl =jm;
+				}else{
+					ju=jm;
+				}
+			}
+			if (numerics::abs(jl-jsav) > dj){
+				cor = false;
+			}else{
+				cor = true;
+			}
+			jsav = jl;
+			return numerics::max(static_cast<int64_t>(0), numerics::min(n-mm, jl-((mm-2)>>1)));
+
+		}
+
+	};
+
+} // namespace numerics
+
@@ -0,0 +1,23 @@
+#pragma once
+
+
+#include "./utils/vector.h"
+#include "./utils/matrix.h"
+
+
+namespace numerics{
+	
+	template <typename T>
+	T max(const T a, const T b){
+
+		if(a < b){
+			return b;
+		}else{
+			return a;
+		}
+	}
+
+
+
+} // namespace numerics
+
@@ -0,0 +1,21 @@
+#pragma once
+
+
+#include "./utils/vector.h"
+#include "./utils/matrix.h"
+
+
+namespace numerics{
+	
+	template <typename T>
+	T min(const T a, const T b){
+
+		if(a < b){
+			return a;
+		}else{
+			return b;
+		}
+	}
+
+} // namespace numerics
+
@@ -5,3 +5,13 @@
 #include "./numerics/inverse.h"
 #include "./numerics/matmul.h"
 #include "./numerics/matvec.h"
+#include "./numerics/min.h"
+#include "./numerics/max.h"
+#include "./numerics/abs.h"
+#include "./numerics/interpolation1d_base.h"                       		// base
+#include "./numerics/interpolation1d/interpolation1d_linear.h"  		// derived
+#include "./numerics/interpolation1d/interpolation1d_polynomial.h"  	// derived
+#include "./numerics/interpolation1d/interpolation1d_cubic_spline.h"  	// derived
+#include "./numerics/interpolation1d/interpolation1d_rational.h"  		// derived 
+#include "./numerics/interpolation1d/interpolation1d_barycentric.h"  		// derived
+
@@ -3,6 +3,11 @@
 CC := g++
 CXXFLAGS := -std=c++14 -Wall -Iinclude -O3 -march=native -fopenmp
 LDFLAGS  := -fopenmp
+# Compiles .h files when there's been a change
+DEPFLAGS := -MMD -MP
+
+CXX ?= g++
+


 # Directories
@@ -12,7 +17,6 @@ OBJ_DIR := obj
 BIN_DIR := bin
 TEST_BIN := $(BIN_DIR)/tests

-# All test sources



@@ -31,6 +35,11 @@ OBJS := $(patsubst $(SRC_DIR)/%.cpp,$(OBJ_DIR)/%.o,$(SRCS))
 # === Test sources ===
 TEST_SRCS := $(shell find test -name 'test_*.cpp')
 TEST_OBJS := $(patsubst test/%.cpp, $(OBJ_DIR)/test/%.o, $(TEST_SRCS))
+
+# Compiles .h files when there's been a change
+DEPS := $(OBJS:.o=.d) $(TEST_OBJS:.o=.d) $(TEST_MAIN:.o=.d)
+-include $(DEPS)
+
 # The single file that defines TEST_MAIN / main()
 TEST_MAIN := $(OBJ_DIR)/test/test_all.o

@@ -53,6 +62,12 @@ export OMP_DYNAMIC
 export OMP_SCHEDULE
 export OMP_DISPLAY_ENV

+# What belongs to "run" vs "test"
+RUN_DEPS := $(OBJS:.o=.d)
+TEST_DEPS := $(TEST_OBJS:.o=.d) $(TEST_MAIN:.o=.d)
+
+RUN_ARTIFACTS  := $(TARGET) $(OBJS) $(RUN_DEPS)
+TEST_ARTIFACTS := $(TEST_BIN) $(TEST_OBJS) $(TEST_MAIN) $(TEST_DEPS)


 # === Default Target ===
@@ -66,11 +81,11 @@ $(TARGET): $(OBJS)
 # === Compiling source files to object files ===
 $(OBJ_DIR)/%.o: $(SRC_DIR)/%.cpp
 	@mkdir -p $(dir $@)
-	$(CXX) $(CXXFLAGS) -c $< -o $@
+	$(CXX) $(CXXFLAGS) $(DEPFLAGS) -c $< -o $@

 # === Run with OpenMP env set only for the run ===
 .PHONY: run
-run: $(TARGET)
+run: clean-test $(TARGET)
 	@echo ">>> OMP_PROC_BIND=$(OMP_PROC_BIND)"
 	@echo ">>> OMP_PLACES=$(OMP_PLACES)"
 	@echo ">>> OMP_MAX_ACTIVE_LEVELS=$(OMP_MAX_LEVELS)"
@@ -96,11 +111,6 @@ run-single: ## Single-level parallel (good default)
 run-nested: ## Two-level nested (outer,inner), adjust to your cores
 	$(MAKE) run OMP_MAX_LEVELS=2 OMP_THREADS="4,8" OMP_PROC_BIND=close OMP_PLACES=cores

-# Clean up
-.PHONY: clean
-clean:
-	rm -rf $(OBJ_DIR) $(BIN_DIR)
-
 # Optional: print debug info
 .PHONY: info
 info:
@@ -111,7 +121,7 @@ info:


 .PHONY: test
-test: $(TEST_BIN)
+test: clean-run $(TEST_BIN)
 	@echo ">>> OMP_PROC_BIND=$(OMP_PROC_BIND)"
 	@echo ">>> OMP_PLACES=$(OMP_PLACES)"
 	@echo ">>> OMP_MAX_ACTIVE_LEVELS=$(OMP_MAX_LEVELS)"
@@ -134,4 +144,18 @@ $(TEST_BIN): $(TEST_OBJS) $(TEST_MAIN)

 $(OBJ_DIR)/test/%.o: test/%.cpp
 	@mkdir -p $(dir $@)
-	$(CXX) $(CXXFLAGS) -c $< -o $@
+	$(CXX) $(CXXFLAGS) $(DEPFLAGS) -c $< -o $@
+
+# Clean up
+.PHONY: clean
+clean:
+	rm -rf $(OBJ_DIR) $(BIN_DIR)
+
+.PHONY: clean-run clean-test
+clean-run:
+	@echo "Cleaning run artifacts..."
+	@rm -f $(RUN_ARTIFACTS)
+
+clean-test:
+	@echo "Cleaning test artifacts..."
+	@rm -f $(TEST_ARTIFACTS)
@@ -0,0 +1,41 @@
+obj/main.o: src/main.cpp include/./utils/utils.h include/./utils/vector.h \
+ include/./utils/matrix.h include/./numerics/numerics.h \
+ include/./numerics/transpose.h include/./numerics/inverse.h \
+ include/./numerics/inverse/inverse_gauss_jordan.h \
+ include/./numerics/inverse/inverse_lu.h include/./decomp/lu.h \
+ include/./numerics/matmul.h include/./core/omp_config.h \
+ include/./numerics/matvec.h include/./numerics/min.h \
+ include/./numerics/max.h include/./numerics/abs.h \
+ include/./numerics/interpolation1d_base.h \
+ include/./numerics/interpolation1d/interpolation1d_linear.h \
+ include/./numerics/interpolation1d/interpolation1d_polynomial.h \
+ include/./numerics/interpolation1d/interpolation1d_cubic_spline.h \
+ include/./numerics/interpolation1d/interpolation1d_rational.h \
+ include/./numerics/interpolation1d/interpolation1d_barycentric.h \
+ include/./decomp/decomp.h include/./modules/grid1d.h \
+ include/utils/vector.h include/./modules/finitedifference1d.h
+include/./utils/utils.h:
+include/./utils/vector.h:
+include/./utils/matrix.h:
+include/./numerics/numerics.h:
+include/./numerics/transpose.h:
+include/./numerics/inverse.h:
+include/./numerics/inverse/inverse_gauss_jordan.h:
+include/./numerics/inverse/inverse_lu.h:
+include/./decomp/lu.h:
+include/./numerics/matmul.h:
+include/./core/omp_config.h:
+include/./numerics/matvec.h:
+include/./numerics/min.h:
+include/./numerics/max.h:
+include/./numerics/abs.h:
+include/./numerics/interpolation1d_base.h:
+include/./numerics/interpolation1d/interpolation1d_linear.h:
+include/./numerics/interpolation1d/interpolation1d_polynomial.h:
+include/./numerics/interpolation1d/interpolation1d_cubic_spline.h:
+include/./numerics/interpolation1d/interpolation1d_rational.h:
+include/./numerics/interpolation1d/interpolation1d_barycentric.h:
+include/./decomp/decomp.h:
+include/./modules/grid1d.h:
+include/utils/vector.h:
+include/./modules/finitedifference1d.h:
@@ -3,9 +3,14 @@
 #include "./decomp/decomp.h"
 #include "./core/omp_config.h"

+#include "./modules/grid1d.h"
+#include "./modules/finitedifference1d.h"
+
+
 #include <iostream>
 #include <stdexcept>

+//#include "./numerics/interpolation/interpolation_linear.h"



@@ -13,11 +18,65 @@
 int main(int argc, char const *argv[])
 {    
    
-    utils::Md A;
-    decomp::LUdcmpd lu(A);

-    // Single-level, 16 threads, runtime may adjust
-    //omp_configure(/*max_levels=*/1, /*dynamic=*/true, /*threads_per_level=*/{16});
+    utils::Vector<double> x(100, 0), y(100,0);
+    for (uint64_t i = 0; i < 100; ++i){
+        x[i] = i+1;
+        y[i] = 1 + i*0.1;
+    }
+    //y[9] = 1.5;
+
+    x.print();
+    y.print();
+
+    double p = 5.5;
+
+
+
+    numerics::interp_linear<double> linear(x,y);
+    numerics::interp_polynomial<double> polynomial(x,y, 3);
+    numerics::interp_cubic_spline<double> cubic_spline(x,y);
+    numerics::interp_rational<double> rational(x,y,2);
+    numerics::interp_barycentric<double> barycentric(x,y, 2);
+
+    std::cout << "=== interpolate: " << p << " ===" << std::endl;
+
+    std::cout << linear.interp(p) << std::endl;
+    std::cout << linear.interp(p) << std::endl;
+    std::cout << polynomial.interp(p) << std::endl; // error = polynomial.dy
+    std::cout << cubic_spline.interp(p) << std::endl;
+    std::cout << rational.interp(p) << std::endl;
+    std::cout << barycentric.interp(p) << std::endl;
+
+    p += 0.01;
+
+    std::cout << "=== interpolate: " << p << " ===" << std::endl;
+
+    std::cout << linear.interp(p) << std::endl;
+    std::cout << polynomial.interp(p) << std::endl;
+    std::cout << cubic_spline.interp(p) << std::endl;
+    std::cout << rational.interp(p) << std::endl;
+    std::cout << barycentric.interp(p) << std::endl;
+
+    p += 0.01;
+
+    std::cout << "=== interpolate: " << p << " ===" << std::endl;
+
+    std::cout << linear.interp(p) << std::endl;
+    std::cout << polynomial.interp(p) << std::endl;
+    std::cout << cubic_spline.interp(p) << std::endl;
+    std::cout << rational.interp(p) << std::endl;
+    std::cout << barycentric.interp(p) << std::endl;
+
+    p += 50.01;
+
+    std::cout << "=== interpolate: " << p << " ===" << std::endl;
+
+    std::cout << linear.interp(p) << std::endl;
+    std::cout << polynomial.interp(p) << std::endl;
+    std::cout << cubic_spline.interp(p) << std::endl;
+    std::cout << rational.interp(p) << std::endl;
+    std::cout << barycentric.interp(p) << std::endl;


    return 0;
@@ -0,0 +1,115 @@
+#include "test_common.h"
+#include "./utils/utils.h"
+#include "./numerics/inverse.h"
+#include "./numerics/matmul.h"
+
+
+TEST_CASE(Inverse_GJ_Basic_3x3) {
+    using T = double;
+    utils::Matrix<T> A(3,3, T{0});
+    // Well-conditioned 3x3
+    A(0,0)=3;   A(0,1)=0.2; A(0,2)=-0.1;
+    A(1,0)=0.1; A(1,1)=2.5; A(1,2)=0.3;
+    A(2,0)=-0.2;A(2,1)=0.4; A(2,2)=2.0;
+
+    auto Ainv = numerics::inverse<T>(A, "Gauss-Jordan");
+    utils::Matrix<T> I;
+    I.eye(3);
+    auto prod = numerics::matmul<T>(A, Ainv);
+
+    CHECK(prod.nearly_equal(I, (T)1e-10), "inverse(GJ): A*A^{-1} ≈ I");
+}
+
+TEST_CASE(Inverse_LU_Basic_3x3) {
+    using T = double;
+    utils::Matrix<T> A(3,3, T{0});
+    A(0,0)=3;   A(0,1)=0.2; A(0,2)=-0.1;
+    A(1,0)=0.1; A(1,1)=2.5; A(1,2)=0.3;
+    A(2,0)=-0.2;A(2,1)=0.4; A(2,2)=2.0;
+
+    auto Ainv = numerics::inverse<T>(A, "LU");
+    utils::Matrix<T> I;
+    I.eye(3);
+    auto prod = numerics::matmul<T>(A, Ainv);
+
+    CHECK(prod.nearly_equal(I, (T)1e-10), "inverse(LU): A*A^{-1} ≈ I");
+}
+
+TEST_CASE(Inverse_GJ_vs_LU_Consistency) {
+    using T = double;
+    utils::Matrix<T> A(3,3, T{0});
+    A(0,0)=4; A(0,1)=1;  A(0,2)=2;
+    A(1,0)=0; A(1,1)=3;  A(1,2)=-1;
+    A(2,0)=0; A(2,1)=0;  A(2,2)=2;
+
+    auto GJ = numerics::inverse<T>(A, "Gauss-Jordan");
+    auto LU = numerics::inverse<T>(A, "LU");
+
+    CHECK(GJ.nearly_equal(LU, (T)1e-12), "inverse: GJ and LU produce nearly the same result");
+}
+
+
+TEST_CASE(Inplace_Inverse_LU) {
+    using T = double;
+    utils::Matrix<T> A(3,3, T{0});
+    A(0,0)=3;   A(0,1)=0.2; A(0,2)=-0.1;
+    A(1,0)=0.1; A(1,1)=2.5; A(1,2)=0.3;
+    A(2,0)=-0.2;A(2,1)=0.4; A(2,2)=2.0;
+
+    auto Ainv_ref = numerics::inverse<T>(A, "LU"); // out-of-place
+    auto A_copy   = A;
+    numerics::inplace_inverse<T>(A_copy, "LU");    // in-place
+
+    CHECK(A_copy.nearly_equal(Ainv_ref, (T)1e-12), "inplace_inverse(LU) equals out-of-place");
+}
+
+TEST_CASE(Inplace_Inverse_GJ) {
+    using T = double;
+    utils::Matrix<T> A(2,2, T{0});
+    A(0,0)=2; A(0,1)=1;
+    A(1,0)=1; A(1,1)=3;
+
+    auto Ainv_ref = numerics::inverse<T>(A, "Gauss-Jordan");
+    auto A_copy   = A;
+    numerics::inplace_inverse<T>(A_copy, "Gauss-Jordan");
+
+    CHECK(A_copy.nearly_equal(Ainv_ref, (T)1e-12), "inplace_inverse(GJ) equals out-of-place");
+}
+
+TEST_CASE(Inverse_Identity) {
+    using T = double;
+    utils::Matrix<T> I;
+    I.eye(3);
+    auto invI = numerics::inverse<T>(I, "LU");
+    CHECK(invI.nearly_equal(I, (T)0), "inverse(I) == I");
+}
+TEST_CASE(Inverse_NonSquare_Throws) {
+    using T = double;
+    utils::Matrix<T> A(2,3, T{0}); // non-square
+    bool threw1=false, threw2=false;
+    try { auto X = numerics::inverse<T>(A, "LU"); (void)X; } catch(...) { threw1=true; }
+    try { numerics::inplace_inverse<T>(A, "Gauss-Jordan"); } catch(...) { threw2=true; }
+    CHECK(threw1 && threw2, "inverse throws on non-square for both methods");
+}
+
+TEST_CASE(Inverse_Singular_Throws) {
+    using T = double;
+    utils::Matrix<T> S(3,3, T{0});
+    S(0,0)=1; S(0,1)=2; S(0,2)=3;
+    S(1,0)=1; S(1,1)=2; S(1,2)=3; // duplicate row -> singular
+    S(2,0)=0; S(2,1)=1; S(2,2)=0;
+
+    bool threw_gj=false, threw_lu=false;
+    try { auto X = numerics::inverse<T>(S, "Gauss-Jordan"); (void)X; } catch(...) { threw_gj=true; }
+    try { auto X = numerics::inverse<T>(S, "LU");           (void)X; } catch(...) { threw_lu=true;  }
+    CHECK(threw_gj && threw_lu, "inverse throws on singular for both methods");
+}
+
+TEST_CASE(Inverse_Unknown_Method_Throws) {
+    using T = double;
+    utils::Matrix<T> A(2,2, T{0});
+    A(0,0)=1; A(1,1)=1;
+    bool threw=false;
+    try { auto X = numerics::inverse<T>(A, "Foobar"); (void)X; } catch(...) { threw=true; }
+    CHECK(threw, "inverse unknown method throws");
+}