FUNCTION ORTHO_RANK,A,TAU=tau,H,R,K
;+
;rankA=Ortho_Rank(A,TAU=tau) is the rank of a matrix A.
;If A is a scalar then rankA=0. tau is the threshold
;for zero values. Default: tau=1e-6.
;
;Uses QR,A,Q,R and counts the number rows of R
;that contain at least one value GE tau.
;
;Arrays H,R,K such that A=H##R##Transpose(K) are
;returned through the arguments if ortho_rank is called
;with rankA=Ortho_Rank(A,TAU=tau,H,R,K)
;
;H. Rhody
;June, 2004
;-


IF N_elements(tau) LE 0 THEN tau=1e-6
sa=Size(A)

IF (sa[0] EQ 0) THEN BEGIN
	H=[1] & R=[A] & K=[1]
	RETURN,0
ENDIF
IF (sa[0] EQ 1) THEN BEGIN
	QR,Transpose(A),Q,R
	H=[1]
	R=Transpose(R)
	K=Q
	return,1
ENDIF

IF sa[0] EQ 2 THEN BEGIN
  IF sa[1] LE sa[2] THEN BEGIN ;n <= m case
  ;Do a SVD to find the rank and sort the columns
  ;in terms of descending singular values
  SVDC,A,s,U,V
  ks=reverse(sort(s))
  rankA=fix(total(s GE tau*max(s)))
  ;Construct a permutation matrix to make the
  ;first rankA columns linearly independent.
  P=(Identity(sa[1]))[ks,*]
  A1=A##P
  QR,A1,Q,R1
  QR,transpose(R1),QQ,RR
  	H=Q
  	R=Transpose(RR)
  	K=P##QQ
  	Return,rankA
  ENDIF ELSE BEGIN ;n>m case
    A0=Transpose(A)
    svdc,A0,s,U,V
    ks=reverse(sort(s))
    rankA=fix(total(s GE tau*max(s)))
    P=(Identity(sa[2]))[ks,*]
    A1=A0##P
    QR,A1,Q,R1
    QR,Transpose(R1),QQ,RR
    H=P##QQ
    R=RR
    K=Q
    return,rankA
    ENDELSE
ENDIF
RETURN,-1 ;Array dimension failure
end