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

program define matorth
        version 5.0

      * scratch

        tempname A At OrthA U Ui V W 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 Format(str) Orth(str) Rank(str) Tol(str)"
        parse "`*'"

      * Stata's singular value decomposition requires "long" matrices !?!
 
        if rowsof(`A') < colsof(`A') {
                mat `At' = `A''
                mat svd `V' `W' `U' = `At'
        } 
        else    mat svd `U' `W' `V' = `A'

      * tolerance for deciding which elements of W are 0

        if "`tol'" == "" {
                * max of singular values
                local nc = colsof(`W')
                scalar `Wmax' = 0
                local i 1
                while `i' <= `nc' {
                        if (`W'[1,`i'] > `Wmax') { scalar `Wmax' = `W'[1,`i'] }
                        local i = `i' + 1
                }
                tempname tol
                scalar `tol' = `Wmax' * `nc' * (2.22E-16) 
        }
        else {
                confirm number `tol'
        }                

      * select columns of U with singular value W(i) > tol
      * note that Stata's SVD returns unsorted singular values
      
        local i 1
        local rnk 0
        while `i' <= `nc' {
                if `W'[1,`i'] > `tol' { 
                        local rnk = `rnk' + 1
                        mat `Ui' = `U'[.,`i']
                        mat `OrthA' = `OrthA' , `Ui'
                }
                local i = `i' + 1
        }                
                
      * display
      
        if "`display'" != "" | "`orth'" == "" {
                di _n in gr "Column space (Image) of `A0'[`nrA',`ncA'] " _c 
                if `rnk' == 0 {
                         di in gr " = " in ye " { 0 } "
                }
                else {
                        di in gr " has dimension " in ye `rnk' in gr " and ortho-basis"
                        if "`format'" != "" {
                                local fmt "format(`format')"
                        }                                
                        mat list `OrthA', noheader `fmt'
                }                        
        }
        
      * return results
      
        if "`orth'" != "" { 
                if `rnk' > 0 { 
                        mat `orth' = `OrthA' 
                }
                else {
                        di in bl "`orth' is -empty- hence left undefined"
                }
        }       
        if "`rank'" != "" { 
                scalar `rank' = `rnk'   
        }
end