#ifdef newc
#include <iostream>
#include <iomanip>
#include <fstream>
#include <cstdlib>
#include <cstring>
#include <cmath>
using namespace std;
#else
#include <iostream.h>
#include <iomanip.h>
#include <fstream.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#endif

// Intel oneMKL LAPACK interface
#include <mkl_lapacke.h>
/* Linear equation solution using Intel oneMKL LAPACK.
a[0..n-1][0..n-1] is the input matrix. b[0..n-1] is input
containing the right-hand side vectors. On output a is
replaced by its matrix inverse, and b is replaced by the
corresponding set of solution vectors.

Mathematical equivalence:
  Solves: A * x = b  =>  x = A^(-1) * b
  Original Gauss-Jordan and LAPACK dgesv/dgetri produce identical results
  within numerical precision. */

int gaussj(double *a, double *b, int n)
{
  // Allocate pivot array and workspace
  lapack_int *ipiv = new lapack_int[n];
  lapack_int info;

  // Make a copy of matrix a for solving (dgesv modifies it to LU form)
  double *a_copy = new double[n * n];
  for (int i = 0; i < n * n; i++) {
    a_copy[i] = a[i];
  }

  // Step 1: Solve linear system A*x = b using LU decomposition
  // LAPACKE_dgesv uses column-major by default, but we use row-major
  info = LAPACKE_dgesv(LAPACK_ROW_MAJOR, n, 1, a_copy, n, ipiv, b, 1);

  if (info != 0) {
    cout << "gaussj: Singular Matrix (dgesv info=" << info << ")" << endl;
    delete[] ipiv;
    delete[] a_copy;
    return 1;
  }

  // Step 2: Compute matrix inverse A^(-1) using LU factorization
  // First do LU factorization of original matrix a
  info = LAPACKE_dgetrf(LAPACK_ROW_MAJOR, n, n, a, n, ipiv);

  if (info != 0) {
    cout << "gaussj: Singular Matrix (dgetrf info=" << info << ")" << endl;
    delete[] ipiv;
    delete[] a_copy;
    return 1;
  }

  // Then compute inverse from LU factorization
  info = LAPACKE_dgetri(LAPACK_ROW_MAJOR, n, a, n, ipiv);

  if (info != 0) {
    cout << "gaussj: Singular Matrix (dgetri info=" << info << ")" << endl;
    delete[] ipiv;
    delete[] a_copy;
    return 1;
  }

  delete[] ipiv;
  delete[] a_copy;

  return 0;
}
// for check usage
/*
int main()
{
  double *A,*b;
  A=new double[9];
  b=new double[3];

  A[0]=0.5; A[1]=1.0/3; A[2]=1;
  A[3]=1;   A[4]=5.0/3; A[5]=3;
  A[6]=2;   A[7]=4.0/3; A[8]=5;

  b[0]=1; b[1]=3; b[2]=2;

  cout<<"initial data:"<<endl;
  for(int i=0;i<3;i++) cout<<A[i*3]<<" "<<A[i*3+1]<<" "<<A[i*3+2]<<" "<<b[i]<<endl;

  gaussj(A, b, 3);

  cout<<"final data:"<<endl;
  for(int i=0;i<3;i++) cout<<A[i*3]<<" "<<A[i*3+1]<<" "<<A[i*3+2]<<" "<<b[i]<<endl;

  delete[] A; delete[] b;
}
*/