SP_MatrixOperations.cpp

Go to the documentation of this file.

/******************************************************************************
SP_MatrixOperations.cpp

Operation on Matrix

SUKITTI PUNAK   (07/14/2007)
UPDATE          (07/14/2007)
******************************************************************************/
#include "SP_MatrixOperations.h"
// Using Inclusion Model (i.e. definitions are included in declarations)
//                       (this name.cpp is included in name.hpp)
// Each friend is defined directly inside its declaration.
//-----------------------------------------------------------------------------
// Member Fn:   HouseHolderMatrix()     ---------------------------------------
// Desc:    Find a HouseHolder matrix of a real square matrix only.
SPtMatrix< double > MatrixOperations::HouseHolderMatrix( const SPtMatrix< double > &A,  // I/P
                                                         SPtMatrix< double > *ptrH,     // O/P
                                                         SPtMatrix< double > *ptrR )    // O/P
{
    // only for square matrix
    if ( A.m_iRows != A.m_iCols )   return NULL;
    if ( ptrR != NULL ) {
        if ( ptrR->m_iRows != A.m_iRows || ptrR->m_iCols != A.m_iCols ) return false;
    }

    SPtMatrix< double > B( A );
    SPtMatrix< double > H( A );

    int i, k;
    double g, s;
    int n = H.m_iRows;
    SPtMatrix< double > W( n, 1 );
    //SPtMatrix< double > V( n, 1 );
    SPtMatrix< double > U( n, 1 );
    
    // Create Identity matrix
    SPtMatrix< double > I( n, n );
    I.SetToIdentity();
    

    for ( k = 0; k < n-1; k++ ) {
        //cout << "\nk = " << k << "\n";

        // Step 1: Set W_i = 0 for i = 0, ... , k-2
        for ( i = 0; i < k; i++ ) {
            W.m_tElement[i][0] = 0.0;
        }

        // Step 2.1: Find g
        g = 0;
        for ( i = k; i < n; i++ ) {
            g += B.m_tElement[i][k] * B.m_tElement[i][k];
        }
        g = sqrt(g);

        // Step 2.2: Find s
        s = sqrt( 2 * g * ( g + fabs(B.m_tElement[k][k]) ) );

        // Step 3.1: Set W_k
        W.m_tElement[k][0] = ( B.m_tElement[k][k] + (B.m_tElement[k][k] >= 0.0 ? +1.0 : -1.0)*g ) / s;

        // Step 3.2: Set W_i
        for ( i = k+1; i < n; i++ ) {
            W.m_tElement[i][0] = B.m_tElement[i][k] / s;
        }

        //cout << "W" << endl;
        //W.DisplayElements();

        U = 2 * B.GetTranspose() * W;

        //cout << "W*U'" << endl;
        //(W*U.GetTranspose()).DisplayElements();

        //H -= W*U.GetTranspose();

        H = I - ( 2 * W * W.GetTranspose() );
        if ( ptrH != NULL )     ptrH[k] = H;

        B = H*B;
    }

    if ( ptrR != NULL )     *ptrR = B;

    return H;
} // END Member Fn: HouseHolderMatrix()
//-----------------------------------------------------------------------------
// Member Fn:   QRFactorization()       ---------------------------------------
// Desc:    Find a QR Factorization of a real square matrix only.
bool MatrixOperations::QRFactorization( const SPtMatrix< double > &A,                       // I/P
                                        SPtMatrix< double > &Q, SPtMatrix< double > &R )    // O/P
{
    // only square matrix
    if ( A.m_iRows != A.m_iCols )   return false;
    if ( Q.m_iRows != A.m_iRows || Q.m_iCols != A.m_iCols ) return false;
    if ( R.m_iRows != A.m_iRows || R.m_iCols != A.m_iCols ) return false;

    int n = A.m_iRows;
    int k;

    // Define R_(0) = A
    //R = A;
    SPtMatrix< double > *arrayH = new SPtMatrix< double >[ n-1 ];

    // For k = 1,...,n-1
    for ( k = 0; k < n-1; k++ ) {
        // Find H_(k) in order to reduce positions k+1,...,n, 
        // in the k_th column of R_(k-1) to zero
        HouseHolderMatrix( A, arrayH, &R );
        //R = ptrH[n-2] * R;
    }

    // Define Q = I;
    // For k = n-1,...,1
    //    Q = H_(k) * Q
    Q.SetToIdentity();
    for ( k = n-2; k >= 0; k-- ) {
        Q = arrayH[k]*Q;
    }

    //cout << "R is\n";
    //R.DisplayElements();
    //cout << "Q is\n";
    //Q.DisplayElements();
    //cout << "Q'*Q is ";
    //(Q.GetTranspose()*Q).DisplayElements();
    //cout << "Q*R is\n";
    //(Q*R).DisplayElements();

    delete [] arrayH;

    return true;
} // END Member Fn: QRFactorization()
//-----------------------------------------------------------------------------
// Member Fn:   EigenValsByQRFactorization()        ---------------------------
// Desc:    Find all eigen values of a real square matrix only.
// I/P:  A is a real square matrix.
// I/P:  (default = 6), times is the number of iterations.
// O/P:  (default = NULL), ptrO is a pointer to a matrix with will be the same 
//       order as I/P Matrix, A.  Its diagonal are the approximated eigen values.
// Return:  A column vector (n-by-1 matrix) contains all approximated eigenvalues.
SPtMatrix< double > MatrixOperations::EigenValsByQRFactorization( const SPtMatrix< double > &A, // I/P
                                                   int times,                       // Number of iterations
                                                   SPtMatrix< double > *ptrO )      // O/P
{
    // only square matrix
    if ( A.m_iRows != A.m_iCols )   return NULL;

    // Given a real n-by-n matrix A:
    // Define A_(1) = A
    // For k = 1 to times
    //      Factor A_(k)   =  Q_(k)R_(k)        (so that Q_(k)'A_(k) = R_(k))
    //      Define A_(k+1) =  R_(k)Q_(k)        (so that A_(k+1) = Q_(k)'A_(k)Q_(k))
    // end.
    double errThreshold = 1.0;
    int limit = 10000;
    int n = A.m_iRows;
    SPtMatrix< double > M( A ), Q( n, n ), R( n, n );
    if ( times < 0 ) {
        errThreshold = pow( 10.0, times );
        times = limit;
    }

    for ( int k = 0; k < times; k++ ) {
        //cout << k << " ";
        
        //M.DisplayElements();
        
        QRFactorization( M, Q, R );
        
        //cout << "R is \n";
        //R.DisplayElements();
        //cout << "Q is \n";
        //Q.DisplayElements();

        M = R*Q;
    
        if ( 1.0 > errThreshold ) {
            bool isEnough = true;
            for ( int r = n/2; r < n; r++ ) {
                for ( int c = 0; c <= n/2; c++ ) {
                    if ( r > c ) {
                        //cout << "(" << r << ", " << c << ") ";
                        if ( fabs(M.m_tElement[r][c]) > errThreshold ) {
                            isEnough = false;
                            break;
                        }
                    }
                }
            }
            if ( isEnough ) {
                cout << "Stopped at Iteration# " << k << "\n";
                break;
            }
        }
    }
    //cout << "M_(" << k << ")\n";
    //M.DisplayElements();

    if ( ptrO != NULL ) *ptrO = M;

    SPtMatrix< double > outMatrix( n, 1 );
    for ( int i = 0; i < n; i++ ) {
        outMatrix.m_tElement[i][0] = M.m_tElement[i][i];
    }

    return M;
    //return outMatrix;
} // END Member Fn: EigenValsByQRFactorization()
//-----------------------------------------------------------------------------
//34567890123456789012345678901234567890123456789012345678901234567890123456789
//--+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----