*! 1.0.4  18Jun1997  Jeroen Weesie/ICS  STB-39 dm49

program define matnull
        version 5.0

      * scratch

        tempname A B Evec Eval Tol

      * parse input

        parse "`*'", parse(",")
        local A0 "`1'"
        capt mat `A' = `A0'
        if _rc {
                di in red "matrix `A0' is undefined/illegal matrix expression"
                exit 198
        }
        mac shift
        local options "Display Format(str) Null(str) Rank(str) Tol(str)"
        parse "`*'"
        
      * matrix to save result
      
        if "`null'" == "" {
                tempname Null
        }
        else {  
                local Null "`null'" 
        }                

      * spectral decomposition

        mat `B' = `A'' * `A'
        mat syme `Evec' `Eval' = `B' 
        local n = colsof(`Eval')
 
      * tolerance for deciding which elements of Eval are 0

        if "`tol'" == "" {
                scalar `Tol' = sqrt(`Eval'[1,1] * `n' * (2.22E-16)) 
        }
        else {
                confirm number `tol'
                scalar `Tol' = sqrt(`tol')        
        }                

      * select columns with eigenvalue Eval(i) < tol^2. Note that Stata's
      * procedure symeigen returns eigenvalues sorted in decreasing order.
      
        local r `n'
        local i 1
        while `i' <= `n' {
                if `Eval'[1,`i'] > `Tol' { 
                        local i = `i'+1 
                }
                else {
                        local r = `i'-1 
                        local i = `n'+2                      
                }                        
        }
        
      * Null space is columns r+1..n
      
        if `r' < `n' {
                local rp = `r'+1
                mat `Null' = `Evec'[.,`rp'...]
        }                

        if "`rank'" != "" { 
                sca `rank' = `r'   
        }      
        
      * display
      
        if "`display'" != "" | "`null'" == "" {
                di in gr "Null space of `A0'[" rowsof(`A') "," colsof(`A') "]" _c 
                if `r' == `n' {
                        di in gr " = " in ye "{ 0 }"
                }
                else {
                        di in gr " has dimension " in ye =(`n'-`r') in gr " and ortho-basis"
                        if "`format'" != "" {
                                local fmt "format(`format')"
                        }                                
                        mat list `Null', noheader `fmt'
                }                        
        }
end