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