matrixn.cpp

basic mathematic classes used for math programming

#include "mathclass.h"

matrixN::matrixN(const matrixN& a)
{
    m=0;
	n=0;
	on=0;
	v=NULL;
    matrixN &c = (*this);

    c.setSize( a.row(), a.column() );

	for( int i=0; i<a.row(); i++ )
    for( int j=0; j<a.column(); j++ )
        c[i][j] = a[i][j];

}
matrixN::matrixN()
{
    on = n = m = 0;
}

matrixN::matrixN( int x, int y )
{
    on = n = x;
    m = y;
    v = new vectorN[n];
	assert(v!=NULL);

    for( int i=0; i<n; i++ )
        v[i].setSize( m );
}


matrixN::matrixN( const vectorN& x, const vectorN& y) // outer product to generate the matrix N. 
{
	on=n = x.getSize();
	m = y.getSize();

	// allocate space
	 v = new vectorN[n];

	 assert(v!=NULL);

    for( int i=0; i<n; i++ )
        v[i].setSize( m );

	for(i=0;i<n;i++)
		for(int j=0;j<m;j++){
		
			v[i][j] = x[i]*y[j];
		}



}

matrixN::matrixN( int x, int y, vectorN *a )
{
    on = n = x;
    m = y;
    v = a;
}

matrixN::~matrixN()
{
    if ( on>0 && m>0 ) delete[] v;
}

void
matrixN::getValue( m_real **d ) const
{
	for( int i=0; i<n; i++ )
	for( int j=0; i<m; j++ )
		v[i][j] = d[i][j];
}

m_real
matrixN::getValue( int row, int column ) const
{
	return v[row][column];
}

void
matrixN::setValue( m_real **d )
{
	for( int i=0; i<n; i++ )
	for( int j=0; i<m; j++ )
		d[i][j] = v[i][j];
}

void
matrixN::setValue( int row, int column, m_real value )
{
	this->v[row].setValue( column, value );
}

void
matrixN::setRow( int x, const vectorN& vec )
{
	for( int i=0; i<m; i++ )
		v[x][i] = vec[i];
}

void      
matrixN::exchangeColumn(int i, int j){

	assert(i<m);
	assert(j<m);

	vectorN coli, colj;
    getColumn(i,coli);
	getColumn(j,colj);
	setColumn(i,colj);
	setColumn(j,coli);
	



}
void      
matrixN::getColumn( int j, vectorN& c){

	//vectorN c;
	assert(j<m);

	c.setExactSize(n);
	
	for(int i=0;i<n;i++){

		c[i] = v[i][j];
	

	}

}

/*
vectorN& matrixN::getColumn(int j){  // Wei Li Trouble

	assert(j<m);
	vectorN c;

	c.setExactSize(n);

	for(int i=0;i<n;i++){

		c[i] = v[i][j];
	

	}

	return c;
	
}*/

void
matrixN::setColumn( int x, const vectorN& vec )
{
	for (int i=0; i<n; i++)
		v[i][x] = vec[i];
}

void
matrixN::setSize( int x, int y )
{

	if( on<x )
	{

		if ( on>0 ) delete[] v;
		v = new vectorN[x];
	
		assert(v!=NULL);
		on=x;
		
	}

	for( int i=0; i<on; i++ )
			v[i].setSize( y );

	n = x;
	m = y;
}

void
matrixN::setExactSize( int rownum, int colnum )
{
	if( on!=rownum )
	{
		if ( on>0 ) delete[] v;
		v = new vectorN[rownum];
		assert(v!=NULL);

		on=rownum;

	}
	for( int i=0; i<rownum; i++ )
	v[i].setExactSize( colnum );
	
	n = rownum;
	m = colnum;
}

matrixN&
matrixN::assign(const matrixN& a )
{
    matrixN &c = (*this);
    c.setSize( a.row(), a.column() );

    for( int i=0; i<a.row(); i++ )
    for( int j=0; j<a.column(); j++ )
        c[i][j] = a[i][j];

    return c;
}

matrixN&
matrixN::operator+=(const matrixN& a )
{
    matrixN &c = (*this);
    c.setSize( a.row(), a.column() );

    for( int i=0; i<a.row(); i++ )
    for( int j=0; j<a.column(); j++ )
        c[i][j] += a[i][j];

    return c;
}

matrixN&
matrixN::operator-=(const matrixN& a )
{
    matrixN &c = (*this);
    c.setSize( a.row(), a.column() );

    for( int i=0; i<a.row(); i++ )
    for( int j=0; j<a.column(); j++ )
        c[i][j] -= a[i][j];

    return c;
}

matrixN&
matrixN::operator = (const matrixN& a)
{
	matrixN &c = (*this);
    c.setSize( a.row(), a.column() );


    for( int i=0; i<a.row(); i++ )
    for( int j=0; j<a.column(); j++ )
        c[i][j] = a[i][j];
    return c;

}

matrixN&
matrixN::operator*=( m_real a )
{
    matrixN &c = (*this);

    for( int i=0; i<c.row(); i++ )
    for( int j=0; j<c.column(); j++ )
        c[i][j] *= a;

    return c;
}

matrixN&
matrixN::operator/=( m_real a )
{
    matrixN &c = (*this);

    for( int i=0; i<c.row(); i++ )
    for( int j=0; j<c.column(); j++ )
        c[i][j] /= a;

    return c;
}

matrixN&
matrixN::addColumn(const vectorN& a){

	 matrixN &c = (*this);
	if(n!=a.getSize()){

		cout<<"the matrix row size doesn't match vector size."<<endl;
	
	}else{

		matrixN old;
		old = c; // add one more function.
		int newcolnum;
		newcolnum = old.column()+1; 
		c.setSize(old.row(),newcolnum);
	    for(int i=0;i<old.row();i++)
			for(int j=0;j<old.column();j++)
				c[i][j] = old[i][j];

		for( i=0;i<c.row();i++){
		
			c[i][newcolnum-1] = a[i];
		}
	}

	return c;
}

matrixN operator*( const matrixN& a, const matrixN& b){ 

    matrixN c ;
    assert( a.column()==b.row() );
    c.setSize( a.row(), b.column() );

    for( int i=0; i<a.row(); i++ )
    for( int j=0; j<b.column(); j++ )
    {
        c[i][j] = 0;
        for( int k=0; k<a.column(); k++ )
            c[i][j] += a[i][k] * b[k][j];
    }

    return c;

};

matrixN&
matrixN::mult(const matrixN& a, const matrixN& b )
{
    matrixN &c = (*this);
    assert( a.column()==b.row() );
    c.setSize( a.row(), b.column() );

    for( int i=0; i<a.row(); i++ )
    for( int j=0; j<b.column(); j++ )
    {
        c[i][j] = 0;
        for( int k=0; k<a.column(); k++ )
            c[i][j] += a[i][k] * b[k][j];
    }

    return c;
}

matrixN&
matrixN::transpose( const matrixN& a )
{
    matrixN &c = (*this);

	matrixN acopy;
	acopy = a;

    c.setSize( a.column(), a.row() );

    for( int i=0; i<a.row(); i++ )
    for( int j=0; j<a.column(); j++ )
        c[j][i] = acopy[i][j];

    return c;
}

ostream& operator<<( ostream& os, const matrixN& a )
{
    for( int i=0; i< a.row(); i++ ) os << a.v[i] << endl;
    return os;
}

istream& operator>>( istream& is, matrixN& a )
{
    for( int i=0; i< a.row(); i++ ) is >> a.v[i];
    return is;
}

void
matrixN::LUdecompose( int* index )
{
	assert( this->row() == this->column() );

	int n = this->row();
	int i, j, k, imax;
	m_real big, dum, sum, temp;

	static vectorN vv; vv.setSize( n );
	matrixN &a = (*this);

	for ( i=0; i<n; i++ )
	{
		big = 0.0f;
		for ( j=0; j<n; j++ )
			if ((temp = fabs(a[i][j])) > big)
				big = temp;

		if (big == 0.0f)
		{
			cerr << "Singular matrix in routine LUdecompose" << endl;
			assert( FALSE );
		}

		vv[i] = 1.0f / big;
	}

	for ( j=0; j<n; j++ )
	{
		for ( i=0; i<j; i++ )
		{
			sum = a[i][j];
			for ( k=0; k<i; k++ )
				sum -= a[i][k] * a[k][j];
			a[i][j] = sum;
		}

		big = 0.0;
		for ( i=j; i<n; i++ )
		{
			sum = a[i][j];
			for ( k=0; k<j; k++ )
				sum -= a[i][k] * a[k][j];
			a[i][j] = sum;
			if ((dum = vv[i] * fabs(sum)) >= big)
			{
				big = dum;
				imax = i;
			}
		}

		if ( j!=imax )
		{
			for ( k=0; k<n; k++ )
			{
				dum = a[imax][k];
				a[imax][k] = a[j][k];
				a[j][k] = dum;
			}
			vv[imax] = vv[j];
		}

		index[j] = imax;
		if (a[j][j] == 0.0f) a[j][j] = EPS;

		if ( j!=n )
		{
			dum = 1.0f / a[j][j];
			for ( i=j+1; i<n; i++ )
				a[i][j] *= dum;
		}
	}
}

void
matrixN::LUsubstitute( int* index, vectorN &b )
{
	assert( this->row() == this->column() );

	int n = this->row();
	matrixN &a = (*this);

	int i, ii = -1, ip, j;
	m_real sum;

	for ( i=0; i<n; i++ )
	{
		ip = index[i];
		sum = b[ip];
		b[ip] = b[i];

		if (ii>-1)
			for ( j=ii; j<i; j++ )
				sum -= a[i][j] * b[j];
		else
			if (sum)
				ii = i;

		b[i] = sum;
	}

	for ( i=n-1; i>=0; i-- )
	{
		sum = b[i];
		for ( j=i+1; j<n; j++ )
			sum -= a[i][j] * b[j];
		b[i] = sum / a[i][i];
	}
}

void
matrixN::LUinverse( matrixN &mat )
{
	assert( this->row() == this->column() );

	int n = this->row();

	static int* index;
	static int index_count = 0;
	if ( index_count<n )
	{
		if ( index_count>0 ) delete[] index;
		index_count = n;
		if ( index_count>0 ) index = new int[index_count];
	}

	static vectorN b;
	b.setSize( n );

	mat.setSize( n, n );

	LUdecompose( index );