Revert "remove C code for QR decomposition which is no longer needed"

This reverts commit 9a98a8b.

src/getF.c

@@ -154,6 +154,119 @@ static void transpose(double *x, double *tx, int nr, int nc)
 	    tx[ij++] = x[j * nr + i];
 }
 
+/* test transpose from an R session */
+
+SEXP test_trans(SEXP x)
+{
+    int nr = nrows(x), nc = ncols(x);
+    SEXP tx = PROTECT(allocMatrix(REALSXP, nc, nr));
+    transpose(REAL(x), REAL(tx), nr, nc);
+    UNPROTECT(1);
+    return tx;
+}
+
+/* Reconstruct data X from weighted QR decomposition. We need this
+   only for CCA, and there the data are weighted by row sums, but we
+   need the original unweighted data. So we do here both qrX and
+   de-weighting. */
+
+static void qrXw(double *qr, int rank, double *qraux, double *X, double *w,
+		 int nr, int nc)
+{
+    int i, j, ij, len = nr*nc, info = 0, qrkind;
+    double dummy = 0, wsqrt;
+    /* Extract  R from qr into upper triangle of X */
+    for(i = 0; i < len; i++)
+	X[i] = 0;
+    for(j = 0; j < nc; j++)
+	for(i = 0; i <= j; i++) {
+	    ij = i + nr*j;
+	    X[ij] = qr[ij];
+	}
+    /* Find data as Qy: if y = R then X = QR. The data will over-write
+       R. No pivoting, and aliased variables will be moved to last
+       columns. Uses Linpack. */
+    qrkind = QY;
+    for(j = 0; j < nc; j++)
+	F77_CALL(dqrsl)(qr, &nr, &nr, &rank, qraux, X + j*nr, X + j*nr,
+			&dummy, &dummy, &dummy, &dummy, &qrkind, &info);
+
+    /* de-weight X */
+    for(i = 0; i < nr; i++) {
+	wsqrt = sqrt(w[i]);
+	for (j = 0; j < nc; j++)
+	    X[i + nr*j] /= wsqrt;
+    }
+}
+
+/* function to test qrX from R. Use with CCA model 'm' as
+   .Call("test_qrXw", m$CCA$QR, weights(m)) */
+
+SEXP test_qrXw(SEXP QR, SEXP w)
+{
+    int nc, nr;
+    double *qr = REAL(VECTOR_ELT(QR, 0));
+    int rank = asInteger(VECTOR_ELT(QR, 1));
+    double *qraux = REAL(VECTOR_ELT(QR, 2));
+    nr = nrows(VECTOR_ELT(QR, 0));
+    nc = ncols(VECTOR_ELT(QR, 0));
+    SEXP X = PROTECT(allocMatrix(REALSXP, nr, nc));
+    qrXw(qr, rank, qraux, REAL(X), REAL(w), nr, nc);
+    UNPROTECT(1);
+    return X;
+}
+
+/* QR decomposition. Actually R does not use this function, but an
+   edited version. From the user point of view, the main difference is
+   that the R function returns rank, but this does not do so
+   directly. So we test here if we can make this compatible with
+   R. This function can be called as .Call("do_QR", x), where x is a
+   (centred) matrix. The function returns an object of class "qr" and
+   similar to R::qr() result object, expect that the order of columns
+   (pivoting) is different than in R.  */
+
+SEXP do_QR(SEXP x)
+{
+    /* set up */
+    int k;
+    int nr = nrows(x), nx = ncols(x);
+    double TOL = 1e-7;
+    SEXP qraux = PROTECT(allocVector(REALSXP, nx));
+    SEXP pivot = PROTECT(allocVector(INTSXP, nx));
+    memset(INTEGER(pivot), 0, nx * sizeof(int));
+    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);
+    /* get rank */
+    for (k = 1; k < nx; k++) {
+	if (fabs(REAL(x)[nr * k + k]) < fabs(TOL * REAL(x)[0]))
+	    break;
+    }
+
+    /* 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(k));
+    SET_VECTOR_ELT(qr, 2, qraux);
+    SET_VECTOR_ELT(qr, 3, pivot);
+    UNPROTECT(4); /* qr, x, pivot, qraux */
+    return qr;
+}
+
 /* Function do_getF is modelled after R function getF embedded in
  * permutest.cca. The do_getF provides a drop-in replacement to the R
  * function, and is called directly the R function */

src/goffactor.c

@@ -47,7 +47,7 @@ static void goffactor(double *ord, int *f, double *w, int *nrow, int *ndim,
 
 #include <math.h> /* sqrt */
 
-static void wcentre(double *x, double *w, int *nr, int *nc)
+void wcentre(double *x, double *w, int *nr, int *nc)
 {
      int i, j, ij;
      double sw, swx;