/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * Geophysical Computational Tools & Library (GCTL)
 *
 * Copyright (c) 2023  Yi Zhang (yizhang-geo@zju.edu.cn)
 *
 * GCTL is distributed under a dual licensing scheme. You can redistribute 
 * it and/or modify it under the terms of the GNU Lesser General Public 
 * License as published by the Free Software Foundation, either version 2 
 * of the License, or (at your option) any later version. You should have 
 * received a copy of the GNU Lesser General Public License along with this 
 * program. If not, see <http://www.gnu.org/licenses/>.
 * 
 * If the terms and conditions of the LGPL v.2. would prevent you from using 
 * the GCTL, please consider the option to obtain a commercial license for a 
 * fee. These licenses are offered by the GCTL's original author. As a rule, 
 * licenses are provided "as-is", unlimited in time for a one time fee. Please 
 * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
 * to include some description of your company and the realm of its activities. 
 * Also add information on how to contact you by electronic and paper mail.
 ******************************************************/

#ifndef _GCTL_MATRIX_H
#define _GCTL_MATRIX_H

#include <typeinfo>
#include "array.h"

namespace gctl
{
	template <typename MatValType> class matrix;

	/*
	 * Here we define some commonly used array types.
	 */
	typedef matrix<int> _2i_matrix; ///< 2-D integer array
	typedef matrix<float> _2f_matrix; ///< 2-D single-precision floating-point array
	typedef matrix<double> _2d_matrix; ///< 2-D double-precision floating-point array
	typedef matrix<std::string> _2s_matrix; ///< 2-D string array
	typedef matrix<std::complex<float> > _2cf_matrix; ///< 1-D complex float array
	typedef matrix<std::complex<double> > _2cd_matrix; ///< 1-D complex double array

	/**
	 * @brief      二维数组模版
	 * 
	 * 二维数组模版是一个简单的二维动态数组实现，包含简单的数组操作与IO操作函数。
	 *
	 * @tparam     MatValType     数组元素模版类型
	 */
	template <typename MatValType>
	class matrix
	{
	protected:
		MatValType **val; ///< 二维数组数据
		size_t row_length, col_length; ///< 二维数组的行列数

	public:
		/**
		 * @brief      默认构造函数
		 */
		matrix();

		/**
		 * @brief      构造函数：初始化二维数组空间
		 *
		 * @param[in]  row_len  二维数组行数
		 * @param[in]  col_len  二维数组列数
		 */
		matrix(size_t row_len, size_t col_len);

		/**
		 * @brief      构造函数：初始化二维数组空间并赋处值
		 *
		 * @param[in]  row_len   二维数组行数
		 * @param[in]  col_len   二维数组列数
		 * @param[in]  init_val  数组元素的初始值
		 */
		matrix(size_t row_len, size_t col_len, MatValType init_val);

		/**
		 * @brief      构造函数：从一维数组初始化二维数组
		 *
		 * @param      init_mat  输入的一维数组
		 * @param[in]  row_len   二维数组行数
		 * @param[in]  col_len   二维数组列数
		 */
		matrix(MatValType *init_mat, size_t row_len, size_t col_len);

		/**
		 * @brief      构造函数：从二维数组初始化二维数组
		 *
		 * @param      init_mat  输入的二维数组
		 * @param[in]  row_len     二维数组行数
		 * @param[in]  col_len     二维数组列数
		 */
		matrix(MatValType **init_mat, size_t row_len, size_t col_len);

		/**
		 * @brief      构造函数：从一维向量初始化二维数组
		 *
		 * @param[in]  b           输入向量的引用
		 * @param[in]  row_len     二维数组行数
		 * @param[in]  col_len     二维数组列数
		 */
		matrix(const std::vector<MatValType> &b, size_t row_len, size_t col_len);

		/**
		 * @brief      构造函数：从二维向量初始化二维数组
		 *
		 * @param[in]  b     输入向量的引用
		 */
		explicit matrix(const std::vector<std::vector<MatValType> > &b);

		/**
		 * @brief      构造函数：从一维数组初始化二维数组
		 *
		 * @param[in]  b           输入数组的引用
		 * @param[in]  row_len     二维数组行数
		 * @param[in]  col_len     二维数组列数
		 */
		matrix(const array<MatValType> &b, size_t row_len, size_t col_len);

		/**
		 * @brief      拷贝构造函数
		 *
		 * @param[in]  b     拷贝数组
		 */
		matrix(const matrix<MatValType> &b);

		/**
		 * @brief      赋值操作符重载
		 *
		 * @param[in]  b     赋值数组
		 *
		 * @return     新建数组
		 */
		matrix<MatValType>& operator= (const matrix<MatValType> &b);

		/**
		 * @brief      析构函数
		 */
		virtual ~matrix();

		/**
		 * @brief      初始化二维数组空间
		 *
		 * @param[in]  row_len  二维数组行数
		 * @param[in]  col_len  二维数组列数
		 */
		void resize(size_t row_len, size_t col_len);

		/**
		 * @brief      构造函数：初始化二维数组空间并赋处值
		 *
		 * @param[in]  row_len   二维数组行数
		 * @param[in]  col_len   二维数组列数
		 * @param[in]  init_val  数组元素的初始值
		 */
		void resize(size_t row_len, size_t col_len, MatValType init_val);

		/**
		 * @brief      构造函数：从一维数组初始化二维数组
		 *
		 * @param      init_mat  输入的一维数组
		 * @param[in]  row_len   二维数组行数
		 * @param[in]  col_len   二维数组列数
		 */
		void resize(MatValType *init_mat, size_t row_len, size_t col_len);

		/**
		 * @brief      构造函数：从二维数组初始化二维数组
		 *
		 * @param      init_mat  输入的二维数组
		 * @param[in]  row_len     二维数组行数
		 * @param[in]  col_len     二维数组列数
		 */
		void resize(MatValType **init_mat, size_t row_len, size_t col_len);

		/**
		 * @brief      构造函数：从一维向量初始化二维数组
		 *
		 * @param[in]  b           输入向量的引用
		 * @param[in]  row_len     二维数组行数
		 * @param[in]  col_len     二维数组列数
		 */
		void resize(const std::vector<MatValType> &b, size_t row_len, size_t col_len);

		/**
		 * @brief      从二维向量初始化二维数组
		 *
		 * @param[in]  b     输入向量的引用
		 */
		void resize(const std::vector<std::vector<MatValType> > &b);

		/**
		 * @brief      构造函数：从一维数组初始化二维数组
		 *
		 * @param[in]  b           输入数组的引用
		 * @param[in]  row_len     二维数组行数
		 * @param[in]  col_len     二维数组列数
		 */
		void resize(const array<MatValType> &b, size_t row_len, size_t col_len);

		/**
		 * @brief      从二维数组初始化二维数组
		 *
		 * @param[in]  b     输入数组
		 */
		void resize(const matrix<MatValType> &b);

		/**
		 * @brief      清除并重置数组空间
		 */
		void clear();

		/**
		 * @brief      Initialize the array with selected random types.
		 *
		 * @param[in]  np1          Mean (Gauss) or low bound value (Even)
		 * @param[in]  np2          Standard deviation (Gauss) or hig bound value (Even).
		 * @param[in]  mode         Random types. 'RdNormal' for Gaussian distributed numbers and 
		 * 'RdUniform' for even distributed numbers.
		 */
		void random_float(MatValType np1, MatValType np2, random_type_e mode = RdNormal, unsigned int seed = 0);

		/**
		 * @brief      Initialize the array with selected random types.
		 *
		 * @param[in]  np1          Mean (Gauss) or low bound value (Even)
		 * @param[in]  np2          Standard deviation (Gauss) or hig bound value (Even).
		 * @param[in]  mode         Random types. 'RdNormal' for Gaussian distributed numbers and 
		 * 'RdUniform' for even distributed numbers.
		 */
		void random_int(MatValType np1, MatValType np2, unsigned int seed = 0);

		/**
		 * @brief      对全体元素进行赋值
		 *
		 * @param[in]  in_val  输入值
		 */
		void assign_all(MatValType in_val);

		/**
		 * @brief      对全体元素进行缩放
		 * 
		 * @param s    缩放的比例
		 */
		void scale(MatValType s);
		
		/**
		 * @brief      按行或列计算元素平均值
		 * 
		 * @param m    均值数组
		 * @param odr  计算行或列的平均值
		 */
		void mean(array<MatValType> &m, matrix_order_e odr = RowMajor);
		
		/**
		 * @brief      计算矩阵的模长
		 * 
		 * @param n_type 模的类型
		 */
		double norm(norm_type_e n_type = L2) const;

		/**
		 * @brief Set elements' value as a sequent.
		 * 
		 * @param st_val  Start value.
		 * @param row_inc Row increasement.
		 * @param col_inc Column increasement.
		 */
		void sequence(MatValType st_val, MatValType row_inc, MatValType col_inc);

		/**
		 * @brief      获取索引位置元素值
		 * 
		 * @note       可用于指针情况
		 *
		 * @param[in]  row_index  行索引值
		 * @param[in]  col_index  列索引值
		 *
		 * @return     索引处元素
		 */
		MatValType &at(size_t row_index, size_t col_index);

		/**
		 * @brief      获取索引位置元素值（常量重载）
		 * 
		 * @note       可用于指针情况
		 *
		 * @param[in]  row_index  行索引值
		 * @param[in]  col_index  列索引值
		 *
		 * @return     索引处元素
		 */
		MatValType &at(size_t row_index, size_t col_index) const;

		/**
		 * @brief      下标符号重载：获取索引位置元素值
		 * 
		 * @note       不可用于指针情况，对于二维数组只需重载一次下标符号就行。
		 *
		 * @param[in]  index  索引值
		 *
		 * @return     索引处元素
		 */
		MatValType *operator[](size_t index);

		/**
		 * @brief      下标符号重载：获取索引位置元素值（常量重载）
		 * 
		 * @note       不可用于指针情况，对于二维数组只需重载一次下标符号就行。
		 *
		 * @param[in]  index  索引值
		 *
		 * @return     索引处元素
		 */
		MatValType *operator[](size_t index) const;

		/**
		 * @brief      获取成员数组裸指针
		 *
		 * @param[in]  row_index  行索引值 row_index = 0 时数组指针为全部矩阵数据的头指针
		 *
		 * @return     单个元素的指针
		 */
		MatValType *get(size_t row_index = 0) const;

		/**
		 * @brief      判断数组是否为空
		 *
		 * @return     数组为空返回真，否则返回假
		 */
		bool empty() const;

		/**
		 * @brief      返回矩阵行大小
		 *
		 * @return     矩阵行大小
		 */
		size_t row_size() const;

		/**
		 * @brief      返回矩阵列大小
		 *
		 * @return     矩阵列大小
		 */
		size_t col_size() const;

		/**
		 * @brief      返回矩阵数组大小
		 * 
		 * @return 数组大小
		 */
		size_t val_size() const;

		/**
		 * @brief      拷贝数组到新的二维数组
		 *
		 * @param      b     拷贝对象数组
		 */
		void copy_to(matrix<MatValType> &b) const;

		/**
		 * @brief      拷贝数组到二维向量
		 *
		 * @param      b     拷贝对象向量
		 */
		void copy_to(std::vector<std::vector<MatValType> > &b) const;

		/**
		 * @brief      提取二维数组的行或者列到输出数组
		 *
		 * @param      b           输出一维数组
		 * @param[in]  index       行或者列的索引值
		 * @param[in]  order       按行或者列输出
		 */
		void extract(array<MatValType> &b, size_t index, matrix_order_e order = RowMajor) const;

		/**
		 * @brief      将二维数组转换为一维数组
		 * 
		 * @warning    此函数会重置输出的一维数组
		 *
		 * @param      b           转换后数组
		 * @param[in]  order       按行或者列优先排列输出
		 */
		void reform(array<MatValType> &b, matrix_order_e order = RowMajor) const;

		/**
		 * @brief element operate function pointer
		 * 
		 */
		typedef void (*foreach_m_ptr)(MatValType &ele_ptr, size_t r_id, size_t c_id);
		
		/**
		 * @brief    Operate on each and every element
		 * 
		 * @param func  operation function
		 */
		void for_each(foreach_m_ptr func);

		/**
		 * @brief   Display the elements.
		 * 
		 */
		void show(std::ostream &os = std::cout, char sep = '\t') const;
	};

	template <typename MatValType>
	matrix<MatValType>::matrix()
	{
		row_length = col_length = 0;
		val = nullptr;
	}

	template <typename MatValType>
	matrix<MatValType>::matrix(size_t row_len, size_t col_len)
	{
		row_length = col_length = 0;
		val = nullptr;
		resize(row_len, col_len);
	}

	template <typename MatValType>
	matrix<MatValType>::matrix(size_t row_len, size_t col_len, MatValType init_val)
	{
		row_length = col_length = 0;
		val = nullptr;
		resize(row_len, col_len, init_val);
	}

	template <typename MatValType>
	matrix<MatValType>::matrix(MatValType *init_mat, size_t row_len, size_t col_len)
	{
		row_length = col_length = 0;
		val = nullptr;
		resize(init_mat, row_len, col_len);
	}

	template <typename MatValType>
	matrix<MatValType>::matrix(MatValType **init_mat, size_t row_len, size_t col_len)
	{
		row_length = col_length = 0;
		val = nullptr;
		resize(init_mat, row_len, col_len);
	}

	template <typename MatValType>
	matrix<MatValType>::matrix(const std::vector<MatValType> &b, size_t row_len, size_t col_len)
	{
		row_length = col_length = 0;
		val = nullptr;
		resize(b, row_len, col_len);
	}

	template <typename MatValType>
	matrix<MatValType>::matrix(const std::vector<std::vector<MatValType> > &b)
	{
		row_length = col_length = 0;
		val = nullptr;
		resize(b);
	}

	template <typename MatValType>
	matrix<MatValType>::matrix(const array<MatValType> &b, size_t row_len, size_t col_len)
	{
		row_length = col_length = 0;
		val = nullptr;
		resize(b, row_len, col_len);
	}

	template <typename MatValType>
	matrix<MatValType>::matrix(const matrix<MatValType> &b)
	{
		row_length = col_length = 0;
		val = nullptr;
		resize(b);
	}

	template <typename MatValType>
	matrix<MatValType>& matrix<MatValType>::operator= (const matrix<MatValType> &b)
	{
		if (b.empty())
		{
			clear();
		}
		else
		{
			resize(b.row_size(), b.col_size());
			for (size_t i = 0; i < row_length; i++)
			{
				for (size_t j = 0; j < col_length; j++)
				{
					val[i][j] = b[i][j];
				}
			}
		}
		return *this;
	}

	template <typename MatValType>
	matrix<MatValType>::~matrix()
	{
		clear();
	}

	template <typename MatValType>
	void matrix<MatValType>::resize(size_t row_len, size_t col_len)
	{
		if (row_len == 0 || col_len == 0)
		{
			throw std::invalid_argument("[gctl::matrix<T>::resize] Invalid matrix size.");
		}

		// nothing needs to be done
		if (row_len == row_length && col_len == col_length)
		{
			return;
		}

		// re-indexing
		if (row_len*col_len == row_length*col_length)
		{
			row_length = row_len;
			col_length = col_len;

			// construct new row-indexing
			MatValType *t_ptr = val[0];
			delete[] val;
			val = new MatValType* [row_len];
			val[0] = t_ptr;

			for (size_t i = 1; i < row_len; i++)
			{
				val[i] = val[i-1] + col_len;
			}
			return;
		}

		// re-allocate
		if (!empty()) clear();

		row_length = row_len;
		col_length = col_len;

		val = new MatValType* [row_len];
		val[0] = new MatValType [row_len*col_len];

		for (size_t i = 1; i < row_len; i++)
		{
			val[i] = val[i-1] + col_len;
		}
		return;
	}

	template <typename MatValType>
	void matrix<MatValType>::resize(size_t row_len, size_t col_len, MatValType init_val)
	{
		resize(row_len, col_len);

		for (size_t i = 0; i < row_len; i++)
		{
			for (size_t j = 0; j < col_len; j++)
			{
				val[i][j] = init_val;
			}
		}
		return;
	}

	template <typename MatValType>
	void matrix<MatValType>::resize(MatValType *init_mat, size_t row_len, size_t col_len)
	{
		if (init_mat == nullptr)
		{
			throw std::domain_error("[gctl::matrix<T>::resize] Invalid pointer.");
		}

		resize(row_len, col_len);

		for (size_t i = 0; i < row_len; i++)
		{
			for (size_t j = 0; j < col_len; j++)
			{
				val[i][j] = init_mat[j + i*col_len];
			}		
		}
		return;
	}

	template <typename MatValType>
	void matrix<MatValType>::resize(MatValType **init_mat, size_t row_len, size_t col_len)
	{
		if (init_mat == nullptr)
		{
			throw std::domain_error("[gctl::matrix<T>::resize] Invalid pointer.");
		}

		resize(row_len, col_len);

		for (size_t i = 0; i < row_len; i++)
		{
			for (size_t j = 0; j < col_len; j++)
			{
				val[i][j] = init_mat[i][j];
			}
		}
		return;
	}

	template <typename MatValType>
	void matrix<MatValType>::resize(const std::vector<MatValType> &b, size_t row_len, size_t col_len)
	{
		if (!b.empty())
		{
			resize(row_len, col_len);

			for (size_t i = 0; i < row_len; i++)
			{
				for (size_t j = 0; j < col_len; j++)
				{
					val[i][j] = b[j + i*col_len];
				}
			}
			return;
		}
		return;
	}

	template <typename MatValType>
	void matrix<MatValType>::resize(const std::vector<std::vector<MatValType> > &b)
	{
		if (!b.empty())
		{
			size_t max_col = 0;
			for (size_t i = 0; i < b.size(); i++)
			{
				max_col = std::max(max_col, b[i].size());
			}

			resize(b.size(), max_col, 0);
			for (size_t i = 0; i < row_length; i++)
			{
				for (size_t j = 0; j < b[i].size(); j++)
				{
					val[i][j] = b[i][j];
				}
			}
			return;
		}
		return;
	}

	template <typename MatValType>
	void matrix<MatValType>::resize(const array<MatValType> &b, size_t row_len, size_t col_len)
	{
		if (!b.empty())
		{
			if (b.size() != row_len*col_len)
			{
				throw std::runtime_error("[gctl::matrix<T>::resize] Invalid matrix parameters.");
			}
			
			resize(row_len, col_len);

			for (size_t i = 0; i < row_len; i++)
			{
				for (size_t j = 0; j < col_len; j++)
				{
					val[i][j] = b[j + i*col_len];
				}
			}
			return;
		}
		return;
	}

	template <typename MatValType>
	void matrix<MatValType>::resize(const matrix<MatValType> &b)
	{
		if (!b.empty())
		{
			resize(b.row_size(), b.col_size());

			for (size_t i = 0; i < row_length; i++)
			{
				for (size_t j = 0; j < col_length; j++)
				{
					val[i][j] = b[i][j];
				}
			}
			return;
		}
		return;
	}

	template <typename MatValType>
	void matrix<MatValType>::clear()
	{
		row_length = col_length = 0;
		
		if (val != nullptr)
		{
			delete[] val[0];
			delete[] val;
			val = nullptr;
		}
		return;
	}

	template <typename MatValType>
	void matrix<MatValType>::random_float(MatValType np1, MatValType np2, random_type_e mode, unsigned int seed)
	{
		static_assert(std::is_floating_point<MatValType>::value, 
			"gctl::matrix<T>::random_float(...) could only be used with a float point type.");

		if (seed == 0) seed = std::chrono::system_clock::now().time_since_epoch().count();
		std::default_random_engine generator(seed);

		if (mode == RdNormal)
		{
			std::normal_distribution<MatValType> dist(np1, np2); // 添加高斯分布的随机值

			for (size_t i = 0; i < row_length; i++)
			{
				for (size_t j = 0; j < col_length; j++)
				{
					val[i][j] = dist(generator);
				}
			}
			return;
		}

		std::uniform_real_distribution<MatValType> dist(np1, np2); // /添加均匀分布的随机值

		for (size_t i = 0; i < row_length; i++)
		{
			for (size_t j = 0; j < col_length; j++)
			{
				val[i][j] = dist(generator);
			}
		}
		return;
	}

	template <typename MatValType>
	void matrix<MatValType>::random_int(MatValType np1, MatValType np2, unsigned int seed)
	{
		static_assert(std::is_integral<MatValType>::value, 
			"gctl::matrix<T>::random_int(...) could only be used with an integral type.");

		if (seed == 0) seed = std::chrono::system_clock::now().time_since_epoch().count();
		std::default_random_engine generator(seed);

		//添加均匀分布的随机值
		std::uniform_int_distribution<MatValType> dist(np1, np2);
		for (size_t i = 0; i < row_length; i++)
		{
			for (size_t j = 0; j < col_length; j++)
			{
				val[i][j] = dist(generator);
			}
		}
		return;
	}

	template <typename MatValType>
	void matrix<MatValType>::assign_all(MatValType in_val)
	{
		for (size_t i = 0; i < row_length; i++)
		{
			for (size_t j = 0; j < col_length; j++)
			{
				val[i][j] = in_val;
			}
		}
		return;
	}

	template <typename MatValType>
	void matrix<MatValType>::scale(MatValType s)
	{
		for (size_t i = 0; i < row_length; i++)
		{
			for (size_t j = 0; j < col_length; j++)
			{
				val[i][j] *= s;
			}
		}
		return;
	}

	template <typename MatValType>
	void matrix<MatValType>::mean(array<MatValType> &m, matrix_order_e odr)
	{
		if (odr == RowMajor)
		{
			m.resize(col_length);
			for (size_t j = 0; j < col_length; j++)
			{
				m[j] = 0.0;
				for (size_t i = 0; i < row_length; i++)
				{
					m[j] += val[i][j];
				}
				m[j] /= row_length;
			}
		}
		else
		{
			m.resize(row_length);
			for (size_t i = 0; i < row_length; i++)
			{
				m[i] = 0.0;
				for (size_t j = 0; j < col_length; j++)
				{
					m[i] += val[i][j];
				}
				m[i] /= col_length;
			}
		}
		return;
	}

	template <typename MatValType>
	double matrix<MatValType>::norm(norm_type_e n_type) const
	{
		MatValType nval = (MatValType) 0;
		if (n_type == L0)
        {
			for (size_t i = 0; i < row_length; i++)
			{
				for (size_t j = 0; j < col_length; j++)
				{
					if (val[i][j] != (MatValType) 0)
					{
						nval += (MatValType) 1;
					}
				}
			}
            return nval;
        }

        if (n_type == L1)
        {
			for (size_t i = 0; i < row_length; i++)
			{
				for (size_t j = 0; j < col_length; j++)
				{
					nval += GCTL_FABS(val[i][j]);
				}
			}
            return nval;
        }

		if (n_type == L2)
		{
			for (size_t i = 0; i < row_length; i++)
			{
				for (size_t j = 0; j < col_length; j++)
				{
					nval += val[i][j]*val[i][j];
				}
			}
			return sqrt(nval);
		}

		nval = val[0][0];
		for (size_t i = 0; i < row_length; i++)
		{
			for (size_t j = 0; j < col_length; j++)
			{
				nval = GCTL_MAX(nval, val[i][j]);
			}
		}
        return sqrt(nval);
	}

	template <typename MatValType>
	void matrix<MatValType>::sequence(MatValType st_val, MatValType row_inc, MatValType col_inc)
	{
		for (size_t i = 0; i < row_length; i++)
		{
			for (size_t j = 0; j < col_length; j++)
			{
				val[i][j] = st_val + i*row_inc + j*col_inc;
			}
		}
		return;
	}

	template <typename MatValType>
	MatValType *matrix<MatValType>::get(size_t row_index) const
	{
#ifdef GCTL_BOUND_CHECK
		if (row_index >= row_length || col_index >= col_length)
			throw std::out_of_range("[gctl::matrix<T>::get] Invalid index.");
#endif // GCTL_BOUND_CHECK

		return val[row_index];
	}

	template <typename MatValType>
	MatValType &matrix<MatValType>::at(size_t row_index, size_t col_index)
	{
#ifdef GCTL_BOUND_CHECK
		if (row_index >= row_length || col_index >= col_length)
			throw std::out_of_range("[gctl::matrix<T>::at] Invalid index.");
#endif // GCTL_BOUND_CHECK

		return val[row_index][col_index];
	}

	template <typename MatValType>
	MatValType &matrix<MatValType>::at(size_t row_index, size_t col_index) const
	{
#ifdef GCTL_BOUND_CHECK
		if (row_index >= row_length || col_index >= col_length)
			throw std::out_of_range("[gctl::matrix<T>::at] Invalid index.");
#endif // GCTL_BOUND_CHECK

		return val[row_index][col_index];
	}

	template <typename MatValType>
	MatValType *matrix<MatValType>::operator[](size_t index)
	{
		return val[index];
	}

	template <typename MatValType>
	MatValType *matrix<MatValType>::operator[](size_t index) const
	{
		return val[index];
	}

	template <typename MatValType>
	bool matrix<MatValType>::empty() const
	{
		if (row_length == 0 && col_length == 0) return true;
		else return false;
	}

	template <typename MatValType>
	size_t matrix<MatValType>::row_size() const
	{
		return row_length;
	}

	template <typename MatValType>
	size_t matrix<MatValType>::col_size() const
	{
		return col_length;
	}

	template <typename MatValType>
	size_t matrix<MatValType>::val_size() const
	{
		return row_length*col_length;
	}

	template <typename MatValType>
	void matrix<MatValType>::copy_to(matrix<MatValType> &b) const
	{
		b.resize(row_length, col_length);

		for (size_t i = 0; i < row_length; i++)
		{
			for (size_t j = 0; j < col_length; j++)
			{
				b[i][j] = val[i][j];
			}
		}
		return;
	}

	template <typename MatValType>
	void matrix<MatValType>::copy_to(std::vector<std::vector<MatValType> > &b) const
	{
		if (!b.empty())
		{
			for (size_t i = 0; i < b.size(); i++)
			{
				b[i].clear();
			}
			b.clear();
		}

		b.resize(row_length);
		for (size_t i = 0; i < row_length; i++)
		{
			b[i].resize(col_length);
		}

		for (size_t i = 0; i < row_length; i++)
		{
			for (size_t j = 0; j < col_length; j++)
			{
				b[i][j] = val[i][j];
			}
		}
		return;
	}

	template <typename MatValType>
	void matrix<MatValType>::extract(array<MatValType> &b, size_t index, matrix_order_e order) const
	{
		if (order == RowMajor)
		{
			if (index >= row_length)
			{
				throw std::out_of_range("[gctl::matrix<T>::extract] Invalid index.");
			}

			b.resize(col_length);
			for (size_t i = 0; i < col_length; i++)
			{
				b[i] = val[index][i];
			}
			return;
		}

		if (index >= col_length)
		{
			throw std::out_of_range("[gctl::matrix<T>::extract] Invalid index.");
		}

		b.resize(row_length);
		for (size_t i = 0; i < row_length; i++)
		{
			b[i] = val[i][index];
		}
		return;
	}

	template <typename MatValType>
	void matrix<MatValType>::reform(array<MatValType> &b, matrix_order_e order) const
	{
		b.resize(row_length*col_length);

		if (order == RowMajor)
		{
			for (size_t i = 0; i < row_length; i++)
			{
				for (size_t j = 0; j < col_length; j++)
				{
					b[j+i*col_length] = val[i][j];
				}
			}
			return;
		}

		for (size_t j = 0; j < col_length; j++)
		{
			for (size_t i = 0; i < row_length; i++)
			{
				b[i+j*row_length] = val[i][j];
			}
		}
		return;
	}

	template <typename MatValType>
	void matrix<MatValType>::for_each(foreach_m_ptr func)
	{
		for (size_t i = 0; i < row_length; i++)
		{
			for (size_t j = 0; j < col_length; j++)
			{
				func(val[i][j], i, j);	
			}
		}
		return;
	}

	template <typename MatValType>
	void matrix<MatValType>::show(std::ostream &os, char sep) const
	{
		for (size_t i = 0; i < row_length; i++)
		{
			os << val[i][0];
			for (size_t j = 1; j < col_length; j++)
			{
				os << sep << val[i][j];
			}
			os << std::endl;	
		}
		return;
	}
};

#endif // _GCTL_MATRIX_H