diff --git a/lapack/cholesky.cpp b/lapack/cholesky.cpp
deleted file mode 100644
index ea3bc12..0000000
--- a/lapack/cholesky.cpp
+++ /dev/null
@@ -1,72 +0,0 @@
-// This file is part of Eigen, a lightweight C++ template library
-// for linear algebra.
-//
-// Copyright (C) 2010-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
-//
-// This Source Code Form is subject to the terms of the Mozilla
-// Public License v. 2.0. If a copy of the MPL was not distributed
-// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
-
-#include "lapack_common.h"
-#include <Eigen/Cholesky>
-
-// POTRF computes the Cholesky factorization of a real symmetric positive definite matrix A.
-EIGEN_LAPACK_FUNC(potrf,(char* uplo, int *n, RealScalar *pa, int *lda, int *info))
-{
-  *info = 0;
-        if(UPLO(*uplo)==INVALID) *info = -1;
-  else  if(*n<0)                 *info = -2;
-  else  if(*lda<std::max(1,*n))  *info = -4;
-  if(*info!=0)
-  {
-    int e = -*info;
-    return xerbla_(SCALAR_SUFFIX_UP"POTRF", &e, 6);
-  }
-
-  Scalar* a = reinterpret_cast<Scalar*>(pa);
-  MatrixType A(a,*n,*n,*lda);
-  int ret;
-  if(UPLO(*uplo)==UP) ret = int(internal::llt_inplace<Scalar, Upper>::blocked(A));
-  else                ret = int(internal::llt_inplace<Scalar, Lower>::blocked(A));
-
-  if(ret>=0)
-    *info = ret+1;
-  
-  return 0;
-}
-
-// POTRS solves a system of linear equations A*X = B with a symmetric
-// positive definite matrix A using the Cholesky factorization
-// A = U**T*U or A = L*L**T computed by DPOTRF.
-EIGEN_LAPACK_FUNC(potrs,(char* uplo, int *n, int *nrhs, RealScalar *pa, int *lda, RealScalar *pb, int *ldb, int *info))
-{
-  *info = 0;
-        if(UPLO(*uplo)==INVALID) *info = -1;
-  else  if(*n<0)                 *info = -2;
-  else  if(*nrhs<0)              *info = -3;
-  else  if(*lda<std::max(1,*n))  *info = -5;
-  else  if(*ldb<std::max(1,*n))  *info = -7;
-  if(*info!=0)
-  {
-    int e = -*info;
-    return xerbla_(SCALAR_SUFFIX_UP"POTRS", &e, 6);
-  }
-
-  Scalar* a = reinterpret_cast<Scalar*>(pa);
-  Scalar* b = reinterpret_cast<Scalar*>(pb);
-  MatrixType A(a,*n,*n,*lda);
-  MatrixType B(b,*n,*nrhs,*ldb);
-
-  if(UPLO(*uplo)==UP)
-  {
-    A.triangularView<Upper>().adjoint().solveInPlace(B);
-    A.triangularView<Upper>().solveInPlace(B);
-  }
-  else
-  {
-    A.triangularView<Lower>().solveInPlace(B);
-    A.triangularView<Lower>().adjoint().solveInPlace(B);
-  }
-
-  return 0;
-}
diff --git a/lapack/eigenvalues.cpp b/lapack/eigenvalues.cpp
deleted file mode 100644
index 921c515..0000000
--- a/lapack/eigenvalues.cpp
+++ /dev/null
@@ -1,62 +0,0 @@
-// This file is part of Eigen, a lightweight C++ template library
-// for linear algebra.
-//
-// Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr>
-//
-// This Source Code Form is subject to the terms of the Mozilla
-// Public License v. 2.0. If a copy of the MPL was not distributed
-// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
-
-#include "lapack_common.h"
-#include <Eigen/Eigenvalues>
-
-// computes eigen values and vectors of a general N-by-N matrix A
-EIGEN_LAPACK_FUNC(syev,(char *jobz, char *uplo, int* n, Scalar* a, int *lda, Scalar* w, Scalar* /*work*/, int* lwork, int *info))
-{
-  // TODO exploit the work buffer
-  bool query_size = *lwork==-1;
-  
-  *info = 0;
-        if(*jobz!='N' && *jobz!='V')                    *info = -1;
-  else  if(UPLO(*uplo)==INVALID)                        *info = -2;
-  else  if(*n<0)                                        *info = -3;
-  else  if(*lda<std::max(1,*n))                         *info = -5;
-  else  if((!query_size) && *lwork<std::max(1,3**n-1))  *info = -8;
-    
-  if(*info!=0)
-  {
-    int e = -*info;
-    return xerbla_(SCALAR_SUFFIX_UP"SYEV ", &e, 6);
-  }
-  
-  if(query_size)
-  {
-    *lwork = 0;
-    return 0;
-  }
-  
-  if(*n==0)
-    return 0;
-  
-  PlainMatrixType mat(*n,*n);
-  if(UPLO(*uplo)==UP) mat = matrix(a,*n,*n,*lda).adjoint();
-  else                mat = matrix(a,*n,*n,*lda);
-  
-  bool computeVectors = *jobz=='V' || *jobz=='v';
-  SelfAdjointEigenSolver<PlainMatrixType> eig(mat,computeVectors?ComputeEigenvectors:EigenvaluesOnly);
-  
-  if(eig.info()==NoConvergence)
-  {
-    make_vector(w,*n).setZero();
-    if(computeVectors)
-      matrix(a,*n,*n,*lda).setIdentity();
-    //*info = 1;
-    return 0;
-  }
-  
-  make_vector(w,*n) = eig.eigenvalues();
-  if(computeVectors)
-    matrix(a,*n,*n,*lda) = eig.eigenvectors();
-  
-  return 0;
-}
diff --git a/lapack/lu.cpp b/lapack/lu.cpp
deleted file mode 100644
index 90cebe0..0000000
--- a/lapack/lu.cpp
+++ /dev/null
@@ -1,89 +0,0 @@
-// This file is part of Eigen, a lightweight C++ template library
-// for linear algebra.
-//
-// Copyright (C) 2010-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
-//
-// This Source Code Form is subject to the terms of the Mozilla
-// Public License v. 2.0. If a copy of the MPL was not distributed
-// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
-
-#include "common.h"
-#include <Eigen/LU>
-
-// computes an LU factorization of a general M-by-N matrix A using partial pivoting with row interchanges
-EIGEN_LAPACK_FUNC(getrf,(int *m, int *n, RealScalar *pa, int *lda, int *ipiv, int *info))
-{
-  *info = 0;
-        if(*m<0)                  *info = -1;
-  else  if(*n<0)                  *info = -2;
-  else  if(*lda<std::max(1,*m))   *info = -4;
-  if(*info!=0)
-  {
-    int e = -*info;
-    return xerbla_(SCALAR_SUFFIX_UP"GETRF", &e, 6);
-  }
-
-  if(*m==0 || *n==0)
-    return 0;
-
-  Scalar* a = reinterpret_cast<Scalar*>(pa);
-  int nb_transpositions;
-  int ret = int(Eigen::internal::partial_lu_impl<Scalar,ColMajor,int>
-                     ::blocked_lu(*m, *n, a, *lda, ipiv, nb_transpositions));
-
-  for(int i=0; i<std::min(*m,*n); ++i)
-    ipiv[i]++;
-
-  if(ret>=0)
-    *info = ret+1;
-
-  return 0;
-}
-
-//GETRS solves a system of linear equations
-//    A * X = B  or  A' * X = B
-//  with a general N-by-N matrix A using the LU factorization computed  by GETRF
-EIGEN_LAPACK_FUNC(getrs,(char *trans, int *n, int *nrhs, RealScalar *pa, int *lda, int *ipiv, RealScalar *pb, int *ldb, int *info))
-{
-  *info = 0;
-        if(OP(*trans)==INVALID)  *info = -1;
-  else  if(*n<0)                 *info = -2;
-  else  if(*nrhs<0)              *info = -3;
-  else  if(*lda<std::max(1,*n))  *info = -5;
-  else  if(*ldb<std::max(1,*n))  *info = -8;
-  if(*info!=0)
-  {
-    int e = -*info;
-    return xerbla_(SCALAR_SUFFIX_UP"GETRS", &e, 6);
-  }
-
-  Scalar* a = reinterpret_cast<Scalar*>(pa);
-  Scalar* b = reinterpret_cast<Scalar*>(pb);
-  MatrixType lu(a,*n,*n,*lda);
-  MatrixType B(b,*n,*nrhs,*ldb);
-
-  for(int i=0; i<*n; ++i)
-    ipiv[i]--;
-  if(OP(*trans)==NOTR)
-  {
-    B = PivotsType(ipiv,*n) * B;
-    lu.triangularView<UnitLower>().solveInPlace(B);
-    lu.triangularView<Upper>().solveInPlace(B);
-  }
-  else if(OP(*trans)==TR)
-  {
-    lu.triangularView<Upper>().transpose().solveInPlace(B);
-    lu.triangularView<UnitLower>().transpose().solveInPlace(B);
-    B = PivotsType(ipiv,*n).transpose() * B;
-  }
-  else if(OP(*trans)==ADJ)
-  {
-    lu.triangularView<Upper>().adjoint().solveInPlace(B);
-    lu.triangularView<UnitLower>().adjoint().solveInPlace(B);
-    B = PivotsType(ipiv,*n).transpose() * B;
-  }
-  for(int i=0; i<*n; ++i)
-    ipiv[i]++;
-
-  return 0;
-}
diff --git a/lapack/svd.cpp b/lapack/svd.cpp
deleted file mode 100644
index 83544cf..0000000
--- a/lapack/svd.cpp
+++ /dev/null
@@ -1,138 +0,0 @@
-// This file is part of Eigen, a lightweight C++ template library
-// for linear algebra.
-//
-// Copyright (C) 2014 Gael Guennebaud <gael.guennebaud@inria.fr>
-//
-// This Source Code Form is subject to the terms of the Mozilla
-// Public License v. 2.0. If a copy of the MPL was not distributed
-// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
-
-#include "lapack_common.h"
-#include <Eigen/SVD>
-
-// computes the singular values/vectors a general M-by-N matrix A using divide-and-conquer
-EIGEN_LAPACK_FUNC(gesdd,(char *jobz, int *m, int* n, Scalar* a, int *lda, RealScalar *s, Scalar *u, int *ldu, Scalar *vt, int *ldvt, Scalar* /*work*/, int* lwork,
-                         EIGEN_LAPACK_ARG_IF_COMPLEX(RealScalar */*rwork*/) int * /*iwork*/, int *info))
-{
-  // TODO exploit the work buffer
-  bool query_size = *lwork==-1;
-  int diag_size = (std::min)(*m,*n);
-  
-  *info = 0;
-        if(*jobz!='A' && *jobz!='S' && *jobz!='O' && *jobz!='N')  *info = -1;
-  else  if(*m<0)                                                  *info = -2;
-  else  if(*n<0)                                                  *info = -3;
-  else  if(*lda<std::max(1,*m))                                   *info = -5;
-  else  if(*lda<std::max(1,*m))                                   *info = -8;
-  else  if(*ldu <1 || (*jobz=='A' && *ldu <*m)
-                   || (*jobz=='O' && *m<*n && *ldu<*m))           *info = -8;
-  else  if(*ldvt<1 || (*jobz=='A' && *ldvt<*n)
-                   || (*jobz=='S' && *ldvt<diag_size)
-                   || (*jobz=='O' && *m>=*n && *ldvt<*n))         *info = -10;
-  
-  if(*info!=0)
-  {
-    int e = -*info;
-    return xerbla_(SCALAR_SUFFIX_UP"GESDD ", &e, 6);
-  }
-  
-  if(query_size)
-  {
-    *lwork = 0;
-    return 0;
-  }
-  
-  if(*n==0 || *m==0)
-    return 0;
-  
-  PlainMatrixType mat(*m,*n);
-  mat = matrix(a,*m,*n,*lda);
-  
-  int option = *jobz=='A' ? ComputeFullU|ComputeFullV
-             : *jobz=='S' ? ComputeThinU|ComputeThinV
-             : *jobz=='O' ? ComputeThinU|ComputeThinV
-             : 0;
-
-  BDCSVD<PlainMatrixType> svd(mat,option);
-  
-  make_vector(s,diag_size) = svd.singularValues().head(diag_size);
-
-  if(*jobz=='A')
-  {
-    matrix(u,*m,*m,*ldu)   = svd.matrixU();
-    matrix(vt,*n,*n,*ldvt) = svd.matrixV().adjoint();
-  }
-  else if(*jobz=='S')
-  {
-    matrix(u,*m,diag_size,*ldu)   = svd.matrixU();
-    matrix(vt,diag_size,*n,*ldvt) = svd.matrixV().adjoint();
-  }
-  else if(*jobz=='O' && *m>=*n)
-  {
-    matrix(a,*m,*n,*lda)   = svd.matrixU();
-    matrix(vt,*n,*n,*ldvt) = svd.matrixV().adjoint();
-  }
-  else if(*jobz=='O')
-  {
-    matrix(u,*m,*m,*ldu)        = svd.matrixU();
-    matrix(a,diag_size,*n,*lda) = svd.matrixV().adjoint();
-  }
-    
-  return 0;
-}
-
-// computes the singular values/vectors a general M-by-N matrix A using two sided jacobi algorithm
-EIGEN_LAPACK_FUNC(gesvd,(char *jobu, char *jobv, int *m, int* n, Scalar* a, int *lda, RealScalar *s, Scalar *u, int *ldu, Scalar *vt, int *ldvt, Scalar* /*work*/, int* lwork,
-                         EIGEN_LAPACK_ARG_IF_COMPLEX(RealScalar */*rwork*/) int *info))
-{
-  // TODO exploit the work buffer
-  bool query_size = *lwork==-1;
-  int diag_size = (std::min)(*m,*n);
-  
-  *info = 0;
-        if( *jobu!='A' && *jobu!='S' && *jobu!='O' && *jobu!='N') *info = -1;
-  else  if((*jobv!='A' && *jobv!='S' && *jobv!='O' && *jobv!='N')
-           || (*jobu=='O' && *jobv=='O'))                         *info = -2;
-  else  if(*m<0)                                                  *info = -3;
-  else  if(*n<0)                                                  *info = -4;
-  else  if(*lda<std::max(1,*m))                                   *info = -6;
-  else  if(*ldu <1 || ((*jobu=='A' || *jobu=='S') && *ldu<*m))    *info = -9;
-  else  if(*ldvt<1 || (*jobv=='A' && *ldvt<*n)
-                   || (*jobv=='S' && *ldvt<diag_size))            *info = -11;
-  
-  if(*info!=0)
-  {
-    int e = -*info;
-    return xerbla_(SCALAR_SUFFIX_UP"GESVD ", &e, 6);
-  }
-  
-  if(query_size)
-  {
-    *lwork = 0;
-    return 0;
-  }
-  
-  if(*n==0 || *m==0)
-    return 0;
-  
-  PlainMatrixType mat(*m,*n);
-  mat = matrix(a,*m,*n,*lda);
-  
-  int option = (*jobu=='A' ? ComputeFullU : *jobu=='S' || *jobu=='O' ? ComputeThinU : 0)
-             | (*jobv=='A' ? ComputeFullV : *jobv=='S' || *jobv=='O' ? ComputeThinV : 0);
-  
-  JacobiSVD<PlainMatrixType> svd(mat,option);
-  
-  make_vector(s,diag_size) = svd.singularValues().head(diag_size);
-  {
-        if(*jobu=='A') matrix(u,*m,*m,*ldu)           = svd.matrixU();
-  else  if(*jobu=='S') matrix(u,*m,diag_size,*ldu)    = svd.matrixU();
-  else  if(*jobu=='O') matrix(a,*m,diag_size,*lda)    = svd.matrixU();
-  }
-  {
-        if(*jobv=='A') matrix(vt,*n,*n,*ldvt)         = svd.matrixV().adjoint();
-  else  if(*jobv=='S') matrix(vt,diag_size,*n,*ldvt)  = svd.matrixV().adjoint();
-  else  if(*jobv=='O') matrix(a,diag_size,*n,*lda)    = svd.matrixV().adjoint();
-  }
-  return 0;
-}
\ No newline at end of file