math/Matrixmxn.cpp

Go to the documentation of this file.

#include <malloc.h>
#include <memory.h>
#include <math.h>
#include <assert.h>
#include "VectorN.h"
#include "Matrixmxn.h"

extern int svd(int m,int n, double **a,double *q,double **u,double **v);

double** AllocateMatrix(int m, int n)
{
        double **matrix=NULL;
    matrix = (double**) malloc(m*sizeof(double*));  // initializes elements to 0
        for (int i=0; i<m; i++)
        {
                matrix[i] = (double*) calloc(n, sizeof(double));
                if (matrix[i] == NULL)
                        return NULL;
                        //Actually, there would be some memory leak here,
                    //    supposedly, all memory allocate should be free before return;
                    //    same problems in the following section.
        }

        return matrix;
}

void FreeMatrix(int m, int n, double **matrix)
{
        for(int i=0; i<m; i++)
        {
                free(matrix[i]);
        }
        free(matrix);
}

bool Matrixmxn::AllocateMatrix()
{
        assert(values == NULL);

        values = ::AllocateMatrix(m_nRows, m_nColumns);

        assert(values != NULL);

        if (values == NULL)
                return false;

        return true;
}

bool Matrixmxn::ReallocateMatrix(int m, int n)
{
        assert(values != NULL);

    if (m == m_nRows)
    {
        for (int i = 0; i < m; i++)
        {
            values[i] = (double*) realloc(values[i], n*sizeof(double));
            if (values[i] == NULL)
                                return false;
        }
    }
    else
    {
		::FreeMatrix(m_nRows, m_nColumns, values);
                values = ::AllocateMatrix(m, n);

                //values = (double**) realloc(values, m*sizeof(double*));
                //if (values == NULL)
                //      return false;
                //for(int i = 0; i < m; i++)
                //{
                //      values[i] = (double*)calloc(n, sizeof(double));
                //      if (values[i]==NULL)
                //              return false;
                //}       
    }

        m_nRows = m;
        m_nColumns = n;
        return true;
}

void Matrixmxn::FreeMatrix()
{
	::FreeMatrix(m_nRows, m_nColumns, values);
        values = NULL;
}

Matrixmxn::Matrixmxn(const int m, const int n)
{
    m_nRows = m;
    m_nColumns = n;

        values = NULL;
        AllocateMatrix();

    if (m == n)
    {
        for (int i = 0; i < m; i++)
        {
            values[i][i] = 1;
        }
    }
}

Matrixmxn::Matrixmxn(const Matrixmxn& right)
{
    m_nRows = right.m_nRows;
    m_nColumns = right.m_nColumns;

        values = NULL;
        AllocateMatrix();

        SetValues(right);

}

Matrixmxn::~Matrixmxn()
{
        FreeMatrix();     //delete values;
}


Matrixmxn Matrixmxn::Transpose()
{
    Matrixmxn returnMe(m_nColumns, m_nRows);
    for( int i = 0; i < m_nColumns; i++ )
    {
        for( int j = 0; j < m_nRows; j++ )
        {
            returnMe.values[i][j] = values[j][i] ;
        }
    }
    return returnMe ;
}

Matrixmxn Matrixmxn::Inverse()
{
        Matrixmxn inv;
        Inverse(inv);
        return inv;
}

bool Matrixmxn::Inverse(Matrixmxn &inv)
{
    int m=m_nRows;
    int n=m_nColumns;

        if (m<n)
        {
                Matrixmxn matrixT = Transpose();
                Matrixmxn matrixTInv;
                matrixT.Inverse(matrixTInv);
                inv = matrixTInv.Transpose();
                return true;
        }

    double *w = (double *)malloc(n*sizeof(double));
    double **u = ::AllocateMatrix(m, m);
    double **v = ::AllocateMatrix(n, n);
    svd(m_nRows, m_nColumns, values, w, u, v);

    Matrixmxn matrixUT(n, m);
    Matrixmxn matrixWinv(n, n);
    Matrixmxn matrixV(n, n);

    matrixV.SetValues(n, n, v);
	::FreeMatrix(n, n, v);

    //Transpose u:
    double **ut = ::AllocateMatrix(n, m);
    for( int i = 0; i < n; i++ )
    {
        for( int j = 0; j < m; j++ )
        {
            ut[i][j] = u[j][i] ;
        }
    }
    matrixUT.SetValues(n, m, ut);
 	::FreeMatrix(m, m, u);
    ::FreeMatrix(n, m, ut);

    // Because of round off, etc, need to zero w[i] values that are
    // negligible
    double wmax = 0.0;
    for (int j = 0; j < n; j++)
        if (fabs(w[j]) > wmax)
            wmax = fabs(w[j]);
    double wmin = wmax*1.0e-6;
    int rank = 0;
    for (j = 0; j < n; j++)
    {
        double winv = 0;
        if ((fabs(w[j]) < wmin)||(fabs(w[j]) < 1e-12))
        {
            w[j] = 0.0;
        }
        else
        {
            rank++;
            winv = 1/w[j];
        }
        matrixWinv.SetValues(j, j, winv);
    }
        free(w);

    // Back substitution - Taken from NR: x = V(W^-1)(U^T*b) (right to left)
        inv.Resize(n, m);
    inv.SetValues(matrixV);
    inv *= matrixWinv;
    inv *= matrixUT;

    return true;
}

double Matrixmxn::Determinant()
{
    if (m_nRows != m_nColumns)
        return 0;
    if (m_nRows == 2)
    {
        return (values[0][0]*values[1][1]-values[1][0]*values[0][1]);
    }
    double det = 0;
    for (int i = 0; i < m_nColumns; i++)
    {
        if (values[0][i] == 0)
            continue;
        Matrixmxn subMatrix(m_nRows-1, m_nColumns-1);
        for (int j = 0, col = 1; j < m_nColumns; j++)
        {
            if (j == i)
                continue;
            for (int k = 1; k < m_nRows; k++)
            {
                subMatrix.SetValues(k, col, values[k][j]);
            }
            col++;
        }
        int sign = (i%2 != 0)?  -1:1;
        det += sign*values[0][i]*subMatrix.Determinant();
    }
    return det;
}


void Matrixmxn::SetValues (const Matrixmxn& right)
{
        Resize(right.m_nRows, right.m_nColumns);
    for( int i = 0; i < m_nRows; i++ )
    {
        for( int j = 0; j < m_nColumns; j++ )
        {
            values[i][j] = right.values[i][j];
        }
    }
}

void Matrixmxn::SetValues(int m, int n, double **val)
{
        Resize(m, n);
    for( int i = 0; i < m_nRows; i++ )
    {
        for( int j = 0; j < m_nColumns; j++ )
        {
            values[i][j] = val[i][j];
        }
    }
}

bool Matrixmxn::Compare (const Matrixmxn& right)
{
    for ( int i = 0; i < m_nRows; i++ )
    {
        for ( int j = 0; j < m_nColumns; j++ )
        {
            double diff = values[i][j] - right.values[i][j];
            double tol = 1e-5;
            if ( ( diff > tol ) || ( diff < -tol ) )
            {
                return false;
            }
        }
    }

    return true;
}


Matrixmxn Matrixmxn::operator + (const Matrixmxn& right) const
{
    Matrixmxn returnMe(m_nRows, m_nColumns);
    for( int i = 0; i < m_nRows; i++ )
    {
        for( int j = 0; j < m_nColumns; j++ )
        {
            returnMe.values[i][j] = values[i][j] + right.values[i][j];
        }
    }
    return returnMe;
}

Matrixmxn Matrixmxn::operator - (const Matrixmxn& right) const
{
    Matrixmxn returnMe(m_nRows, m_nColumns);
    for( int i = 0; i < m_nRows; i++ )
    {
        for( int j = 0; j < m_nColumns; j++ )
        {
            returnMe.values[i][j] = values[i][j] - right.values[i][j];
        }
    }
    return returnMe;
}

Matrixmxn Matrixmxn::operator * (const Matrixmxn& right) const
{
    int rightCols = right.GetColumns();
    int rightRows = right.GetRows();
    assert(m_nColumns == rightRows);
    Matrixmxn returnMe(m_nRows, rightCols);
    for( int row = 0; row < m_nRows; row++ )
    {
        for( int col = 0; col < rightCols; col++ )
        {
            returnMe.values[ row ][ col ] = 0 ;
            for( int index = 0; index < m_nColumns ; index++ )
            {
                returnMe.values[ row ][ col ] += values[ row ][ index ] * right.values[ index ][ col ] ;
            }
        }
    }
    return returnMe ;
}

Matrixmxn Matrixmxn::operator * (const double& right) const
{
    Matrixmxn returnMe(m_nRows, m_nColumns);
    for( int i = 0; i < m_nRows; i++ )
    {
        for( int j = 0; j < m_nColumns; j++ )
        {
            returnMe = right*values[i][j];
        }
    }
    return returnMe;
}

VectorN Matrixmxn::operator *(const VectorN& right) const
{
    int vectorLen = right.Length();
    assert(vectorLen == m_nColumns);

    VectorN returnMe;
    returnMe.SetLength(m_nRows);
    for( int i = 0; i < m_nRows; i++ )
    {
        returnMe[i] = 0;
        for( int j = 0; j < m_nColumns; j++ )
        {
            returnMe[i] += values[i][j] * right[j];
        }
    }
    return returnMe;
}

Matrixmxn& Matrixmxn::operator = (const Matrixmxn& right)
{
    Resize(right.GetRows(), right.GetColumns());
    SetValues(right);
    return *this;
}

Matrixmxn& Matrixmxn::operator *= (const Matrixmxn& right)
{
    int rightCols = right.GetColumns();
    double **result = ::AllocateMatrix(m_nRows, rightCols);

    for( int row = 0; row < m_nRows; row++ )
    {
        for( int col = 0; col < rightCols; col++ )
        {
            result[ row ][ col ] = 0 ;
            for( int index = 0; index < m_nColumns ; index++ )
            {
                result[ row ][ col ] += values[ row ][ index ] * right.values[ index ][ col ] ;
            }
        }
    }

    SetValues(m_nRows, rightCols, result);
    ::FreeMatrix(m_nRows, rightCols, result);
    return *this;
}

Matrixmxn& Matrixmxn::operator *= (const double& right)
{
    for( int row = 0; row < m_nRows; row++ )
    {
        for( int col = 0; col < m_nColumns; col++ )
        {
            values[ row ][ col ] *= right ;
        }
    }

    return *this;
}

Matrixmxn& Matrixmxn::operator += (const Matrixmxn& right)
{
    for( int i = 0; i < m_nRows; i++ )
    {
        for( int j = 0; j < m_nColumns; j++ )
        {
            values[i][j] += right.values[i][j];
        }
    }
    return *this;
}

Matrixmxn& Matrixmxn::operator -= (const Matrixmxn& right)
{
    for( int i = 0; i < m_nRows; i++ )
    {
        for( int j = 0; j < m_nColumns; j++ )
        {
            values[i][j] -= right.values[i][j];
        }
    }
    return *this;
}

bool Matrixmxn::operator == (const Matrixmxn& right)
{
    return Compare(right);
}

double& Matrixmxn::operator ()(const unsigned int row, const unsigned int col)
{
    return values[row][col];
}

double Matrixmxn::GetValues(const int m, const int n) const
{
    return values[m][n];
}

// Set data in entry row m, column n
void Matrixmxn::SetValues(const int m, const int n, const double data)
{
    values[m][n] = data;
}

// Resize the matrix to m rows, n columns
bool Matrixmxn::Resize(const int m, const int n)
{
    if(m*n == 0)            // cannot allocate 0 size matrix
    {
        return false;
    }
    if ((m == m_nRows) && (n == m_nColumns))
    {
        return true;
    }

    if(values == NULL)
    {
        return false;
    }

        return ReallocateMatrix(m, n);
}

VectorN Matrixmxn::GetColumnVector(int col) const
{
        VectorN returnMe;
        returnMe.SetLength(m_nRows);
        for (int i=0; i<m_nRows; i++)
        {
                returnMe[i] = values[i][col];
        }
        return returnMe;
}

VectorN Matrixmxn::GetRowVector(int row) const
{
        VectorN returnMe;
        returnMe.SetLength(m_nColumns);
        for (int i=0; i<m_nColumns; i++)
        {
                returnMe[i] = values[row][i];
        }
        return returnMe;
}

---