matrix.cpp

一个矩阵类的源码

//---------------------------------------------------------------------------
// 16.jul.2002 13:24 GMT
//   Matrix module implementation of TVector and TMatrix classes for 
// transparent vector and matrix operation.
//   Find more information at:
//                http://www.fuzzy.ru/source/matrix.phtml
//   Contact authors:
//       Boris Izyumov    [e-mail]
//       Olivier Riff     [e-mail]
// Bug fixes:
//   09.apr.2002  Fixed initialization error in TMatrix constructor
//                [thanx to: Theera Lilitwarangkul <[e-mail]>]
//   11.may.2002  Fixed memory leak in multiplication and addition
//                Previously module ate memory because (God knows why)
//                on every operator function call matrix was created
//                with 'new' and never deleted
//   12.jul.2002  Wiped windows.h inclusion from module (was nonsense)
//   16.jul.2002      In cooperation with Olivier Riff, Projet IMAG-PRIMA
//                    [http://www-prima.imag.fr]
//                Improved performance of Determinant and Invert functions
//                by using LU-decomposition, and added two functions, 
//                LU_decomp and LU_backsub taken from "Numerical Recipes in C",
//                http://www.ulib.org/webRoot/Books/Numerical_Recipes/
//---------------------------------------------------------------------------
#ifdef __BORLANDC__
#pragma hdrstop
#endif //__BORLANDC__
//---------------------------------------------------------------------------
#include "matrix.h"

#ifndef min
#define min(a,b) ((a) < (b) ? (a) : (b))
#endif // min

#define TINY 1.0e-16

TVector::TVector ()
{
    m_nBufCount = 0;
    m_nRows = 0;
    m_pcBuf = NULL;
    m_dVals = NULL;
    m_bVertical = true;
}

TVector::TVector (const TVector& m)
{
    int i;
    m_nBufCount = 0;
    m_nRows = 0;
    m_pcBuf = NULL;
    m_dVals = NULL;
    m_bVertical = m.m_bVertical;
    SetSize (m.m_nRows);
    for (i = 0; i<m_nRows; i++)
        m_dVals[i] = m.m_dVals[i];
}

TVector::TVector (int rows, double* val)
{
    m_nBufCount = 0;
    m_nRows = 0;
    m_pcBuf = NULL;
    m_dVals = NULL;
    m_bVertical = true;
    SetSize (rows);
    if (val) for (int i = 0; i < m_nRows; i ++) m_dVals[i] = val[i];
}

TVector::~TVector ()
{
    if (m_pcBuf) delete[] m_pcBuf;
    GrowBy(-m_nRows);
}

TVector TVector::operator- ()        //      -(*this)
{
    TVector tmp(*this);
    for(int i = 0; i < tmp.m_nRows; i ++) tmp.m_dVals[i] = -tmp.m_dVals[i];
    return tmp;
}

TVector& TVector::operator= (const TVector& m)
{
    int i;
    SetSize (m.m_nRows);
    m_bVertical = m.m_bVertical;
    for (i = 0; i<m_nRows; i++)
        m_dVals[i] = m.m_dVals[i];
    return *this;
}

TVector& TVector::operator= (const double* v)
{
    for (int i = 0; i<m_nRows; i++) m_dVals[i] = v[i];
    return *this;
}

TMatrix TVector::operator* (TMatrix& m)        //      (*this)*m
{
    TMatrix tmp;
    Multiply (&m, &tmp);
    return tmp;
}

TVector TVector::operator+ (TVector& m)        //      (*this) + m
{
    TVector tmp;
    Add (&m, &tmp);
    return tmp;
}

TVector TVector::operator- (TVector& m)        //      (*this) - m
{
    TVector tmp;
    Add (&(-m), &tmp);
    return tmp;
}

void TVector::GrowBy(int inc)
{
    double* vals = m_dVals;
    int oldcnt = m_nBufCount ? (1 << m_nBufCount) : 0;
    int newcnt;
    if (inc > 0) {
        if (m_nRows+inc > oldcnt) {
            if (m_nBufCount) m_nBufCount++;
            else m_nBufCount = 5;
            while ((newcnt = (1 << m_nBufCount)) < m_nRows+inc) m_nBufCount ++;
            m_dVals = new double[newcnt];
            if (m_nRows) {
                memcpy((void*)m_dVals,(void*)vals,sizeof(double)*m_nRows);
                delete[] vals;
            }
        }
        m_nRows += inc;
    }
    else if (inc < 0) {
        if (m_nRows+inc > 0) {
            m_nRows += inc;
            return; // can only empty
        }
        m_nBufCount = 0;
        m_nRows = 0;
        if (m_dVals) delete[] m_dVals;
        m_dVals = NULL;
    }
}

void TVector::Reset ()
{
    GrowBy(-m_nRows);
}

void TVector::SetSize (int i)
{
    if (i == m_nRows) return;
    GrowBy(i-m_nRows);
}

void TVector::Transpose ()
{
    m_bVertical = !m_bVertical;
}

void TVector::MakeIdentity ()
{
    int i;
    for (i = 0; i<m_nRows; i++) m_dVals[i] = 1.0;
}

void TVector::MakeZero ()
{
    int i;
    for (i = 0; i<m_nRows; i++) m_dVals[i] = 0.0;
}

bool TVector::Multiply (TMatrix* R, TMatrix* Dest)
{
    int i, j, k;
    double val;
    TMatrix tmp;

//  1. Vector is VERTICAL, i.e. [N * 1]
    if (m_bVertical && (R->m_nRows == 1))
    {
        tmp.SetSize (m_nRows, R->m_nCols);
        for (i = 0; i<m_nRows; i++)
            for (j = 0; j<R->m_nCols; j++)
                tmp.m_dVals[i][j] = m_dVals[i]*R->m_dVals[0][j];
    }
//  2. Vector is HORIZONTAL, i.e. [1 * N] * [R * C] = [1 * C]
    else if (!m_bVertical && (m_nRows == R->m_nRows)) {
        tmp.SetSize (1, R->m_nCols);
        for (j = 0; j<R->m_nCols; j++)
        {
            val = 0.0;
            for (k = 0; k<m_nRows; k++)
                val += m_dVals[k]*R->m_dVals[k][j];
            tmp.m_dVals[0][j] = val;
        }
    }
    else
        return false;
    *Dest = tmp;
    return true;
//  this->Multiply (Some, WhereTo); WhereTo = this*Some
}

bool TVector::Add (TVector* R, TVector* Dest)
{
    if (m_nRows != R->m_nRows)
        return false;
    int i;
    TVector* tmp = new TVector (m_nRows);

    for (i = 0; i<m_nRows; i++)
        tmp->m_dVals[i] = m_dVals[i] + R->m_dVals[i];
    *Dest = *tmp;
    return true;
}

bool TVector::Vertical ()
{
    return m_bVertical;
}

int TVector::Rows()
{
    return m_nRows;
}

void TVector::Add(double val)
{
    GrowBy(1);
    m_dVals[m_nRows-1] = val;
}

void TVector::Insert(int idx, double val)
{
    if (!((idx >= 0) && (idx <= m_nRows))) return;
    GrowBy(1);
    for (int i = m_nRows-1; i > idx; i --)
        m_dVals[i] = m_dVals[i-1];
    m_dVals[idx] = val;
}

void TVector::Delete(int idx)
{
    for (int i = idx; i < m_nRows-1; i ++) m_dVals[i] = m_dVals[i+1];
    GrowBy(-1);
}

void TVector::Swap(int idx1, int idx2)
{
    double val = m_dVals[idx1];
    m_dVals[idx1] = m_dVals[idx2];
    m_dVals[idx2] = val;
}

double TVector::Max()
{
    double d;
    int i;
    if(!m_nRows) return 0.0;
    d = m_dVals[0];
    for (i = 1; i < m_nRows; i ++) if (d < m_dVals[i]) d = m_dVals[i];
    return d;
}

double TVector::Min()
{
    double d;
    int i;
    if(!m_nRows) return 0.0;
    d = m_dVals[0];
    for (i = 1; i < m_nRows; i ++) if (d > m_dVals[i]) d = m_dVals[i];
    return d;
}

double& TVector::Val(int i)
{
    if ((i < m_nRows) && (i >= 0)) return m_dVals[i];
    return m_dVals[0];
}

double& TVector::operator[] (int idx)
{
    if ((idx < m_nRows) && (idx >= 0)) return m_dVals[idx];
    return m_dVals[0];
}

const char* TVector::Packed()
{
    if (m_pcBuf) delete[] m_pcBuf;
    m_pcBuf = new char[m_nRows*24];
    char tmp[32];

    strcpy(m_pcBuf,"{");
    for (int i = 0; i<m_nRows; i++) {
        if (i) strcat(m_pcBuf,"; ");
        sprintf(tmp,"%g",m_dVals[i]);
        strcat(m_pcBuf,tmp);
    }
    strcat(m_pcBuf,"}");
    return m_pcBuf;
}
//-- TMatrix ------------------------------------------------------
TMatrix::TMatrix ()
{
    m_nRows = 0;
    m_nCols = 0;
    m_dVals = NULL;
}

TMatrix::TMatrix (const TMatrix& m)
{
    int i, j;
    m_nRows = 0;
    m_nCols = 0;
    m_dVals = NULL;
    SetSize (m.m_nRows, m.m_nCols);
    for (i = 0; i<m_nRows; i++)
        for (j = 0; j<m_nCols; j++) m_dVals[i][j] = m.m_dVals[i][j];
}

TMatrix::TMatrix (int rows)
{
    m_nRows = 0;
    m_nCols = 0;
    m_dVals = NULL;
    SetSize (rows, 1);
}

TMatrix::TMatrix (int rows, int cols)
{
    m_nRows = 0;
    m_nCols = 0;
    m_dVals = NULL;
    SetSize (rows, cols);
}

TMatrix::~TMatrix ()
{
    FreeMem ();
}

TMatrix& TMatrix::operator= (const TMatrix& m)
{
    int i, j;
    SetSize (m.m_nRows, m.m_nCols);
    for (i = 0; i<m_nRows; i++)
        for (j = 0; j<m_nCols; j++) m_dVals[i][j] = m.m_dVals[i][j];
    return *this;
}

TMatrix TMatrix::operator* (TMatrix& m)        //      (*this)*m
{
    TMatrix tmp;
    Multiply (&m, &tmp);
    return tmp;
}

TMatrix TMatrix::operator* (TVector& v)        //      (*this)*m
{
    TMatrix tmp;
    Multiply (&v, &tmp);
    return tmp;
}

TMatrix TMatrix::operator+ (TMatrix& m)        //      (*this)*m
{
    TMatrix tmp;
    Add (&m, &tmp);
    return tmp;
}

void TMatrix::FreeMem ()
{
    int i;