PRO QR,A,Q,R
;+
;QR,A,Q,R computes matrices Q and R such that A=QR where R is an
;upper triangular matrix. Q is constructed by a sequence of
;Householder transformations.
;
;Harvey Rhody
;Laboratory for Imaging Algorithms and Systems
;Rochester Institute of Technology
;March, 2004
;-

;Find the size of A
sa=size(A,/dim)
IF n_elements(sa) EQ 1 THEN BEGIN
	;Handle the scalar case
	Q=[1.0]
	r=[A]
	RETURN
ENDIF

m=sa[1] & n=sa[0]
I=Identity(m)
Q=Identity(m)
R=DOUBLE(A)
FOR k=0,min([m-1,n])-1 DO BEGIN
	v=R[k,*]
	sigma=(v[k] GE 0) ? 1.0 : -1.0
	s=-sigma*sqrt(total(v[k:*]^2))
	u=v
	IF k GT 0 THEN u[0:k-1]=0
   u[k]=v[k]-s
   b=s*u[k]
   IF b NE 0 THEN P=(I+u##transpose(u)/b) ELSE P=I
   R=P##R
   Q=P##Q

ENDFOR
	Q=transpose(Q)
END