tidy svd pp_def

Author: mohawk2

lib/PDL/MatrixOps.pd

@@ -610,36 +610,27 @@ To compare with L<PDL::LinearAlgebra>:
 EOF
 );
 
-######################################################################
-### svd
 pp_def(
        "svd",
        HandleBad => 0,
-       Pars => 'a(n,m); [t]w(wsize=CALC($SIZE(n) * ($SIZE(m) + $SIZE(n)))); [o]u(n,m); [o,phys]z(n); [o]v(n,n);',
+       Pars => 'a(n,m); [o]u(n,m); [o,phys]s(n); [o]v(n,n);
+         [t]w(wsize=CALC($SIZE(n) * ($SIZE(m) + $SIZE(n))));',
        GenericTypes => ['D'],
        RedoDimsCode => '
         if ($SIZE(m)<$SIZE(n))
             $CROAK("svd requires input ndarrays to have m >= n; you have m=%td and n=%td. Try inputting the transpose.  See the docs for svd.",$SIZE(m),$SIZE(n));
        ',
-       Code => '
-              extern void SVD( double *W, double *Z, int nRow, int nCol );
-              double *t = $P(w);
-              loop (m,n) %{
-                *t++ = $a();
-              %}
-              SVD($P(w), $P(z), $SIZE(m), $SIZE(n));
-              loop (n) %{
-                $z() = sqrt($z());
-              %}
-              t = $P(w);
-              loop (m,n) %{
-                $u() = *t++/$z();
-              %}
-              loop (n1,n0) %{
-                $v() = *t++;
-              %}
-',
-      , Doc => q{
+       CHeader => 'void SVD( double *W, double *Z, int nRow, int nCol );',
+       Code => <<'EOF',
+  double *t = $P(w);
+    loop (m,n) %{ *t++ = $a(); %}
+  SVD($P(w), $P(s), $SIZE(m), $SIZE(n));
+  loop (n) %{ $s() = sqrt($s()); %}
+  t = $P(w);
+    loop (m,n) %{ $u() = *t++/$s(); %}
+    loop (n1,n0) %{ $v() = *t++; %}
+EOF
+      Doc => q{
 =for ref
 
 Singular value decomposition of a matrix.
@@ -654,9 +645,7 @@ x n, C<$v> is n x n, and there are <=n singular values. As long as m
 >= n, the original matrix can be reconstructed as follows:
 
     ($u,$s,$v) = svd($x);
-    $ess = zeroes($x->dim(0),$x->dim(0));
-    $ess->slice("$_","$_").=$s->slice("$_") foreach (0..$x->dim(0)-1); #generic diagonal
-    $a_copy = $u x $ess x $v->transpose;
+    $a_copy = $u x stretcher($s) x $v->transpose;
 
 If m==n, C<$u> and C<$v> can be thought of as rotation matrices that
 convert from the original matrix's singular coordinates to final