return temp;
}


// 求解聯立方程組(solve simultaneous equation)
inline CMatrix CMatrix::Solve (const CMatrix& v) const 
{
	unsigned int i,j,k;
	double a1;
	
	VERIFY (Row == Col && Col == v.Row);
	
	CMatrix temp(Row,Col+v.Col);
	for (i=0; i < Row; i++)
	{
		for (j=0; j < Col; j++)
			temp.Val[i][j] = Val[i][j];
		for (k=0; k < v.Col; k++)
			temp.Val[i][Col+k] = v.Val[i][k];
	}
	for (k=0; k < Row; k++)
	{
		int indx = temp.pivot(k);
		VERIFY (indx != -1);
		
		a1 = temp.Val[k][k];
		for (j=k; j < temp.Col; j++)
			temp.Val[k][j] /= a1;
		
		for (i=k+1; i < Row; i++)
		{
			a1 = temp.Val[i][k];
			for (j=k; j < temp.Col; j++)
				temp.Val[i][j] -= a1 * temp.Val[k][j];
		}
	}
	CMatrix s(v.Row,v.Col);
	for (k=0; k < v.Col; k++)
		for (int m=int(Row)-1; m >= 0; m--)
		{
			s.Val[m][k] = temp.Val[m][Col+k];
			for (j=m+1; j < Col; j++)
				s.Val[m][k] -= temp.Val[m][j] * s.Val[j][k];
		}
		return s;
}


// 矩陣所有元素置零，有可能改變大小(Null)
inline void CMatrix::Null (const unsigned int& row, const unsigned int& col) 
{
	if (row != Row || col != Col)
		realloc( row,col);
	
	for (unsigned int i=0; i < Row; i++)
		for (unsigned int j=0; j < Col; j++)
			Val[i][j] = double(0);
		return;
}

// 矩陣所有元素置一，有可能改變大小(Null)
inline void CMatrix::Ones (const unsigned int& row, const unsigned int& col) 
{
	if (row != Row || col != Col)
		realloc( row,col);
	
	for (unsigned int i=0; i < Row; i++)
		for (unsigned int j=0; j < Col; j++)
			Val[i][j] = double(1);
		return;
}

// 矩陣所有元素置零(Null)
inline void CMatrix::Null() 
{
	for (unsigned int i=0; i < Row; i++)
		for (unsigned int j=0; j < Col; j++)
			Val[i][j] = double(0);
		return;
}

// 矩陣所有元素置一
inline void CMatrix::Ones() 
{
	for (unsigned int i=0; i < Row; i++)
		for (unsigned int j=0; j < Col; j++)
			Val[i][j] = double(1);
		return;
}


// 置為單位矩陣(unity)
inline void CMatrix::Unit (const unsigned int& row) 
{
	if (row != Row || row != Col)
		realloc( row, row);
	
	for (unsigned int i=0; i < Row; i++)
		for (unsigned int j=0; j < Col; j++)
			Val[i][j] = i == j ? double(1) : double(0);
		return;
}

// 置為單位矩陣(unity)
inline void CMatrix::Unit () 
{
	unsigned int row = min(Row,Col);
	Row = Col = row;
	
	for (unsigned int i=0; i < Row; i++)
		for (unsigned int j=0; j < Col; j++)
			Val[i][j] = i == j ? double(1) : double(0);
		return;
}


// 主元素法換行(partial pivoting method)
inline int CMatrix::pivot (unsigned int row)
{
	int k = int(row);
	double amax,temp;
	
	amax = -1;
	for (unsigned int i=row; i < Row; i++)
		if ( (temp = fabs( Val[i][row])) > amax && temp != 0.0)
		{
			amax = temp;
			k = i;
		}
		if (Val[k][row] == double(0))
			return -1;
		if (k != int(row))
		{
			double* rowptr = Val[k];
			Val[k] = Val[row];
			Val[row] = rowptr;
			return k;
		}
		return 0;
}

// 計算行列式(determinant)
inline double CMatrix::Det () 
{
	unsigned int i,j,k;
	double piv,detVal = double(1);
	
	VERIFY (Row == Col);
	
	CMatrix temp(*this);
	
	for (k=0; k < Row; k++)
	{
		int indx = temp.pivot(k);
		if (indx == -1)
			return 0;
		if (indx != 0)
			detVal = - detVal;
		detVal = detVal * temp.Val[k][k];
		for (i=k+1; i < Row; i++)
		{
			piv = temp.Val[i][k] / temp.Val[k][k];
			for (j=k+1; j < Row; j++)
				temp.Val[i][j] -= piv * temp.Val[k][j];
		}
	}
	return detVal;
}

// 計算2范數(norm)
inline double CMatrix::Norm () 
{
	double retVal = double(0);
	
	for (unsigned int i=0; i < Row; i++)
		for (unsigned int j=0; j < Col; j++)
			retVal += Val[i][j] * Val[i][j];
		retVal = sqrt( retVal);
		
		return retVal;
}


// 計算條件數(condition number)
inline double CMatrix::Cond () 
{
	CMatrix inv(Row,Col);
	
	inv = ! (*this);
	double retVal = Norm() * inv.Norm();
	
	return retVal;
}


// 計算余子式(cofactor)
inline double CMatrix::Cofact (unsigned int row, unsigned int col) 
{
	unsigned int i,i1,j,j1;
	
	VERIFY(Row == Col);
	
	VERIFY(row <= Row && col <= Col);
	
	
	CMatrix temp (Row-1,Col-1);
	
	for (i=i1=0; i < Row; i++)
	{
		if (i == row)
			continue;
		for (j=j1=0; j < Col; j++)
		{
			if (j == col)
				continue;
			temp.Val[i1][j1] = Val[i][j];
			j1++;
		}
		i1++;
	}
	double  cof = temp.Det();
	if ((row+col)%2 == 1)
		cof = -cof;
	
	return cof;
}


// 計算伴隨矩陣(adjoin)
inline CMatrix CMatrix::Adj () 
{
	VERIFY (Row == Col);
	
	CMatrix temp(Row,Col);
	
	for (unsigned int i=0; i < Row; i++)
		for (unsigned int j=0; j < Col; j++)
			temp.Val[i][j] = Cofact(i,j);
		
		temp = ~temp;
		return temp;
}

// 是否奇異(singular)
inline bool CMatrix::IsSingular () 
{
	if (Row != Col)
		return false;
	return (Det() == double(0));
}

// 是對角陣(diagonal)
inline bool CMatrix::IsDiagonal () 
{
	if (Row != Col)
		return false;
	for (unsigned int i=0; i < Row; i++)
		for (unsigned int j=0; j < Col; j++)
			if (i != j && Val[i][j] != double(0))
				return false;
			return true;
}

// 是標量(scalar)
inline bool CMatrix::IsScalar () 
{
	if (!IsDiagonal())
		return false;
	double v = Val[0][0];
	for (unsigned int i=1; i < Row; i++)
		if (Val[i][i] != v)
			return false;
		return true;
}

// 是單位陣(unit matrix)
inline bool CMatrix::IsUnit () 
{
	if (IsScalar() && Val[0][0] == double(1))
		return true;
	return false;
}

// 是零矩陣(Null)
inline bool CMatrix::IsNull () 
{
	for (unsigned int i=0; i < Row; i++)
		for (unsigned int j=0; j < Col; j++)
			if (Val[i][j] != double(0))
				return false;
			return true;
}

// 是對稱矩陣(symmetric)
inline bool CMatrix::IsSymmetric () 
{
	if (Row != Col)
		return false;
	for (unsigned int i=0; i < Row; i++)
		for (unsigned int j=0; j < Col; j++)
			if (Val[i][j] != Val[j][i])
				return false;
			return true;
}

// 反對稱矩陣(skew-symmetric)
inline bool CMatrix::IsSkewSymmetric () 
{
	if (Row != Col)
		return false;
	for (unsigned int i=0; i < Row; i++)
		for (unsigned int j=0; j < Col; j++)
			if (Val[i][j] != -Val[j][i])
				return false;
			return true;
}

// 是上三角陣(upper triangular)
inline bool CMatrix::IsUpperTiangular () 
{
	if (Row != Col)
		return false;
	for (unsigned int i=1; i < Row; i++)
		for (unsigned int j=0; j < i-1; j++)
			if (Val[i][j] != double(0))
				return false;
			return true;
}

// 是下三角陣(lower triangular)
inline bool CMatrix::IsLowerTiangular () 
{
	if (Row != Col)
		return false;
	
	for (unsigned int j=1; j < Col; j++)
		for (unsigned int i=0; i < j-1; i++)
			if (Val[i][j] != double(0))
				return false;
			
			return true;
}
inline CMatrix Exp(const CMatrix & m)
{
	CMatrix temp(m.Row,m.Col);
	
	for (unsigned int i=0; i < m.Row; i++)
		for (unsigned int j=0; j < m.Col; j++)
		{
			temp.Val[i][j] = exp(m.Val[i][j]);
			if(temp.Val[i][j]<0.0)temp.Val[i][j]=0.0;
		}

		
		return temp;
}

//點乘
inline CMatrix CMatrix::Mul (const CMatrix &m)
{	
	VERIFY(m.Col == Col);
	VERIFY(m.Row == Row);
	
	CMatrix temp(Row,Col);
	
	for (unsigned int i=0; i < Row; i++)
		for (unsigned int j=0; j < Col; j++)
		
			temp.Val[i][j] = Val[i][j]*m.Val[i][j];
			
		return temp;
}

// 點除(A./B)
inline CMatrix CMatrix::Div (const CMatrix& m)
{	
	VERIFY(m.Col == Col);
	VERIFY(m.Row == Row);
	
	CMatrix temp(Row,Col);
	
	for (unsigned int i=0; i < Row; i++)
		for (unsigned int j=0; j < Col; j++)
		{
			VERIFY(m.Val[i][j] != 0);
			temp.Val[i][j] = Val[i][j]/m.Val[i][j];
		}
			
		return temp;
}

inline CMatrix CMatrix::Div ()
{	

	CMatrix temp(Row,Col);
	
	for (unsigned int i=0; i < Row; i++)
		for (unsigned int j=0; j < Col; j++)
		{
			VERIFY(Val[i][j] != 0);
			temp.Val[i][j] = 1.0/Val[i][j];
		}
			
		return temp;
}