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 <[email protected]> on behalf of Martin Maechler
<[email protected]>
Sent: 25 October 2016 11:08
To: Wojciech Musial (Voitek)
Cc: [email protected]
Subject: Re: [Rd] typo or stale info in qr man
>>>>> Wojciech Musial (Voitek) <[email protected]>
>>>>> 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
______________________________________________
[email protected] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
______________________________________________
[email protected] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
