##-*- Mode: Perl -*- ##====================================================================== ## Header Administrivia ##====================================================================== our $VERSION = '0.06'; pp_setversion($VERSION); ##------------------------------------------------------ ## Header: pods pp_addpm({At=>'Top'},<<'EOPM'); =pod =head1 NAME PDL::SVDLIBC - PDL interface to Doug Rohde's SVD C Library =head1 SYNOPSIS use PDL; use PDL::SVDLIBC; ##--------------------------------------------------------------------- ## Input matrix (dense) ##--------------------------------------------------------------------- $n = 100; ##-- number of columns $m = 50; ##-- number of rows $a = random(double,$n,$m); ##-- random matrix ##--------------------------------------------------------------------- ## Output pdls ##--------------------------------------------------------------------- $d = $n; ##-- max number of output dimensions $ut = zeroes(double,$m,$d); ##-- left singular components $s = zeroes(double,$d); ##-- singular values (diagnonal vector) $vt = zeroes(double,$n,$n); ##-- right singular components ##--------------------------------------------------------------------- ## Singular Value Decomposition (dense) ##--------------------------------------------------------------------- svdlas2d($a, $maxiters, $end, $kappa, $ut, $s, $vt); ##--------------------------------------------------------------------- ## Singular Value Decomposition (sparse) ##--------------------------------------------------------------------- use PDL::CCS; ($ptr,$rowids,$nzvals) = ccsencode($a); $ptr->reshape($ptr->nelem+1); $ptr->set(-1, $rowids->nelem); svdlas2($ptr, $rowids, $nzvals, $m, $maxiters, $end, $kappa, $ut, $s, $vt); =head1 DESCRIPTION PDL::SVDLIBC provides a PDL interface to the SVDLIBC routines for singular value decomposition of large sparse matrices. SVDLIBC is available from http://tedlab.mit.edu/~dr/SVDLIBC/ =cut EOPM ## /pm additions ##------------------------------------------------------ ##------------------------------------------------------ ## Exports: None #pp_export_nothing(); ##------------------------------------------------------ ## Includes / defines pp_addhdr(<<'EOH'); #include /*#define CDEBUG 1*/ /*#define DEBUG_CODE 1*/ #if defined(CDEBUG) || defined(DEBUG_CODE) # define SVD_VERBOSITY_DEFAULT 1 #else # define SVD_VERBOSITY_DEFAULT 0 #endif EOH ##====================================================================== ## C Utilities ##====================================================================== pp_addhdr(<<'EOH'); /*--------------------------------------------------------------*/ double **p2pp_dbl(int nrows, int ncols, double *p) { int i; #ifdef CDEBUG int j; #endif double **matrix; if (!(p && nrows && ncols)) return NULL; matrix = malloc(nrows*sizeof(double**)); //New(0,matrix,nrows,double*); for (i=0; i < nrows; i++) { matrix[i] = p + (i*ncols); #ifdef CDEBUG printf("p2pp_dbl(nr=%d,nc=%d,p=%p) : (p+%ld*%ld)=%p\n", nrows,ncols,p, i,ncols,matrix[i]); for (j=0; j < ncols; j++) { printf("p2pp_dbl: [i=%ld][j=%ld] : val=%g\n", i, j, matrix[i][j]); } #endif } return matrix; } /*--------------------------------------------------------------*/ static void pp2pdl(int nrows, int ncols, double **pp, double *p) { int i,j; for (i=0; iUt->rows, svdr->Ut->cols, svdr->Ut->value, up); pp2pdl(1, svdr->d, &(svdr->S), sp); pp2pdl(svdr->Vt->rows, svdr->Vt->cols, svdr->Vt->value, vp); } /*--------------------------------------------------------------*/ #ifdef DEBUG_CODE static void showpp_int(const char *name, int nrows, int ncols, int **pp) { int i,j; printf("\n"); for (i=0; i= 1) SVDVerbosity = SvIV(ST(0)); RETVAL = SVDVerbosity; OUTPUT: RETVAL '); ##-------------------------------------------------------------- ## SVD Globals: version pp_addpm(<<'EOPM'); =pod =head2 PDL::SVDLIBC::svdVersion() Returns a string representing the SVDLIBC version this module was compiled with. =cut EOPM pp_addxs('',' char * svdVersion() CODE: RETVAL = SVDVersion; OUTPUT: RETVAL '); ##====================================================================== ## SVD Utilities pp_addpm(<<'EOPM'); =pod =head1 SVD Utilities =cut EOPM pp_def ('_svdccsencode', #Pars=> q(double a(n,m); longlong [o]ptr(n1); longlong [o]rowids(nnz); double [o]nzvals(nnz)), Pars => q(double a(n,m); int [o]ptr(n1); int [o]rowids(nnz); double [o]nzvals(nnz)), Code => (' struct dmat dm; SMat sm; dm.rows = $SIZE(m); dm.cols = $SIZE(n); /*printf("p2pp_dbl()\n");*/ dm.value = p2pp_dbl($SIZE(m), $SIZE(n), $P(a)); /*printf("svdConvertDtoS()\n");*/ sm = svdConvertDtoS(&dm); /*printf("pp2pdl()\n");*/ pp2pdl_int(1, $SIZE(n)+1, &sm->pointr, (int *) $P(ptr)); pp2pdl_int(1, sm->vals, &sm->rowind, (int *) $P(rowids)); pp2pdl (1, sm->vals, &sm->value, (double *)$P(nzvals)); /*-- cleanup --*/ /*printf("cleanup()\n");*/ if (dm.value) free(dm.value); if (sm) svdFreeSMat(sm); '), ); #pp_add_exported('', 'svdccsencode'); ##------------------------------------------------------ ## svdlas2a() : singular value decomposition (convenience) pp_add_exported('','svdlas2a'); pp_addpm(<<'EOPM'); =pod =head2 svdlas2a =for sig int ptr(nplus1); int rowids(nnz); double nzvals(nnz); int nrows(); ##-- default: max($rowids)+1 int d(); ##-- default: nplus1-1 int iterations(); ##-- default: 2*$d double end(2); ##-- default: [-1e-30,1e-30] double kappa(); ##-- default: 1e-6 double [o]ut(m,d); ##-- default: new double [o] s(d); ##-- default: new double [o]vt(n,d); ##-- default: new Uses a variant of the single-vector Lanczos method (Lanczos, 1950) to compute the singular value decomposition of a sparse matrix with $nrows() rows and data encoded in Harwell-Boeing sparse format in the input parameters $ptr(), $rowids(), and $nzvals(). See L<"PDL::CCS"> for a way to acquire these parameters from a dense input matrix, but note that for svdlas2(), the column pointer $ptr() is of size ($n+1) for a dense matrix $a with $n columns, where $ptr($n)==$nnz is the total number of nonzero values in $a. $iterations() is the maximum number of Lanczos iterations to perform. $end() specifies two endpoints of an interval within which all unwanted eigenvalues lie. $kappa() is a double containing the relative accuracy of Ritz values acceptable as eigenvalues. The left singular components are returned in the matrix $ut(), the singular values themselved in the vector $s(), and the right singular components in the matrix $vt(). Note that $ut() and $vt() are transposed, and must be specified explicitly in the call, so that the degree of reduction (the size parameter $d) can be determined. If $d==$n, then a full decomposition will be computed, and on return, $ut() and $vt() should be transposed instances of the matrices $u() and $v() as returned by PDL::MatrixOps::svd(). The Lanczos method as used here seems to be consistently the fastest. This algorithm has the drawback that the low order singular values may be relatively imprecise, but that is not a problem for most users who only want the higher-order values or who can tolerate some imprecision. See also: svdlas2d() =cut ## ($iters,$end,$kappa,$ut,$s,$vt) = svddefaults($nrows,$ncols,$d, $iters,...) ## + returns default values sub svddefaults { my ($nrows,$ncols,$d, $iters,$end,$kappa,$ut,$s,$vt) = @_; $nrows = $nrows->at(0) if (UNIVERSAL::isa($nrows,'PDL')); $ncols = $ncols->at(0) if (UNIVERSAL::isa($ncols,'PDL')); $d = $ncols if (!defined($d)); $iters = 2*$d if (!defined($iters)); $end = pdl(double,[-1e-30,1e-30]) if (!defined($end)); $kappa = pdl(double,1e-6) if (!defined($kappa)); $ut = PDL->zeroes(double,$nrows,$d) if (!defined($ut)); $s = PDL->zeroes(double,$d) if (!defined($s)); $vt = PDL->zeroes(double,$ncols,$d) if (!defined($vt)); return ($iters,$end,$kappa,$ut,$s,$vt); } sub svdlas2a { my ($ptr,$rowids,$nzvals, $nrows,$d, @args) = @_; $nrows = $rowids->flat->max+1 if (!defined($nrows)); @args = svddefaults($ptr->dim(0)-1,$nrows,$d,@args); svdlas2($ptr,$ropwids,$nzvals,$nrows,@args); return @args[3..5]; } EOPM ##------------------------------------------------------ ## svdLAS2() : singular value decomposition (low-level) pp_def ('svdlas2', Pars => join("\n ", '', q(int ptr(nplus1);), ##-- longlong q(int rowids(nnz);), ##-- longlong q(double nzvals(nnz);), q(int nrows();), q(int iterations();), ##-- longlong q(double end(2);), q(double kappa();), q(double [o]ut(m,d);), q(double [o] s(d);), q(double [o]vt(n,d);), '', ), Code => (' struct smat sm; SVDRec svdr; /*-- setup sparse matrix --*/ sm.rows = (int) $nrows(); sm.cols = (int) $SIZE(n); sm.vals = (int) $SIZE(nnz); sm.pointr = (int *)$P(ptr); sm.rowind = (int *)$P(rowids); sm.value = (double *)$P(nzvals); /*-- view decoded stuff --*/ #ifdef DEBUG_CODE printf("sm: rows=%ld; cols=%ld; vals=%ld\n", sm.rows, sm.cols, sm.vals); showpp_int("sm.pointr", 1, sm.cols+1, &sm.pointr); showpp_int("sm.rowind", 1, sm.vals, &sm.rowind); showpp_dbl("sm.value", 1, sm.vals, &sm.value); printf("--\n"); #endif /*-- run SVD --*/ svdr = svdLAS2(&sm, $SIZE(d), $iterations(), $P(end), $kappa()); svdrec2pdls(svdr, $P(ut), $P(s), $P(vt)); /*-- output some stats (debug) --*/ #ifdef DEBUG_CODE printf("svdrec.d=%d ~ d=%d\n", svdr->d, $SIZE(d)); printf("svdrec.Ut.rows=%ld ~ d=%d\n", svdr->Ut->rows, $SIZE(d)); printf("svdrec.Ut.cols=%ld ~ m=%d\n", svdr->Ut->cols, $SIZE(m)); printf("svdrec.Vt.rows=%ld ~ n=%d\n", svdr->Vt->cols, $SIZE(n)); printf("svdrec.Vt.cols=%ld ~ d=%d !~ n=%d\n", svdr->Vt->cols, $SIZE(d), $SIZE(n)); #endif /*-- pre-cleanup --*/ sm.pointr = NULL; sm.rowind = NULL; sm.value = NULL; /*-- cleanup --*/ if (svdr) svdFreeSVDRec(svdr); '), Doc => q( Guts for svdlas2a(). No default instantiation, and slightly different calling conventions. ), ); #pp_add_exported('', 'svdlas2'); ##------------------------------------------------------ ## svdlas2ad() : singular value decomposition (dense): convenience pp_add_exported('','svdlas2ad'); pp_addpm(<<'EOPM'); =pod =head2 svdlas2ad =for sig double a(n,m); int d(); ##-- default: $n int iterations(); ##-- default: 2*$d double end(2); ##-- default: [-1e-30,1e-30] double kappa(); ##-- default: 1e-6 double [o]ut(m,d); ##-- default: new double [o] s(d); ##-- default: new double [o]vt(n,d); ##-- default: new As for svdlas2(), but implicitly converts the dense input matrix $a() to sparse format before computing the decomposition. =cut sub svdlas2ad { my ($a,$d, @args) = @_; @args = svddefaults($a->dim(1),$a->dim(0),$d,@args); svdlas2d($a,@args); return @args[3..5]; } EOPM ##------------------------------------------------------ ## svdlas2d() : singular value decomposition (dense) pp_def ('svdlas2d', Pars => join("\n ", '', q(double a(n,m);), q(int iterations();), q(double end(2);), q(double kappa();), q(double [o]ut(m,d);), q(double [o] s(d);), q(double [o]vt(n,d);), '', ), Code => (' struct dmat dm; SMat smp; SVDRec svdr; /*-- setup dense matrix --*/ dm.rows = $SIZE(m); dm.cols = $SIZE(n); dm.value = p2pp_dbl($SIZE(m), $SIZE(n), $P(a)); /*-- generate sparse matrix --*/ smp = svdConvertDtoS(&dm); /*-- compute svd --*/ svdr = svdLAS2(smp, $SIZE(d), $iterations(), $P(end), $kappa()); svdrec2pdls(svdr, $P(ut), $P(s), $P(vt)); /*-- cleanup --*/ if (svdr) svdFreeSVDRec(svdr); if (smp) svdFreeSMat(smp); if (dm.value) free(dm.value); '), Doc => q( Guts for _svdlas2d(). ), ); #pp_add_exported('', 'svdlas2d'); ##====================================================================== ## Footer Administrivia ##====================================================================== ##------------------------------------------------------ ## footer: pm additions pp_addpm(<<'EOPM'); ##--------------------------------------------------------------------- =pod =head1 ACKNOWLEDGEMENTS Perl by Larry Wall. PDL by Karl Glazebrook, Tuomas J. Lukka, Christian Soeller, and others. SVDLIBC by Dough Rohde. SVDPACKC by Michael Berry, Theresa Do, Gavin O'Brien, Vijay Krishna and Sowmini Varadhan. =cut ##---------------------------------------------------------------------- =pod =head1 KNOWN BUGS Globals still lurk in the depths of SVDLIBC. =cut ##--------------------------------------------------------------------- =pod =head1 AUTHOR Bryan Jurish Emoocow@cpan.orgE =head2 Copyright Policy Copyright (C) 2005, Bryan Jurish. All rights reserved. This package is free software, and entirely without warranty. You may redistribute it and/or modify it under the same terms as Perl itself. =head1 SEE ALSO perl(1), PDL(3perl), PDL::CCS(3perl), SVDLIBC documentation. =cut EOPM # Always make sure that you finish your PP declarations with # pp_done pp_done(); ##----------------------------------------------------------------------