tsqr_elemental.cpp

#include <El.hpp>
#include <stdio.h>
#include <iostream>
#include <vector>
#include <stdlib.h>
#include <cmath>
#include <time.h>

#include <boost/program_options.hpp>
#include <omp.h>

using namespace El;

std::tuple<double, double> validation(DistMatrix<double,VC,  STAR>& Q, DistMatrix<double,STAR,STAR>& R, DistMatrix<double,VC,  STAR>& A){

    const Grid& g = A.Grid();
    const Int m = A.Height();
    const Int n = A.Width();
    const Int maxDim = Max(m,n);
    const double eps = 5*10e-16;
    const double oneNormA = OneNorm( A );

    // Form I - Q^H Q
    DistMatrix<double> Z(g);

    Identity( Z, n, n );

    Herk( UPPER, ADJOINT, -1.0, Q, 1.0, Z );

    const double infOrthogError = HermitianInfinityNorm( UPPER, Z );

    const double relOrthogError = infOrthogError / (eps*maxDim); 

    // Form A - Q R
    LocalGemm( NORMAL, NORMAL, -1.0, Q, R, 1.0, A );

    const double infError = InfinityNorm( A );

    const double relError = infError / (eps*maxDim*oneNormA);

    return std::make_tuple(relOrthogError, relError);

}

DistMatrix<double,STAR,STAR> ConstructR(AbstractDistMatrix<double>& A, qr::TreeData<double>& treeData )
{
    if( A.RowDist() != STAR ){

        LogicError("Invalid row distribution for TSQR");

    }

    const Grid& g = A.Grid();

    DistMatrix<double,CIRC,CIRC> RRoot(g);

    if( A.ColRank() == 0 )
    {
        const Int n = A.Width();

        auto R = qr::ts::RootQR(A,treeData);

        auto RTop = R( IR(0,n), IR(0,n) );

        RRoot.CopyFromRoot( RTop );

        MakeTrapezoidal( UPPER, RRoot );

    }
    else