#include <string>
#include <iostream>
#include <complex>

typedef std::complex<double> dcomplex;

extern "C"{
  void dgemm_(const char* transa, const char* transb, const int* m, const int* n, const int* k, const double* alpha, const double* A, const int* lda, const double* B, const int* ldb, const double* beta, double* C, const int* ldc);
  void zgemm_(const char* transa, const char* transb, const int* m, const int* n, const int* k, const dcomplex* alpha, const dcomplex* A, const int* lda, const dcomplex* B, const int* ldb, const dcomplex* beta, dcomplex* C, const int* ldc);
  void dgeev_(const char* jobvl, const char* jobvr,  const int* n,  double* A, const int* lda, double* wr, double* wi, double* vl, const int* ldvl, double* vr, const int* ldvr, double* work, const int* lwork, int* info);
  void zgeev_(const char* jobvl, const char* jobvr,  const int* n,  dcomplex* A,const int* lda, dcomplex* w, dcomplex* vl, const int* ldvl, dcomplex* vr, const int* ldvr, dcomplex* work, const int* lwork, double* rwork, int* info);
}


// Implementation for multiplication of real matrices
inline void xgemm(const std::string& transa, const std::string& transb, const int m, const int n,
                  const int k, const double alpha, const double* A,
                  const int lda, const double* B, const int ldb, const double beta,
                  double* C, const int ldc)
{
  dgemm_(transa.c_str(), transb.c_str(), &m, &n, &k, &alpha, A, &lda, B, &ldb, &beta, C, &ldc);
}

// Implementation for multiplication of complex matrices
inline void xgemm(const std::string& transa, const std::string& transb, const int m, const int n,
                  const int k, const dcomplex& alpha, const dcomplex* A,
                  const int lda, const dcomplex* B, const int ldb, const dcomplex& beta,
                  dcomplex* C, const int ldc)
{
  zgemm_(transa.c_str(), transb.c_str(), &m, &n, &k, &alpha, A, &lda, B, &ldb, &beta, C, &ldc);
}

// Eigenvalues of real matrix
inline int xgeev_(int N, double* A, int lda, dcomplex* w, double* work, int lwork)
{
  int one = 1, info = 0;
  double* wr = new double[N];
  double* wi = new double[N];  
  dgeev_("N", "N",  &N,  A, &lda, wr, wi, NULL, &one, NULL, &one, work, &lwork, &info);
  for (int i=0; i<N; i++) w[i] = dcomplex(wr[i], wi[i]);
  delete[] wr;
  delete[] wi;
  if (info)
    std::cerr << "Can't compute eigenvalues! " << info << std::endl;
  return info;
}
// Eigenvalues of a complex matrix
inline int xgeev_(int N, dcomplex* A, int lda, dcomplex* w, dcomplex* work, int lwork)
{
  int info = 0, one = 1;
  double* rwork = new double[2*N];
  zgeev_("N", "N",  &N,  A, &lda, w, NULL, &one, NULL, &one, work, &lwork, rwork, &info);
  if (info) std::cerr << "Can't compute eigenvalues and eigenvectors! " << info << std::endl;
  delete[] rwork;
  return info;
}