--- a/scripts/linear-algebra/krylov.m
+++ b/scripts/linear-algebra/krylov.m
@@ -1,162 +1,188 @@
-# Copyright (C) 1996,1998 A. Scottedward Hodel 
+# Copyright (C) 1993, 1998, 1999 A. Scottedward Hodel
+#
+# This file is part of Octave.
 #
-# This file is part of Octave. 
+# Octave is free software; you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published by the
+# Free Software Foundation; either version 2, or (at your option) any
+# later version.
+#
+# Octave is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+# for more details.
 #
-# Octave is free software; you can redistribute it and/or modify it 
-# under the terms of the GNU General Public License as published by the 
-# Free Software Foundation; either version 2, or (at your option) any 
-# later version. 
-# 
-# Octave is distributed in the hope that it will be useful, but WITHOUT 
-# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 
-# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License 
-# for more details.
-# 
-# You should have received a copy of the GNU General Public License 
-# along with Octave; see the file COPYING.  If not, write to the Free 
-# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. 
+# You should have received a copy of the GNU General Public License
+# along with Octave; see the file COPYING.  If not, write to the Free
+# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+
+function [Uret,H,nu] = krylov(A,V,k,eps1,pflg);
+  # function [U,H,nu] = krylov(A,V,k{,eps1,pflg});
+  # construct orthogonal basis U of block Krylov subspace;
+  # 	[V AV A^2*V ... A^(k+1)*V];
+  # method used: householder reflections to guard against loss of
+  # orthogonality
+  # eps1: threshhold for 0 (default: 1e-12)
+  # pflg: flag to use row pivoting  (improves numerical behavior)
+  #   0 [default]: no pivoting; prints a warning message if trivial
+  #                null space is corrupted
+  #   1          : pivoting performed
+  #
+  # outputs:
+  #   Uret: orthogonal basis of block krylov subspace
+  #   H: Hessenberg matrix; if V is a vector then A U = U H 
+  #      otherwise H is meaningless
+  # nu: dimension of span of krylov subspace (based on eps1)
+  # if B is a vector and k > m-1, krylov returns H = the Hessenberg 
+  # decompostion of A.
+  #
+  # Reference: Hodel and Misra, "Partial Pivoting in the Computation of 
+  #     Krylov Subspaces", In preparation.
+
+  defeps = 1e-12;
+  if(nargin < 3 | nargin > 5)
+    usage("[U,nu] = krylov(A,V,k{,eps1,pflg})")
+  elseif(nargin < 5)
+    pflg = 0;        # default permutation flag
+  endif
+  if(nargin < 4)
+    eps1 = defeps;    # default tolerance parameter
+  endif
+  if(isempty(eps1)) eps1 = defeps; endif
+  
+  na = is_square(A);
+  if( !na ) error("A(%d x %d) must be square",rows(A),columns(A)); endif
  
-function [U,H,k1] = krylov(A,v,k,eps1,pflg);
-  # function [U,H,k1] = krylov(A,v,k{,eps1});
-  # construct orthogonal basis U of Krylov subspace;
-  #     span([v Av A^2*v ... A^(k+1)*v]);
-  # via householder reflections; reflections are multiplied to obtain U
-  # inputs:
-  #   A: square matrix
-  #   v: column vector
-  #   k: desired krylov subspace dimension (as above)
-  #   eps1: threshhold for 0 relative to norm of current column (default: 1e-12)
-  #   pflg: permutation flag (default 0): avoid using zero rows of
-  #      [A,v] as householder pivots; this avoids spurious non-zero entries
-  #      in returned orthogonal matrix U, but destroys the Householder matrix
-  #      structure of H.
-  # outputs:
-  #   U: (n x k+1) orthogonal basis of Krylov subspace. A*U(:,1:k) = U*H
-  #   H: (pflg=0): Hessenberg matrix satisfying A*U(:,1:k) = U*H
-  #      (pflg=1): Workspace; does not satisfy above equation.
-  # k1: dimension of span of krylov subspace (based on eps1)
-  # if k > m-1, krylov returns the Hessenberg decompostion of A.
+  [m,kb] = size(V); 
+  if(m != na);
+    error("A(%d x %d), V(%d x %d): argument dimensions do not match", ...
+      na,na,m,kb)
+  endif
 
-  # Written by A. S. Hodel 1992
-  # $Revision: 1.2 $
-  # $Log$
-
-  save_empty_list_elements_ok = empty_list_elements_ok;
-  empty_list_elements_ok = 1;
+  if( !is_scalar(k) )
+    error("krylov: third argument must be a scalar integer")
+  endif
 
-  if    (nargin > 5)   usage("[U,H,k1] = krylov(A,v,k{,eps1,pflg})");
-  elseif(nargin < 5)   pflg = 0;
-  elseif(nargin < 4)   eps1 = 1e-12; endif
-  na = is_square(A);
-  if(!na) error("krylov: A(%dx%d) must be square",na,na); endif
-  [m,n] = size(v);
-  if(!is_vector(v)) error("krylov: v(%dx%d) must be a vector",m,n);
-  elseif(length(v) != na)
-    error("krylov: A(%dx%d), v(%dx1); mismatch",na,na,length(v));
-  endif
-  v = vec(v);    # make sure it's a column vector
-  if(!is_scalar(k))
-    error("krylov: k(%dx%d) must be a scalar",rows(k), columns(k));
-  elseif( k > m)
-    warning("krylov: k is too large; reducing to %d",m-1);
-    k = m-1;
+  Vnrm = norm(V,Inf);
+  
+  # check for trivial solution
+  if(Vnrm == 0)
+    Uret = []; nu = 0;  return;
   endif
 
-  # check for zero input vector
-  if(norm(v) == 0) U = []; H = []; k1 = 0; return endif
+  # identify trivial null space
+  abm = max(abs([A,V]'));  nzidx = find(abm != 0);  zidx = find(abm == 0);
+
+  # set up vector of pivot points
+  pivot_vec = 1:na;
 
-  # indices of legal pivot points (identifies trivial null space).
-  abm = max(abs([A,v]')); nzidx = find(abm != 0); zidx = find(abm == 0);
+  iter = 0;
+  alpha = [];  
+  nh = 0;
+  while (length(alpha) < na) & (columns(V) > 0) & (iter < k)
+    iter++;
 
-  # check permutation flag
-  if(pflg)
-    # permute zero rows of [A,v] to trailing rows
-    permvec = [vec(nzidx); vec(zidx)];
-    pmat = eye(na); pmat = pmat(permvec,:);
-    ipermvec = pmat'*vec(1:na);
-    A = A(permvec,permvec);
-    v = v(permvec);
-  endif
+    # get orthogonal basis of V
+    jj = 1;
+    while(jj <= columns(V) & length(alpha) < na)
+      nu = length(alpha)+1;   # index of next Householder reflection
+     
+      short_pv = pivot_vec(nu:na);
+      q = V(:,jj);
+      short_q = q(short_pv);
 
-  k1 = k+1;		# Assume subspace basis has full rank
-  [m,n] = size(A);
-  [hv,alpha(1),z] = housh(v,1,0);
-
-  # initial orthogonal vector
-  q = zeros(n,1);
-  q(1) = 1;
-  q = q - alpha*hv*hv'*q;
-  normq = norm(q);
-  normres = normq;
+      if(norm(short_q) < eps1)
+        # insignificant column; delete
+        nv = columns(V);
+        if(jj != nv) 
+          [V(:,jj),V(:,nv)] = swap(V(:,jj),V(:,nv)); 
+          # FIX ME: H columns should be swapped too.  Not done since
+          # Block Hessenberg structure is lost anyway.
+        endif
+        V = V(:,1:(nv-1));
+        nu = nu - 1;    # one less reflection
 
-  U(:,1) = hv;
-  j = 1;
-  while(j <= k & normres > eps1*normq)
-    # multiply to get new vector;
-    q = A*q;
-    normq = norm(q);
+      else
+        # new householder reflection
+        if(pflg)
+          # locate max magnitude element in short_q
+          asq = abs(short_q);
+	  maxv = max(asq);
+          maxidx = find(asq == maxv);
+          pivot_idx = short_pv(maxidx(1));
+           
+          # see if need to change the pivot list
+          if(pivot_idx != pivot_vec(nu))
+            swapidx = maxidx(1) + (nu-1);
+            [pivot_vec(nu),pivot_vec(swapidx)] = ...
+              swap(pivot_vec(nu),pivot_vec(swapidx));
+          endif
+        endif
+        
+        # isolate portion of vector for reflection
+        idx = pivot_vec(nu:na);
+        jdx = pivot_vec(1:nu);
+
+        [hv,av,z] = housh(q(idx),1,0);
+        alpha(nu) = av;
+        U(idx,nu) = hv;
 
-    # multiply by householder reflections to obtain projected vector and the
-    # next column of H
-    for i=1:j
-     hv = U(i:n,i);
-     av = alpha(i);
-     q(i:n,1) = q(i:n,1)-av*hv*(hv'*q(i:n,1));
-    endfor
+        # reduce V per the reflection
+        V(idx,:) = V(idx,:) - av*hv*(hv' * V(idx,:));
+        if(iter > 1)
+          # FIX ME: not done correctly for block case
+          H(nu,nu-1) = V(pivot_vec(nu),jj);
+        endif
+
+        # advance to next column of V
+        jj=jj+1;
+      endif
+    endwhile
 
-    i = j+1;
-    # compute and apply next householder vector;
-    if(i <= n)
-      [hv,av,z] = housh(q(i:n,1),1,0);
-      alpha(i) = av;
-      q(i:n,1) = q(i:n,1)-av*hv*(hv'*q(i:n,1));
-      U(i:n,i) = hv;
-      H(1:i,j) = q(1:i);
-    else
-      av = 0;
-      H(:,j) = q;	# complete Hessenberg decomposition
+    if(columns(V) > 0)
+      # construct next Q and multiply
+      Q = zeros(size(V));
+      for kk=1:columns(Q)
+        Q(pivot_vec(nu-columns(Q)+kk),kk) = 1;
+      endfor
+      
+      # apply Householder reflections
+      for ii = nu:-1:1
+        idx = pivot_vec(ii:na);
+        hv = U(idx,ii);
+        av = alpha(ii);
+        Q(idx,:) = Q(idx,:) - av*hv*(hv'*Q(idx,:));
+      endfor
     endif
 
-    # see if done yet
-    normres = norm(q(i:n));
-    if(normres > eps1*normq)
-      j = j+1;
-    else   
-      k1 = min(k1,j);	# time to quit; norm of residual is small
-    endif
-  
-    # back out new q vector for next pass;
-    j1 = columns(U);
-    q = zeros(n,1);
-    q(j1) = 1;
-    for i=j1:-1:1;
-      hv = U(i:n,i);
-      av = alpha(i);
-      q(i:n,1) = q(i:n,1)-av*hv*(hv'*q(i:n,1));
-    endfor
+    # multiply to get new vector;
+    V = A*Q;
+    # project off of previous vectors
+    nu = length(alpha);
+    for i=1:nu
+      hv = U(:,i);  av = alpha(i);
+      V = V - av*hv*(hv'*V);
+      H(i,nu-columns(V)+(1:columns(V))) = V(pivot_vec(i),:);
+    end
+
   endwhile
 
+  # Back out complete U matrix
   # back out U matrix ;
   j1 = columns(U);
   for i=j1:-1:1;
-   hv = U(i:n,i);
-   av = alpha(i);
-   U(:,i) = zeros(n,1);
-   U(i,i) = 1;
-   U(i:n,i:j1) = U(i:n,i:j1)-av*hv*(hv'*U(i:n,i:j1));
+    idx = pivot_vec(i:na);
+    hv = U(idx,i);
+    av = alpha(i);
+    U(:,i) = zeros(na,1);
+    U(idx(1),i) = 1;
+    U(idx,i:j1) = U(idx,i:j1)-av*hv*(hv'*U(idx,i:j1));
   endfor
-
-  # check permutation flag
-  if(pflg)
-    # permute rows of U back to original coordinates
-    U = U(ipermvec,:);
-  endif
-
-  # check for spurious nonzero entries
-  if( max(max( abs(U(zidx,:)) )) > eps1 )
+  
+  nu = length(alpha);
+  Uret = U;
+  if( max(max( abs(Uret(zidx,:)) )) > 0)
     warning("krylov: trivial null space corrupted; set pflg=1 or eps1>%e",eps1);
   endif
 
-  empty_list_elements_ok = save_empty_list_elements_ok;
-
 endfunction