[Rd] typo or stale info in qr man

[Jari Oksanen]

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-bounces at r-project.org> on behalf of Martin Maechler <maechler at stat.math.ethz.ch>
Sent: 25 October 2016 11:08
To: Wojciech Musial (Voitek)
Cc: R-devel at r-project.org
Subject: Re: [Rd] typo or stale info in qr man

>>>>> Wojciech Musial (Voitek) <wojciech.musial at 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 at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel