*! 1.0.4  18Jun1997  Jeroen Weesie/ICS  STB-39 dm49
*         code was modeled after -rank- in Matlab 4.0

program define matrank
        version 5.0

      * scratch

        tempname A U w V wmax 

      * parse input     

        parse "`*'", parse(",")
        local A0 "`1'" 
        capt mat `A' = `A0'
        if _rc {
            di in red "`A0' is undefined or an illegal matrix expression"
            exit 198    
        }
        local nrA = rowsof(`A')
        local ncA = colsof(`A')        

        mac shift
        local options "Display Rank(str) Tol(str)"
        parse "`*'"
      
      * Stata's "mat svd" requires #rows >= #cols  !!!!
        
        if `nrA' < `ncA' { 
                mat `A' = `A'' 
        }
        mat svd `U' `w' `V' = `A'
        local nc = colsof(`w')

      * tolerance for small singular values (see Matlab 4., rank.m)

        if "`tol'" == "" {
                tempname tol
                scalar `wmax' = 0
                local i 1
                while `i' <= `nc' {
                        if (`w'[1,`i'] > `wmax') { 
                                scalar `wmax' = `w'[1,`i'] 
                        }
                        local i = `i' + 1
                }
                scalar `tol' = `wmax' * `nc' * (2.22E-16) 
        } 
        else {
                confirm number `tol'
        }
        
      * rank = # sv's > tol 

        local r 0
        local i 1
        while `i' <= `nc' {
                if `w'[1,`i'] > `tol' { 
                        local r = `r' + 1 
                }
                local i = `i' + 1
        }

      * output

        if "`display'" != "" | "`rank'" == "" {
                di in gr "rank `A0'[`nrA',`ncA'] = " in ye `r'
        }
        else {
                scalar `rank' = `r'
        }
end