File Coverage

blib/lib/PDL/SVDLIBC.pm
Criterion Covered Total %
statement 40 40 100.0
branch 19 28 67.8
condition n/a
subroutine 9 9 100.0
pod 4 5 80.0
total 72 82 87.8


line stmt bran cond sub pod time code
1              
2             #
3             # GENERATED WITH PDL::PP! Don't modify!
4             #
5             package PDL::SVDLIBC;
6              
7             @EXPORT_OK = qw( PDL::PP _svdccsencode svdlas2a PDL::PP svdlas2 svdlas2aw PDL::PP svdlas2w svdlas2ad PDL::PP svdlas2d PDL::PP svdindexND svdindexNDt PDL::PP svdindexccs PDL::PP svderror );
8             %EXPORT_TAGS = (Func=>[@EXPORT_OK]);
9              
10 2     2   561269 use PDL::Core;
  2         9  
  2         11  
11 2     2   408 use PDL::Exporter;
  2         4  
  2         10  
12 2     2   48 use DynaLoader;
  2         3  
  2         155  
13              
14              
15              
16             $PDL::SVDLIBC::VERSION = 0.19;
17             @ISA = ( 'PDL::Exporter','DynaLoader' );
18             push @PDL::Core::PP, __PACKAGE__;
19             bootstrap PDL::SVDLIBC $VERSION;
20              
21              
22              
23              
24             =pod
25              
26             =head1 NAME
27              
28             PDL::SVDLIBC - PDL interface to Doug Rohde's SVD C Library
29              
30             =head1 SYNOPSIS
31              
32             use PDL;
33             use PDL::SVDLIBC;
34              
35             ##---------------------------------------------------------------------
36             ## Input matrix (dense)
37             ##---------------------------------------------------------------------
38             $n = 100; ##-- number of columns
39             $m = 50; ##-- number of rows
40             $a = random(double,$n,$m); ##-- random matrix
41              
42             ##---------------------------------------------------------------------
43             ## Output pdls
44             ##---------------------------------------------------------------------
45             $d = $n; ##-- max number of output dimensions
46             $ut = zeroes(double,$m,$d); ##-- left singular components
47             $s = zeroes(double,$d); ##-- singular values (diagnonal vector)
48             $vt = zeroes(double,$n,$d); ##-- right singular components
49              
50             ##---------------------------------------------------------------------
51             ## Singular Value Decomposition (dense)
52             ##---------------------------------------------------------------------
53             svdlas2d($a, $maxiters, $end, $kappa, $ut, $s, $vt);
54              
55             ##---------------------------------------------------------------------
56             ## Singular Value Decomposition (sparse, using direct whichND()-encoding)
57             ##---------------------------------------------------------------------
58             $which = whichND($a)->qsortvec();
59             $nzvals = indexND($a,$which);
60              
61             svdlas2w($which, $nzvals, $n, $m, $maxiters, $end, $kappa, $ut, $s, $vt);
62              
63             ##---------------------------------------------------------------------
64             ## Singular Value Decomposition (sparse, using PDL::CCS encoding)
65             ##---------------------------------------------------------------------
66             use PDL::CCS;
67             ($ptr,$rowids,$nzvals) = ccsencode($a);
68             $ptr->reshape($ptr->nelem+1);
69             $ptr->set(-1, $rowids->nelem);
70              
71             svdlas2($ptr, $rowids, $nzvals, $m, $maxiters, $end, $kappa, $ut, $s, $vt);
72              
73             ##---------------------------------------------------------------------
74             ## SVD decoding (lookup)
75             ##---------------------------------------------------------------------
76             $vals = svdindexND ($u, $s, $v, $which);
77             $vals = svdindexNDt($ut,$s,$vt, $which);
78             $vals = svdindexccs($u, $s, $v, $ptr,$rowids);
79             $err = svderror ($u, $s, $v, $ptr,$rowids,$nzvals);
80              
81             =head1 DESCRIPTION
82              
83             PDL::SVDLIBC provides a PDL interface to the SVDLIBC routines
84             for singular value decomposition of large sparse matrices.
85             SVDLIBC is available from http://tedlab.mit.edu/~dr/SVDLIBC/
86              
87             =cut
88              
89              
90              
91              
92              
93              
94              
95             =head1 FUNCTIONS
96              
97              
98              
99             =cut
100              
101              
102              
103              
104              
105 2     2   10 use strict;
  2         4  
  2         1206  
106              
107             =pod
108              
109             =head1 SVDLIBC Globals
110              
111             There are several global data structures still lurking in the
112             SVDLIBC code, so expect problems if you are trying to run more
113             than one 'las2' procedure at once (even in different processes).
114              
115             PDL::SVDLIBC provides access to (some of) the SVDLIBC globals
116             through the following functions, which are not exported.
117              
118             =cut
119              
120              
121              
122             =pod
123              
124             =head2 PDL::SVDLIBC::verbosity()
125              
126             =head2 PDL::SVDLIBC::verbosity($level)
127              
128             Get/set the current SVDLIBC verbosity level.
129             Valid values for $level are between 0 (no messages) and
130             2 (many messages).
131              
132             =cut
133              
134              
135              
136              
137             =pod
138              
139             =head2 PDL::SVDLIBC::svdVersion()
140              
141             Returns a string representing the SVDLIBC version
142             this module was compiled with.
143              
144             =cut
145              
146              
147              
148              
149             =pod
150              
151             =head1 SVD Utilities
152              
153             =cut
154              
155              
156              
157              
158              
159             =head2 _svdccsencode
160              
161             =for sig
162              
163             Signature: (double a(n,m); indx [o]ptr(n1); indx [o]rowids(nnz); double [o]nzvals(nnz))
164              
165              
166             =for ref
167              
168             info not available
169              
170              
171             =for bad
172              
173             _svdccsencode does not process bad values.
174             It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.
175              
176              
177             =cut
178              
179              
180              
181              
182              
183              
184             *_svdccsencode = \&PDL::_svdccsencode;
185              
186              
187              
188              
189             =pod
190              
191             =head2 svdlas2a
192              
193             =for sig
194              
195             indx ptr(nplus1);
196             indx rowids(nnz);
197             double nzvals(nnz);
198             indx m(); ##-- default: max($rowids)+1
199             int d(); ##-- default: max(nplus1-1,m)
200             int iterations(); ##-- default: 2*$d
201             double end(2); ##-- default: [-1e-30,1e-30]
202             double kappa(); ##-- default: 1e-6
203             double [o]ut(m,d); ##-- default: new
204             double [o] s(d); ##-- default: new
205             double [o]vt(n,d); ##-- default: new
206              
207             Uses a variant of the single-vector Lanczos method (Lanczos, 1950)
208             to compute the singular value decomposition of a sparse matrix with
209             $m() rows and data encoded
210             in Harwell-Boeing sparse format in the input parameters $ptr(), $rowids(),
211             and $nzvals(). See L<"PDL::CCS"> for a way to acquire these parameters
212             from a dense input matrix, but note that for svdlas2(), the
213             column pointer $ptr() is of size ($n+1) for a dense matrix $a with
214             $n columns, where $ptr($n)==$nnz is the total number of nonzero
215             values in $a.
216              
217             $iterations() is the maximum number of Lanczos iterations to perform.
218              
219             $end() specifies two endpoints of an interval within which all unwanted
220             eigenvalues lie.
221              
222             $kappa() is a double containing the relative accuracy of Ritz
223             values acceptable as eigenvalues.
224              
225             The left singular components are returned in the matrix $ut(),
226             the singular values themselved in the vector $s(), and the
227             right singular components in the matrix $vt(). Note that
228             $ut() and $vt() are transposed, and must be specified explicitly
229             in the call, so that the degree of reduction (the size parameter $d)
230             can be determined. If $d==$n, then a full decomposition
231             will be computed, and on return, $ut() and $vt() should be transposed
232             instances of the matrices $u() and $v() as returned by PDL::MatrixOps::svd().
233              
234             The Lanczos method as used here seems to be consistently the
235             fastest. This algorithm has the drawback that the low order singular
236             values may be relatively imprecise, but that is not a problem for most
237             users who only want the higher-order values or who can tolerate some
238             imprecision.
239              
240             See also: svdlas2aw(), svdlas2d()
241              
242             =cut
243              
244             ## ($iters,$end,$kappa,$ut,$s,$vt) = svddefaults($n=$nrows,$m=$ncols,$d, $iters,...)
245             ## + returns default values
246             ## + changed calling conventions in v0.14
247             ## - WAS: svddefaults($nrows,$cols, $d,$iters,...) ##-- SVDLIBC-style (col-primary)
248             ## ~ svddefaults($m, $n, $d,$iters,...) ##-- SVDLIBC-style (for dense $a(n,m))
249             ## - NOW: svddefaults($n, $m, $d,$iters,...) ##-- pdl-style
250             sub svddefaults {
251 7     7 0 19 my ($n,$m,$d, $iters,$end,$kappa,$ut,$s,$vt) = @_;
252 7 50       23 $n = $n->at(0) if (UNIVERSAL::isa($n,'PDL'));
253 7 50       23 $m = $m->at(0) if (UNIVERSAL::isa($m,'PDL'));
254 7 50       19 $d = ($n >= $m ? $n : $m) if (!defined($d));
    100          
255 7 50       20 $iters = 2*$d if (!defined($iters));
256 7 50       25 $end = pdl(double,[-1e-30,1e-30]) if (!defined($end));
257 7 50       184 $kappa = pdl(double,1e-6) if (!defined($kappa));
258 7 50       371 $ut = PDL->zeroes(double,$m,$d) if (!defined($ut));
259 7 50       416 $s = PDL->zeroes(double, $d) if (!defined($s));
260 7 50       290 $vt = PDL->zeroes(double,$n,$d) if (!defined($vt));
261 7         311 return ($iters,$end,$kappa,$ut,$s,$vt);
262             }
263              
264             sub svdlas2a {
265 2     2 1 113795 my ($ptr,$rowids,$nzvals, $m,$d, @args) = @_;
266 2 100       12 $m = $rowids->flat->max+1 if (!defined($m));
267 2         87 @args = svddefaults($ptr->dim(0)-1,$m,$d,@args);
268 2         128 svdlas2($ptr,$rowids,$nzvals,$m,@args);
269 2         22 return @args[3..5];
270             }
271              
272              
273              
274              
275              
276             =head2 svdlas2
277              
278             =for sig
279              
280             Signature: (
281             indx ptr(nplus1);
282             indx rowids(nnz);
283             double nzvals(nnz);
284             indx m();
285             int iterations();
286             double end(2);
287             double kappa();
288             double [o]ut(m,d);
289             double [o] s(d);
290             double [o]vt(n,d);
291             )
292              
293              
294             Guts for svdlas2a().
295             No default instantiation, and slightly different calling conventions.
296              
297              
298             =for bad
299              
300             svdlas2 does not process bad values.
301             It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.
302              
303              
304             =cut
305              
306              
307              
308              
309              
310              
311             *svdlas2 = \&PDL::svdlas2;
312              
313              
314              
315              
316             =pod
317              
318             =head2 svdlas2aw
319              
320             =for sig
321              
322             indx which(nnz,2); ##-- sorted indices of non-zero values
323             double nzvals(nnz); ##-- non-zero values
324             indx n(); ##-- default: max($indx(0,:))+1
325             indx m(); ##-- default: max($indx(1,:))+1
326             int d(); ##-- default: max(n,m)
327             int iterations(); ##-- default: 2*$d
328             double end(2); ##-- default: [-1e-30,1e-30]
329             double kappa(); ##-- default: 1e-6
330             double [o]ut(m,d); ##-- default: new
331             double [o] s(d); ##-- default: new
332             double [o]vt(n,d); ##-- default: new
333              
334             As for svdlas2a(), but implicitly converts the index-encoded matrix
335             ($which(),$nzvals()) to an internal CCS-like sparse format
336             before computing the decomposition.
337             Should be slightly more efficient than using PDL::CCS::ccsencode()
338             or similar if you already have $which() and $nzvals() available.
339             These can be attained for a dense matric $a() e.g. by:
340              
341             $which = $a->whichND->qsortvec->xchg(0,1);
342             $nzvals = $a->indexND($which->xchg(0,1));
343              
344             For convenience, $which() will be implicitly transposed if it is passed
345             as a list-of-vectors C<$whichND(2,nnz)> such as returned by L,
346             but it must still be lexicographically sorted.
347              
348             See also: svdlas2a(), svdlas2d()
349              
350             =cut
351              
352             sub svdlas2aw {
353 3     3 1 5898 my ($which,$nzvals, $n,$m,$d, @args) = @_;
354 3 100       22 $which = $which->xchg(0,1) if ($which->dim(1) > $which->dim(0));
355 3 100       14 $n = $which->slice(":,0")->max+1 if (!defined($n));
356 3 100       135 $m = $which->slice(":,1")->max+1 if (!defined($m));
357 3         92 @args = svddefaults($n,$m,$d,@args);
358 3         225 svdlas2w($which,$nzvals,$n,$m,@args);
359 3         34 return @args[3..5];
360             }
361              
362              
363              
364              
365              
366             =head2 svdlas2w
367              
368             =for sig
369              
370             Signature: (
371             indx whichi(nnz,Two);
372             double nzvals(nnz);
373             indx n();
374             indx m();
375             int iterations();
376             double end(2);
377             double kappa();
378             double [o]ut(m,d);
379             double [o] s(d);
380             double [o]vt(n,d);
381             )
382              
383              
384             Guts for svdlas2a().
385             No default instantiation, and slightly different calling conventions.
386              
387              
388             =for bad
389              
390             svdlas2w does not process bad values.
391             It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.
392              
393              
394             =cut
395              
396              
397              
398              
399              
400              
401             *svdlas2w = \&PDL::svdlas2w;
402              
403              
404              
405              
406             =pod
407              
408             =head2 svdlas2ad
409              
410             =for sig
411              
412             double a(n,m);
413             int d(); ##-- default: max($n,$m)
414             int iterations(); ##-- default: 2*$d
415             double end(2); ##-- default: [-1e-30,1e-30]
416             double kappa(); ##-- default: 1e-6
417             double [o]ut(m,d); ##-- default: new
418             double [o] s(d); ##-- default: new
419             double [o]vt(n,d); ##-- default: new
420              
421             As for svdlas2(), but implicitly converts the dense input matrix
422             $a() to sparse format before computing the decomposition.
423              
424             =cut
425              
426             sub svdlas2ad {
427 2     2 1 3534 my ($a,$d, @args) = @_;
428 2         12 @args = svddefaults($a->dim(0),$a->dim(1),$d,@args);
429 2         117 svdlas2d($a,@args);
430 2         17 return @args[3..5];
431             }
432              
433              
434              
435              
436              
437             =head2 svdlas2d
438              
439             =for sig
440              
441             Signature: (
442             double a(n,m);
443             int iterations();
444             double end(2);
445             double kappa();
446             double [o]ut(m,d);
447             double [o] s(d);
448             double [o]vt(n,d);
449             )
450              
451              
452             Guts for _svdlas2d().
453              
454              
455             =for bad
456              
457             svdlas2d does not process bad values.
458             It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.
459              
460              
461             =cut
462              
463              
464              
465              
466              
467              
468             *svdlas2d = \&PDL::svdlas2d;
469              
470              
471              
472              
473              
474             =head2 svdindexND
475              
476             =for sig
477              
478             Signature: (
479             u(d,m);
480             s(d);
481             v(d,n);
482             indx which(Two,nnz);
483             [o] vals(nnz);
484             )
485              
486              
487             Lookup selected values in an SVD-encoded matrix, L-style.
488             Should be equivalent to:
489              
490             ($u x stretcher($s) x $v->xchg(0,1))->indexND($which)
491              
492             or its PDL-friendlier variant:
493              
494             ($u * $s)->matmult($v->xchg(0,1))->indexND($which)
495              
496             ... but only computes the specified values $which(), avoiding
497             memory bottlenecks for large sparse matrices.
498             This is a pure PDL::PP method, so you can use e.g.
499             C for the SVD-encoded matrix if you wish.
500              
501              
502              
503             =for bad
504              
505             svdindexND does not process bad values.
506             It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.
507              
508              
509             =cut
510              
511              
512              
513              
514              
515              
516             *svdindexND = \&PDL::svdindexND;
517              
518              
519              
520              
521             =pod
522              
523             =head2 svdindexNDt
524              
525             =for sig
526              
527             ut(m,d); s(d); vt(n,d); indx which(Two,nnz); [o] vals(nnz);
528              
529             Wrapper for L accepting transposed singular components
530             $ut() and $vt() as returned by e.g. L.
531              
532             =cut
533              
534             sub svdindexNDt {
535 2     2 1 2720 return svdindexND($_[0]->xchg(0,1),$_[1],$_[2]->xchg(0,1),@_[3..$#_]);
536             }
537              
538              
539              
540              
541              
542             =head2 svdindexccs
543              
544             =for sig
545              
546             Signature: (
547             u(d,m);
548             s(d);
549             v(d,n);
550             indx ptr(nplus1);
551             indx rowids(nnz);
552             [o] vals(nnz);
553             )
554              
555              
556             Lookup selected values in an SVD-encoded matrix using L-style indexing
557             as for L.
558              
559              
560              
561             =for bad
562              
563             svdindexccs does not process bad values.
564             It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.
565              
566              
567             =cut
568              
569              
570              
571              
572              
573              
574             *svdindexccs = \&PDL::svdindexccs;
575              
576              
577              
578              
579              
580             =head2 svderror
581              
582             =for sig
583              
584             Signature: (
585             u(d,m);
586             s(d);
587             v(d,n);
588             indx ptr(nplus1);
589             indx rowids(nnz);
590             nzvals(nnz);
591             [o]err();
592             )
593              
594              
595             Compute sum of squared errors for a sparse SVD-encoded matrix.
596             Should be equivalent to:
597              
598             sum( ($a - ($u x stretcher($s) x $v->xchg(0,1)))**2 )
599              
600             ... but computes all values on-the-fly, avoiding
601             memory bottlenecks for large sparse matrices.
602             This is a pure PDL::PP method, so you can use e.g.
603             C for the SVD-encoded matrix if you wish.
604              
605             Error contributions are computed even for "missing" (zero) values,
606             so running time is O(n*m).
607             Consider using L or L
608             to compute error rates
609             only for non-missing values if you have a large sparse matrix, e.g.:
610              
611             $svdvals = svdindexccs($u,$s,$v, $ptr,$rowids);
612             $err_nz = ($nzvals-$svdvals)->pow(2)->sumover;
613              
614              
615              
616             =for bad
617              
618             svderror does not process bad values.
619             It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.
620              
621              
622             =cut
623              
624              
625              
626              
627              
628              
629             *svderror = \&PDL::svderror;
630              
631              
632              
633              
634             ##---------------------------------------------------------------------
635             =pod
636              
637             =head1 ACKNOWLEDGEMENTS
638              
639             Perl by Larry Wall.
640              
641             PDL by Karl Glazebrook, Tuomas J. Lukka, Christian Soeller, and others.
642              
643             SVDLIBC by Dough Rohde.
644              
645             SVDPACKC by Michael Berry, Theresa Do, Gavin O'Brien, Vijay Krishna and Sowmini Varadhan.
646              
647             =cut
648              
649             ##----------------------------------------------------------------------
650             =pod
651              
652             =head1 KNOWN BUGS
653              
654             Globals still lurk in the depths of SVDLIBC.
655              
656             =cut
657              
658              
659             ##---------------------------------------------------------------------
660             =pod
661              
662             =head1 AUTHOR
663              
664             Bryan Jurish Emoocow@cpan.orgE
665              
666             =head1 COPYRIGHT AND LICENSE
667              
668             Copyright (c) 2005-2015, Bryan Jurish. All rights reserved.
669              
670             This package is free software, and entirely without warranty.
671             You may redistribute it and/or modify it under the same terms
672             as Perl itself, either version 5.20.2 or any newer version of Perl 5
673             you have available.
674              
675             The SVDLIBC sources included in this distribution are themselves
676             released under a BSD-like license. See the file
677             F in the PDL-SVDLIBC source distribution
678             for details.
679              
680             =head1 SEE ALSO
681              
682             perl(1), PDL(3perl), PDL::CCS(3perl), SVDLIBC documentation.
683              
684             =cut
685              
686              
687              
688             ;
689              
690              
691              
692             # Exit with OK status
693              
694             1;
695              
696