// Mat.cpp: implementation of the Mat class.
// 矩陣運算 2006年1月由本人編制,如有任何問題,請聯系xiaofc4395@sina.com
// 所有涉及矩陣運算的首先必須初始化，即要調用 zero函數
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "RanMethod.h"
#include "mat.h"
#include "math.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

mat::mat()
{
	row=1;
	col=1;
	p=NULL;
}

mat::~mat()
{
	Freethis();
}

void mat::Freethis()
{
	if (p)
	{
		delete[] p;
		p = NULL;
	}
}

mat::mat(const mat& other)
{
	p=NULL;
	row = other.row;
	col = other.col;
	zero(row,col);
	memcpy(p, other.p, sizeof(double)*row*col);
}

mat::mat(int m,int n,double a[])
{
	Freethis();
	row=int(fabs(m));	col=int(fabs(n));
	int nSize = row*col;
	ASSERT(nSize!=0);
	p = new double[nSize];
	ASSERT(p != NULL);				        	  // 內存分配失敗
	BOOL bad=true;
	if (IsBadReadPtr(p, sizeof(double) * nSize)) bad=false;
	ASSERT(bad);
	for (int i=0;i<nSize;i++)
	{
		p[i]=a[i];
	}
}

void mat::zero(int m, int n)
{
	Freethis();
	row=int(fabs(m));	col=int(fabs(n));
	int nSize = row*col;
	ASSERT(nSize!=0);	
	p = new double[nSize];	
	ASSERT(p != NULL);
	BOOL bad=true;
	if (IsBadReadPtr(p, sizeof(double) * nSize)) bad=false;
	ASSERT(bad);
	memset(p,0, sizeof(double) * nSize);   // 將各元素值置0
}

//------------顯示某個元素-------------
double mat::r(int nRow, int nCol) const
{
	if (nCol <=0 || nCol > col || nRow <= 0 || nRow > row)
		return 0.0;
	if (p != NULL)
		return p[ (nCol-1) +(nRow-1)* col];
	else
		return 0.0;
}

//------------設置某個元素-------------
void mat::set(int nRow, int nCol, double value)
{
	BOOL a1=true;
	if (nCol <= 0 || nCol > col || nRow <= 0 || nRow > row) 
		a1=false;
	ASSERT(a1);
	if (p == NULL) 
		ASSERT(false);	
	p[(nCol-1) +(nRow-1)* col] = value;
}

//------------設置某個元素-------------
void mat::set(int nRow, int nCol, int value1)
{
	double value=double(value1);
	set(nRow, nCol,value);
}

//------------顯示某一行-------------
mat mat::cpr(int m)  const
{
	BOOL a1=true;
	if( (m <= 0) || (m > row) ) a1=false;
	ASSERT(a1);
	if (p == NULL) ASSERT(false);
	mat out;out.zero(1,col);
	for (int i=1;i<=col;i++)
	{
		out.set(1,i,this->r(m,i));
	}
	return out;
}

//------------顯示某一列-------------
mat mat::cpc(int n) const
{
	BOOL a1=true;
	if( (n <= 0) && ( n> col) )  a1=false;
	ASSERT(a1);
	if (p == NULL) 
		ASSERT(false);
	mat out;out.zero(row,1);
	for (int i=1;i<=row;i++)
	{
		out.set(i,1,this->r(i,n));
	}
	return out;
}

//------從第m列開始，顯示從m-n列之間的數據---
mat mat::cpc(int m,int n) const
{
	ASSERT( (m>0 ) && ( m <=n) );
	ASSERT( (n>=m) && ( n <=col) );
	mat out;out.zero(row,n-m+1);
	for (int i=1;i<=(n-m+1);i++)
	{
		mat temp=cpc(i+m-1);
		out.setc(i,temp);
	}
	return out;
}

//--從k列開始顯示一個m*n的矩陣
mat mat::cpc(int k,int m,int n) const
{
	ASSERT( (k>0 ) && ( k<=col));
	ASSERT( (m<=row ) && ((k+n-1)<=col));
	mat out;out.zero(m,n);
	for(int i=1;i<=m;i++)
	{
		for(int j=1;j<=n;j++)
		{
			out.set(i,j,r(i,j+k-1));
		}
	}
	return out;
}

//--從k行開始顯示一個m*n的矩陣
mat mat::cpr(int k,int m,int n) const
{
	ASSERT( (k>0 ) && ( k<=row));
	ASSERT( (n>0 ) && ( n<=col));
	ASSERT( (k<=row ) && ((k+m-1)<=row));
	mat out;out.zero(m,n);
	for(int i=1;i<=m;i++)
	{
		for(int j=1;j<=n;j++)
		{
			out.set(i,j,r(i+k-1,j));
		}
	}
	return out;
}

//--從k行,j列開始顯示一個m*n的矩陣
mat mat::cprc(int k,int j, int m,int n) const
{
	ASSERT( (k>0 ) && ( k<=row));
	ASSERT( (m>0 ) && ( m<=row));
	ASSERT( (j>0 ) && ( j<=col));
	ASSERT( (n>0 ) && ( n<=col));
	ASSERT( (k<=row ) && ((k+m-1)<=row));
	ASSERT( (j<=row ) && ((j+n-1)<=col));
	mat out;out.zero(m,n);
	for(int i=1;i<=m;i++)
	{
		for(int rrr=1;rrr<=n;rrr++)
		{
			out.set(i,rrr,r(i+k-1,rrr+j-1));
		}
	}
	return out;
}

//------------設置某行元素-------------
void mat::setr(int nRow,const mat& b)
{
	BOOL a1=true;
	if (nRow <= 0 || nRow > row) a1=false;
	ASSERT(a1);
	if (p == NULL) ASSERT(false);
	for (int i=1;i<=col;i++)
	{
		set(nRow,i,b.r(1,i));
	}
}

//------------設置某列元素-------------
void mat::setc(int m,const mat& b)
{
	BOOL a1=true;
	if ((m <=0 || m > col)) a1=false;
	ASSERT(a1);
	if (p == NULL) ASSERT(false);
	for (int i=1;i<=row;i++)
	{
		set(i,m,b.r(i,1));
	}
}

//----從第k列開始，設置一個m*n的矩陣
void mat::setm(int k,const mat& b)
{
	int m=b.row;
	int n=b.col;
	ASSERT( k>0);
	ASSERT( (k+n-1)<=col);
	mat pre,after,out;
	if (k==1)
	{
		after=cpc(n+1,col);
		out=b/=after;
	}
	else
	{
		if (k+n<=col)
		{
			pre=cpc(1,k-1);
			after=cpc(k+n,col);
			out=pre/=b/=after;
		}
		else
		{
			pre=cpc(1,k-1);
			out=pre/=b;
		}
	}
	zero(out.row,out.col);
	memcpy(p, out.p, sizeof(double)*row*col);
}

//-------------轉置-------------------
mat mat::T() const
{
	mat out; out.zero(col,row);
	for (int i=1;i<=col;i++)
	{
		for (int j=1;j<=row;j++)
		{
			out.set(i,j,this->r(j,i));			
		}
	}
	return out;
}

//-------------單位陣----------------
void mat::eyemat(int n)
{
	zero(n,n);
	for (int i=1;i<=n;i++)
	{ 
		for (int j=1;j<=n;j++)
		{
			if (i==j)
				set(i,j,1.0);				
			else
				set(i,j,0.0);
		}
	}
}

//-------------付值----------------
mat& mat::operator=(const mat& other)
{
	if(other.p !=NULL)
	{
		if (&other != this)
		{
			zero(other.row,other.col);
			memcpy(p, other.p, sizeof(double)*row*col);
		}
	}
	return *this;
}

//-------------乘一個實數----------------
mat mat::operator*(double value) const
{
	mat	result(*this) ;	
	for (int i = 1 ; i <= row ; ++i)
	{
		for (int j = 1 ; j <= col; ++j)
			 result.set(i,j,result.r(i,j)*value);
	}
	return result ;
}

//-------------乘一個整數----------------
mat mat::operator*(int value) const
{
	mat	result(*this) ;	
	for (int i = 1 ; i <= row ; i++)
	{
		for (int j = 1 ; j <= col; j++)
			 result.set(i,j,result.r(i,j)*value);
	}
	return result ;
}

//-------------矩陣乘----------------
mat	mat::operator*(const mat& other) const
{
	// 首先檢查行列數是否符合要求
	ASSERT (col == other.row);
	mat result;result.zero(row,other.col);
	double	value ;
	for (int i = 1 ; i <=row ; i++)
	{
		for (int j = 1 ; j <= other.col ; j++)
		{
			value = 0.0 ;
			for (int k = 1 ; k <= col ; k++)
			{
				value += this->r(i, k) * other.r(k, j);
			}
			result.set(i,j,value) ;
		}
	}
	return result ;
}

//-------------矩陣+----------------
mat	mat::operator+(const mat& other) const
{
	ASSERT (row == other.row && col == other.col); // 首先檢查行列數是否相等
	mat	result;result.zero(row,col);
	for (int i = 1; i <= row ; i++)
	{
		for (int j = 1; j <= col; j++)
		{
			result.set(i,j,r(i,j)+other.r(i,j));
		}
	}
	return result ;
}

//-------------矩陣-  ----------------
mat	mat::operator-(const mat& other) const
{
	ASSERT (row == other.row && col == other.col); // 首先檢查行列數是否相等
	mat	result;result.zero(row,col);
	for (int i = 1 ; i <= row ; ++i)
	{
		for (int j = 1 ; j <= col; ++j)
			result.set(i,j,r(i,j)-other.r(i,j));
	}
	return result;
}

//-------------矩陣/----------------
mat mat::operator/(double a) const
{
	mat	result;result.zero(row,col);	
	for (int i = 1 ; i <= row ; ++i)
	{
		for (int j = 1 ; j <= col; ++j)
			 result.set(i,j,this->r(i,j)/a);
	}
	return result ;
}

//-------------矩陣/----------------
mat mat::operator/(int a) const
{
	mat	result;result.zero(row,col);
	for (int i = 1 ; i <= row ; ++i)
	{
		for (int j = 1 ; j <= col; ++j)
			 result.set(i,j,this->r(i,j)/a);
	}
	return result ;
}

//----------------矢量除-----------------
mat mat::operator/(const mat& other) const
{
	mat out;
	int m=other.row;
	int n=other.col;
	BOOL a1=true;
	if (m !=n) 	
	{
		a1=false;
		ASSERT(a1);
	}
	else
	{
		mat thismat;thismat.zero(row,col);
		for (int i=1;i<=row;i++)
		{
			for (int j=1;j<=col;j++)
			{
				thismat.set(i,j,r(i,j));
			}
		}
		mat inva=other.invmat();
		out=thismat*inva;
	}
	return out;
}

//-----------------------矢量差乘------------------
mat mat::operator^(const mat& c) const
{
	mat out;	
	int m=c.row;
	int n=c.col;
	BOOL a1;
	if ( (m == row ) && ( n == col)  && (m*n == 3) )
	{
		a1=true;
		if (m =3)
		{
			out.zero(3,1);
			out.set(1,1,r(2,1)*c.r(3,1)-r(3,1)*c.r(2,1));
			out.set(2,1,r(3,1)*c.r(1,1)-r(1,1)*c.r(3,1));
			out.set(3,1,r(1,1)*c.r(2,1)-r(2,1)*c.r(1,1));
		}
		else
		{
			out.zero(1,3);
			out.set(1,1,r(1,2)*c.r(1,3)-r(1,3)*c.r(1,2));
			out.set(2,1,r(1,3)*c.r(1,1)-r(1,1)*c.r(1,3));
			out.set(3,1,r(1,1)*c.r(1,2)-r(1,2)*c.r(1,1));
		}
	}
	else
	{
		a1=false;
		ASSERT(a1);
	}
	return out;
}

//---------------矢量點乘--------------------
double mat::operator|(const mat& a) const
{
	//無論是矩陣還是向量，只對第一列有效/列
	double out=0.0;
	int m=a.row;
	int n=a.col;
	if (( m==row ) && (n ==col) )
	{
		double temp=0.0;
		if (m==1)
		{
			for (int i=1;i<=n;i++)
			{
				temp+=a.r(1,i)*r(1,i);
			}
			out=temp;
		}
		else
		{
			if (n==1)
			{
				for (int i=1;i<=m;i++)
				{
					temp+=a.r(i,1)*r(i,1);
				}
				out=temp;
			}	
		}
	}
	return out;
}

//--------------- 2范數 -----------------------
double mat::normmat() const
{
	double out=0.0;
	int m=row;
	int n=col;
	BOOL a1=true;
	if ((m !=1) && (n !=1)) a1=false; 
	ASSERT(a1);
	if (m==1)
	{
		double temp=0.0;
		for (int i=1;i<=n;i++)
		{
			temp+=r(1,i)*r(1,i);
		}
		out=sqrt(temp);
	}
	else
	{
		if (n==1)
		{
			double temp=0.0;
			for (int i=1;i<=m;i++)
			{
				temp+=r(i,1)*r(i,1);
			}
			out=sqrt(temp);
		}
	}
	return out;
}


//---------------伴隨陣-------------------------
mat mat::adj() const
{
	mat out;out.zero(row,col);
	for(int i=1;i<=row;i++)
	{
		for(int j=1;j<=col;j++)
		{
			if (((i+j)%2)!=0)
			{
				out.set(i,j,(-1.0)*(sub_mat(i,j).detmat()));
			}
			else
			{
				out.set(i,j,       (sub_mat(i,j).detmat()));
			}
		}
	}
	mat reuslt=out.T();
	return reuslt;
}

//---------------行列式的值-----------------------
double mat::detmat() const
{
	double out=0.0;
	BOOL a1=true;
	if ( row !=col) 
		a1=false; 
	ASSERT(a1);

	if ( (col==1) &&(row==1) )
	{
		out=r(1,1);
	}
	else
	{
		if ((col==2)&&(row==2))
		{
			out=r(1,1)*r(2,2)-r(1,2)*r(2,1);
		}
		else
		{
			int i=1;
			for(int j=1;j<=col;j++)
			{
				if ( ((i+j)%2 )!=0)
				{
					out+=(-1.0)*sub_mat(i,j).detmat()*r(i,j);
				}
				else
				{
					out+=       sub_mat(i,j).detmat()*r(i,j);
				}
			}
		}
	}
	return out;
}

//--------------------余子式------------------------
mat mat::sub_mat(int rr,int cc) const
{
	mat out;out.zero(row-1,col-1);
	int m=out.row;
	int n=out.col;
	for (int i=1;i<=m;i++)
	{
		for (int j=1;j<=n;j++)
		{
			if ( i< rr)
			{
				if ( j< cc)
				{
					out.set(i,j,r(i,j));
				}
				else
				{
					out.set(i,j,r(i,j+1));
				}
			}
			else
			{
				if (i>=rr)
				{
					if ( j< cc)
					{
						out.set(i,j,r(i+1,j));
					}
					else
					{
						out.set(i,j,r(i+1,j+1));
					}
				}
			}
		}
	}
	return out;
}

//--------------------逆矩陣--------------------------------
mat mat::invmat() const
{
	mat out;out.zero(row,col);

	BOOL a1=true;
	if (col!=row)
	{
		a1=false; 
		ASSERT(a1);
	}
	else
	{
		double d=this->detmat();
		if (d !=0.0)
		{
			double d1=1/d;
			mat rusult=adj();
			int m=rusult.row;
			int n=rusult.col;
			for (int i=1;i<=m;i++)
			{
				for (int j=1;j<=n;j++)
				{
					out.set(i,j,rusult.r(i,j)*d1);			
				}
			}
		}
	}
	return out;
}

//------2個線性空間的結合(return this)----------------
void mat::operator||(mat a)
{	
	int m1=a.row;	int n1=a.col;
	int m=row;      int n=col;
	BOOL a1=true;
	if (m1 != m)
	{
		a1=false;
		ASSERT(a1);
	}
	else
	{
		mat temp;temp.zero(m,n);
		for (int i=1;i<=m;i++)
		{
			for (int j=1;j<=n;j++)
			{
				temp.set(i,j,r(i,j));
			}
		}
		zero(m1,n+n1);		
		for (int ii=1;ii<=row;ii++)
		{
			for (int jj=1;jj<=col;jj++)
			{
				if (jj<=n)
					this->set(ii,jj,temp.r(ii,jj));
				else
					this->set(ii,jj,a.r(ii,jj-n));
			}
		}
	} 
}

//------2個線性空間的結合(return ohers)----------------
mat mat::operator|=(const mat& a) const
{
	mat out;out.zero(row,col);
	memcpy(out.p, p, sizeof(double)*row*col);
	out||a;
	return out;
}

//------2個線性行與列均不相同空間的結合(return ohers)----------------
mat mat::operator /=(const mat& a) const
{
	int m1=a.row,n=a.col;
	int m;
	if (m1> row)
		m=m1;
	else
		m=row;
	mat out;out.zero(m,col+n);
	for (int i=1;i<=out.row;i++)
	{
		for (int j=1;j<=out.col;j++)
		{
			if (j<=col)
				out.set(i,j,r(i,j));
			else
				out.set(i,j,a.r(i,j-col));
		}
	}
	return out;
}

//------行向量與一個實數的結合----------------
void mat::operator||(double a)
{
	int m=row;      int n=col;
	ASSERT(m < 1.5);
	mat temp;temp.zero(1,n);
	for (int j=1;j<=n;j++)
	{
		temp.set(1,j,r(1,j));
	}
	zero(m,n+1);
	for (int jj=1;jj<=col;jj++)
	{
		if (jj<=n)
			this->set(1,jj,temp.r(1,jj));
		else
			this->set(1,jj,a);
	} 
}


int mat::length()
{
	return col*row;
}


mat mat::Delc(int n) const
{
	ASSERT( 1<=n  && n<=col);
	mat out;
	if(col==1)
	{
		out.zero(row,1);
	}
	else
	{
		out.zero(row,col-1);
		for (int i=1;i<=row;i++)
		{
			for (int j=1;j<=col-1;j++)
			{
				if (j < n)
					out.set(i,j,r(i,j));
				else
					out.set(i,j,r(i,j+1));
			}
		}		
	}
	return out;
}

mat mat::Delr(int n) const
{
	mat out;out.zero(row-1,col);
	if ( n <= row)
	{
		for (int j=1;j<=col;j++)
		{
			for (int i=1;i<=row-1;i++)
			{
				if (i < n)
					out.set(i,j,r(i,j));
				else
					out.set(i,j,r(i+1,j));
			}
		}		
	}
	return out;
}

void mat::setmat(int m, int n, double a[])
{
	Freethis();
	row=m;	col=n;
	int nSize = row*col;
	ASSERT (nSize > 0);
	p = new double[nSize];	
	if (p == NULL)
		exit(0);				        	// 內存分配失敗
	if (IsBadReadPtr(p, sizeof(double) * nSize))
		exit(0);
	for (int i=0;i<nSize;i++)
	{
		p[i]=a[i];
	}
}