TDecompQRH does not work for m>n matrices

[JSNielsen]

ROOT Version: ROOT 6.28/06
Platform: Ubuntu 24
Compiler: gcc

It states in the User’s guide section 14.6.6 that it should "Decompose a (m xn ) - matrix (A) with m >= n.
In this piece of code I am unable to do so:

#include "TMatrix.h"
#include "TDecompQRH.h"
#include "iostream"

using namespace std;

int main(){
    TMatrixD A(3, 2);
    for(int i = 0; i < A.GetNrows(); i++){
        for(int j = 0; j < A.GetNcols(); j++){
            A[i][j] = 1.0-2.0*rand() / RAND_MAX;
        }
    }
    cout << "A: " << endl;
    A.Print();
    TDecompQRH decomp(A);
    bool succes = decomp.Decompose();
    cout << succes << endl;
    auto R = decomp.GetTriangularMatrix();
    auto Q = decomp.GetOrthogonalMatrix();
    cout << "A: " << endl;
    A.Print();
    cout << "QR: " << endl;
    (Q*R).Print();
}

It crashes with error messages:

Error in <TMatrixTColumn_const(const TMatrixT &,Int_t)>: column index out of bounds

 *** Break *** segmentation violation

It works fine if A has m=n (that is, the decomposeQR.c example works). Am I doing something wrong?

[dastudillo]

I think the problem is actually on this line:

    auto Q = decomp.GetOrthogonalMatrix();

Looking at the ROOT code for GetOrthogonalMatrix, maybe this line

         const TVectorD qc_k = TMatrixDColumn_const(fQ, k);

should be this instead:

         const TVectorD qc_k = TMatrixDRow_const(fQ, k);

because fQ in your case is 3x2 (rows x cols), so k=2 in the first iteration of that loop (since nQrow=3) and fQ only has column indices (k) 0 and 1, thus the complaint.

[JSNielsen]

Hi,

Thanks for the answer.
It seems likely that that would be the required change, but when I tried recompiling ROOT with the change, I still got the error message.