matrix.cpp

一个矩阵类的源码

    if (!m_dVals) return;
    for (i = 0; i<m_nRows; i++)
        delete[] m_dVals[i];
    delete[] m_dVals;
    m_dVals = NULL;
}

void TMatrix::GetMem ()
{
    int i;
    if (m_dVals) FreeMem ();
    if ((m_nRows>0) && (m_nCols>0)) {
        m_dVals = new double*[m_nRows];
        for (i = 0; i<m_nRows; i++) {
            m_dVals[i] = new double[m_nCols];
            memset(m_dVals[i],0,sizeof(double)*m_nCols);
        }
    }
    else {
        m_nRows = m_nCols = 0;
        m_dVals = NULL;
    }
}

void TMatrix::SetSize (int i, int j)
{
    if ((i == m_nRows) && (j == m_nCols)) return;
    FreeMem ();
    m_nRows = i;
    m_nCols = j;
    GetMem ();
}

void TMatrix::Reset ()
{
    SetSize (0, 0);
}

void TMatrix::Transpose ()
{
    int i, j;
    TMatrix tmp (*this);
    if (m_nCols != m_nRows) SetSize (m_nCols, m_nRows);
    for (i = 0; i<m_nRows; i++)
        for (j = 0; j<m_nCols; j++)
            m_dVals[i][j] = tmp.m_dVals[j][i];
}

void TMatrix::MakeIdentity ()
{
    int i, j;
    if (m_nRows != m_nCols)
        SetSize (m_nRows, m_nRows);
    for (i = 0; i<m_nRows; i++)
        for (j = 0; j<m_nCols; j++)
            m_dVals[i][j] = ((i == j) ? 1.0 : 0.0);
}

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

bool TMatrix::Multiply (TMatrix* R, TMatrix* Dest)
{
    if (m_nCols != R->m_nRows)
        return false;
    int i, j, k;
    double val;
    TMatrix tmp(m_nRows,R->m_nCols);

    for (i = 0; i<m_nRows; i++)
    {
        for (j = 0; j<R->m_nCols; j++)
        {
            val = 0.0;
            for (k = 0; k<m_nCols; k++)
                val += m_dVals[i][k]*R->m_dVals[k][j];
            tmp.m_dVals[i][j] = val;
        }
    }
    *Dest = tmp;
    return true;
//  this->Multiply (Some, WhereTo); WhereTo = this*Some
}

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

bool TMatrix::Add (TMatrix* R, TMatrix* Dest)
{
    if ((m_nCols != R->m_nCols) || (m_nRows != R->m_nRows))
        return false;
    int i, j;
    TMatrix tmp(m_nRows,m_nCols);
/*  Dest->FreeMem ();
    Dest->m_nRows = m_nRows;
    Dest->m_nCols = R->m_nCols;
    Dest->GetMem ();*/

    for (i = 0; i<m_nRows; i++)
        for (j = 0; j<m_nCols; j++)
            tmp.m_dVals[i][j] = m_dVals[i][j] + R->m_dVals[i][j];
    *Dest = tmp;
    return true;
//  this->Multiply (Some, WhereTo); WhereTo = this*Some
}

bool TMatrix::Square ()
{
    if (m_nRows == m_nCols)
        return true;
    return false;
}

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

int TMatrix::Cols ()
{
    return m_nCols;
}

int TMatrix::AddLine (TVector& vctr)
{
    double** tmp;
    int i;
    if (m_dVals)
    {
        tmp = new double*[m_nRows+1];
        for (i = 0; i<m_nRows; i++)
            tmp[i] = m_dVals[i];
        tmp[i] = new double[m_nCols];
        delete[] m_dVals;
        m_dVals = tmp;
        m_nRows ++;
    }
    else
        SetSize (1, vctr.m_nRows);
    for (i = 0; i<m_nCols; i++)
        m_dVals[m_nRows-1][i] = vctr[i];
    return m_nRows-1;
}

void TMatrix::GrowBy (int rows, int cols)
{
    double** ppd;
    int i, j, c = m_nCols, r = m_nRows;

    ppd = m_dVals;
    m_dVals = NULL;
    SetSize (m_nRows + rows, m_nCols + cols);
    for (i = 0; i<min(r,m_nRows); i++)
        for (j = 0; j<min(c,m_nCols); j++)
            m_dVals[i][j] = ppd[i][j];
    for (i = 0; i<r; i++)
        delete[] ppd[i];
    delete[] ppd;
}

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

double* TMatrix::operator[] (int idx)
{
    int i;
    if ((idx >= m_nRows) || (idx < 0))
        i = 0;
    else
        i = idx;
    return m_dVals[i];
}

double TMatrix::Determinant()
{
    if(!Square()) return 0.0;
    TMatrix tmp(*this);
    int* indx = new int[m_nRows];
    double d = 0.0,f,sum = 0.0;
    if(LU_decomp(tmp,indx,d)) {
        for(int i = 0; i < m_nRows; i ++) {
            if((f = fabs(tmp.m_dVals[i][i])) < TINY) {
                d = 0.0;
                break;
            }
            sum += log(f);
            if(f != tmp.m_dVals[i][i]) d *= -1.0;
        }
    }
    delete[] indx;
    if(d == 0.0) return 0.0;
    return d*exp(sum);
}

double TMatrix::AlgAdd (int row, int col)
{
    char* rows = new char[m_nRows];
    char* cols = new char[m_nCols];
    memset (rows, 1, m_nRows*sizeof(char));
    memset (cols, 1, m_nCols*sizeof(char));
    rows[row] = cols[col] = 0;
    double d = (((row+col)%2) ? -1.0 : 1.0)*Minor (rows, cols);
    delete[] rows;
    delete[] cols;
    return d;
}

double TMatrix::Minor (char* rowsused, char* colsused)
{
    int i, j;
    int cj = 0;
    double sum = 0.0;
    for (i = 0; i<m_nRows; i++)
        if (rowsused[i])
            break;
    if (i == m_nRows)
        return 1.0;
    for (j = 0; j<m_nCols; j++)
    {
        if (!colsused[j]) continue;
        if (m_dVals[i][j] != 0.0)
        {
            rowsused[i] = colsused[j] = 0;
            sum += ((cj%2) ? -1.0 : 1.0)*m_dVals[i][j]*Minor (rowsused, colsused);
            rowsused[i] = colsused[j] = 1;
        }
        cj ++;
    }
    return sum;
}

bool TMatrix::Invert()
{
    if(!Square()) return false;
    TMatrix tmp(*this);
    double* col = new double[m_nRows];
    int* indx = new int[m_nRows];
    int i, j;
    double d = 0.0;
    if(!LU_decomp(tmp,indx,d)) return false; // Decompose the matrix just once.
    for(j = 0; j < m_nCols; j++) { // Find inverse by columns. 
        for(i = 0; i < m_nRows; i++) col[i] = (i == j) ? 1.0 : 0.0;
        LU_backsub(tmp,indx,col);
        for(i = 0; i < m_nRows; i++) m_dVals[i][j] = col[i];
    }
    delete[] col;
    delete[] indx;
    return true;
}

bool TMatrix::LU_decomp(TMatrix& m,int* indx,double& d)
{
    int n = m.m_nRows;
    double* vv = new double[n]; //vv stores the implicit scaling of each row. 
    double big,dum,sum,temp;
    double** a = m.m_dVals;
    int imax = 0,i,j,k;
  
    d = 1.0; // No row interchanges yet. 
    for(i = 0; i < n; i++) { // Loop over rows to get the implicit scaling information.
        big = 0.0; 
        for(j = 0; j < n; j++)
            if((temp = fabs(a[i][j])) > big) big = temp; 
        if(big == 0.0) // No nonzero largest element. 
            return false;
        vv[i] = 1.0/big; // Save the scaling. 
    } 
    for (j = 0; j < n; j++) { // This is the loop over columns of Crout's method. 
        for(i = 0; i < j; i++) { // This is equation (2.3.12) except for i = j. 
            sum = a[i][j]; 
            for(k = 0; k < i; k++)
                sum -= a[i][k]*a[k][j];
            a[i][j] = sum;
        } 
        big = 0.0; // Initialize for the search for largest pivot element. 
        for(i = j; i < n; i++) { // This is i = j of equation (2.3.12) and 
                                 // i=j+1: ::N of equation (2.3.13). 
            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) { //Is the figure of merit 
                // for the pivot better than the best so far?
                big = dum; 
                imax = i; 
            } 
        } 
        if(j != imax) { // Do we need to interchange rows? 
            for (k = 0; k < n; k++) { // Yes, do so...
                dum = a[imax][k];
                a[imax][k] = a[j][k];
                a[j][k] = dum; 
            } 
            d = -d; // ...and change the parity of d.
            vv[imax] = vv[j]; // Also interchange the scale factor. 
        }
        indx[j] = imax; 

        // If the pivot element is zero the matrix is singular 
        // (at least to the precision of the algorithm). For some 
        // applications on singular matrices, it is desirable 
        // to substitute TINY for zero.
        if(fabs(a[j][j]) < TINY) a[j][j] = TINY;
        if(j != n) { // Now, finally, divide by the pivot element. 
            dum = 1.0 / a[j][j]; 
            for(i = j+1; i < n; i++) a[i][j] *= dum;
        } 
    }
    delete[] vv;
    return true;
}

void TMatrix::LU_backsub(TMatrix& m,int* indx,double* b)
{
    double** a = m.m_dVals;
    double sum = 0.0;
    int ii = -1,ip = 0,i,j,n = m.m_nRows;

    for(i = 0; i < n; i++) { 
        // When ii is set to a positive value, it will become the 
        // index of the first nonvanishing element of b. We now do 
        // the forward substitution, equation (2.3.6). The only 
        // new wrinkle is to unscramble the permutation as we go. 
        ip = indx[i]; 
        sum = b[ip]; 
        b[ip] = b[i]; 
        if(ii != -1) for(j = ii; j <= i-1; j++)
            sum -= a[i][j]*b[j];
        else if(sum) ii=i;
        // A nonzero element was encountered, so from now on we will 
        // have to do the sums in the loop above. 
        b[i] = sum; 
    } 
    for(i = n-1; i >= 0; i--) { 
      // Now we do the backsubstitution, equation (2.3.7). 
      sum = b[i]; 
      for(j = i+1; j < n; j++) sum -= a[i][j]*b[j];
      b[i] = sum/a[i][i]; // Store a component of the solution vector X. 
    } // All done!
}