Skip to content

Commit

Permalink
Merge branch 'ubsan'
Browse files Browse the repository at this point in the history
Should fix "ubsan issue" of CRAN: argument resid in Linpack call
could be a null pointer. In that case resid was not needed nor
referenced in the Linpack function, but it was a null pointer any way.
  • Loading branch information
jarioksa committed Dec 22, 2024
2 parents 1499ab9 + e9e2cb1 commit 7f38987
Showing 1 changed file with 14 additions and 17 deletions.
31 changes: 14 additions & 17 deletions src/getF.c
Original file line number Diff line number Diff line change
Expand Up @@ -180,7 +180,7 @@ static void qrXw(double *qr, int rank, double *qraux, int *pivot, double *X,
double *w, int nr, int nc, int discard)
{
int i, j, ij, len = nr*nc, info = 0, qrkind;
double dummy = 0, wsqrt;
double dummy[1] = {0.0}, wsqrt;
double *xwork = (double *) R_alloc(len, sizeof(double));
/* Extract R from qr into upper triangle of X */
for(i = 0; i < len; i++)
Expand All @@ -201,7 +201,7 @@ static void qrXw(double *qr, int rank, double *qraux, int *pivot, double *X,
for(j = 0; j < nc; j++) {
if (pivot[j] >= 0)
F77_CALL(dqrsl)(qr, &nr, &nr, &rank, qraux, xwork + j*nr,
X + pivot[j]*nr, &dummy, &dummy, &dummy, &dummy,
X + pivot[j]*nr, dummy, dummy, dummy, dummy,
&qrkind, &info);
}

Expand Down Expand Up @@ -331,13 +331,10 @@ SEXP do_getF(SEXP perms, SEXP E, SEXP QR, SEXP QZ, SEXP effects,
}

double *fitted = (double *) R_alloc(nr * nc, sizeof(double));
/* separate resid needed only in some cases */
double *resid;
if (PARTIAL || FIRST)
resid = (double *) R_alloc(nr * nc, sizeof(double));
double *resid = (double *) R_alloc(nr * nc, sizeof(double));
/* work array and variables for QR decomposition */
double *qty = (double *) R_alloc(nr, sizeof(double));
double dummy = 0.0;
double dummy[1] = {0.0};
int info, qrkind;
/* Weighted methods currently need re-evaluation of QR
decomposition (probably changed in the future, but now for the
Expand Down Expand Up @@ -414,16 +411,16 @@ SEXP do_getF(SEXP perms, SEXP E, SEXP QR, SEXP QZ, SEXP effects,
qrkind = RESID;
for(i = 0; i < nc; i++)
F77_CALL(dqrsl)(Zqr, &nr, &nr, &Zqrank, Zqraux, rY + i*nr,
&dummy, qty, &dummy, rY + i*nr, &dummy,
dummy, qty, dummy, rY + i*nr, dummy,
&qrkind, &info);
/* distances need symmetric residuals */
if (DISTBASED) {
transpose(rY, transY, nr, nr);
qrkind = RESID;
for(i = 0; i < nc; i++)
F77_CALL(dqrsl)(Zqr, &nr, &nr, &Zqrank, Zqraux,
transY + i*nr, &dummy, qty, &dummy,
rY + i*nr, &dummy, &qrkind, &info);
transY + i*nr, dummy, qty, dummy,
rY + i*nr, dummy, &qrkind, &info);
}

}
Expand All @@ -438,8 +435,8 @@ SEXP do_getF(SEXP perms, SEXP E, SEXP QR, SEXP QZ, SEXP effects,
qrkind = RESID;
for (i = 0; i < nx; i++)
F77_CALL(dqrsl)(Zqr, &nr, &nr, &Zqrank, Zqraux,
qr + i*nr, &dummy, qty, &dummy, qr + i*nr,
&dummy, &qrkind, &info);
qr + i*nr, dummy, qty, dummy, qr + i*nr,
dummy, &qrkind, &info);
}
for(i = 0; i < nx; i++)
pivot[i] = i + 1;
Expand All @@ -456,7 +453,7 @@ SEXP do_getF(SEXP perms, SEXP E, SEXP QR, SEXP QZ, SEXP effects,
for (p = 0; p < (nterms - 1); p++) {
for (i = 0; i < nc; i++)
F77_CALL(dqrsl)(qr, &nr, &nr, term + p, qraux, rY + i*nr,
&dummy, qty, &dummy, &dummy, fitted + i*nr,
dummy, qty, dummy, dummy, fitted + i*nr,
&qrkind, &info);
ev = getEV(fitted, nr, nc, DISTBASED);
rans[k + p*nperm] = ev - ev0;
Expand All @@ -469,8 +466,8 @@ SEXP do_getF(SEXP perms, SEXP E, SEXP QR, SEXP QZ, SEXP effects,
else
qrkind = FIT;
for (i = 0; i < nc; i++)
F77_CALL(dqrsl)(qr, &nr, &nr, &qrank, qraux, rY + i*nr, &dummy,
qty, &dummy, resid + i*nr, fitted + i*nr,
F77_CALL(dqrsl)(qr, &nr, &nr, &qrank, qraux, rY + i*nr, dummy,
qty, dummy, resid + i*nr, fitted + i*nr,
&qrkind, &info);

/* Eigenvalues: either sum of all or the first If the sum of
Expand All @@ -485,8 +482,8 @@ SEXP do_getF(SEXP perms, SEXP E, SEXP QR, SEXP QZ, SEXP effects,
qrkind = FIT;
for(i = 0; i < nc; i++)
F77_CALL(dqrsl)(qr, &nr, &nr, &qrank, qraux,
transY + i*nr, &dummy, qty, &dummy,
&dummy, fitted + i*nr, &qrkind, &info);
transY + i*nr, dummy, qty, dummy,
dummy, fitted + i*nr, &qrkind, &info);
ev1 = eigenfirst(fitted, nr);
} else {
ev1 = svdfirst(fitted, nr, nc);
Expand Down

0 comments on commit 7f38987

Please sign in to comment.