Shot in the dark: Anyone ever use CLAPACK routines?
Jerry Feldman
gaf at blu.org
Wed May 19 15:36:20 EDT 2010
On 05/19/2010 02:34 PM, bruce.labitt at autoliv.com wrote:
> 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!
>
So do I. I'm dubious about using -ffloat-store, but if it works, then
that's fine. I took a look at our libraries, and we use libf77.a and
libif77.a.
--
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
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 253 bytes
Desc: OpenPGP digital signature
Url : http://mail.gnhlug.org/mailman/private/gnhlug-discuss/attachments/20100519/1adeab3e/attachment-0001.bin
More information about the gnhlug-discuss
mailing list