And that missing functionality is that Linpack/Lapack routines do not return 
rank and have a different style of pivoting? For other aspects, the 
user-interface is very similar in dqrdc2 in R and in dqrdc in Linpack. Another 
difference seems to be that the final pivoting reported to the user is 
different: R keeps the original order except for aliased variables, but Linpack 
makes either wild shuffling or no pivoting at all. I haven't looked at dqpq3 in 
Lapack, but it appears to return no rank either (don't know about shuffling the 
columns). It seems that using Linpack dqrdc directly is not always compatible 
with dqrdc2 of R although it returns similar objects. That is, when packing up 
the Linpack function to produce an object with same items as qr.default (qr, 
rank, qraux, pivot, class "qr"), the result object may not yield similar 
results in base::qr.fitted, base::qr.resid etc as base::qr.default result (but 
I haven't had time for thorough testing).

This is how I tried to do the packing (apologies for clumsy coding):

SEXP do_QR(SEXP x, SEXP dopivot)
{
    /* set up */
    int i;
    int nr = nrows(x), nx = ncols(x);
    int pivoting = asInteger(dopivot);
    SEXP qraux = PROTECT(allocVector(REALSXP, nx));
    SEXP pivot = PROTECT(allocVector(INTSXP, nx));
    /* do pivoting or keep the order of columns? */
    if (pivoting)
        memset(INTEGER(pivot), 0, nx * sizeof(int));
    else
        for(i = 0; i < nx; i++)
            INTEGER(pivot)[i] = i+1;
    double *work = (double *) R_alloc(nx, sizeof(double));
    int job = 1;
    x = PROTECT(duplicate(x));

    /* QR decomposition with Linpack */
    F77_CALL(dqrdc)(REAL(x), &nr, &nr, &nx, REAL(qraux),
                    INTEGER(pivot), work, &job);

    /* pack up */
    SEXP qr = PROTECT(allocVector(VECSXP, 4));
    SEXP labs = PROTECT(allocVector(STRSXP, 4));
    SET_STRING_ELT(labs, 0, mkChar("qr"));
    SET_STRING_ELT(labs, 1, mkChar("rank"));
    SET_STRING_ELT(labs, 2, mkChar("qraux"));
    SET_STRING_ELT(labs, 3, mkChar("pivot"));
    setAttrib(qr, R_NamesSymbol, labs);
    SEXP cl = PROTECT(allocVector(STRSXP, 1));
    SET_STRING_ELT(cl, 0, mkChar("qr"));
    classgets(qr, cl);
    UNPROTECT(2); /* cl, labs */
    SET_VECTOR_ELT(qr, 0, x);
    SET_VECTOR_ELT(qr, 1, ScalarInteger(nx)); /* not really the rank,
                                                 but no. of columns */
    SET_VECTOR_ELT(qr, 2, qraux);
    SET_VECTOR_ELT(qr, 3, pivot);
    UNPROTECT(4); /* qr, x, pivot, qraux */
    return qr;
}


cheers, Jari Oksanen
________________________________________
From: R-devel <r-devel-boun...@r-project.org> on behalf of Martin Maechler 
<maech...@stat.math.ethz.ch>
Sent: 25 October 2016 11:08
To: Wojciech Musial (Voitek)
Cc: R-devel@r-project.org
Subject: Re: [Rd] typo or stale info in qr man

>>>>> Wojciech Musial (Voitek) <wojciech.mus...@gmail.com>
>>>>>     on Mon, 24 Oct 2016 15:07:55 -0700 writes:

    > man for `qr` says that the function uses LINPACK's DQRDC, while it in
    > fact uses DQRDC2.

which is a modification of LINPACK's DQRDC.

But you are right, and I have added to the help file (and a tiny
bit to the comments in the Fortran source).

When this change was done > 20 years ago, it was still hoped
that the numerical linear algebra community or more specifically
those behind LAPACK would eventually provide this functionality
with LAPACK (and we would then use that),
but that has never happened according to my knowledge.

Thank you for the 'heads up'.

Martin Maechler
ETH Zurich

______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Reply via email to