function QRdata = recqr(QR)
  [M, N] = size(QR);
  MN = min(M,N);
  tau = zeros(1,N);

  T = zeros(MN);
  Y = zeros(M, MN);

  ttall = 0;
  tthouse = 0;
  ttupd = 0;

  ttall_start=cputime();
  for j = 1:MN,
    #j
    kahead = bitand(j-1, 1+bitcmp(j-1));
    k = min(1, kahead);
    while k < kahead,
      # Form T
      Ti = j-2*k;
      Tj = j-k;

      Y1cols = Ti:Ti+k-1;
      Y2cols = Tj:Tj+k-1;
      
      YY = Y(Tj:M, Y1cols)' * Y(Tj:M, Y2cols);
      T(Y1cols, Y2cols) = -T(Y1cols, Y1cols) * YY * T(Y2cols, Y2cols);
      
      k *= 2;
    endwhile

    kahead = min(kahead, MN-j+1);

    if kahead >= 1,
      Ycols = j-k:j-1;
      update_cols = j:j+kahead-1;
      update_rows = j-k:M;
      ttupd_start = cputime();
      # parentheses for performance...
      QR(update_rows,update_cols) -= \
 	  (Y(update_rows,Ycols) * T(Ycols, Ycols)') * \
 	  (Y(update_rows,Ycols)' * QR(update_rows,update_cols));
      ttupd += cputime() - ttupd_start;
    endif

    tthouse_start = cputime();
    [Y(j:M,j), tau(j), alpha] = make_house(QR(j:M,j));
    tthouse += cputime() - tthouse_start;
    QR(j+1:M,j) = Y(j+1:M,j);
    QR(j,j) = alpha;
    T(j,j) = tau(j);
  endfor

  if M < N,
    QR(:,M+1:N) -= (Y*T')*(Y'*QR(:,M+1:N));
  endif

#  ttall = cputime()-ttall_start
#  tthouse
#  ttupd

  QRdata.QR = QR;
  QRdata.tau = tau;

endfunction