qr_decompose.h

Go to the documentation of this file.

 
#ifndef NUMERICAL_METHODS_QR_DECOMPOSE_H_
#define NUMERICAL_METHODS_QR_DECOMPOSE_H_
 
#include <cmath>
#include <cstdlib>
#include <iomanip>
#include <iostream>
#include <limits>
#include <numeric>
#include <valarray>
#ifdef _OPENMP
#include <omp.h>
#endif

namespace qr_algorithm {
template <typename T>
std::ostream &operator<<(std::ostream &out,
                         std::valarray<std::valarray<T>> const &v) {
    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;
}

template <typename T>
std::ostream &operator<<(std::ostream &out, std::valarray<T> const &v) {
    const int width = 10;
    const char separator = ' ';
 
    out.precision(4);
    for (size_t row = 0; row < v.size(); row++) {
        out << std::right << std::setw(width) << std::setfill(separator)
            << v[row];
    }
 
    return out;
}

template <typename T>
inline double vector_dot(const std::valarray<T> &a, const std::valarray<T> &b) {
    return (a * b).sum();
    // could also use following
    // return std::inner_product(std::begin(a), std::end(a), std::begin(b),
    // 0.f);
}

template <typename T>
inline double vector_mag(const std::valarray<T> &a) {
    double dot = vector_dot(a, a);
    return std::sqrt(dot);
}

template <typename T>
std::valarray<T> vector_proj(const std::valarray<T> &a,
                             const std::valarray<T> &b) {
    double num = vector_dot(a, b);
    double deno = vector_dot(b, b);

    if (deno <= std::numeric_limits<double>::epsilon()) {
        std::cerr << "[" << __func__ << "] Possible division by zero\n";
        return a;  // return vector a back
    }
 
    double scalar = num / deno;
 
    return b * scalar;
}

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        
) {
    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();
        }
    }
}
}  // namespace qr_algorithm
 
#endif  // NUMERICAL_METHODS_QR_DECOMPOSE_H_