// test_pardiso.cpp : This file contains the 'main' function. Program execution begins and ends there.
//

#include <iostream>

#include "mkl_pardiso.h"

#include "mkl_blas.h"
#include "mkl_rci.h"
#include "mkl_spblas.h"

#define IFORMAT "%i"
#define MKL int


void test_pardiso()
{
    // From example pardiso_unsym.c
    
    /* Matrix data. */
    MKL_INT n = 5;
    MKL_INT ia[6] = { 1, 4, 6, 9, 12, 14 };
    MKL_INT ja[13] =
    { 1, 2,    4,
      1, 2,
            3, 4, 5,
      1,    3, 4,
         2,       5
    };
    double a[13] =
    { 1.0,-1.0,     -3.0,
     -2.0, 5.0,
                4.0, 6.0, 4.0,
     -4.0,      2.0, 7.0,
           8.0,          -5.0
    };
    MKL_INT mtype = 11;       /* Real unsymmetric matrix */
    // Descriptor of main sparse matrix properties
    struct matrix_descr descrA;
    // Structure with sparse matrix stored in CSR format
    sparse_matrix_t       csrA;
    sparse_operation_t    transA;
    /* RHS and solution vectors. */
    double b[5], x[5], bs[5], 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. */
    /* -------------------------------------------------------------------- */
    /* .. 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 */
    
    iparm[2] = 1; //# procs
    
    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[12] = 0;
    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 = 1;           /* Print statistical information  */
    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: " IFORMAT, error);
        exit(1);
    }
    printf("\nReordering completed ... ");
    printf("\nNumber of nonzeros in factors = " IFORMAT, iparm[17]);
    printf("\nNumber of factorization MFLOPS = " IFORMAT, 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: " IFORMAT, error);
        exit(2);
    }
    printf("\nFactorization completed ... ");
    /* -------------------------------------------------------------------- */
    /* .. Back substitution and iterative refinement. */
    /* -------------------------------------------------------------------- */
    phase = 33;
    
    descrA.type = SPARSE_MATRIX_TYPE_GENERAL;
    descrA.mode = SPARSE_FILL_MODE_UPPER;
    descrA.diag = SPARSE_DIAG_NON_UNIT;
    mkl_sparse_d_create_csr(&csrA, SPARSE_INDEX_BASE_ONE, n, n, ia, ia + 1, ja, a);
    
    /* Set right hand side to one. */
    for ( i = 0; i < n; i++ )
    {
        b[i] = 1;
    }
    //  Loop over 3 solving steps: Ax=b, AHx=b and ATx=b
    for ( i = 0; i < 3; i++ )
    {
        iparm[11] = i;        /* Conjugate transposed/transpose solve */
        if ( i == 0 )
            transA = SPARSE_OPERATION_NON_TRANSPOSE;
        else if ( i == 1 )
            transA = SPARSE_OPERATION_CONJUGATE_TRANSPOSE;
        else
            transA = SPARSE_OPERATION_TRANSPOSE;
    
        printf("\n\nSolving system with iparm[11] = " IFORMAT " ...\n", iparm[11]);
        PARDISO(pt, &maxfct, &mnum, &mtype, &phase,
            &n, a, ia, ja, &idum, &nrhs, iparm, &msglvl, b, x, &error);
        if ( error != 0 )
        {
            printf("\nERROR during solution: " IFORMAT, error);
            exit(3);
        }
    
        printf("\nThe solution of the system is: ");
        for ( j = 0; j < n; j++ )
        {
            printf("\n x [" IFORMAT "] = % f", j, x[j]);
        }
        printf("\n");
        // Compute residual
        mkl_sparse_d_mv(transA, 1.0, csrA, descrA, x, 0.0, bs);
        res = 0.0;
        res0 = 0.0;
        for ( j = 1; j <= n; j++ )
        {
            res += (bs[j - 1] - b[j - 1]) * (bs[j - 1] - b[j - 1]);
            res0 += b[j - 1] * b[j - 1];
        }
        res = sqrt(res) / sqrt(res0);
        printf("\nRelative residual = %e", res);
        // Check residual
        if ( res > 1e-10 )
        {
            printf("Error: residual is too high!\n");
            exit(10 + (int)i);
        }
    }
    mkl_sparse_destroy(csrA);
    
    /* -------------------------------------------------------------------- */
    /* .. 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);
    
    int iFF = 0;

}

// Run program: Ctrl + F5 or Debug > Start Without Debugging menu
// Debug program: F5 or Debug > Start Debugging menu

// Tips for Getting Started: 
//   1. Use the Solution Explorer window to add/manage files
//   2. Use the Team Explorer window to connect to source control
//   3. Use the Output window to see build output and other messages
//   4. Use the Error List window to view errors
//   5. Go to Project > Add New Item to create new code files, or Project > Add Existing Item to add existing code files to the project
//   6. In the future, to open this project again, go to File > Open > Project and select the .sln file



int main()
{
    std::cout << "Hello World!\n";

    //test_pardiso();
    test_matrix_check();
};