change C like we changed R: no re-weighting in CCA permutations

similar change was made in the R code in commit e78bbee, and
now both C and R give the same results which are compatible
with Canoco 3.12.

R/permutest.cca.R

@@ -78,8 +78,7 @@ permutest.default <- function(x, ...)
     {
         if (!is.matrix(indx))
             indx <- matrix(indx, nrow=1)
-        out <- .Call("do_getF", indx, E, Q, QZ, w, first, isPartial,
-                     isCCA, isDB)
+        out <- .Call("do_getF", indx, E, Q, QZ, first, isPartial, isDB)
         if (!isPartial && !first)
             out[,2] <- Chi.tot - out[,1]
         out <- cbind(out, (out[,1]/q)/(out[,2]/r))

src/getF.c

@@ -9,16 +9,10 @@
 #include <Rinternals.h>
 #include <R_ext/Linpack.h> /* QR */
 #include <R_ext/Lapack.h>  /* SVD, eigen */
-#include <R_ext/Applic.h> /* R version of QR decomposition dqrdc2: not
-			     in public API */
 
 #include <math.h> /* sqrt */
 #include <string.h> /* memcpy */
 
-/* The following file is in goffactor.c file in vegan */
-extern
-void wcentre(double *, double *, int *, int *);
-
 /* LINPACK uses the same function (dqrsl) to find derived results from
  * the QR decomposition. It uses decimal coding to define the kind of
  * output with following alternatives (although we will not use all of
@@ -268,13 +262,13 @@ SEXP do_QR(SEXP x)
  * permutest.cca. The do_getF provides a drop-in replacement to the R
  * function, and is called directly the R function */
 
-SEXP do_getF(SEXP perms, SEXP E, SEXP QR, SEXP QZ, SEXP w, SEXP first,
-	     SEXP isPartial, SEXP isCCA, SEXP isDB)
+SEXP do_getF(SEXP perms, SEXP E, SEXP QR, SEXP QZ, SEXP first,
+	     SEXP isPartial, SEXP isDB)
 {
     int i, j, k, ki,
 	nperm = nrows(perms), nr = nrows(E), nc = ncols(E),
 	FIRST = asInteger(first), PARTIAL = asInteger(isPartial),
-	WEIGHTED = asInteger(isCCA), DISTBASED = asInteger(isDB);
+	DISTBASED = asInteger(isDB);
     double ev1;
     SEXP ans = PROTECT(allocMatrix(REALSXP, nperm, 2));
     double *rans = REAL(ans);
@@ -303,27 +297,6 @@ SEXP do_getF(SEXP perms, SEXP E, SEXP QR, SEXP QZ, SEXP w, SEXP first,
     double *qty = (double *) R_alloc(nr, sizeof(double));
     double dummy;
     int info, qrkind;
-    /* Weighted methods currently need re-evaluation of QR
-       decomposition (probably changed in the future, but now for the
-       compatibility with the current code). For this we need to
-       reconstruct constraints and conditions. */
-    int nx = ncols(VECTOR_ELT(QR, 0)), nz = 0, *pivot, *zpivot;
-    double *wperm, *Xorig, *Zorig, *qrwork, *zqrwork, qrtol=1e-7; 
-    if (WEIGHTED) {
-	if (PARTIAL) {
-	    nz = ncols(VECTOR_ELT(QZ, 0));
-	    Zorig = (double *) R_alloc(nr * nz, sizeof(double));
-	    qrXw(Zqr, Zqrank, Zqraux, Zorig, REAL(w), nr, nz);
-	    zpivot = (int *) R_alloc(nz, sizeof(int));
-	    zqrwork = (double *) R_alloc(2 * nz, sizeof(double));
-	}
-	pivot = (int *) R_alloc(nx > nz ? nx : nz, sizeof(int));
-	wperm = (double *) R_alloc(nr, sizeof(double));
-	Xorig = (double *) R_alloc(nr * nx, sizeof(double));
-	qrXw(qr, qrank, qraux, Xorig, REAL(w), nr, nx);
-	pivot = (int *) R_alloc(nx, sizeof(int));
-	qrwork = (double *) R_alloc(2 * nx, sizeof(double));
-    }
 
     /* distance-based methods need to transpose data */
     double *transY;
@@ -349,21 +322,10 @@ SEXP do_getF(SEXP perms, SEXP E, SEXP QR, SEXP QZ, SEXP w, SEXP first,
 		else   /* shuffle rows */
 		    rY[i + nr*j] = REAL(E)[ki + nr*j];
 	    }
-	    if (WEIGHTED)
-		wperm[i] = REAL(w)[ki];
 	}
 
 	/* Partial model: qr.resid(QZ, Y) with LINPACK */
 	if (PARTIAL) {
-	    /* Re-do QR decomposition with changed weights */
-	    if (WEIGHTED) {
-	        memcpy(Zqr, Zorig, nr * nz * sizeof(double));
-		/* wcentre is in goffactor.c file */
-		wcentre(Zqr, wperm, &nr, &nz);
-		/* dqrdc2 is not in R API */
-		F77_CALL(dqrdc2)(Zqr, &nr, &nr, &nz, &qrtol, &Zqrank,
-	                         Zqraux, zpivot, zqrwork);
-	    }
 	    qrkind = RESID;
 	    for(i = 0; i < nc; i++)
 		F77_CALL(dqrsl)(Zqr, &nr, &nr, &Zqrank, Zqraux, rY + i*nr,
@@ -382,14 +344,6 @@ SEXP do_getF(SEXP perms, SEXP E, SEXP QR, SEXP QZ, SEXP w, SEXP first,
 
 	/* CONSTRAINED COMPONENT */
 
-	/* Re-weight constraints are re-do QR */
-	if (WEIGHTED) {
-	    memcpy(qr, Xorig, nr * nx * sizeof(double));
-	    wcentre(qr, wperm, &nr, &nx);
-	    F77_CALL(dqrdc2)(qr, &nr, &nr, &nx, &qrtol, &qrank,
-			     qraux, pivot, qrwork);  
-	}
-
 	/* qr.fitted(QR, Y) + qr.resid(QR, Y) with LINPACK */
 	if (PARTIAL || FIRST)
 	    qrkind = FIT + RESID;