Thanks Chris, On Sun, May 21, 2017 at 01:23:29PM -0400, Chris Marshall wrote: > PDL::LinearAlgebra is bindings to LAPACK which > is a fortran library for numeric computation. > > I would guess that the transpose is because fortran > uses column major ordering while PDL uses row major > ordering.
I'm starting to understand this. The fastest dimension is the first, both in Fortran and in PDL, but PDL interprets it as the column dimension for matrices and fortran as the row dimension. Thus, I guess all calls to the PDL::LinearAlgebra package would require transposition of all matrices. Would this be correct? I made a few tests for this idea and they seemed to work. > > I would expect the optimized numeric library (LAPACK) > performs better than the simple implementation of svd > in PDL::MatrixOps. Much better. For my particular calculations, with 600x300 matrices, I got a ~200 times speedup. Using gesdd I get a further 40% increase in speed. > > Documentation patches and better examples are > welcome for PDL::LinearAlgebra and others. Got it. I guess a short warning about column vs. row major for each routine would be adequate. An alternative would be to wrap the library routines with perl routines that perform the transpositions. As lapack expects that consecutive data along columns are adjacent, I guess the transpositions should be done first and then the data copied to the arrays that would actually be fed to lapack. This is what the [phys] flag in the signature does, right? Best regards, Luis > > --Chris > > On 5/21/2017 00:57, Luis Mochan wrote: > > I'm making a calculation using singular value decomposition SVD to > > solve 2N equations with N unknowns, with N= a few hundreds, where I > > know the equations are not all linearly independent and I expect the > > system to have a unique well defined solution . I noticed that the > > time taken to compute the solution using the routine 'svd' from > > PDL::MatrixOps seems to scale with the size of the problem to the > > fourth power. Is this correct? Anyway, I realized there are much faster > > routines, such as gesvd, available through PDL::LinearAlgebra::Real, which > > interfaces > > to lapack/blas. Nevertheless, I found some strange behavior. It seems > > gesvd cannot create its output parameters so I have to pass them > > explicitly. Furthermore, if I pass a > > rectangular matrix with more rows than columns I get an error > > *** Error in `perl': double free or corruption (!prev): > > 0x000055d3420976f0 *** > > or a bus error or a segmentation fault. > > Other aproaches I tried produce the error > > ** On entry to DGESVD parameter number 9 (or 11) had an illegal value > > Nevertheless, if I transpose some of the arguments it seems I get the > > correct result without changing the order of the matrices. The program > > below should illustrate these points: > > > > ################# > > use feature 'say'; > > use PDL; > > use PDL::NiceSlice; > > use PDL::LinearAlgebra::Real; > > $N=5; > > $a=sequence($N,2*$N); > > $u=zeroes($N,2*$N); > > $v=zeroes($N,$N); > > $s=zeroes($N); > > $i=pdl(long,0); > > > > #The slow solution: > > ($u,$s,$v)=svd($a); > > say $u x ($s->(*1,:)*$v->transpose); > > #-> [ > > # [0,1,2,3,4] > > # [5,6,7,8,9] > > # ... > > # [45,46,47,48,49] > > # ] > > > > #gesvd($a,2,2, $s,$u,$v,$i); produces > > #*** Error in `perl': double free or corruption (!prev): 0x000055d3420976f0 > > *** > > > > #My fast solution: > > gesvd($a->transpose,2,2, $s,$u->transpose,$v,$i); > > say $u x ($s->(*1,:) * $v->transpose); > > #-> [ > > # [0,1,2,3,4] > > # [5,6,7,8,9] > > # ... > > # [45,46,47,48,49] > > # ] > > ##################### > > > > Calling gesvd this way seems to run much faster than svd and produces > > the same results, so far, but I would like to understand if this is the > > correct way to use this routine. > > > > Best regards, > > Luis > > > > > ------------------------------------------------------------------------------ > Check out the vibrant tech community on one of the world's most > engaging tech sites, Slashdot.org! http://sdm.link/slashdot > _______________________________________________ > pdl-general mailing list > [email protected] > https://lists.sourceforge.net/lists/listinfo/pdl-general > -- o W. Luis Mochán, | tel:(52)(777)329-1734 /<(*) Instituto de Ciencias Físicas, UNAM | fax:(52)(777)317-5388 `>/ /\ Apdo. Postal 48-3, 62251 | (*)/\/ \ Cuernavaca, Morelos, México | [email protected] /\_/\__/ GPG: 791EB9EB, C949 3F81 6D9B 1191 9A16 C2DF 5F0A C52B 791E B9EB ------------------------------------------------------------------------------ Check out the vibrant tech community on one of the world's most engaging tech sites, Slashdot.org! http://sdm.link/slashdot _______________________________________________ pdl-general mailing list [email protected] https://lists.sourceforge.net/lists/listinfo/pdl-general
