qr_decomposition.cpp File Reference

Program to compute the QR decomposition of a given matrix. More...

#include <array>
#include <cmath>
#include <cstdlib>
#include <ctime>
#include <iostream>
#include "./qr_decompose.h"

Include dependency graph for qr_decomposition.cpp:

Go to the source code of this file.

Functions
int	main (void)
template<typename T>
void	qr_decompose (const std::valarray< std::valarray< T > > &A, std::valarray< std::valarray< T > > *Q, std::valarray< std::valarray< T > > *R)
template<typename T>
std::ostream &	operator<< (std::ostream &out, std::valarray< std::valarray< T > > const &v)

Detailed Description

Program to compute the QR decomposition of a given matrix.

Author

Krishna Vedala

Definition in file qr_decomposition.cpp.

Function Documentation

main()

int main

(

void

)

main function

Definition at line 23 of file qr_decomposition.cpp.

               {
    unsigned int ROWS, COLUMNS;
 
    std::cout << "Enter the number of rows and columns: ";
    std::cin >> ROWS >> COLUMNS;
 
    std::cout << "Enter matrix elements row-wise:\n";
 
    std::valarray<std::valarray<double>> A(ROWS);
    std::valarray<std::valarray<double>> Q(ROWS);
    std::valarray<std::valarray<double>> R(COLUMNS);
    for (int i = 0; i < std::max(ROWS, COLUMNS); i++) {
        if (i < ROWS) {
            A[i] = std::valarray<double>(COLUMNS);
            Q[i] = std::valarray<double>(COLUMNS);
        }
        if (i < COLUMNS) {
            R[i] = std::valarray<double>(COLUMNS);
        }
    }
 
    for (int i = 0; i < ROWS; i++)
        for (int j = 0; j < COLUMNS; j++) std::cin >> A[i][j];
 
    std::cout << A << "\n";
 
    clock_t t1 = clock();
    qr_decompose(A, &Q, &R);
    double dtime = static_cast<double>(clock() - t1) / CLOCKS_PER_SEC;
 
    std::cout << Q << "\n";
    std::cout << R << "\n";
    std::cout << "Time taken to compute: " << dtime << " sec\n ";
 
    return 0;
}

operator<<()

template<typename T>

std::ostream & qr_algorithm::operator<<	(	std::ostream &	out,
		std::valarray< std::valarray< T > > const &	v )

operator to print a matrix

Definition at line 33 of file qr_decompose.h.

                                                               {
    const int width = 12;
    const char separator = ' ';
 
    out.precision(4);
    for (size_t row = 0; row < v.size(); row++) {
        for (size_t col = 0; col < v[row].size(); col++)
            out << std::right << std::setw(width) << std::setfill(separator)
                << v[row][col];
        out << std::endl;
    }
 
    return out;
}

qr_decompose()

template<typename T>

void qr_algorithm::qr_decompose	(	const std::valarray< std::valarray< T > > &	A,
		std::valarray< std::valarray< T > > *	Q,
		std::valarray< std::valarray< T > > *	R )

Decompose matrix \(A\) using Gram-Schmidt process.

\begin{eqnarray*} \text{given that}\quad A &=& *\left[\mathbf{a}_1,\mathbf{a}_2,\ldots,\mathbf{a}_{N-1},\right]\\ \text{where}\quad\mathbf{a}_i &=& \left[a_{0i},a_{1i},a_{2i},\ldots,a_{(M-1)i}\right]^T\quad\ldots\mbox{(column vectors)}\\ \text{then}\quad\mathbf{u}_i &=& \mathbf{a}_i *-\sum_{j=0}^{i-1}\text{proj}_{\mathbf{u}_j}\mathbf{a}_i\\ \mathbf{e}_i &=&\frac{\mathbf{u}_i}{\left|\mathbf{u}_i\right|}\\ Q &=& \begin{bmatrix}\mathbf{e}_0 & \mathbf{e}_1 & \mathbf{e}_2 & \dots & \mathbf{e}_{N-1}\end{bmatrix}\\ R &=& \begin{bmatrix}\langle\mathbf{e}_0\,,\mathbf{a}_0\rangle & \langle\mathbf{e}_1\,,\mathbf{a}_1\rangle & \langle\mathbf{e}_2\,,\mathbf{a}_2\rangle & \dots \\ 0 & \langle\mathbf{e}_1\,,\mathbf{a}_1\rangle & \langle\mathbf{e}_2\,,\mathbf{a}_2\rangle & \dots\\ 0 & 0 & \langle\mathbf{e}_2\,,\mathbf{a}_2\rangle & \dots\\ \vdots & \vdots & \vdots & \ddots \end{bmatrix}\\ \end{eqnarray*}

Parameters

A	input matrix to decompose
Q	output decomposed matrix
R	output decomposed matrix

Definition at line 146 of file qr_decompose.h.

  {
    std::size_t ROWS = A.size();        // number of rows of A
    std::size_t COLUMNS = A[0].size();  // number of columns of A
    std::valarray<T> col_vector(ROWS);
    std::valarray<T> col_vector2(ROWS);
    std::valarray<T> tmp_vector(ROWS);
 
    for (int i = 0; i < COLUMNS; i++) {
        /* for each column => R is a square matrix of NxN */
        int j;
        R[0][i] = 0.; /* make R upper triangular */
 
        /* get corresponding Q vector */
#ifdef _OPENMP
// parallelize on threads
#pragma omp for
#endif
        for (j = 0; j < ROWS; j++) {
            tmp_vector[j] = A[j][i]; /* accumulator for uk */
            col_vector[j] = A[j][i];
        }
        for (j = 0; j < i; j++) {
            for (int k = 0; k < ROWS; k++) {
                col_vector2[k] = Q[0][k][j];
            }
            col_vector2 = vector_proj(col_vector, col_vector2);
            tmp_vector -= col_vector2;
        }
 
        double mag = vector_mag(tmp_vector);
 
#ifdef _OPENMP
// parallelize on threads
#pragma omp for
#endif
        for (j = 0; j < ROWS; j++) Q[0][j][i] = tmp_vector[j] / mag;
 
            /* compute upper triangular values of R */
#ifdef _OPENMP
// parallelize on threads
#pragma omp for
#endif
        for (int kk = 0; kk < ROWS; kk++) {
            col_vector[kk] = Q[0][kk][i];
        }
 
#ifdef _OPENMP
// parallelize on threads
#pragma omp for
#endif
        for (int k = i; k < COLUMNS; k++) {
            for (int kk = 0; kk < ROWS; kk++) {
                col_vector2[kk] = A[kk][k];
            }
            R[0][i][k] = (col_vector * col_vector2).sum();
        }
    }
}