#include "optimization/AdditionalIceUtils.h"
#include "vectort.h" // ICE
#include "ludecomp.h" //ICE
#include "mateigen.h" //ICE
#include <cassert>
using namespace opt;
using namespace ice;

matrix_type & MatrixAddVal ( matrix_type & mat, double val )
{
	for(int i=0; i< mat.rows(); ++i)
	{
		for(int j =0; j<mat.cols(); ++j)
		{
			mat[i][j] += val;
		}
	}

	return mat;
}


double MatrixSum(const matrix_type &mat)
{
	double sum=0.0;
	for(int i=0; i< mat.rows(); ++i)
	{
		for(int j =0; j<mat.cols(); ++j)
		{
			sum += mat[i][j];
		}
	}
	return sum;
}



void linSolve(const opt::matrix_type &A, opt::matrix_type &x, const opt::matrix_type b)
{
	int n=A.cols();

	assert(x.rows() == n && b.rows() == n );

	matrix_type LU;
	IVector indx;
	Vector xv(n);
	Vector iv(n);

	for(int i=0; i < n;++i)
	{
		iv[i] = b[i][0];
	}

	// LU-Zerlegung mit Pivotisierung
	LUDecompositionPacked(A,LU,indx,true);
	
	// Lösung M*x1=i1
	xv = LUSolve(LU,indx,iv);
	
	x=matrix_type(n,1);
	for(int i=0; i < n;++i)
	{
		x[i][0] = xv[i];
	}
	
}

void getEig(const opt::matrix_type &A, opt::matrix_type &L, opt::matrix_type &V)
{
	int n= A.cols();
	L = matrix_type(n,1);
	V = matrix_type(n,n);

	Vector l(n);

	Eigenvalue(A,l,V);
	for(int i=0; i < n;++i)
	{
		L[i][0] = l[i];
	}

}