#include <iostream>
#include <math.h>

using namespace std;

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

//sparse rhs iparm[30]=1
int main(void){

  void *pt[64];
  int mtype = 2; //spd
  int solver = 0;
  int iparm[64];

  for(int i=0; i<64; i++) iparm[i] = 0;
  iparm[0] = 1;
  iparm[24] = 0;
  iparm[30] = 1;
  for(int i=0; i<64; i++) pt[i] = 0;

  iparm[2] = 1;
  iparm[51] = 1;
  int error = 0;

  if(error)
    cerr<<endl<<"ERROR: IN pardisoinit. error = "<<error<<endl;
  
  // factorization and solve
  // Identity matrix
  //              1 2 3 4 5 6 7 8 9 10 11 12 13 14
  //         1  / 1 0               0             \
  //         2 |    1             0    0           |
  //         3 |      1 0         0       0        |
  //         4 |        1                    0     |
  //         5 |          1 0       0              |
  //         6 |            1          0        0  |
  //         7 |              1 0         0     0  |
  //         8 |                1            0     |
  //         9 |                  1                |
  //         10|                    1              |
  //         11|                       1           |
  //         12|                          1        |
  //         13|                             1     |
  //         14\                                1  /

  int maxfct = 1;
  int mnum = 1;
  int phase = 12;

  int n = 14;
  int nz = 30;
  double a[30] = {1,0,              0,
		    1,            0,   0,           
		      1,0,        0,      0,        
		        1,                   0,     
		          1,0,      0,             
		            1,         0,       0,  
		              1,0,        0,    0,  
                                1,           0,     
		                  1,                
		                    1,             
		                       1,           
		                          1,        
		                             1,     
		                                1 };
  

  int ia[15] = {1,4,7,11,13,16,19,23,25,26,27,28,29,30,31}; //fortran notation

  int ja[30]    = {1,2,              10,
		     2,           9,   11,           
		      3,4,        9,      12,        
		        4,                   13,     
		          5,6,      10,             
		            6,         11,       14,  
		              7,8,        12,    14,  
                                8,           13,     
		                  9,                
		                    10,             
		                       11,           
		                          12,        
		                             13,     
		                                14 }; //fortran notation
  
  int perm[n];
  int nrhs = 1;
  
  iparm[0] = 1;
  iparm[2] = 1;
  iparm[51] = 1;
  int msglvl = 0;
  double b[n];
  for(int i=0; i<n; i++)
    b[i] = 0;

  double x[n];
  for(int i=0; i<n; i++)
    x[i] = 0;
  
  iparm[24] = 0;
  iparm[30] = 1; // sparse rhs and sparse sol
  for(int i=0; i<n; i++)
    perm[i] = 0;
  perm[8] = 1;
  perm[9] = 1;
  perm[10] = 1;
  perm[11] = 1;
  perm[12] = 1;
  //perm[13] = 1; comment out perm[13] it will be OK

  cout<<endl<<"perm before reordering"<<endl;
  for(int i=0; i<n; i++) cout<<perm[i]<<" ";

  pardiso(pt, &maxfct, &mnum, &mtype, &phase,
	  &n, a, ia, ja, perm, &nrhs,
	  iparm, &msglvl, b, x, &error);

  cout<<endl<<"perm after reordering"<<endl;
  for(int i=0; i<n; i++) cout<<perm[i]<<" ";

  if(error){
    cerr<<endl<<"ERROR: WHEN ANALYSIS AND FACTORIZATION IN pardiso"<<endl
	<<"error = "<<error<<endl<<endl;
    return error;
  }

  // SOLVE
  phase = 33;
  pardiso(pt, &maxfct, &mnum, &mtype, &phase,
	  &n, a, ia, ja, perm, &nrhs,
	  iparm, &msglvl, b, x, &error);
  if(error){
    cerr<<endl<<"ERROR: WHEN SOLVE IN pardiso"<<endl
	<<"error = "<<error<<endl;
    return error;
  }
  
  // RELEASE MEMORY
  phase = -1; //release memory
  pardiso(pt, &maxfct, &mnum, &mtype, &phase,
	  &n, a, ia, ja, perm, &nrhs,
	  iparm, &msglvl, b, x, &error);
  if(error)
    cerr<<endl<<"ERROR: WHEN RELEASE MEMORY IN pardiso"<<endl
	<<"error = "<<error<<endl;
  
  return error;

}