matrix.cpp

斯坦福Energy211/CME211课《c++编程——地球科学科学家和工程师》的课件

// ENERGY211/CME211
// 
// matrix.cpp - Implementation file for Project 2
//

#include "matrix.h"	// Must include corresponding header file
#include "vector.h"

#include <iostream>	// Need this to use cin, cout, etc.
#include <cmath>

using namespace std;

// Text of error messages, used by Matrix::ReportError
// Do not change this!
char *Matrix::ErrorMessages[] = {
	"invalid size specification", // ERROR_INVALID_SIZE
	"dimension mismatch", // ERROR_SIZE_MISMATCH
	"index out of range", // ERROR_INVALID_INDEX
	"invalid range specification", // ERROR_INVALID_RANGE
	"reshape to inconsistent size", // ERROR_INVALID_RESHAPE
	"invalid assignment to submatrix" // ERROR_INVALID_ASSIGNMENT
} ;

////////////////////////////////////////////////////////////

// Output operator for Range class
// Do not change this either!
ostream& operator<<( ostream& out, const Range& r )
{
	out << "(" << r.get_row1() << ":" << r.get_row2() << "," << r.get_col1() << ":" << r.get_col2() << ")\n";
	return out;
}

////////////////////////////////////////////////////////////
// Implementation of Matrix class - public operations
//
// NOW you can start changing things!
////////////////////////////////////////////////////////////

Matrix::Matrix( int rows, int cols, double *data )
{
	m_issubmatrix = false;
	m_data = NULL;
	m_rowdata = NULL;
	Set( rows, cols, data );
}

Matrix::Matrix( const Matrix& A )
{
	m_issubmatrix = false;
	m_data = NULL;
	m_rowdata = NULL;
	Set(A);
}

Matrix::~Matrix()
{
	if ( m_issubmatrix == false )
		deallocate();
}
	
void Matrix::Zeros( int rows, int cols )
{
	if ( rows <= 0 || cols < 0 )
		ReportError( ERROR_INVALID_SIZE );
	if ( cols == 0 )
		cols = rows;
	Set( rows, cols );
}

void Matrix::Ones( int rows, int cols )
{
	if ( rows <= 0 || cols < 0 )
		ReportError( ERROR_INVALID_SIZE );
	if ( cols == 0 )
		cols = rows;
	Set( rows, cols );
	for ( int i = 0; i < rows; i++ )
		for ( int j = 0; j < cols; j++ )
			m_rowdata[i][j] = 1.0;
}

void Matrix::Identity( int rows, int cols )
{
	if ( rows <= 0 || cols < 0 )
		ReportError( ERROR_INVALID_SIZE );
	if ( cols == 0 )
		cols = rows;
	Set( rows, cols );
	int minrc = ( rows < cols ? rows : cols );
	for ( int i = 0; i < minrc; i++ )
		m_rowdata[i][i] = 1.0;
}

Matrix Matrix::operator+( const Matrix& B ) const
{
	if ( m_rows != B.get_rows() || m_cols != B.get_cols() )
		ReportError( ERROR_SIZE_MISMATCH );
	Matrix C( m_rows, m_cols );
	for ( int i = 0; i < m_rows; i++ )
		for ( int j = 0; j < m_cols; j++ )
			C[i][j] = m_rowdata[i][j] + B[i][j];
	return C;
}

Matrix Matrix::operator-( const Matrix& B ) const
{
	if ( m_rows != B.get_rows() || m_cols != B.get_cols() )
		ReportError( ERROR_SIZE_MISMATCH );
	Matrix C( m_rows, m_cols );
	for ( int i = 0; i < m_rows; i++ )
		for ( int j = 0; j < m_cols; j++ )
			C[i][j] = m_rowdata[i][j] - B[i][j];
	return C;
}

Matrix Matrix::operator*( const Matrix& B ) const
{
	if ( m_cols != B.get_rows() )
		ReportError( ERROR_SIZE_MISMATCH );
	Matrix C( m_rows, B.m_cols );
	for ( int i = 0; i < C.m_rows; i++ )
		for ( int j = 0; j < C.m_cols; j++ )
			for ( int k = 0; k < m_cols; k++ )
				C[i][j] += m_rowdata[i][k] * B[k][j];
	return C;
}

Matrix Matrix::operator*( double s ) const
{
	Matrix B( m_rows, m_cols );
	for ( int i = 0; i < m_rows; i++ )
		for ( int j = 0; j < m_cols; j++ )
			B[i][j] = m_rowdata[i][j] * s;
	return B;
}

Matrix& Matrix::operator=( const Matrix& A )
{
	Set(A);
	return *this;
}

Matrix& Matrix::operator+=( const Matrix& A )
{
	if ( m_rows != A.get_rows() || m_cols != A.get_cols() )
		ReportError( ERROR_SIZE_MISMATCH );
	for ( int i = 0; i < m_rows; i++ )
		for ( int j = 0; j < m_cols; j++ )
			m_rowdata[i][j] += A[i][j];
	return *this;
}

Matrix& Matrix::operator-=( const Matrix& A )
{
	if ( m_rows != A.get_rows() || m_cols != A.get_cols() )
		ReportError( ERROR_SIZE_MISMATCH );
	for ( int i = 0; i < m_rows; i++ )
		for ( int j = 0; j < m_cols; j++ )
			m_rowdata[i][j] -= A[i][j];
	return *this;
}

Matrix& Matrix::operator*=( double s )
{
	for ( int i = 0; i < m_rows; i++ )
		for ( int j = 0; j < m_cols; j++ )
			m_rowdata[i][j] *= s;
	return *this;
}

Matrix& Matrix::operator/=( double s )
{
	for ( int i = 0; i < m_rows; i++ )
		for ( int j = 0; j < m_cols; j++ )
			m_rowdata[i][j] /= s;
	return *this;
}

double *Matrix::operator[]( int i ) const
{
	if ( i < 0 || i >= m_rows )
		ReportError( ERROR_INVALID_INDEX );
	if ( m_rowdata == NULL )
		return NULL;
	return m_rowdata[i];
}

Matrix Matrix::operator[]( Range r ) const
{
	if ( r.get_row2() == -1 )
		r.set_row2(m_rows - 1);
	if ( r.get_col2() == -1 )
		r.set_col2(m_cols - 1);
	
	if ( r.get_row1() < 0 || r.get_row1() >= m_rows )
		ReportError( ERROR_INVALID_INDEX );
	if ( r.get_row2() < 0 || r.get_row2() >= m_rows )
		ReportError( ERROR_INVALID_INDEX );
	if ( r.get_col1() < 0 || r.get_col1() >= m_cols )
		ReportError( ERROR_INVALID_INDEX );
	if ( r.get_col2() < 0 || r.get_col2() >= m_cols )
		ReportError( ERROR_INVALID_INDEX );
	
	if ( r.get_row1() > r.get_row2() )
		ReportError( ERROR_INVALID_RANGE );
	if ( r.get_col1() > r.get_col2() )
		ReportError( ERROR_INVALID_RANGE );
	
	Matrix S;
	S.m_rows = r.get_row2() - r.get_row1() + 1;
	S.m_cols = r.get_col2() - r.get_col1() + 1;
	S.m_data = m_data;
	S.m_rowdata = new double *[S.m_rows];
	for ( int i = 0; i < S.m_rows; i++ )
		S.m_rowdata[i] = (double *) (m_data + ( i + r.get_row1() ) * m_cols + r.get_col1() );
	S.m_issubmatrix = true;
	return S;
}

void Matrix::Reshape( const Matrix& A, int rows, int cols )
{
	if ( rows * cols != A.get_rows() * A.get_cols() )
		ReportError( ERROR_INVALID_RESHAPE );
	Set( rows, cols );
	for ( int i = 0; i < rows; i++ )
	{
		for ( int j = 0; j < cols; j++ )
		{
			int index = j * rows + i;
			int jA = index / A.m_rows;
			int iA = index % A.m_rows;
			m_rowdata[i][j] = A[iA][jA];
		}
	}
}

////////////////////////////////////////////////////////////
// Matrix class - private member functions
////////////////////////////////////////////////////////////

void Matrix::allocate()
{
	deallocate();
	int n_elts = m_rows * m_cols;
	if ( n_elts > 0 )
	{
		m_data = new double[n_elts];
		for ( int i = 0; i < n_elts; i++ )
			m_data[i] = 0;
		m_rowdata = new double *[m_rows];
		for ( int i = 0; i < m_rows; i++ )
			m_rowdata[i] = m_data + (i * m_cols);
	}
}

void Matrix::deallocate()
{
	delete [] m_data;
	delete [] m_rowdata;
	m_data = NULL;
	m_rowdata = NULL;
}

void Matrix::Set( int rows, int cols, double *data )
{
	m_rows = rows;
	m_cols = cols;
	allocate();
	if ( data != NULL )
	{
		for ( int i = 0; i < m_rows; i++ )
		{
			for ( int j = 0; j < m_cols; j++ )
			{
				int index = i * m_cols + j;
				m_rowdata[i][j] = data[index];
			}
		}
	}
}

void Matrix::Set( const Matrix& A )
{
	if ( m_issubmatrix == false )
		Set( A.get_rows(), A.get_cols() );
	else if ( m_rows != A.get_rows() || m_cols != A.get_cols() )
		ReportError( ERROR_INVALID_ASSIGNMENT );
	for ( int i = 0; i < m_rows; i++ )
		for ( int j = 0; j < m_cols; j++ )
			m_rowdata[i][j] = A[i][j];
}

void Matrix::ReportError( ErrorCode code )
{
	string prefix( "Error: " );
	throw std::runtime_error( prefix + ErrorMessages[code] );
}

////////////////////////////////////////////////////////////
// Non-class operator overloads for Matrix objects
////////////////////////////////////////////////////////////

ostream& operator<<( ostream& out, const Matrix& A )
{
	for ( int i = 0; i < A.get_rows(); i++ )
	{
		for ( int j = 0; j < A.get_cols(); j++ )
		{
			out.width(4);
			out << A[i][j];
		}
		out << endl;
	}
	return out;
}

istream& operator>>( istream& in, Matrix& A )
{
	int rows;
	in >> rows;
	if ( in == false )
		return in;
	int cols;
	in >> cols;
	if ( in == false )
		return in;
	
	A.Zeros( rows, cols );
	
	for ( int i = 0; i < rows; i++ )
	{
		for ( int j = 0; j < cols; j++ )
		{
			double value;
			in >> value;
			if ( in == false )
				break;
			A[i][j] = value;
		}
	}
	return in;
}

Matrix operator*( double s, const Matrix& A )
{
	return A * s;
}

Matrix Matrix::TriU( int diag ) const {
	Matrix U( m_rows, m_cols );
	for ( int i = 0; i < m_rows; i++ )
		for ( int j = 0; j < m_cols; j++ )
			if ( i <= j - diag )
				U[i][j] = m_rowdata[i][j];
	return U;
}

Matrix Matrix::TriL( int diag ) const {
	Matrix L( m_rows, m_cols );
	for ( int i = 0; i < m_rows; i++ )
		for ( int j = 0; j < m_cols; j++ )
			if ( i >= j - diag )
				L[i][j] = m_rowdata[i][j];
	return L;
}

ColVector Matrix::GetDiag( int diag ) const {
	int len1;
	int len2;
	if ( diag >= 0 )
	{
		len1 = m_rows;
		len2 = m_cols - diag;
	}
	else
	{
		len1 = m_rows + diag;
		len2 = m_cols;
	}
	int len = ( len1 < len2 ? len1 : len2 );
	ColVector v( len );
	for ( int i = 0; i < len; i++ )
		if ( diag >= 0 )
			v[i] = m_rowdata[i][i + diag];
		else
			v[i] = m_rowdata[i + diag][i];
	return v;
}

void Matrix::SetDiag( const ColVector& d, int diag ) {
	int rows, cols;
	if ( diag >= 0 )
		rows = d.length() + diag;
	else if ( diag < 0 )
		rows = d.length() - diag;
	cols = rows;
	if ( m_rows * m_cols == 0 )
		Set( rows, cols );
	else if ( m_rows != rows || m_cols != cols )
		ReportError( ERROR_SIZE_MISMATCH );
	for ( int i = 0; i < d.length(); i++ )
	{
		if ( diag >= 0 )
			m_rowdata[i][i+diag] = d[i];
		else
			m_rowdata[i-diag][i] = d[i];	
	}
}

double Matrix::norm1( const Matrix& A ) {
	double maxnorm = 0.0;
	for ( int j = 0; j < A.get_cols(); j++ ) {
		double normj = 0.0;
		for ( int i = 0; i < A.get_rows(); i++ )
			normj += abs( A[i][j] );
		if ( normj > maxnorm )
			maxnorm = normj;
	}
	return maxnorm;
}