// Evgenii B. Rudnyi, http://Evgenii.Rudnyi.Ru

#ifdef USECLOCK
#include <time.h>
#else
#include <sys/time.h>
#endif
#include <vector>

class Matrix : public std::vector<double>
{
  size_t m;
  size_t n;

public:
  Matrix() : m(0), n(0) {}
  Matrix(size_t m_, size_t n_) : std::vector<double>(m_*n_), m(m_), n(n_) {}
  Matrix(size_t m_, size_t n_, double val) : std::vector<double>(m_*n_, val), m(m_), n(n_) {}

  void resize(size_t m_, size_t n_)
    {m = m_; n = n_; std::vector<double>::resize(m_*n_);}
  void reserve(size_t m_, size_t n_)
    {std::vector<double>::reserve(m_*n_);}
  void clear()
    {m = n = 0; std::vector<double>::clear();}

  size_t NRows() const {return m;}
  size_t NCols() const {return n;}

  double& operator()(size_t i, size_t j)
  {
    return operator[](i + j*m);
  }
  const double& operator()(size_t i, size_t j) const
  {
    return operator[](i + j*m);
  }
	void swap(Matrix &y)
	{
		std::vector<double>::swap(y);
		std::swap(n, y.n);
		std::swap(m, y.m);
	}
	void clearMemory()
	{
		Matrix empty;
	  swap(empty);
	}
};

/// Simple timing
#ifdef USECLOCK
class Timing
{
	double st;
	double gettime() {return clock();}
public:
	Timing() {st = gettime();}
	void reset() {st = gettime();}
	double time() {return (gettime() - st)/CLK_TCK;}
};
#else
class Timing
{
	timeval tim;
//Sun is different with gettimeofday from others
#ifndef __sun__
	struct timezone zone;
#endif	
	double st;
	double gettime()
	{
#ifdef __sun__
		gettimeofday(&tim, (void*)0);
#else
		gettimeofday(&tim, &zone);
#endif	
		return double(tim.tv_sec) + tim.tv_usec*1e-6;
	}
public:
	Timing() {st = gettime();}
	void reset() {st = gettime();}
	double time() {return (gettime() - st);}
};
#endif

#ifdef MKL
#define dgetf2_ DGETF2
#define dgetfr_ DGETFR
#define dgetrs_ DGETRS
#endif

extern "C" {
  int dgetf2_(int *m, int *n, double *a, int *lda, int *ipiv, int *info);
  int dgetrf_(int *m, int *n, double *a, int *lda, int *ipiv, int *info);
// you may need add int at the end to specify the lenght of trans
	int dgetrs_(const char *trans, int *n, int *nrhs, double *a, int *lda, int *ipiv, double *b, int *ldb, int *info);
}

int dgetf2(Matrix &A, std::vector<int> &ipiv)
{
	int m = A.NRows();
	int n = A.NCols();
	int info;
	dgetf2_(&m, &n, &*A.begin(), &m, &*ipiv.begin(), &info);
	return info;
}

int dgetrf(Matrix &A, std::vector<int> &ipiv)
{
	int m = A.NRows();
	int n = A.NCols();
	int info;
	dgetrf_(&m, &n, &*A.begin(), &m, &*ipiv.begin(), &info);
	return info;
}

int dgetrs(Matrix &A, Matrix &B, std::vector<int> &ipiv)
{
	int m = A.NRows();
	int n = A.NCols();
	int nrhs = B.NCols();
	int info;
	dgetrs_("N", &m, &nrhs, &*A.begin(), &m, &*ipiv.begin(), &*B.begin(), &m, &info);
	return info;
}