MatrixOperations Class Reference

#include <SP_MatrixOperations.h>

List of all members.

Static Public Member Functions
static SPtMatrix< double >	EigenValsByQRFactorization (const SPtMatrix< double > &, int=-6, SPtMatrix< double > *=NULL)
static SPtMatrix< double >	HouseHolderMatrix (const SPtMatrix< double > &, SPtMatrix< double > *=NULL, SPtMatrix< double > *=NULL)
static bool	QRFactorization (const SPtMatrix< double > &, SPtMatrix< double > &, SPtMatrix< double > &)

---

Detailed Description

Definition at line 23 of file SP_MatrixOperations.h.

---

Member Function Documentation

SPtMatrix< double > MatrixOperations::EigenValsByQRFactorization	(	const SPtMatrix< double > &	A,
		int	times = -6,
		SPtMatrix< double > *	ptrO = NULL
	)			[static]

Definition at line 143 of file SP_MatrixOperations.cpp.

{
    // 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()

SPtMatrix< double > MatrixOperations::HouseHolderMatrix	(	const SPtMatrix< double > &	A,
		SPtMatrix< double > *	ptrH = NULL,
		SPtMatrix< double > *	ptrR = NULL
	)			[static]

Definition at line 17 of file SP_MatrixOperations.cpp.

{
    // 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()

bool MatrixOperations::QRFactorization	(	const SPtMatrix< double > &	A,
		SPtMatrix< double > &	Q,
		SPtMatrix< double > &	R
	)			[static]

Definition at line 91 of file SP_MatrixOperations.cpp.

{
    // 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()

---

The documentation for this class was generated from the following files:

• G:/spunak_xp/Works/TAPs/Include/TAPs/External/SP_MatrixOperations.h
• G:/spunak_xp/Works/TAPs/Include/TAPs/External/SP_MatrixOperations.cpp