Shot in the dark: Anyone ever use CLAPACK routines?
bruce.labitt at autoliv.com
bruce.labitt at autoliv.com
Wed May 19 14:34:19 EDT 2010
gnhlug-discuss-bounces at mail.gnhlug.org wrote on 05/19/2010 01:36:03 PM:
> gnhlug-discuss-bounces at mail.gnhlug.org wrote on 05/19/2010 09:19:37 AM:
> > On 05/18/2010 08:13 PM, Bruce Labitt wrote:
> > > As the subject line indicates - a total shot in the dark...
> > >
> > > Prototyping Platform: Ubuntu 10.04 x86-64
> > > Libraries: BLAS from ATLAS, CLAPACK
> > >
> > > I'm trying to use some CLAPACK routines to perform matrix
> manipulation,
> > > in particular, the zgesvd routine to do a singular value
decomposition
>
> > > (SVD). My code is working for a 2x2 matrix, but it does not work
for
> a
> > > 9x9. I posted the code at the lapack-forum
> > >
> > > http://icl.cs.utk.edu/lapack-forum/viewtopic.php?
> > f=2&t=1839&sid=a44c7f5bb3f4836d77568664db0e1c89
> > >
> > >
> > > which works for a 2x2 and fails for a 9x9, with a Segmentation
Fault.
> > >
> > > I'm suspicious that it is 99% operator (me) error. (Fair guess :-P
)
> > >
> > > In particular, I'm worried about stuff like declaring:
> > >
> > > doublecomplex A[m][m];
> > > where doublecomplex is defined in f2c.h as struct{ double r; double
i;
> }
> > >
> > > Is is better in general (more portable) to use something like
> > >
> > > doublecomplex A[m*m] instead?
> > >
> > > For those who may not know, CLAPACK is a C version of LAPACK, which
> was
> > > originally written in FORTRAN (gasp). It is the Linear Algebra
> library
> > > that both OSS and closed source use. I know that Numpy & Scipy use
> > > LAPACK, as well as MATLAB. I'm using CLAPACK because it can be
built
> > > entirely in C. (FORTRAN is not available on the 'final' platform)
> > >
> > > If anyone has a few spare moments, I'd appreciate a quick look and
any
>
> > > helpful comments you may have. FWIW, I used valgrind and saw that
> even
> > > when I got the correct answer, there were tons of warnings and
errors
> > > reported. (These errors were DEEP inside of the CLAPACK library.)
> > >
> > > Note: if anyone is adventurous enough to try this at home (or
anywhere
>
> > > else) you will need to change the libs and includes to point to your
> > > blas and lapack libs. I built my BLAS using ATLAS with the nof77
> option
> > > (C only!), and linked to CLAPACK (which I also built). Most folks
> would
> > > just install the libs from their repositories...
> > >
> > >
> >
> > <Jerry replied>
> > I've used CLAPACK, but as an underlying library for other things. My
> > company uses CLAPACK 3.0 in our product and I am unaware of any
issues.
> > We do complex things like Stochastic systems. In addition about 99% of
> > our code is pure C++ with a very small amount of C Flex and Bison.
> > Additionally, the Python Numpy and Scipy packages use CLAPACK, and we
do
> > use a bit of the lapack stuff in our Python modules. After looking at
> > your code, the first thing I looked for is initializations, and it
does
> > appear that your initializations of matricies A and I are ok. I also
> > think that defining 'doublecomplex A[m][m]' is correct. I'm wondering
if
> > your LDx values are correct. What happens if you define 'doublecomplex
> > A[m+1][m+1];' as well as I and C.
>
> > Jerry Feldman <gaf at blu.org>
> > Boston Linux and Unix
> > PGP key id: 537C5846
> > PGP Key fingerprint: 3D1B 8377 A3C0 A5F2 ECBB CA3B 4607 4319 537C
5846
>
======================================================================================
>
> Jerry, thank you for responding to my shot in the dark. I'm going to
edit
> your response and put it in order. The subject is arcane enough without
> trying to unravel top posting. Indeed, I'm having trouble with the
subject
> matter itself. Please indulge me this time.
>
>
======================================================================================
> <Jerry asked>
> "I'm wondering if your LDx values are correct."
>
> <Bruce responded>
> Here is the documentation on zgesvd I found, which tell one how to set
up
> LDx (and everything else:
>
> /* Fortran documentation from zgesvd.c */
> /* ZGESVD computes the singular value decomposition (SVD) of a complex
*/
> /* M-by-N matrix A, optionally computing the left and/or right singular
> */
> /* vectors. The SVD is written */
>
> /* A = U * SIGMA * conjugate-transpose(V) */
>
> /* where SIGMA is an M-by-N matrix which is zero except for its */
> /* min(m,n) diagonal elements, U is an M-by-M unitary matrix, and */
> /* V is an N-by-N unitary matrix. The diagonal elements of SIGMA */
> /* are the singular values of A; they are real and non-negative, and */
> /* are returned in descending order. The first min(m,n) columns of */
> /* U and V are the left and right singular vectors of A. */
>
> /* Note that the routine returns V**H, not V. */
>
> /* Arguments */
> /* ========= */
>
> /* JOBU (input) CHARACTER*1 */
> /* Specifies options for computing all or part of the matrix U:
> */
> /* = 'A': all M columns of U are returned in array U: */
> /* = 'S': the first min(m,n) columns of U (the left singular
*/
> /* vectors) are returned in the array U; */
> /* = 'O': the first min(m,n) columns of U (the left singular
*/
> /* vectors) are overwritten on the array A; */
> /* = 'N': no columns of U (no left singular vectors) are */
> /* computed. */
>
> /* JOBVT (input) CHARACTER*1 */
> /* Specifies options for computing all or part of the matrix */
> /* V**H: */
> /* = 'A': all N rows of V**H are returned in the array VT; */
> /* = 'S': the first min(m,n) rows of V**H (the right singular
*/
> /* vectors) are returned in the array VT; */
> /* = 'O': the first min(m,n) rows of V**H (the right singular
*/
> /* vectors) are overwritten on the array A; */
> /* = 'N': no rows of V**H (no right singular vectors) are */
> /* computed. */
>
> /* JOBVT and JOBU cannot both be 'O'. */
>
> /* M (input) INTEGER */
> /* The number of rows of the input matrix A. M >= 0. */
>
> /* N (input) INTEGER */
> /* The number of columns of the input matrix A. N >= 0. */
>
> /* A (input/output) COMPLEX*16 array, dimension (Lif ( mysvd(
JOBU,
> JOBVT, m, n, (doublecomplex)*C, lda, *S, (doublecomplex)*UU, ldu,
> (doublecomplex)*VV, ldvt, (doublecomplex)*wk, lwork) == 0)DA,N) */
> /* On entry, the M-by-N matrix A. */
> /* On exit, */
> /* if JOBU = 'O', A is overwritten with the first min(m,n) */
> /* columns of U (the left singular vectors, */
> /* stored columnwise); */
> /* if JOBVT = 'O', A is overwritten with the first min(m,n) */
> /* rows of V**H (the right singular vectors, */
> /* stored rowwise); */
> /* if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A */
> /* are destroyed. */
>
> /* LDA (input) INTEGER */
> /* The leading dimension of the array A. LDA >= max(1,M). */
>
> /* S (output) DOUBLE PRECISION array, dimension (min(M,N)) */
> /* The singular values of A, sorted so that S(i) >= S(i+1). */
>
> /* U (output) COMPLEX*16 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 unitary 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. */
>
> /* VT (output) COMPLEX*16 array, dimension (LDVT,N) */
> /* If JOBVT = 'A', VT contains the N-by-N unitary matrix */
> /* V**H; */
> /* if JOBVT = 'S', VT contains the first min(m,n) rows of */
> /* V**H (the right singular vectors, stored rowwise); */
> /* if JOBVT = 'N' or 'O', VT is not referenced. */
>
> /* LDVT (input) INTEGER */
> /* The leading dimension of the array VT. LDVT >= 1; if */
> /* JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N). */
>
> /* WORK (workspace/output) COMPLEX*16 array, dimension
(MAX(1,LWORK))
> */
> /* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
>
> /* LWORK (input) INTEGER */
> /* The dimension of the array WORK. */
> /* LWORK >= MAX(1,2*MIN(M,N)+MAX(M,N)). */
> /* For good performance, LWORK should generally be larger. */
>
> /* If LWORK = -1, then a workspace query is assumed; the
routine
> */
> /* only calculates the optimal size of the WORK array, returns
*/
> /* this value as the first entry of the WORK array, and no
error
> */
> /* message related to LWORK is issued by XERBLA. */
>
> /* RWORK (workspace) DOUBLE PRECISION array, dimension (5*min(M,N))
*/
> /* On exit, if INFO > 0, RWORK(1:MIN(M,N)-1) contains the */
> /* unconverged superdiagonal elements of an upper bidiagonal */
> /* matrix B whose diagonal is in S (not necessarily sorted). */
> /* B satisfies A = U * B * VT, so it has the same singular */
> /* values as A, and singular vectors related by U and VT. */
>
> /* INFO (output) INTEGER */
> /* = 0: successful exit. */
> /* < 0: if INFO = -i, the i-th argument had an illegal value.
*/
> /* > 0: if ZBDSQR did not converge, INFO specifies how many */
> /* superdiagonals of an intermediate bidiagonal form B */
> /* did not converge to zero. See the description of RWORK
> */
> /* above for details. */
>
> Since I'm using an mxm matrix, LDA = LDU = LDVT = m, right? That is how
I
> read it.
>
> <then Jerry commented>
> "What happens if you define 'doublecomplex A[m+1][m+1];' as well as I
and
> C.
>
> <finally, Bruce responds to Jerry's comment>
> I have not tried that, though I'm doubtful... (based on the above
> documentation.)
> Nonetheless, I'll try it.
>
> I just wrote a script in Numpy/Scipy which uses the same data set which
I
> posted at the LAPACK-forum. It solves the svd problem for both the 2x2
> and 9x9 case. So the dataset itself is not ill-behaved. Jeesh, I like
> python sooooooo much better. Too bad this HAS to be C.
>
> Any comments are welcome.
>
> -Bruce
>
Comment on my own post. Probably nothing is wrong with MY CODE.
http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=2&t=745
The claim is that zgesvd_ function hangs. The fix is allegedly to modify
the make.inc file for CLAPACK.
The change is adding -ffloat-store to the end of NOOPT.
"NOOPT = -O0 -I$(TOPDIR)/INCLUDE -fPIC -ffloat-store"
Hope that fixes it. I'll find out in 4+ hours!
-Bruce
******************************
Neither the footer nor anything else in this E-mail is intended to or constitutes an <br>electronic signature and/or legally binding agreement in the absence of an <br>express statement or Autoliv policy and/or procedure to the contrary.<br>This E-mail and any attachments hereto are Autoliv property and may contain legally <br>privileged, confidential and/or proprietary information.<br>The recipient of this E-mail is prohibited from distributing, copying, forwarding or in any way <br>disseminating any material contained within this E-mail without prior written <br>permission from the author. If you receive this E-mail in error, please <br>immediately notify the author and delete this E-mail. Autoliv disclaims all <br>responsibility and liability for the consequences of any person who fails to <br>abide by the terms herein. <br>
******************************
More information about the gnhlug-discuss
mailing list