#define EIGEN_NO_DEBUG_SMALL_PRODUCT_BLOCKS
#include "sparse_solver.h"

#if defined(DEBUG)
#undef DEBUG
#endif

#include <Eigen/AccelerateSupport>

template<typename MatrixType,typename DenseMat>
int generate_sparse_rectangular_problem(MatrixType& A, DenseMat& dA, int maxRows = 300, int maxCols = 300)
{
  typedef typename MatrixType::Scalar Scalar;
  int rows = internal::random<int>(1, maxRows);
  int cols = internal::random<int>(1, maxCols);
  double density = (std::max)(8.0 / (rows * cols), 0.01);
  
  A.resize(rows,cols);
  dA.resize(rows,cols);
  initSparse<Scalar>(density, dA, A, ForceNonZeroDiag);
  A.makeCompressed();
  return rows;
}

template<typename MatrixType,typename DenseMat>
int generate_sparse_square_symmetric_problem(MatrixType& A, DenseMat& dA, int maxSize = 300)
{
  typedef typename MatrixType::Scalar Scalar;
  int rows = internal::random<int>(1, maxSize);
  int cols = rows;
  double density = (std::max)(8.0 / (rows * cols), 0.01);
  
  A.resize(rows,cols);
  dA.resize(rows,cols);
  initSparse<Scalar>(density, dA, A, ForceNonZeroDiag);
  dA = dA * dA.transpose();
  A  = A * A.transpose();
  A.makeCompressed();
  return rows;
}

template<typename Scalar, typename Solver> void test_accelerate_ldlt()
{
  typedef SparseMatrix<Scalar> MatrixType; 
  typedef Matrix<Scalar,Dynamic,1> DenseVector;
 
  MatrixType A;
  Matrix<Scalar,Dynamic,Dynamic> dA;

  generate_sparse_square_symmetric_problem(A, dA);

  DenseVector b = DenseVector::Random(A.rows());

  Solver solver;
  solver.compute(A);

  if (solver.info() != Success)
  {
    std::cerr << "sparse LDLT factorization failed\n";
    exit(0);
    return;
  }

  DenseVector x = solver.solve(b);

  if (solver.info() != Success)
  {
    std::cerr << "sparse LDLT factorization failed\n";
    exit(0);
    return;
  }

  //Compare with a dense solver
  DenseVector refX = dA.ldlt().solve(b);
  VERIFY((A * x).isApprox(A * refX, test_precision<Scalar>()));
}

template<typename Scalar, typename Solver> void test_accelerate_llt()
{
  typedef SparseMatrix<Scalar> MatrixType; 
  typedef Matrix<Scalar,Dynamic,1> DenseVector;
 
  MatrixType A;
  Matrix<Scalar,Dynamic,Dynamic> dA;

  generate_sparse_square_symmetric_problem(A, dA);

  DenseVector b = DenseVector::Random(A.rows());

  Solver solver;
  solver.compute(A);

  if (solver.info() != Success)
  {
    std::cerr << "sparse LLT factorization failed\n";
    exit(0);
    return;
  }

  DenseVector x = solver.solve(b);

  if (solver.info() != Success)
  {
    std::cerr << "sparse LLT factorization failed\n";
    exit(0);
    return;
  }

  //Compare with a dense solver
  DenseVector refX = dA.llt().solve(b);
  VERIFY((A * x).isApprox(A * refX, test_precision<Scalar>()));
}

template<typename Scalar, typename Solver> void test_accelerate_qr()
{
  typedef SparseMatrix<Scalar> MatrixType; 
  typedef Matrix<Scalar,Dynamic,1> DenseVector;
 
  MatrixType A;
  Matrix<Scalar,Dynamic,Dynamic> dA;

  generate_sparse_rectangular_problem(A, dA);

  DenseVector b = DenseVector::Random(A.rows());

  Solver solver;
  solver.compute(A);

  if (solver.info() != Success)
  {
    std::cerr << "sparse QR factorization failed\n";
    exit(0);
    return;
  }

  DenseVector x = solver.solve(b);

  if (solver.info() != Success)
  {
    std::cerr << "sparse QR factorization failed\n";
    exit(0);
    return;
  }

  //Compare with a dense solver
  DenseVector refX = dA.colPivHouseholderQr().solve(b);
  VERIFY((A * x).isApprox(A * refX, test_precision<Scalar>()));
}

template<typename Scalar> 
void run_tests()
{
  typedef SparseMatrix<Scalar> MatrixType; 

  test_accelerate_ldlt<Scalar, AccelerateLDLT<MatrixType, Lower | Symmetric> >();
  test_accelerate_ldlt<Scalar, AccelerateLDLTUnpivoted<MatrixType, Lower | Symmetric> >();
  test_accelerate_ldlt<Scalar, AccelerateLDLTSBK<MatrixType, Lower | Symmetric> >();
  test_accelerate_ldlt<Scalar, AccelerateLDLTTPP<MatrixType, Lower | Symmetric> >();

  test_accelerate_ldlt<Scalar, AccelerateLDLT<MatrixType, Upper | Symmetric> >();
  test_accelerate_ldlt<Scalar, AccelerateLDLTUnpivoted<MatrixType, Upper | Symmetric> >();
  test_accelerate_ldlt<Scalar, AccelerateLDLTSBK<MatrixType, Upper | Symmetric> >();
  test_accelerate_ldlt<Scalar, AccelerateLDLTTPP<MatrixType, Upper | Symmetric> >();

  test_accelerate_llt<Scalar, AccelerateLLT<MatrixType, Lower | Symmetric> >();

  test_accelerate_llt<Scalar, AccelerateLLT<MatrixType, Upper | Symmetric> >();

  test_accelerate_qr<Scalar, AccelerateQR<MatrixType, 0> >();
}

EIGEN_DECLARE_TEST(accelerate_support)
{
  CALL_SUBTEST_1(run_tests<float>());
  CALL_SUBTEST_2(run_tests<double>());
}