|
|
Ublas : |
From: Nico Galoppo (ngaloppo_at_[hidden])
Date: 2008-03-22 00:15:11
I got an answer on the lapack list. Apparently, it's a convention issue:
Hello,
I will try to answer. I hope one of our Fortran guru will back me up if I
am wrong.
Well I think that if you enter the routine with LDU=0 then the Fortran
statement:
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.
Julien.
2008/3/15 Dima Sorkin <dsorkin_at_[hidden]>:
> Hi.
> In such cases I usually post the question also to
> lapack_at_[hidden]
> They typically give very helpful responses.
> Maybe such checking is Lapack's bug ...
>
> Regards,
> Dima.
>
> ãÉÔÉÒÕÀ Nico Galoppo <ngaloppo_at_[hidden]>:
>
>
>
> > Fres,
> >
> > thanks for your explanation. However, I am confused with the way I should
> > call the bindings then. Why is it necessary to allocate memory for U, if it
> > is not going to be referenced (as stated in the comments of the fortran
> > code).
> >
> > Currently, I've 'circumvented' that requirement by calling gesvd as such:
> >
> > gesvd('O', 'O', 'S', A, s, A, VT);
> >
> > However, I am not sure if this is safe to do.
> >
> > --nico
> >
> > On Sat, Mar 15, 2008 at 4:03 AM, Kresimir Fresl <fresl_at_[hidden]> wrote:
> >
> > > Nico Galoppo wrote:
> > >
> > > > Upon calling the lapack gesvd in numeric/bindings/lapack/gesvd.cpp,
> > > > with JobU = 'O', I'm getting an assertion error on line 418:
> > > >
> > > > 417 assert ((jobu == 'N' && traits::leading_dimension (u) >= 1)
> > > > 418 || (jobu == 'O' && traits::leading_dimension (u) >=
> > > 1)
> > > > 419 || (jobu == 'A' && traits::leading_dimension (u) >=
> > > m)
> > > > 420 || (jobu == 'S' && traits::leading_dimension (u) >=
> > > m));
> > > >
> > > > I don't understand why this assertion is correct. jobu == 'O' means
> > > > that the left-hand vectors should be written to input matrix A, not to
> > > > u, hence, the u input to gesvd should essentially be able to be a null
> > > > pointer (or zero-size matrix).
> > > >
> > > > Am I missing something?
> > >
> > > This is the requirement in the "original" LAPACK xgesvd.
> > >
> > > DGESVD subroutine is in http://www.netlib.org/lapack/double/dgesvd.f
> > > declared as:
> > >
> > > 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
> > >
> > > Regards,
> > >
> > > fres
> > >
> > > _______________________________________________
> > > ublas mailing list
> > > ublas_at_[hidden]
> > > http://lists.boost.org/mailman/listinfo.cgi/ublas
> > >
> >
> >
> >
> > --
> > Nico Galoppo :: http://www.ngaloppo.org
> >
>
>
>
-- Nico Galoppo :: http://www.ngaloppo.org