[Rd] Proposal for diagonal-only chol2inv option

[Ravi Varadhan]

Hi Martin,

You are correct. I just realized that because the output matrix ans is allocated as an sz x sz matrix (rather than retaining the original n x n dimensions), using n as the leading dimension in my edited code causes out-of-bounds memory access when size < ncol(x).

Adam's original code correctly uses sz for both the LAPACK leading dimension and the pointer arithmetic, so it handles the submatrix case perfectly without any issues.

Thanks for double-checking and catching my mistake! 

Ravi

-----Original Message-----
From: Martin Maechler <maechler using stat.math.ethz.ch> 
Sent: Friday, September 25, 2026 3:25 AM
To: Ravi Varadhan <ravi.varadhan using jhu.edu>
Cc: J C Nash <profjcnash using gmail.com>; Martin Maechler <maechler using stat.math.ethz.ch>; Adam marstrand <adamboegh using gmail.com>; r-devel using r-project.org
Subject: Re: [Rd] Proposal for diagonal-only chol2inv option


      External Email - Use Caution



>>>>> Ravi Varadhan
>>>>>     on Mon, 21 Sep 2026 22:09:17 +0000 writes:

    > The code snippet provided by Adam works correctly when size = ncol(x), the default.  However, it would not work if the user wanted the diagonal of Cholesky inverse of a submatrix (size < ncol(x)).  I edited it slightly to handle this more general case.
    > if (only_diag) {
    > F77_CALL(dtrtri)("U", "N", &sz, REAL(ans), &n, &info
    > FCONE FCONE);
    > } else {
    > F77_CALL(dpotri)("U", &sz, REAL(ans), &n, &info FCONE);
    > }
    > if (only_diag) {
    > SEXP d = PROTECT(allocVector(REALSXP, sz)); nprot++;
    > BLAS_INT inc = (BLAS_INT) n;
    > for (int i = 0; i < sz; i++) {
    > BLAS_INT len = (BLAS_INT) (sz - i);
    > double *x = REAL(ans) + i + (size_t) n* i;
    > REAL(d)[i] = F77_CALL(ddot)(&len, x, &inc, x, &inc);
    > }
    > ans = d;
    > } else {
    > for (int j = 0; j < sz; j++)
    > for (int i = j + 1; i < sz; i++)
    > REAL(ans)[i + (size_t) n * j] = REAL(ans)[j + (size_t) n * i];
    > }

    > Ravi

Thank you, Ravi, for mentioning the  size = <n> option, also in conjunction with the newly proposed  `diag.only = TRUE` option (which I had already added to the R sources ("R-devel") a couple of hours before your e-mail came in [r90579 | maechler |
2026-09-21 18:08:09 +0200].

But I've checked now somewhat systematically, and not seen any reason to change, as it seems the code we got from Adam is also working fine for  size < ncol(x).
(and now actually I'm even guessing your code would *not* work correctly  in all cases .. or if it would it would allocate more memory  than necessary (in the case size < n)).

I've committed these checks as equality assertions using
 stopifnot() at the end of  help(chol2inv) now,  i.e., file src/library/base/man/chol2inv.Rd  in svn rev 90588.

Martin


  [.........................]
  [.........................]