FUNCTION gramschmidt,A,R
;+
;Q=gramSchmidt(A,R) returns Q containing the Gram-Schmidt
;orthogonalization of the columns of A. An upper triangular
;matrix R such that A=Q##R is returned in the second argument.
;
;H. Rhody July, 2000
;-
SA=SIZE(A) & m=SA[2] & n=SA[1]
Q=FLTARR(n,m) & R=FLTARR(n,n)
e1=A[0,*] ;First column of A
R[0,0]=norm(e1)
Q[0,*]=e1/R[0,0]

FOR k=1,n-1 DO BEGIN
	R[k,0:k-1]=TRANSPOSE(Q[0:k-1,*])##A[k,*]
	ek=A[k,*]-Q[0:k-1,*]##R[k,0:k-1]
	R[k,k]=norm(ek)
	IF R[k,k] GT 5E-7 THEN ek=ek/R[k,k] ELSE ek=0*ek
	Q[k,*]=ek
ENDFOR
RETURN,Q
END