Mercurial > hg > octave-nkf
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 +++ b/scripts/special-matrix/hankel.m @@ -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); -