matrix.cpp

Ho-Kashyap算法

//Matrix.cpp

#include "stdafx.h"
#include "Matrix.h"

//成员函数实现
Matrix::Matrix(int row,int col,double value)
{
	if (row<0 || col<0)
		throw exception("Matrix::Matrix(int,int):illegal input");
	m_row = row;												//没说的，矩阵初始化，各元素置零
	m_col = col;
	m_c.resize(row);
	for (int i=0; i<row; i++)
	{
		m_c[i].resize(col,value);
	}
}

Matrix::Matrix(const double* first,const double* last,int row,int col)
{
	if (row<1 || col<1)
		throw exception("Matrix::Matrix(double*,double*,int,int):illegal input");
	if (row*col!=(last-first))
		throw exception("Matrix::Matrix(double*,double*,int,int):arguements not match");

	const double *p=first;										//工作指针初始化
	m_row = row;												//行列初始化
	m_col = col;
	m_c.resize(row);											//向量大小初始化
	for (int i=0; i<row; i++)									//矩阵各元素初始化
	{
		m_c[i].resize(col);
		for (int j=0; j<col; j++)
		{
		 	if (fabs(*p)>=MATRIX_MIN)
				m_c[i][j] = *p++;
		    else
		    	m_c[i][j] = 0.0;
		}
	}
}

void Matrix::SetRow(int row,double value)
{
	if (row<0)
		throw exception("Matrix::SetRow():illegal input");
	m_row = row;
	m_c.resize(row,m_vector<double>(m_col,value));
}


void Matrix::SetCol(int col,double value)
{
	if (col<0)
		throw exception("Matrix::SetCol():illegal input");
	m_col = col;
	for (int i=0; i<m_row; i++)
		m_c[i].resize(col,value);
}

bool Matrix::operator == (const Matrix& s) const
{
	if (m_row!=s.m_row || m_col!=s.m_col) return false;		//行列数不等，矩阵不等
	for (int i=0; i<m_row; i++)
	    for (int j=0; j<m_col; j++)
	        if (fabs(m_c[i][j]-s.m_c[i][j])>=MATRIX_MIN)
				return false;								//有元素值不对应相等，矩阵不等
	return true;
}

Matrix Matrix::operator + (const Matrix& s) const
{
	if (m_row!=s.m_row || m_col!=s.m_col)
		throw exception("Matrix::operator + ():illegal input");
	Matrix sum(m_row,m_col);
	for (int i=0; i<m_row; i++)
		for (int j=0; j<m_col; j++)
		{
			sum.m_c[i][j] = m_c[i][j] + s.m_c[i][j];
			if (fabs(sum.m_c[i][j])<MATRIX_MIN)
			   sum.m_c[i][j] = 0.0;
		}
	return sum;
}

Matrix& Matrix::operator += (const Matrix& s)
{
	if (m_row!=s.m_row || m_col!=s.m_col)
		throw exception("Matrix::operator += ():illegal input");
	for (int i=0; i<m_row; i++)
		for (int j=0; j<m_col; j++)
		{
			m_c[i][j] += s.m_c[i][j];
			if (fabs(m_c[i][j])<MATRIX_MIN)
			    m_c[i][j] = 0.0;
		}
	return *this;
}

Matrix Matrix::operator - (const Matrix& s) const
{
	if (m_row!=s.m_row || m_col!=s.m_col)
		throw exception("Matrix::operator - ():illegal input");
	Matrix dif(m_row,m_col);
	for (int i=0; i<m_row; i++)
		for (int j=0; j<m_col; j++)
		{
			dif.m_c[i][j] = m_c[i][j] - s.m_c[i][j];
			if (fabs(dif.m_c[i][j])<MATRIX_MIN)
			   dif.m_c[i][j] = 0.0;
		}
	return dif;
}

Matrix& Matrix::operator -= (const Matrix& s)
{
	if (m_row!=s.m_row || m_col!=s.m_col)
		throw exception("Matrix::operator -= ():illegal input");
	for (int i=0; i<m_row; i++)
		for (int j=0; j<m_col; j++)
		{
			m_c[i][j] -= s.m_c[i][j];
			if (fabs(m_c[i][j])<MATRIX_MIN)
			   m_c[i][j] = 0.0;
		}
	return *this;
}

Matrix Matrix::operator * (double time) const
{
	Matrix mul=*this;
	for (int i=0; i<m_row; i++)
		for (int j=0; j<m_col; j++)
		{
			mul.m_c[i][j] = m_c[i][j]*time;
			if (fabs(mul.m_c[i][j])<MATRIX_MIN)
			   mul.m_c[i][j] = 0.0;
		}
	return mul;
}

Matrix& Matrix::operator *= (double time)
{
	for (int i=0; i<m_row; i++)
		for (int j=0; j<m_col; j++)
		{
			m_c[i][j] *= time;
			if (fabs(m_c[i][j])<MATRIX_MIN)
			   m_c[i][j] = 0.0;
		}
	return *this;
}

Matrix Matrix::operator * (const Matrix& s) const
{
	if (m_col!=s.m_row)
		throw exception("Matrix::operator * (Matrix&):illegal input");
	Matrix mul(m_row,s.m_col);
	for (int i=0; i<m_row; i++)
	{
		for (int j=0; j<s.m_col; j++)
		{
			mul.m_c[i][j] = 0.0;
			for (int k=0; k<m_col; k++)
				mul.m_c[i][j] += m_c[i][k]*s.m_c[k][j];
			if (fabs(mul.m_c[i][j])<MATRIX_MIN)
			    mul.m_c[i][j] = 0.0;
		}
	}
	return mul;
}

Matrix& Matrix::operator *= (const Matrix& s)
{
	if (m_col!=s.m_row)
		throw exception("Matrix::operator *= (Matrix&):illegal input");
	Matrix mul(m_row,s.m_col,0.0);
	for (int i=0; i<m_row; i++)
	{
		for (int j=0; j<s.m_col; j++)
		{
			for (int k=0; k<m_col; k++)
				mul.m_c[i][j] += m_c[i][k]*s.m_c[k][j];
			if (fabs(mul.m_c[i][j])<MATRIX_MIN)
			    mul.m_c[i][j] = 0.0;
		}
	}
	return *this=mul;
}

Matrix Matrix::operator | (const Matrix& B) const
{
	//if (m_row!=m_col || m_row!=B.m_row || B.m_col!=1)
	if (m_row!=B.m_row)
		throw exception("Matrix::operator | ():illegal input");
	Matrix extend(*this);
	for (int i=0; i<m_row; i++)
	{
		extend.m_c[i].resize(extend.m_col+B.m_col);
		for (int j=0; j<B.m_col; j++)
			extend.m_c[i][extend.m_col+j] = B.m_c[i][j];
	}
	extend.m_col += B.m_col;
	return extend;
}

Matrix Matrix::Transpose() const
{
	Matrix trans(m_col,m_row);
	for (int i=0; i<m_row; i++)
		for (int j=0; j<m_col; j++)
			trans.m_c[j][i] = m_c[i][j];
	return trans;
}

Matrix Matrix::Converse() const
{
	//检查矩阵，抛掷异常
	if (m_row<=0 || m_col<=0)
		throw exception("Matrix::Converse():矩阵参数非法");
	if (m_row != m_col)
		throw exception("Matrix::Converse():矩阵参数非法");
	int n=m_row;

	//用单位矩阵生成增广矩阵
	Matrix work(*this|this->Unit());

	// 高斯消元
	vector<int> msg=work.m_gause();
	
	//检查是否可逆
	if (msg.size()>0)
		throw exception("Matrix::Converse():矩阵不可逆");

	//返回结果
	Matrix Con(n,n);
	for (int i=0; i<n; i++)
		for (int j=0; j<n; j++)
			Con[i][j] = work[i][n+j];
	return Con;
}

Matrix Matrix::Adjoint() const
{
	if (m_row!=m_col)
		throw "Matrix::Adjoint():非方阵";
	double det=Determinate();
	if (fabs(det)>=MATRIX_MIN)
	{
		return det*Converse();
	}
	else
	{
		int n=m_row;
		Matrix adj(n,n);
		for (int i=0; i<n; i++)
			for (int j=0; j<n; j++)
				adj[i][j] = Cofactor(j+1,i+1);
		return adj;
	}
}

double Matrix::Cofactor(int i,int j) const
{
	i--,j--;                        	//适应数组习惯
	if (m_row!=m_col || m_row<2)
	    throw exception("Matrix::Cofactor():矩阵不满足该运算所需条件");
	if (i<0 || i>=m_row || j<0 || j>=m_col)
	    throw exception("Matrix::Cofactor():illegal input");
	Matrix work(m_row-1,m_col-1);
	for (int row=0,r=0; row<m_row; row++)
	{
		if (row==i)//Pay attention!
			continue;
		for (int col=0,c=0; col<m_col; col++)
		{
			if (col==j)
			    continue;
			work.m_c[r][c] = m_c[row][col];
			c++;
		}
		r++;
	}
	int sign=1;
	if ((i+j)&1) sign*=-1;          	//（i+j）为奇数时,代数余子式符号为-1
	/*cout<<"*this:"<<endl<<*this;
	cout<<"work with i="<<i+1<<";j="<<j+1<<";sign="<<sign<<endl<<work;
	system("pause");*/
	return sign*work.Determinate();
}

double Matrix::Determinate() const
{
 	if (m_col!=m_row)
	    throw exception("Matrix::Determinate():非方阵");
	if (m_row<1)
	    throw exception("Matrix::Determinate():矩阵不满足该运算所需条件");
    //通过转换为上三角矩阵求行列式的值
	Matrix work(*this);
	int n=m_row,sign=1;
	for (int i=0; i<n; i++)                                 //初等行变换为上三角矩阵
	{
		double maxi=fabs(work.m_c[i][i]);
		int maxIndex = i,j;