Shot in the dark: Anyone ever use CLAPACK routines?
bruce.labitt at autoliv.com
bruce.labitt at autoliv.com
Wed May 19 13:36:03 EDT 2010
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
******************************
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