Possible bug in many of the QR/LQ LAPACK calls
Hi, thanks for putting out this software. I've been using it to solve high dimensional PDEs this summer. It really is an excellent package. Your paper on the topic is also extremely informative.
I think that there may be an error which is repeated many times throughout the code. It has to do with the block sizes used when calling the QR and LQ factorizations of LAPACK. For example, line 362 of attac.c reads
dgeqrt_(&n[0], &rx[1], &rx[1], x[0], &n[0], tau, &rx[1], work, &ierr);
In an application where the rank is low, it is typical that rx[1]
is smaller than n[0]
. When adding a list of tensors with the formal sum routine, before a compression this may be violated. That makes the 3rd argument of the LAPACK dgeqrt_()
call incorrect, as LAPACK expects a block size no large than the minimum of rows and columns. This leads to a crash on after the stderr
printout. I believe this can be fixed by introducing a new variable into every function which calls a LAPACK function with a blocksize argument. Then set blocksize to be the minimum of rows and columns.
A related error may be propagated into the the various calls to dgeqrt_()
and dgelqt_()
, but it is a bit harder to trace. I think a few of the arguments here would have to be changed to this proposed new blocksize integer.