#include "optimization/AdditionalIceUtils.h"
//#include "vectort.h" // ICE
// #include "ludecomp.h" //ICE
// #include "mateigen.h" //ICE
#include <cassert>

#include <core/algebra/LUDecomposition.h>
#include <core/algebra/EigValues.h>
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 OPTIMIZATION::matrix_type &A, OPTIMIZATION::matrix_type &x, const OPTIMIZATION::matrix_type b)
{
  int n=A.cols();

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

  matrix_type LU;
  NICE::VectorT<int> indx;
  NICE::Vector xv(n);
  NICE::Vector iv(n);

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

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

void getEig(const OPTIMIZATION::matrix_type &A, OPTIMIZATION::matrix_type &L, OPTIMIZATION::matrix_type &V)
{
  std::cerr << "This feature is not supported currently" << std::endl;
//   int n= A.cols();
//   L = matrix_type(n,1);
//   V = matrix_type(n,n);
// 
//   NICE::Vector l(n);
// 
//   NICE::EVArnoldi eig;
//   eig.getEigenvalues (A, l, V, n);
//   for(int i=0; i < n;++i)
//   {
//     L(i,0) = l[i];
//   }

}