FUNCTION ORTHO_FACTOR,Ain,TAU=tau,H,R,K,s,U,V
;+
;rankA=Ortho_Factor(A,H,R,K,TAU=tau) returns the rank
;of a matrix A and provides orthogonal matrices
;H and K and an upper-triangular matrix R such that
;A=H##R##Transpose(K). H,R and K are computed by the
;QR algorithm.
;
;The SVDC decomposition is available through arguments
;s,U,V by calling
;rankA=Ortho_Factor(A,H,R,K,s,U,V)
;
;KEYWORD
;TAU is the threshold for rank determination. Rank is
;calculated by examining the singular values. If
;s[k] is the last first singular value such that
;s[k] LT s[0]*tau then the array has rank=k. The default
;value is TAU=1e-6.
;
;If A is a scalar then rankA=0. The SVD arrays are not
;defined for this case.
;
;If A is a row vector then rankA=1. The SVD arrays are
;not defined for this case.
;
;Returns -1 in cases not handled, such as arrays with
;dimensions greater than 2D.
;
;USES algorithm QR and IDL procedure SVDC.
;
;H. Rhody
;June, 2004
;-

A=double(Ain)

IF N_elements(tau) LE 0 THEN tau=1e-10
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