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 would expect the optimized numeric library (LAPACK) performs better than the simple implementation of svd in PDL::MatrixOps. Documentation patches and better examples are welcome for PDL::LinearAlgebra and others. --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
