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

program define matginv
        version 5.0
      
      * scratch

        tempname A At InvA 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 Format(str) Ginv(str) Rank(str) Tol(str)"
        parse "`*'"

      * 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 inverting the diagonal elements of W

        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'
        }                

      * invert the elements of W

        local i 1
        local rnk 0
        while `i' <= colsof(`W') {
                if `W'[1,`i'] > `tol' { 
                        mat `W'[1,`i'] = 1/(`W'[1,`i'])
                        local rnk = `rnk' + 1
                }
                else { 
                        mat `W'[1,`i'] = 0
                }
                local i = `i' + 1
        }           

      * produce InvA 

        if `rnk' == 0 {
                mat `InvA' = J(`ncA',`nrA',0)
        }
        else {      
                mat `W' = diag(`W')
                mat `InvA' = `V' * `W'
                mat `InvA' = `InvA' * `U''
        }

      * display results

        if "`ginv'" == "" | "`display'" != "" {
                if "`format'" != "" {       
                        local fmt "format(`format')" 
                }  
                di _n in gr "Moore-Penrose inverse of `A0'" /* 
                  */ "[`nrA',`ncA'] with rank = " in ye `rnk' 
                mat list `InvA', noheader `fmt'
        }
 
      * save results

        if "`rank'" != "" { scalar `rank' = `rnk' }
        if "`ginv'" != "" { matrix `ginv' = `InvA' } 
end
exit

matginv A, [ ginv(Ai) rank(r) tol(#)]

computes the Moore-Penrose inverse Ai of A, i.e., the unqiue 
matrix that satisfies the 4 conditions
 (1)  Ai*A*Ai = Ai
 (2)  A*Ai*A  = A
 (3)  Ai*Ai'  is a projection on the null space of A
 (4)  Ai'*Ai  is a projection on the image of A

if A is invertible, Ai = inverse(A)

if A is n by m, Ai = m by n

In accordance with Matlab, if A is a matrix of zeros, the 
MP-inverse is all-zeroes also.