#include <stdio.h>
	#include <stdlib.h>
	#include <malloc.h>
	#include <mkl.h>
	#include <math.h>
	#include <time.h>
	#include <string.h>

	#include <windows.h>
	int main(int argc, char* argv[])
	{

	int nnz=4746;
	int n=459;
	int i,j;

	MKL_INT mtype = 11; /* Real unsymmetric matrix */
	MKL_INT nrhs = 1; /* Number of right hand sides. */
	/* Pardiso control parameters. */
	MKL_INT maxfct, mnum, phase, error, msglvl;
	/* Auxiliary variables. */
	double ddum; /* Double dummy */
	MKL_INT idum; /* Integer dummy. */


	FILE *valf=fopen("a.txt","r");
	FILE *colIndf=fopen("ja.txt","r");
	FILE *rowPtrf=fopen("ia.txt","r");
	if (valf==NULL) {
	printf("file open failed");
	return 0;
	}


	double *a=(double*)malloc(nnz*sizeof(double));
	int *ja=(int*)malloc(nnz*sizeof(int));
	int *ia=(int*)malloc((n+1)*sizeof(int));

	for (i=0;i<nnz;i++) {
	fscanf(valf,"%lf",&a[i]);
	}

	for (i=0;i<nnz;i++) {
	fscanf(colIndf,"%d",&ja[i]);

	}
   
	for (i=0;i<n+1;i++) {
	fscanf(rowPtrf,"%d",&ia[i]);
	}

	fclose(valf);
	fclose(colIndf);
	fclose(rowPtrf);

	double *b=(double*)malloc(nrhs*n*sizeof(double));
	double *x=(double*)malloc(nrhs*n*sizeof(double));
	for (i=0;i<nrhs*n;i++) {
	x[i]=0;
	}
	printf("val: %lf\n",a[0]);
	printf("ia[n]: %d /t ia[n+1]: %d\n",ia[n-1], ia[n]);
	ia[n]=4746;

	for (i=0;i<nrhs*n;i++) {
	b[i]=0;
	}

	/* 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 */
	int *pt[64];
	MKL_INT *perm, *iparm;
	perm = (MKL_INT*) malloc (sizeof (MKL_INT)*n);
	iparm = (MKL_INT*) malloc (sizeof (MKL_INT)*64);
	/* -------------------------------------------------------------------- */
	/* .. Initialize the internal solver memory pointer. This is only */
	/* necessary for the FIRST call of the PARDISO solver. */
	/* -------------------------------------------------------------------- */
	for (i = 0; i < 64; i++) 
	{
		pt[i] = 0;
	}
	/* -------------------------------------------------------------------- */
	/* .. Initialize the internal solver memory pointer. This is only */
	/* necessary for the FIRST call of the PARDISO solver. */
	/* -------------------------------------------------------------------- */
	for (i = 0; i < n; i++) 
	{
		perm[i] = 0;
	}
	/* -------------------------------------------------------------------- */
	/* .. Setup Pardiso control parameters. */
	/* -------------------------------------------------------------------- */
	for (i = 0; i < 64; i++) 
	{
		iparm[i] = 0;
	}
	iparm[0] = 1;	/*	Default Values	*/
	iparm[26] = 1; /* mtype 11 double precision */
	iparm[34] = 1; /* = 1 for zero-based indexing; = 0 for one-based indexing */
	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 */

	phase = 11;
	pardiso (pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, perm, &nrhs, iparm, &msglvl, b, x, &error);
	if (error != 0) {
		printf("\nERROR during symbolic factorization: %d", error);
		exit(1);
	}
	printf("\nReordering completed ... ");
	printf("\nNumber of nonzeros in factors = %d", iparm[17]);
	printf("\nNumber of factorization MFLOPS = %d", iparm[18]);

	phase = 22;
	pardiso (pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, perm, &nrhs, iparm, &msglvl, b, x, &error);
		if (error != 0)
		{
			printf ("\nERROR during numerical factorization: %lld\n", error);
		    exit (2);
		}
		printf ("\nFactorization completed ... ");

		 phase = 331;
		pardiso (pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, perm, &nrhs, iparm, &msglvl, b, x, &error);
		if (error != 0)
	{
		printf ("\nERROR during numerical factorization: %lld\n", error);
		exit (3);
	}
	printf ("\nFactorization completed ... ");

	phase = -1;
	pardiso (pt, &maxfct, &mnum, &mtype, &phase, &n, a, ia, ja, perm, &nrhs, iparm, &msglvl, b, x, &error);
		if (error != 0)
	{
		printf ("\nERROR during release: %lld\n", error);
		exit (-1);
	}
	
	system("pause");

	return 0;
	}