QR-PerfTest.cpp

Go to the documentation of this file.

 
#include <chrono>
#include <iomanip>
#include <iostream>
 
#include <linalg/HouseholderQR.hpp>
#include <Eigen/QR>
using EigenMat =
    Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>;
 
int main(int argc, const char *argv[]) {
    // First (optional) command line argument is the matrix size
    size_t size = 500;
    if (argc > 1) {
        size = std::stoi(argv[1]);
    }
 
    // Generate a random matrix of the given size
    Matrix A = Matrix::random(size, size, -1, +1);
    // Convert it to an Eigen matrix
    EigenMat A_eigen = EigenMat::Zero(size, size);
    std::copy(A.begin(), A.end(), A_eigen.data());
 
    {
        // Compute the QR factorization using our naive implementation:
        HouseholderQR qr;
        Matrix A_cpy = A; // make a copy here to prevent allocation while timing
        auto start   = std::chrono::high_resolution_clock::now();
        qr.compute(std::move(A_cpy));
        auto finish = std::chrono::high_resolution_clock::now();
 
        std::chrono::duration<double> elapsed = finish - start;
        Matrix QR      = qr.apply_Q(qr.get_R()); // Multiply Q×R
        double err_fro = (A - QR).normFro();
        std::cout << "Elapsed time HouseholderQR: " << elapsed.count() << " s\n"
                  << "Error QR - A in Frobenius norm: " << std::scientific
                  << err_fro << std::defaultfloat << std::endl;
    }
    {
        // Now perform the factorization using Eigen:
        Eigen::HouseholderQR<EigenMat> qr_eigen;
        EigenMat A_eigen_cpy = A_eigen;
        auto start           = std::chrono::high_resolution_clock::now();
        qr_eigen.compute(std::move(A_eigen_cpy));
        auto finish = std::chrono::high_resolution_clock::now();
 
        std::chrono::duration<double> elapsed = finish - start;
        EigenMat QR = qr_eigen.matrixQR().triangularView<Eigen::Upper>();
        qr_eigen.householderQ().applyThisOnTheLeft(QR); // Multiply Q×R
        auto err_fro = (A_eigen - QR).norm();
        std::cout << "Elapsed time Eigen: " << elapsed.count() << " s\n"
                  << "Error QR - A in Frobenius norm: " << std::scientific
                  << err_fro << std::defaultfloat << std::endl;
    }
}