--- triu.m.orig	2009-10-22 17:36:28.000000000 +0200
+++ triu.m	2009-10-22 17:49:10.000000000 +0200
@@ -28,25 +28,70 @@
 
   if (nargin > 0)
     if (isstruct (x))
-       error ("tril: structure arrays not supported");
+       error ("triu: structure arrays not supported");
      endif 
     [nr, nc] = size (x);
   endif
   if (nargin == 1)
     k = 0;
   elseif (nargin == 2)
-    if ((k > 0 && k > nc) || (k < 0 && k < -nr))
+    if ((k > 0 && k > nc) || (k < 0 && k < -nr+1))
       error ("triu: requested diagonal out of range");
     endif
   else
     print_usage ();
   endif
-
-  retval = resize (resize (x, 0), nr, nc);
-  for j = max (1, k+1) : nc
-    nr_limit = min (nr, j-k);
-    retval (1:nr_limit, j) = x (1:nr_limit, j);
-  endfor
+  if (nr == nc)
+    if (k <= 0)
+      retval = x;      
+      for j = 1 : nr-1+k
+	retval((j-1)*nr+j+1-k:j*nr) = 0;
+      endfor	
+    else
+% k > 0
+      retval = zeros(size(x));
+      for j = nc : -1 : k+1
+	idx = (j-1)*nr+1:(j-1)*nr-k+j;
+	retval(idx) = x(idx);
+      endfor	
+    endif
+  elseif (nr > nc)
+    if (k <= floor((1+nc-nr)/2))
+      retval = x;      
+      for j = 1 : min(nc,nr-1+k)
+	retval((j-1)*nr+j+1-k:j*nr) = 0;
+      endfor	
+    else
+% k > floor((1+nc-nr)/2)
+      retval = zeros(size(x));
+      for j = nc : -1 : max(1,k+1)
+	idx = (j-1)*nr+1:(j-1)*nr-k+j;
+	retval(idx) = x(idx);
+      endfor	
+    endif
+  else
+% nc > nr
+    if (k <= floor((1+nc-nr)/2))
+      retval = x;      
+      for j = 1 : k-1
+	retval((j-1)*nr+1:j*nr) = 0;
+      endfor	
+      for j = max(1,k) : nr-1+k
+	retval((j-1)*nr+j+1-k:j*nr) = 0;
+      endfor	
+    else
+% k > floor((1+nc-nr)/2)
+      retval = zeros(size(x));
+      for j = nc : -1 : k+nr+1
+	idx = (j-1)*nr+1:j*nr;
+	retval(idx) = x(idx);
+      end	
+      for j = min(nc,k+nr) : -1 : k+1
+	idx = (j-1)*nr+1:(j-1)*nr-k+j;
+	retval(idx) = x(idx);
+      endfor	
+    endif
+  endif
 
 endfunction