//
// test pardiso solver
//
#include <fstream>
#include <complex>
#include "mkl.h"
#include <iomanip>
#include <iostream>
using namespace std;

void debug_pardiso_loaddump_all(int &maxfct, int &mnum, int &mtype, int &phase, int &nTotalUnknown, complex<double> * &a, int *&ia, 
				int *&ja, int *&perm, int &nrhs, int *&iparm, int &msglvl, complex<double> *& vectorB, complex<double> *&vectorX, int &error)
{
	int NUM64 = 64;

	ifstream debugFile;
	debugFile.open("pardiso.dat", std::ios_base::in);

	debugFile >> nTotalUnknown;
	ia = new int[nTotalUnknown + 1];
	for (int i=0; i <= nTotalUnknown; ++i)
		debugFile >> ia[i];

        ja = new int [ia[nTotalUnknown]-1];
	for ( int i=0; i < ia[nTotalUnknown]-1; ++i )
		debugFile >> ja[i];

	a = new complex<double> [ia[nTotalUnknown]-1];
	double dr, di;
	for ( int i=0; i < ia[nTotalUnknown]-1; ++i )
	{
	        debugFile >> dr >> di; 
            a[i] = complex<double>(dr,di); 
        }
	
	perm = new int [nTotalUnknown];

	debugFile.close();
}

int main()
{
    int *ia;
    int *ja;
    complex<double> * a;

    int  *iparm = new int[64];
    void *pt[64];
    for ( int i = 0; i < 64; ++i )
    {
        iparm[i] = 0; 
        pt[i] = 0;
    }

    // Symbolic factorization
    int phase = 11;
    int maxfct = 1;
    int mnum   = 1;
    int mtype = 6;
    int *perm;
    int nrhs = 1;
    int msglvl = 1;
    int	 PARDISOerror = 0;
    complex<double> * dumm  = new complex<double>;
    int nUnknowns;

    debug_pardiso_loaddump_all(maxfct, mnum, mtype, phase, nUnknowns,a, ia, 
				ja, perm, nrhs, iparm, msglvl, dumm, dumm, PARDISOerror);

    iparm[0] = 1; 
    iparm[1] = 2; 
    iparm[7] = 1; 
    iparm[9] = 13; 


    /////////////////////////////////////////////////////////
    // Change this parameter to produce different error flags
    iparm[23] = 1; 
    /////////////////////////////////////////////////////////

    iparm[24] = 1; 
    iparm[30] = 1; // partial solve
    iparm[59] = 1; // disable OOC

    for (int i = 0; i< 64; ++i)
	std::cout << i << "= " << iparm[i]  << std::endl;

    PARDISO(pt, &maxfct, &mnum, &mtype, &phase, &nUnknowns, a, ia, 
        ja, perm, &nrhs, iparm, &msglvl, dumm, dumm, &PARDISOerror);

    std::cout << "pardiso error code: " << PARDISOerror << std::endl;

    delete [] iparm;
    return 0;
}
