#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sys/mman.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <assert.h>
#include "mkl.h"

#define alloc(TYPE, SIZE) (TYPE *)mkl_malloc(sizeof(TYPE) * (SIZE), 64)

int main (int argc, char ** argv)
{
	/* Matrix data. */
	srand(1);
	assert(argc == 3);
	MKL_INT n = atoi(argv[1]);
	MKL_INT nz = atoi(argv[2]);
	MKL_INT mtype = 11;		/* Real unsymmetric matrix */
	/* RHS and solution vectors. */

	FILE * data_file = fopen("data.csr", "rb");
	double * data = alloc(double, nz);;
	fread(data, sizeof(double), nz, data_file);

	FILE * rows_file = fopen("rows.csr", "rb");
	MKL_INT * irn = alloc(MKL_INT, n + 1);
	fread(irn, sizeof(MKL_INT), n + 1, rows_file);

	FILE * cols_file = fopen("cols.csr", "rb");
	MKL_INT * jcn = alloc(MKL_INT, nz);
	fread(jcn, sizeof(MKL_INT), nz, cols_file);

	double * b = alloc(double, n);
	double * x = alloc(double, n);
	double * bs = alloc(double, n);
	double res, res0;
	MKL_INT nrhs = 1;		/* Number of right hand sides. */
	/* Internal solver memory pointer pt, */
	/* 32-bit: int pt[64]; 64-bit: long int pt[64] */
	/* or void *pt[64] should be OK on both architectures */
	void *pt[64];
	/* Pardiso control parameters. */
	MKL_INT iparm[64];
	MKL_INT maxfct, mnum, phase, error, msglvl;
	/* Auxiliary variables. */
	MKL_INT i, j;
	double ddum;			/* Double dummy */
	MKL_INT idum;			/* Integer dummy. */
	char *uplo;
	/* -------------------------------------------------------------------- */
	/* .. Setup Pardiso control parameters. */
	/* -------------------------------------------------------------------- */
	for (i = 0; i < 64; i++)
	{
		iparm[i] = 0;
	}
	iparm[0] = 1;			/* No solver default */
	iparm[1] = 3;			/* Fill-in reordering from METIS */
	/* Numbers of processors, value of OMP_NUM_THREADS */
	iparm[2] = 1;
	iparm[3] = 0;			/* No iterative-direct algorithm */
	iparm[4] = 0;			/* No user fill-in reducing permutation */
	iparm[5] = 0;			/* Write solution into x */
	iparm[6] = 0;			/* Not in use */
	iparm[7] = 2;			/* Max numbers of iterative refinement steps */
	iparm[8] = 0;			/* Not in use */
	iparm[9] = 13;		/* Perturb the pivot elements with 1E-13 */
	iparm[10] = 1;		/* Use nonsymmetric permutation and scaling MPS */
	iparm[11] = 0;		/* Conjugate transposed/transpose solve */
	iparm[12] = 1;		/* Maximum weighted matching algorithm is switched-on (default for non-symmetric) */
	iparm[13] = 0;		/* Output: Number of perturbed pivots */
	iparm[14] = 0;		/* Not in use */
	iparm[15] = 0;		/* Not in use */
	iparm[16] = 0;		/* Not in use */
	iparm[17] = -1;		/* Output: Number of nonzeros in the factor LU */
	iparm[18] = -1;		/* Output: Mflops for LU factorization */
	iparm[19] = 0;		/* Output: Numbers of CG Iterations */
	iparm[26] = 1;
	iparm[34] = 1;
	iparm[59] = 1;
	maxfct = 1;			/* Maximum number of numerical factorizations. */
	mnum = 1;			/* Which factorization to use. */
	msglvl = 1;			/* Print statistical information in file */
	error = 0;			/* Initialize error flag */
	for (i = 0; i < 64; i++)
	{
		pt[i] = 0;
	}
	phase = 11;
	pardiso_64 (pt, &maxfct, &mnum, &mtype, &phase,
			&n, data, irn, jcn, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);
	if (error != 0)
	{
		printf ("\nERROR during symbolic factorization: %lld", error);
		exit (1);
	}
	phase = 22;
	pardiso_64 (pt, &maxfct, &mnum, &mtype, &phase,
			&n, data, irn, jcn, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);
	if (error != 0)
	{
		printf ("\nERROR during numerical factorization: %lld", error);
		exit (2);
	}
	phase = 33;
	b[0] = 1;
	for (i = 1; i < n; i++)
	{
		b[i] = 0;
	}
	iparm[11] = 2;
	pardiso_64 (pt, &maxfct, &mnum, &mtype, &phase,
			&n, data, irn, jcn, &idum, &nrhs, iparm, &msglvl, b, x, &error);
	if (error != 0)
	{
		printf ("\nERROR during solution: %lld", error);
		exit (3);
	} else {
		for (i = 0; i < 10; i++) {
			printf("%f\n", x[i]);
		}

	}

	phase = -1;			/* Release internal memory. */
	pardiso_64 (pt, &maxfct, &mnum, &mtype, &phase,
			&n, &ddum, irn, jcn, &idum, &nrhs,
			iparm, &msglvl, &ddum, &ddum, &error);
	return 0;
}