I will try to answer. I hope one of our Fortran guru will back me up if I
Well I think that if you enter the routine with LDU=0 then the Fortran
DOUBLE PRECISION U(LDU,*)
is not Fortran-correct. So the descrition of LDU says that LDU needs to be
LDU > 0, no matter whether U is referenced or not. So since the
description says it, we check it.
It can be a little bit weird because your code can either crash or run
correctly. This will depend on your compiler. If it crashes then so
be it. If it runs correctly then why do we stop the execution?
The philosophy is the following: If your code runs correctly with a
statement like DOUBLE PRECISION U(0,*) we prefer to exit the routine to
tell you that this is not Fortran-correct.
This avoids users to have problems when they change compilers and this
does not cost much.
On Sat, 15 Mar 2008, Nico Galoppo wrote:
there seems to be an odd requirement in dgesvd. If JOBU = 'O' (overwrite A
with the left eigenvectors), then U is never referenced, but there is a
check in the code that checks for LDU>=1 and aborts if it's not (see code
excerpts below). Can anyone explain? Is this a bug?
DGESVD subroutine is in http://www.netlib.org/lapack/double/dgesvd.f
SUBROUTINE DGESVD( JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT,
$ WORK, LWORK, INFO )
* .. Scalar Arguments ..
CHARACTER JOBU, JOBVT
INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N
* .. Array Arguments ..
DOUBLE PRECISION A( LDA, * ), S( * ), U( LDU, * ),
$ VT( LDVT, * ), WORK( * )
Description of the arguments says:
* U (output) DOUBLE PRECISION array, dimension (LDU,UCOL)
* (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'.
* If JOBU = 'A', U contains the M-by-M orthogonal matrix U;
* if JOBU = 'S', U contains the first min(m,n) columns of U
* (the left singular vectors, stored columnwise);
* if JOBU = 'N' or 'O', U is not referenced.
* LDU (input) INTEGER
* The leading dimension of the array U. LDU >= 1; if
* JOBU = 'S' or 'A', LDU >= M.
And there is a check later in the code:
ELSE IF( LDU.LT.1 .OR. ( WNTUAS .AND. LDU.LT.M ) ) THEN
INFO = -9