[Rd] Proposal for diagonal-only chol2inv option

[Martin Maechler]

>>>>> 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
 

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