changeset 13888:c78ac846fcbc

hankel.m: Recode for 3.5X speedup * hankel.m: Recode for 3.5X speedup
author Rik <octave@nomad.inbox5.com>
date Fri, 18 Nov 2011 12:28:20 -0800
parents e0fb702a62a4
children aaefd6b28188
files scripts/special-matrix/hankel.m
diffstat 1 files changed, 31 insertions(+), 50 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/special-matrix/hankel.m	Fri Nov 18 12:25:19 2011 -0800
+++ b/scripts/special-matrix/hankel.m	Fri Nov 18 12:28:20 2011 -0800
@@ -49,69 +49,50 @@
 
 function retval = hankel (c, r)
 
-  if (nargin == 1)
-    r = resize (resize (c, 0), size(c));
-  elseif (nargin != 2)
+  if (nargin < 1 || nargin > 2)
     print_usage ();
   endif
 
-  [c_nr, c_nc] = size (c);
-  [r_nr, r_nc] = size (r);
-
-  if ((c_nr != 1 && c_nc != 1) || (r_nr != 1 && r_nc != 1))
-    error ("hankel: expecting vector arguments");
-  endif
-
   if (nargin == 1)
-    r (1) = c (length (c));
-  endif
+
+    if (! isvector (c))
+      error ("hankel: C must be a vector");
+    endif
 
-  if (c_nc != 1)
-    c = c.';
-  endif
+    nr = length (c);
+    nc = nr;
+    data = [c(:) ; zeros(nr, 1)];
+
+  else
 
-  if (r_nr != 1)
-    r = r.';
-  endif
+    if (! (isvector (c) && isvector (r))) 
+      error ("hankel: C and R must be vectors");
+    elseif (r(1) != c(end))
+      warning ("hankel: column wins anti-diagonal conflict");
+    endif
 
-  nc = length (r);
-  nr = length (c);
+    nr = length (c);
+    nc = length (r);
+    data = [c(:) ; r(2:end)(:)];
 
-  if (r (1) != c (nr))
-    warning ("hankel: column wins anti-diagonal conflict");
   endif
-
-  ## This should probably be done with the colon operator...
-
-  retval = resize (resize (c, 0), nr, nc);
-
-  for i = 1:min (nr, nc)
-    retval (1:nr-i+1, i) = c (i:nr);
-  endfor
-
-  tmp = 1;
-  if (nc <= nr)
-    tmp = nr - nc + 2;
-  endif
-
-  for i = nr:-1:tmp
-    retval (i, 2+nr-i:nc) = r (2:nc-nr+i);
-  endfor
+   
+  slices = cellslices (data, 1:nc, nr:1:nc+nr-1);
+  retval = horzcat (slices{:});
 
 endfunction
 
-%!assert(hankel(1:3),[1,2,3;2,3,0;3,0,0])
-%!assert(hankel(1),[1]);
-%!assert(hankel(1:3,3:6),[1,2,3,4;2,3,4,5;3,4,5,6]);
-%!assert(hankel(1:3,3:4),[1,2;2,3;3,4]);
-%!assert(hankel(1:3,4:6),[1,2,3;2,3,5;3,5,6]);
 
-%!assert((hankel (1) == 1 && hankel ([1, 2]) == [1, 2; 2, 0]
-%! && hankel ([1, 2], [2; -1; -3]) == [1, 2, -1; 2, -1, -3]));
-
-%!error hankel ([1, 2; 3, 4], [1, 2; 3, 4]);
+%!assert (hankel (1), [1])
+%!assert (hankel ([1, 2]), [1, 2; 2, 0])
+%!assert (hankel ([1, 2], [2; -1; -3]), [1, 2, -1; 2, -1, -3])
+%!assert (hankel (1:3), [1,2,3;2,3,0;3,0,0])
+%!assert (hankel (1:3,3:6), [1,2,3,4;2,3,4,5;3,4,5,6])
+%!assert (hankel (1:3,3:4), [1,2;2,3;3,4])
+%!assert (hankel (1:3,4:6), [1,2,3;2,3,5;3,5,6])
 
 %!error hankel ();
+%!error hankel (1, 2, 3);
+%!error <C must be a vector> hankel ([1, 2; 3, 4])
+%!error <C and R must be vectors> hankel (1:4, [1, 2; 3, 4])
 
-%!error hankel (1, 2, 3);
-