|
|
Ublas : |
From: Dima Sorkin (dsorkin_at_[hidden])
Date: 2008-03-15 09:40:28
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
>