/*
  CLAPACK demo to illustrate the solution of a linear system Ax=b
  using LU decomposition.
  You may explore many other variants at www.netlib.org/lapack
  J.K. Johnstone, 2004
*/

#include <iostream>
using namespace std;

// to avoid C++ incompatibilities (after all, it is expecting C)
extern "C" {               // make sure these are visible in directory lookup
#include "f2c.h"           // our solution: copy them to /usr/include
#include "clapack.h"
}

int main (int argc, char **argv)
{
  // goal is to solve Ax=b; A will also be overwritten with its LU decomposition
  float *A;  // matrix as a 1d array in column-major order (a la Fortran)
  A = new float[9];
  float *B;  // on entry, b; on exit, x
  B = new float[3];

  // example from p. 99 Golub and van Loan
  int i;
  for (i=0; i<8; i++) A[i] = i+1;  A[8] = 10;
  for (i=0; i<3; i++) B[i] = 1;
  cout << "A = (";
  for (i=0; i<8; i++) cout << A[i] << ",";
  cout << A[8] << ")" << endl;

  // use call-by-reference parameters, in another adaptation to Fortran
  integer N=3, NRHS=1, LDA=3, LDB=3, INFO;  
  integer *IPIV; IPIV = new integer[3];

  sgesv_ (&N, &NRHS, A, &LDA, IPIV, B, &LDB, &INFO); 
                   // don't forget the extra _ at the end of the function name
  cout << "INFO = " << INFO << endl;
  if (INFO>=0)
    {
      cout << "Solution x = (" << B[0] << "," 
	                       << B[1] << "," 
	                       << B[2] << ")" << endl;
      cout << "IPIV = (" << IPIV[0] << "," << IPIV[1] << "," << IPIV[2] 
	   << ")" << endl;
      cout << "LU decomposition: (";
      for (i=0; i<8; i++) cout << A[i] << ",";
      cout << A[8] << ")" << endl;
    }
  delete [] A;
  delete [] B;
  delete [] IPIV;
  return 0;
}