|
|
Ublas : |
From: Karl Meerbergen (Karl.Meerbergen_at_[hidden])
Date: 2006-01-26 09:36:20
All I can say is that the bindings are written for column major matrices.
Your matrix is row_major. I do not believe it is sufficient to transpose it,
since for a row_major matrix the leading dimension corresponds to the number
of columns. Lapack assumes it corresponds to the number of rows.
You could pass trans(M) as an argument to lapack. trans(M) creates an object
of type matrix_unary2<M,scalar_identity<typename M::value_type> >. There are
no bindings for this kind of object. These should be written.
This would be the clean way of doing it.
Karl
On Thursday 26 January 2006 15:17, Peter Melchior wrote:
> Hello,
>
> I want to do a SVD for non-square real matrices.
> When I have a Matrix with more columns than rows, I get this runtime
> error:
> ** On entry to DGESVD parameter number 6 had an illegal value
>
> When I have a matrix with more rows than columns I get a wrong result
> (checked against Mathematica and IDL).
> But when using a square matrix, everything works fine.
> I'm definitely sure, that the U, S and V matrix have the correct sizes.
>
> This is the code I use:
>
> template <class T>
> inline void svd(boost::numeric::ublas::matrix<T>& M,
> boost::numeric::ublas::matrix<T>& U, boost::numeric::ublas::matrix<T>&
> S, boost::numeric::ublas::matrix<T>& Vt) {
> // dimensions of matrices
> int n = matrix.size1();
> int m = matrix.size2();
> int length = min(n,m);
>
> // since LAPACK works on column-major matrices, transpose it.
> boost::numeric::ublas::matrix<T> transpose =
> boost::numeric::ublas::trans(M);
> boost::numeric::ublas::vector<T> s(length);
> boost::numeric::bindings::lapack::gesvd('A','A',transpose,s,U,Vt);
>
> // transpose U and Vt
> U = boost::numeric::ublas::trans(U);
> Vt = boost::numeric::ublas::trans(Vt);
>
> // addressing the diagonal elements of S to store singular values
>
> boost::numeric::ublas::matrix_vector_range<boost::numeric::ublas::matrix<T>
> > (S, boost::numeric::ublas::range (0, length),
> boost::numeric::ublas::range (0, length)) = s; }
>
>
>
> Does anybody have a clue what is going on here?
>
>
> Best regards,
-- ============================================== Look at our unique training program and Register on-line at http://www.fft.be/?id=35 ---------------------------------------------- Karl Meerbergen Free Field Technologies 16 place de l'Université B-1348 Louvain-la-Neuve - BELGIUM Company Phone: +32 10 45 12 26 Company Fax: +32 10 45 46 26 Mobile Phone: +32 474 26 66 59 Home Phone: +32 2 306 38 10 mailto:Karl.Meerbergen_at_[hidden] http://www.fft.be ========================================================================================== NEW ADDRESS FROM FEBRUARY 1st ONWARD: Axis Park Louvain-la-Neuve rue Emile Francqui, 1 B-1435 Mont-Saint Guibert - BELGIUM Same telephone, same fax, same mail ============================================