tnt_lapack.h

#ifndef TNT_LAPACK_HPP
#define TNT_LAPACK_HPP
 
extern "C"
{
#include "f2c.h"
#include "clapack.h"
}
 
 
namespace TNT
{
 
  Fortran_Array2D<double> Lapack_LinearSolve( const Fortran_Array2D<double> &A,  
      const Fortran_Array2D<double> &B)
  {
      
      integer m = A.dim1();
      integer n = A.dim2();
 
      if (m != n)
        return Fortran_Array2D<double>();
 
      integer nrhs = B.dim2();
 
      if (n != B.dim1())
        return Fortran_Array2D<double>();
 
 
      Fortran_Array2D<double> X(B.copy());
      Fortran_Array2D<double> lu(A.copy());
      Fortran_Array1D<integer>  ipiv(m, integer(0));
 
      integer  info=0;
 
      dgesv_( &n, &nrhs, &lu(1,1), &m, &ipiv(1), &X(1,1), &n, &info);
 
      return X;
  }
 
 
 
}
 
 
#endif
// TNT_LAPACK_HPP