// MKLlibTest.cpp : This file contains the 'main' function. Program execution begins and ends there.
//
//#define EIGEN_USE_MKL_ALL
//#define MKL_DIRECT_CALL
//#define MKL_NUM_THREADS 10

#include <omp.h>
#include <set>
#include <map>
#include <queue>
#include <fstream>
#include <iostream>
#include <iomanip>
#include <sstream>
#include <thread>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <Eigen/SparseCore>
#include <Eigen/Eigenvalues>
#include <unsupported/Eigen/KroneckerProduct>
#include <unsupported/Eigen/ArpackSupport>
#include <unsupported/Eigen/CXX11/Tensor>
#include "unsupported/Eigen/ArpackSupport"
#include <unsupported/Eigen/SparseExtra>
#include "mkl.h"
#include "performanceCounter.h"
using vec2i = Eigen::Vector2i;
using vec2 = Eigen::Vector2d;
using vec3 = Eigen::Vector3d;
using vec4 = Eigen::Vector4d;
using vec3f = Eigen::Vector3f;
using mat3 = Eigen::Matrix3d;
using DMatrix = Eigen::MatrixXd;
using IMatrix = Eigen::MatrixXi;
using DSparse = Eigen::SparseMatrix<double>;
using DVector = Eigen::VectorXd;
using IVector = Eigen::VectorXi;

static void EigenSparseToCuSparseTranspose(
	//const Eigen::SparseMatrix<double, Eigen::RowMajor> &mat, int *row, int *col, double *val)
	const DSparse &mat, int *row, int *col, double *val)
{
	const int num_non0 = mat.nonZeros();
	const int num_outer = mat.cols() + 1;

	memcpy(row, mat.outerIndexPtr(), sizeof(int) * num_outer);

	memcpy(col, mat.innerIndexPtr(), sizeof(int) * num_non0);

	memcpy(val, mat.valuePtr(), sizeof(double) * num_non0);

	for (int el = 0; el < num_outer; el++) {
		row[el] = row[el] + 1;
	}
	for (int el = 0; el < num_non0; el++) {
		col[el] = col[el] + 1;
	}
}

int main()
{
	// Define the sparse matrix object
	Eigen::SparseMatrix<double> Lap;
	Eigen::SparseMatrix<double> Vol;
	// Read the matrix from the text file
	Eigen::loadMarket(Lap,"D:/csource/InvMat_vegaTestMyL2L0_Temp - DiscreteL0ADMMSubspaceRevised/Lap.txt");
	Eigen::loadMarket(Vol, "D:/csource/InvMat_vegaTestMyL2L0_Temp - DiscreteL0ADMMSubspaceRevised/Vol.txt");
	PerformanceCounter overallT,substep0, substep1, substep2, substep3;
	overallT.StartCounter();
	substep0.StartCounter();
	// Print the matrix to the console
	MKL_INT  N = Lap.rows();
	int num_non0 = Lap.nonZeros();
	int num_outer = Lap.cols() + 1;
	MKL_INT* ia = new MKL_INT[num_outer];
	MKL_INT* ja = new MKL_INT[num_non0];
	double* Vala = new double[num_non0];
	EigenSparseToCuSparseTranspose(Lap, ia, ja, Vala);
	int num_non0b = Vol.nonZeros();
	int num_outerb = Vol.cols() + 1;
	MKL_INT* ib = new MKL_INT[num_outerb];
	MKL_INT* jb = new MKL_INT[num_non0b];
	double* Valb = new double[num_non0b];
	EigenSparseToCuSparseTranspose(Vol, ib, jb, Valb);
	substep0.StopCounter();
	substep1.StartCounter();
	int nr = 3;
	////exact for a 3*3 matrix, not my matrix
	double   Eig[3] = { 1.0, 9.0, 27.0 };      //  

	/* mkl_sparse_d_gv input parameters */
	char         which = 'S'; /* Which eigenvalues to calculate. ('L' - largest (algebraic) eigenvalues, 'S' - smallest (algebraic) eigenvalues) */
	MKL_INT      pm[128];     /* This array is used to pass various parameters to Extended Eigensolver Extensions routines. */
	MKL_INT      k0 = nr;     /* Desired number of max/min eigenvalues */

	/* mkl_sparse_d_gv output parameters */
	MKL_INT      k;           /* Number of eigenvalues found (might be less than k0). */
	double*		 E = new double[k0];
	//double       E[3];        /* Eigenvalues */
	double*		 X = new double[k0*N];
	//double       X[3 * 3];      /* Eigenvectors */
	double*		 res = new double[k0];
	//double       res[k0];      /* Residual */

	/* Local variables */
	MKL_INT      info;               /* Errors */
	MKL_INT      compute_vectors = 0;/* Flag to compute eigenvectors */
	MKL_INT      tol = 7;            /* Tolerance */
	MKL_INT      i;

	//Eigen::PardisoLDLT <Eigen::SparseMatrix<double> slo;
	/* Sparse BLAS IE variables */
	sparse_matrix_t A = NULL, B = NULL; /* Handle containing sparse matrix in internal data structure */
	struct matrix_descr descrA, descrB; /* Structure specifying sparse matrix properties */
	/////
	//SPARSE_MATRIX_TYPE_GENERAL = 20,   /*    General case                    */
	//	SPARSE_MATRIX_TYPE_SYMMETRIC = 21,   /*    Triangular part of              */
	//	SPARSE_MATRIX_TYPE_HERMITIAN = 22,   /*    the matrix is to be processed   */
	//	SPARSE_MATRIX_TYPE_TRIANGULAR = 23,
	//	SPARSE_MATRIX_TYPE_DIAGONAL = 24,   /* diagonal matrix; only diagonal elements will be processed */
	//	SPARSE_MATRIX_TYPE_BLOCK_TRIANGULAR = 25,
	//	SPARSE_MATRIX_TYPE_BLOCK_DIAGONAL = 26    /* block-diagonal matrix; only diagonal blocks will be processed */

	////////////

	/* Create handle for matrix A stored in CSR format */
	//descrA.type = SPARSE_MATRIX_TYPE_GENERAL; /* Full matrix is stored */
	descrA.type = SPARSE_MATRIX_TYPE_SYMMETRIC;
	descrA.mode = SPARSE_FILL_MODE_UPPER;
	descrA.diag = SPARSE_DIAG_NON_UNIT;
	//descrA.type = SPARSE_MATRIX_TYPE_DIAGONAL;
	//	struct matrix_descr {
	//	sparse_matrix_type_t  type;       /* matrix type: general, diagonal or triangular / symmetric / hermitian */
	//	sparse_fill_mode_t    mode;       /* upper or lower triangular part of the matrix ( for triangular / symmetric / hermitian case) */
	//	sparse_diag_type_t    diag;       /* unit or non-unit diagonal ( for triangular / symmetric / hermitian case) */
	//};

	//mkl_sparse_d_create_csr(&A, SPARSE_INDEX_BASE_ONE, N, N, ia, ia + 1, ja, EigenSPTest.valuePtr());
	sparse_status_t Astate = mkl_sparse_d_create_csr(&A, SPARSE_INDEX_BASE_ONE, N, N, ia, ia + 1, ja, Vala);
	//mkl_ccsrcsc();
	std::cout << "Astate:" << Astate << std::endl;
	//mkl_sparse_d_create_csr(&A, SPARSE_INDEX_BASE_ONE, N, N, ia, ia + 1, ja, a);
	//descrB.type = SPARSE_MATRIX_TYPE_GENERAL; /* Full matrix is stored */
	descrB.type = SPARSE_MATRIX_TYPE_SYMMETRIC;
	descrB.mode = SPARSE_FILL_MODE_UPPER;
	descrB.diag = SPARSE_DIAG_NON_UNIT;
	//descrB.type = SPARSE_MATRIX_TYPE_DIAGONAL;
	/*mkl_sparse_d_create_csr(&B, SPARSE_INDEX_BASE_ONE, N, N, ib, ib + 1, jb, EigenSPTestB.valuePtr());*/
	sparse_status_t Bstate = mkl_sparse_d_create_csr(&B, SPARSE_INDEX_BASE_ONE, N, N, ib, ib + 1, jb, Valb);
	std::cout << "Bstate:" << Bstate << std::endl;
	/////////////////
	//SPARSE_STATUS_SUCCESS = 0,    /* the operation was successful */
	//	SPARSE_STATUS_NOT_INITIALIZED = 1,    /* empty handle or matrix arrays */
	//	SPARSE_STATUS_ALLOC_FAILED = 2,    /* internal error: memory allocation failed */
	//	SPARSE_STATUS_INVALID_VALUE = 3,    /* invalid input value */
	//	SPARSE_STATUS_EXECUTION_FAILED = 4,    /* e.g. 0-diagonal element for triangular solver, etc. */
	//	SPARSE_STATUS_INTERNAL_ERROR = 5,    /* internal error */
	//	SPARSE_STATUS_NOT_SUPPORTED = 6     /* e.g. operation for double precision doesn't support other types */
	////////////////////////////////////
	/* Step 2. Call mkl_sparse_ee_init to define default input values */
	mkl_sparse_ee_init(pm);
	//dfeast_scsrgv();
	pm[1] = 6; /* Set tolerance */
	pm[2] = 0;
	pm[6] = 1;
	/* Step 3. Solve the standard Ax = ex eigenvalue problem. */
	for (i = 0; i < 10; i++)
	{
		std::cout << "pm_" << i << ":" << pm[i] << std::endl;
	}
	substep1.StopCounter();
	substep2.StartCounter();
	info = mkl_sparse_d_gv(&which, pm, A, descrA, B, descrB, k0, &k, E, X, res);
	substep2.StopCounter();
	substep3.StartCounter();
	printf("mkl_sparse_d_gv output info after \n");
	for (i = 0; i < 10; i++)
	{
		std::cout << "pm_" << i << ":" << pm[i] << std::endl;
	}
	//info = mkl_sparse_d_ev(&which, pm, A, descrA, k0, &k, E, X, res);
	printf("mkl_sparse_d_gv output info %d \n", info);
	if (info != 0)
	{
		printf("Routine mkl_sparse_d_gv returns code of ERROR: %i", (int)info);
	}

	printf("*************************************************\n");
	printf("************** REPORT ***************************\n");
	printf("*************************************************\n");
	printf("#mode found/subspace %d %d \n", k, k0);
	printf("Index/Exact Eigenvalues/Estimated Eigenvalues/Residuals\n");
	for (i = 0; i < k; i++)
	{
		printf("   %d  %.15e %.15e %.15e \n", i, Eig[i], E[i], res[i]);
	}
	mkl_sparse_destroy(A);
	mkl_sparse_destroy(B);
	substep3.StopCounter();
	overallT.StopCounter();
	std::cout << "substep0:" << substep0.GetElapsedTime() << std::endl;
	std::cout << "substep1:" << substep1.GetElapsedTime() << std::endl;
	std::cout << "substep2:" << substep2.GetElapsedTime() << std::endl;
	std::cout << "substep3:" << substep3.GetElapsedTime() << std::endl;
	std::cout << "overallT:" << overallT.GetElapsedTime() << std::endl;
	std::cout << "End!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!:" << std::endl;
	getchar();
}

// 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