#include <iostream>
#include <fstream>
#include <complex>

#ifndef LAPACK_COMPLEX_CUSTOM
#define LAPACK_COMPLEX_CUSTOM
#define lapack_complex_float std::complex<float>
#define lapack_complex_double std::complex<double>
#include "lapacke.h"
#endif

namespace tira {

	template <typename type>
	class Matrix {
		type* m_ptr;
		int m_rows, m_cols;

		void allocate(int rows, int cols) {
			if (m_ptr != NULL)
				free(m_ptr);
			m_rows = rows;
			m_cols = cols;
			m_ptr = (type*)malloc(m_rows * m_cols * sizeof(type));
		}

	public:
		/// Default constructor, creates an empty matrix (generally filled as a parameter)
		Matrix() {
			m_rows = m_cols = 0;
			m_ptr = NULL;
		}

		/// Contructor, creates a rows x cols matrix with undefined values
		Matrix(int rows, int cols) : Matrix() {
			allocate(rows, cols);
		}

		/// De-allocate data and clear the matrix for garbage collection
		void del() {
			m_rows = 0;
			m_cols = 0;
			free(m_ptr);
			m_ptr = NULL;
		}

		void create(int rows, int cols) {
			allocate(rows, cols);
		}

		int rows() { return m_rows; }
		int cols() { return m_cols; }
		type* get_ptr() { return m_ptr; }

		void force_square() {
			if (m_rows != m_cols)
				throw "Matrix must be square!";
			if (m_rows == 0 || m_cols == 0)
				throw "Matrix undefined!";
		}

		inline void set(int row, int col, type val) {
			m_ptr[col * m_rows + row] = val;
		}

		inline type get(int row, int col) {
			return m_ptr[col * m_rows + row];
		}

		/// Return the number of bytes in the matrix array
		size_t bytes() { return m_rows * m_cols * sizeof(type); }

		/// Create a deep copy of the current matrix and return it
		Matrix copy() {
			Matrix result(m_rows, m_cols);
			memcpy(result.m_ptr, m_ptr, bytes());
			return result;
		}

		/// Calculate the determinant of the matrix
		type det() {
			force_square();										//throw an exception if the matrix isn't square
			Matrix tmp = copy();								//copy the current matrix to create space for LUP decomposition
			Matrix<int> idx = getrf(tmp);						//perform LUP decomposition

			// multiply all elements along the diagonal
			type determinant = get(0, 0);						//initialize the determinant
			for (int i = 1; i < m_rows; i++)					//for each remaining element of the diagonal
				determinant *= tmp.get(i, i);					//calculate the running factor
			tmp.del();											//delete memory for the temporary matrix

			// calculate the number of row swaps
			int swaps = 0;										//initialize the number of row swaps
			for (int i = 0; i < m_rows; i++)					//iterate through each element of the swap vector
				if (idx.get(i, 0) != i + 1)						//if the element does not equal the row number
					swaps++;									//a row was swapped, so increment the swap counter
			if (swaps % 2 != 0)									//if there was an odd number of row swaps
				determinant *= (-1);							//multiply the determinant by -1
			return determinant;									//return the determinant
		}

		/// Calculate the trace of the matrix
		type trace() {
			force_square();										//throw an exception if the matrix isn't square
			type trace = 0;										//initialize the trace summation
			for (int i = 0; i < m_rows; i++)					//for each element along the diagonal
				trace += get(i, i);								//sum the trace
			return trace;										//return the matrix trace result
		}

		/// Output the matrix to a CSV file
		void csv(std::string filename) {
			std::ofstream outfile(filename.c_str());			//open the file for writing
			for (int r = 0; r < m_rows; r++) {					//for each row
				for (int c = 0; c < m_cols; c++) {				//for each column
					outfile << get(r, c);						//output the value to a file
					if (c != m_cols - 1) outfile << ",";		//if this isn't the last column, add a comma
				}
				if (r != m_rows - 1) outfile << std::endl;		//if this isn't the last row, end the line
			}
		}

		/* Static functions used to generate specialized matrices */

		/// Static function for generating an identity matrix
		static Matrix I(int dim) {
			Matrix result(dim, dim);							//create and allocate space for the matrix
			memset(result.m_ptr, 0, dim * dim * sizeof(type));	//set all elements to zero
			for (int i = 0; i < dim; i++)						//for each element on the diagonal
				result.set(i, i, (type)1);						//set that element to 1
			return result;										//return the completed matrix
		}

		/// Static function for generating a geometric matrix for robustness testing
		static Matrix geometric_matrix(int dim) {
			Matrix result(dim, dim);							//create and allocate space for the matrix
			for (int r = 0; r < dim; r++)						//for each row
				for (int c = 0; c < dim; c++)					//for each column
					result.set(r, c, pow((type)(2 + r), c));	//set the geometric coefficient
			return result;										//return the completed matrix
		}

		/// Static function analytically generates a solution (right hand side) vector for a geometric matrix
		static Matrix geometric_vector(int dim) {
			Matrix result(dim, 1);								//generate and allocate space for the vector
			for (int r = 0; r < dim; r++)
				result.set(r, 0, (pow((type)(2 + r), dim) - 1) / (type)(r + 1));
			return result;										//return the completed vector
		}

		/* Static functions used for LAPACK and BLAS */

		/// Calls the appropriate LAPACK GETRF function for performing LUP decomposition
		static Matrix<int> getrf(Matrix<float> M) {
			Matrix<int> idx(M.m_rows, 1);
			LAPACKE_sgetrf(LAPACK_COL_MAJOR, (int)M.m_rows, (int)M.m_cols, M.m_ptr, (int)M.m_rows, idx.get_ptr());
			return idx;
		}
		static Matrix<int> getrf(Matrix<double> M) {
			Matrix<int> idx(M.m_rows, 1);
			LAPACKE_dgetrf(LAPACK_COL_MAJOR, (int)M.m_rows, (int)M.m_cols, M.m_ptr, (int)M.m_rows, idx.get_ptr());
			return idx;
		}
	};
}