#include "intflapack.h"

extern "C" {
void sgesv_(int_4* n, int_4* nrhs, r_4* a, int_4* lda, 
            int_4* ipiv, r_4* b, int_4* ldb, int_4* info);
void dgesv_(int_4* n, int_4* nrhs, r_8* a, int_4* lda, 
            int_4* ipiv, r_8* b, int_4* ldb, int_4* info);
void cgesv_(int_4* n, int_4* nrhs, complex<r_4>* a, int_4* lda, 
            int_4* ipiv, complex<r_4>* b, int_4* ldb, int_4* info);
void zgesv_(int_4* n, int_4* nrhs, complex<r_8>* a, int_4* lda, 
            int_4* ipiv, complex<r_8>* b, int_4* ldb, int_4* info);
}

void rztest_lapack(TArray<r_4>& aa, TArray<r_4>& bb)
{
  if ( aa.NbDimensions() != 2 ) throw(SzMismatchError("rztest_lapack(TMatrix<r_4> A Not a Matrix"));
  if ( aa.SizeX() !=  aa.SizeY()) throw(SzMismatchError("rztest_lapack(TMatrix<r_4> A Not a square Matrix"));
  if ( bb.NbDimensions() != 2 ) throw(SzMismatchError("rztest_lapack(TMatrix<r_4> A Not a Matrix"));
  if ( bb.SizeX() !=  aa.SizeX() ) throw(SzMismatchError("rztest_lapack(TMatrix<r_4> A <> B "));
  if ( !bb.IsPacked() || !bb.IsPacked() ) 
    throw(SzMismatchError("rztest_lapack(TMatrix<r_4> Not packed A or B "));
  
  int_4 n = aa.SizeX();
  int_4 nrhs = bb.SizeY();
  int_4 lda = n;
  int_4 ldb = bb.SizeX();
  int_4 info;
  int_4* ipiv = new int_4[n];
  sgesv_(&n, &nrhs, aa.Data(), &lda, ipiv, bb.Data(), &ldb, &info);
  cout << "rztest_lapack/Info= " << info << endl;
  cout << aa << "\n" << bb << endl;
  return;
}