#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "mkl_pardiso.h"
#include "mkl_types.h"

MKL_INT
main (int argc,char *argv[]){
  /* Matrix data. */
  MKL_INT n,nnz,*ia,*ja;
  double *a,*b,*bs,*x,res,res0;

  MKL_INT mtype = 11;		/* Real unsymmetric matrix */
  /* RHS and solution vectors. */
  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;
  FILE *fil; char vers[198]; int nbad;

    if(argc < 2){
       fprintf(stderr,"Usage: unsym <Fortran unformatted file>\n");
       exit(1);
       }
    if(!(fil=fopen(argv[1],"rb"))){
       fprintf(stderr,"Cannot open %s for reading\n",argv[1]);
       exit(1);
       }
    fread(&idum,sizeof(int),1,fil); //header
    fread(&n,sizeof(int),1,fil); fread(&nnz,sizeof(int),1,fil);
    fread(&idum,sizeof(int),1,fil); //trailer

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

    fread(&idum,sizeof(int),1,fil); //header
    fread(a,sizeof(double),nnz,fil);
    fread(ja,sizeof(int),nnz,fil);
    fread(ia,sizeof(int),n+1,fil);
    fread(b,sizeof(double),n,fil);
    fread(&idum,sizeof(int),1,fil); //trailer

    fclose(fil);

    mkl_get_version_string (vers, 198);
    puts(vers);
/* -------------------------------------------------------------------- */
/* .. Setup Pardiso control parameters. */
/* -------------------------------------------------------------------- */
  for (i = 0; i < 64; i++)iparm[i] = 0;
  iparm[0] = 1;			/* No solver default */
  iparm[1] = 2;			/* 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] = 5;			/* 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 */
  maxfct = 1;			/* Maximum number of numerical factorizations. */
  mnum = 1;			/* Which factorization to use. */
  msglvl = 0;			/* Print statistical information in file */
  error = 0;			/* Initialize error flag */
/* -------------------------------------------------------------------- */
/* .. 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;
    }
/* -------------------------------------------------------------------- */
/* .. Reordering and Symbolic Factorization. This step also allocates */
/* all memory that is necessary for the factorization. */
/* -------------------------------------------------------------------- */
  phase = 11;
  PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
	   &n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &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]);
/* -------------------------------------------------------------------- */
/* .. Numerical factorization. */
/* -------------------------------------------------------------------- */
  phase = 22;
  PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
	   &n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, &ddum, &ddum, &error);
  if (error != 0)
    {
      printf ("\nERROR during numerical factorization: %d", error);
      exit (2);
    }
  printf ("\nFactorization completed ... ");
/* -------------------------------------------------------------------- */
/* .. Back substitution and iterative refinement. */
/* -------------------------------------------------------------------- */
  phase = 33;
  iparm[11] = 0;  // Solve Ax = b
  uplo = "non-transposed";
  printf ("\n\nSolving %s system...\n", uplo);
  PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
       &n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, b, x, &error);
  if (error != 0)	{
	  printf ("\nERROR during solution: %d", error);
	  exit (3);
	}
/*
      printf("\nScanning for Indefs. in x[]..."); nbad=0;
      for(j=0; j<n; j++){
         union {double x; int i[2];} xi;
         xi.x=x[j];
         if(xi.i[1]==0xFFF80000){printf("%4d ",j+1); nbad++;}
         }
      printf("\n%d indefinite x[] values\n",nbad);
      printf("\nScanning for zeroes in x[]..."); nbad=0;
      for(j=0; j<n; j++){
         if(x[j]==0.0){printf("%4d ",j+1); nbad++;}
         }
      printf("\n%d zero x[] values\n",nbad); */
      printf ("\nThe solution of the system is (first ten shown) : ");
      for (j = 0; j < 10; j++){
         printf ("\n x [%4d] = %12.4e",j+1, x[j]);
         }
      printf ("\n . . .");
      for (j = n-11; j < n; j++){
         printf ("\n x [%4d] = %12.4e",j+1, x[j]);
         }

// Compute residual
      printf("\nComputing residuals");
      mkl_dcsrgemv (uplo, &n, a, ia, ja, x, bs);
      res = 0.0;
      res0 = 0.0;
      for (j = 0; j < n; j++){
	  res  += (bs[j] - b[j]) * (bs[j] - b[j]);
	  res0 += b[j] * b[j];
	}
      res0 = sqrt (res) / sqrt (res0);
      printf ("\nResid_norm, relative residual = %e %e\n", res,res0);
// Check residual
      if (res > 1e-10)	{
	  printf ("Error: residual is too high!\n");
	  exit (10);
    }

/* -------------------------------------------------------------------- */
/* .. Termination and release of memory. */
/* -------------------------------------------------------------------- */
  phase = -1;			/* Release internal memory. */
  PARDISO (pt, &maxfct, &mnum, &mtype, &phase,
     &n, &ddum, ia, ja, &idum, &nrhs,
     iparm, &msglvl, &ddum, &ddum, &error);
  return 0;
}