*! version 1.2.0  18JUN2011
capture program drop dea
program define dea, rclass
    version 10.0

/** Terms Description:
 * -----------------------------------------------------------------------------
 * RTS: Return To Scale
 * CRS: Constant Return to Scale
 * VRS: Variant Return to Scale
 * IRS: Increasing Return to Scale
 * DRS: Decreasing Retrn to Scale
 * -----------------------------------------------------------------------------
 */

// syntax checking and validation-----------------------------------------------
// rts - return to scale, ort - orientation
// input varlist = output varlist
// example:
//     dea Employee Area = Sales Profits, rts(CRS) ort(IN) tol1(1e-14)
//     dea Employee Area = Sales Profits, sav
// -----------------------------------------------------------------------------
    // returns 1 if the first nonblank character of local macro `0' is a comma,
    // or if `0' is empty.
	if replay() {
        dis as err "ivars and ovars required."
        exit 198
    }

    // get and check invarnames
    gettoken word 0 : 0, parse(" =:,")
    while `"`word'"' != ":" & `"`word'"' != "=" {
        if `"`word'"' == "," | `"`word'"'=="" {
                error 198
        }
        local invarnames `invarnames' `word'
        gettoken word 0 : 0, parse(" :,")
    }
    unab invarnames : `invarnames'

    #del ;
    syntax varlist(min=1) [if] [in] [using/]
    [,
        RTS(string)         // ignore case sensitive,[{CRS|CCR}|{BCC|VRS}|DRS|IRS]
        ORT(string)         // ignore case sensitive,[{IN|INPUT}|{OUT|OUTPUT}]
        STAGE(integer 2)    // stage 1 or 2
        TOL1(real 1e-14)    // entering or leaving value tolerance
        TOL2(real 1e-8) 	// B inverse tolerance: 2.22e-12
        MINSUBScript        // minimal subscript pivot
        TRACE               // Whether or not to do the log
        SAVing(string)      // log file name
        REPLACE             // Whether or not to replace the log file
    ];
    #del cr

    local num_invar : word count `invarnames'
    local i 1
    while (`i'<=`num_invar') {
        local invarn : word `i' of `invarnames'
        local junk : subinstr local invarnames "`invarn'" "", ///
            word all count(local j)
        if `j' > 1 {
            di as error ///
                "cannot specify the same input variable more than once"
            error 498
        }
        local i = `i' + 1
    }

    // default model - CRS(Constant Return to Scale)
    if ("`rts'" == "") local rts = "CRS"
    else {
        local rts = upper("`rts'")
        if ("`rts'" == "CCR") local rts = "CRS"
        else if ("`rts'" == "BCC") local rts = "VRS"
        else if (~("`rts'" == "CRS" | "`rts'" == "VRS" | ///
                   "`rts'" == "DRS" | "`rts'" == "IRS")) {
            di as err "option rts allow for case-insensitive " _c
            di as err "CRS (eq CCR) or VRS (eq BCC) or DRS or IRS or nothing."
            exit 198
        }
    }

    // default orientation - Input Oriented
    if ("`ort'" == "") local ort = "IN"
    else {
        local ort = upper("`ort'")
        if ("`ort'" == "I" | "`ort'" == "IN" | "`ort'" == "INPUT") {
            local ort = "IN"
        }
        else if ("`ort'" == "O" | "`ort'" == "OUT" | "`ort'" == "OUTPUT") {
            local ort = "OUT"
        }
        else {
            di as err "option ort allowed for case-insensitive " _c
            di as err "(i|in|input|o|out|output) or nothing."
            exit 198
        }
    }

    // default stage - 1
    if (~("`stage'" == "1" | "`stage'" == "2")) {
        dis as err "option stage allows 1 or 2 or nothing exclusively."
        exit 198
    }

    if ("`using'" != "") use "`using'", clear
    if (~(`c(N)' > 0 & `c(k)' > 0)) {
        dis as err "dataset required!"
        exit 198
    }

// end of syntax checking and validation ---------------------------------------

    set more off
    capture log close dealog
    log using "dea.log", replace text name(dealog)
    preserve

    if ("`if'" != "" | "`in'" != "") {
        qui keep `in' `if'  // filtering : keep in range [if exp]
    }

    deanormal, ivars(`invarnames') ovars(`varlist') ///
               rts(`rts') ort(`ort') stage(`stage') ///
               tol1(`tol1') tol2(`tol2') `minsubscript' ///
               `trace' saving(`saving') `replace'
    return add

    restore, preserve
    log close dealog
end

********************************************************************************
* DEA Normal - Data Envelopment Analysis Normal
********************************************************************************
program define deanormal, rclass
    #del ;
    syntax , IVARS(string) OVARS(string) RTS(string) ORT(string)
    [
        STAGE(integer 2) TOL1(real 1e-3) TOL2(real 1e-12) MINSUBScript
        TRACE SAVing(string) REPLACE
    ];
    #del cr

    preserve
    // di _n as input ///
    //     "options: RTS(`rts') ORT(`ort') STAGE(`stage') SAVing(`saving')"
    // di as input "Input Data:"
    // list

    // -------------------------------------------------------------------------

    tempname dmuIn dmuOut frameMat deamainrslt dearslt vrsfrontier crslambda
    mkDmuMat `ivars', dmumat(`dmuIn') sprefix("i")
    mkDmuMat `ovars', dmumat(`dmuOut') sprefix("o")
    local dmuCount = colsof(`dmuIn')
    if ("`minsubscript'" == "") local minsubscript = "nominsubscript"
    if ("`ort'" == "OUT") local minrank = 1
    else local minrank = 0

    mata: _mkframemat("`frameMat'", "`dmuIn'", "`dmuOut'", ///
            "`rts'", "`ort'", `stage')
    deamain `dmuIn' `dmuOut' `frameMat' `rts' `ort' `stage' ///
            `tol1' `tol2' `minsubscript' `trace'
    matrix `deamainrslt' = r(deamainrslt)

    mata: _dmurank("`deamainrslt'", ///
            `=rowsof(`dmuIn')', `=rowsof(`dmuOut')', `minrank', `tol2')
    matrix `dearslt' = r(rank), `deamainrslt'

    // use mata function '_setup_dearslt_names' because the maximum string
    // variable length needs to be kept under the 244 for all the time
    mata: _setup_dearslt_names("`dearslt'", "`dmuIn'", "`dmuOut'")

    if ("`rts'" == "VRS") {
        // caution : join order! (CRS -> VRS -> DRS)
        // 1. VRS TE(VRS Technical Efficiency)
        matrix `deamainrslt' = r(deamainrslt)
        matrix `vrsfrontier' = `deamainrslt'[1...,1]

        // 2. CRS TE(CRS Technical Efficiency)
        local stageForVRS = 2
        mata: _mkframemat("`frameMat'", "`dmuIn'", "`dmuOut'", ///
                "CRS", "`ort'", `stageForVRS')
        deamain `dmuIn' `dmuOut' `frameMat' "CRS" `ort' `stageForVRS'  ///
                `tol1' `tol2' `minsubscript' `trace'
        matrix `deamainrslt' = r(deamainrslt)
        matrix `vrsfrontier' = `deamainrslt'[1...,1], `vrsfrontier'
        matrix `crslambda' = `deamainrslt'[1...,2..(`dmuCount' + 1)]

        // 3. DRS TE(DRS Technical Efficiency)
        local stageForVRS = 1
        mata: _mkframemat("`frameMat'", "`dmuIn'", "`dmuOut'", ///
                "DRS", "`ort'", `stageForVRS')
        deamain `dmuIn' `dmuOut' `frameMat' "DRS" `ort' `stageForVRS' ///
                `tol1' `tol2' `minsubscript' `trace'
        matrix `deamainrslt' = r(deamainrslt)
        matrix `vrsfrontier' = `vrsfrontier', `deamainrslt'[1...,1] // VRS
        mata: _roundmat("`vrsfrontier'", 1e-14) // round off

        //
        matrix `vrsfrontier' = `vrsfrontier', J(rowsof(`vrsfrontier'), 2, 0)
        matrix rownames `vrsfrontier' = `: colfullnames `dmuIn''
        matrix colnames `vrsfrontier' = CRS_TE VRS_TE DRS_TE SCALE RTS

        if ("`trace'" == "trace") matrix list `crslambda'
        forvalues i = 1/`=rowsof(`vrsfrontier')' {
            matrix `vrsfrontier'[`i',4] = /*
                */ float(`vrsfrontier'[`i',1]/`vrsfrontier'[`i',2])

            /*******************************************************************
             * if CRS(CCR) == VRS(BCC) then CRS
             * else
             *     if the sum of ¥ë equals to 1, then mark CRS
             *     if sum of ¥ë is greater than 1, then mark DRS
             *     if sum of ¥ë is less than 1, then mark IRS
             ******************************************************************/
            // (-1:drs, 0:crs, 1:irs)
            if (`vrsfrontier'[`i',1] == `vrsfrontier'[`i',2]) {
                matrix `vrsfrontier'[`i',5] = 0
            }
            else {
                local sumlambda = 0
                forvalues j = 1/`dmuCount' {
                    if (`crslambda'[`i', `j'] < .) {
                        local sumlambda = `sumlambda' + `crslambda'[`i', `j']
                    }
                }
                if (`sumlambda' < 1) {
                    matrix `vrsfrontier'[`i',5] = 1 // irs mark
                }
                else { // if (`sumlambda' >= 1) {
                    matrix `vrsfrontier'[`i',5] = -1 // drs mark
                }
            }
        }
    }

    // -------------------------------------------------------------------------
    // REPORT
    // -------------------------------------------------------------------------
    di as result ""
    di as input "options: RTS(`rts') ORT(`ort') STAGE(`stage')"
    di as result "`rts'-`ort'PUT Oriented DEA Efficiency Results:"
    matrix list `dearslt', noblank nohalf noheader f(%9.6g)

    if ("`saving'" != "") {
        // if the saved file exists and doesn't replace, make the backup file.
        if ("`replace'" == "") {
            local dotpos = strpos("`saving'",".")
            if (`dotpos' > 0) {
                mata: _file_exists("`saving'")
            }
            else {
                mata: _file_exists("`saving'.dta")
            }
            if r(fileexists) {
                local curdt = subinstr("`c(current_date)'", " ", "", .) + /*
                    */ subinstr("`c(current_time)'", ":", "", .)
                if (`dotpos' > 0) {
                    #del ;
                    local savefn = substr("`saving'", 1, `dotpos' - 1) +
                                   "_bak_`curdt'" +
                                   substr("`saving'",`dotpos', .);
                    #del cr
                    qui copy "`saving'" "`savefn'", replace
                }
                else {
                    local savefn = "`saving'_bak_`curdt'" + ".dta"
                    qui copy "`saving'.dta" "`savefn'", replace
                }
            }
        }

        if ("`rts'" != "VRS") {
            restore, preserve
            svmat `dearslt', names(eqcol)
            capture {
                renpfix _
                renpfix ref ref_
                renpfix islack is_
                renpfix oslack os_
            }
            capture save `saving', replace

            // di as result ""
            // di as result "DEA Result file:"
            // list
        }
    }
    return matrix dearslt = `dearslt'

    if ("`rts'" == "VRS") {
        if ("`trace'" == "trace") {
            di _n(2) as result "CRS lambda(¥ë)"
            matrix list `crslambda' , noblank nohalf noheader f(%9.6f)
        }

        // to verify: if you don't have to verify, you can comment this sentence.
        di _n(2) as result "VRS Frontier(-1:drs, 0:crs, 1:irs)"
        matrix list `vrsfrontier', noblank nohalf noheader f(%9.6f)

        di _n(2) as result "VRS Frontier:"
        restore, preserve
        svmat `vrsfrontier', names(col)
        rename RTS RTSNUM
        qui generate RTS = "drs" if RTSNUM == -1
        qui replace RTS = "-" if RTSNUM == 0
        qui replace RTS = "irs" if RTSNUM == 1
        drop DRS_TE RTSNUM
        format CRS_TE VRS_TE SCALE %9.6f
        list
        if ("`saving'"!="") {
            capture save `saving', replace
        }
    }
    restore, preserve
end

********************************************************************************
* DEA Main - Data Envelopment Analysis Main
********************************************************************************
program define deamain, rclass
    args dmuIn dmuOut frameMat rts ort stage tol1 tol2 minsubscript trace

    tempname efficientVec deamainrslt

    // stage step 1.
    if ("`trace'" == "trace") {
        di _n(2) as txt "RTS(`rts') ORT(`ort') 1st stage."
    }
    mata: _dealp("`frameMat'", "`dmuIn'", "`dmuOut'", "`rts'", "`ort'", ///
        1, `tol1', `tol2', "`minsubscript'", "`efficientVec'", "`trace'")
    matrix `deamainrslt' = r(dealprslt)

    // stage step 2.
    if ("`stage'" == "2") {
        if ("`trace'" == "trace") {
            di _n(2) as txt "RTS(`rts') ORT(`ort') 2nd stage."
        }
        matrix `efficientVec' = `deamainrslt'[1...,1]

        mata: _dealp("`frameMat'", "`dmuIn'", "`dmuOut'", "`rts'", "`ort'", ///
            2, `tol1', `tol2', "`minsubscript'", "`efficientVec'", "`trace'")
        matrix `deamainrslt' = r(dealprslt)
    }

    return matrix deamainrslt = `deamainrslt'
end

********************************************************************************
* Data Import and Conversion
********************************************************************************

// Make DMU Matrix -------------------------------------------------------------
program define mkDmuMat
    #del ;
    syntax varlist(numeric) [if] [in], DMUmat(name)
    [
        SPREfix(string)
    ];
    #del cr

    qui ds
    // variable not found error
    if ("`varlist'" == "") {
        di as err "variable not found"
        exit 111
    }

    // make matrix
    mkmat `varlist' `if' `in', matrix(`dmumat') rownames(dmu)
    matrix roweq `dmumat' = "dmu"
    matrix coleq `dmumat' = `=lower("`sprefix'") + "slack"'
    matrix `dmumat' = `dmumat''
end

// Start of the MATA Definition Area -------------------------------------------

version 10
mata:
mata set matastrict on

/**
 * Declare the variable's boundary condition structure.
 */
struct BoundCond {
    real scalar val, lower, upper, free
}

/**
 * Declare the LP's parameter structure.
 */
struct LpParamStruct {
    string scalar rts       // return to scale(CRS|VRS|IRS|DRS)
	real scalar isin        // if 1 then 'in', other then 'out'
	real scalar stagestep   // stage step. 1 or 2

    real scalar dmus        // number of dmus
    real scalar slacks      // number of slacks
	real scalar artificials // number of artificials

	real scalar tol1        // tolerance 1
    real scalar tol2        // tolerance 2
    real scalar isminsubscript // whether min subscript or not.
    string scalar trace     // whether trace or not.
}

/**
 * Declare the LP's result structure.
 */
struct LpResultStruct {
	real scalar xVal	// value of objective.
	real matrix XB		// basic feasible solution.
	real scalar rc		// return code(zero means success)
	string scalar rmsg	// return message
}

/**
 * make frame matrix and set matrix value at the param frameMat
 * rts - return to scale, ort - orientation
 */
function _mkframemat( string scalar frameMat,
                      string scalar dmuIn,
                      string scalar dmuOut,
                      string scalar rts,
                      string scalar ort,
                      real scalar stage )
{
    real matrix F, DI, DO
    real scalar row, col, sig, isin
    real scalar dmus, slackins, slackouts, slacks
    real scalar frows, fcols

    DI = st_matrix(dmuIn)
    DO = st_matrix(dmuOut)
    if (cols(DI) != cols(DO)) {
        _error(3200, "in and out count of dmu is not match!")
    }

    // basic value setting for artificial variabels
    isin = (ort == "IN"); sig = (isin ? -1 : 1)

    dmus = cols(DI) // or cols(DO), because cols(DI) == cols(DO)
    slackins = rows(DI); slackouts = rows(DO)
    if (rts == "CRS") {
		slacks = slackins + slackouts
	    // target coefficient\slackins\slackouts
        frows = 1 + slackins + slackouts
		// target coefficient,theta,dmus,slackins,slackouts,rhs
        fcols = 1 + 1 + dmus + slackins + slackouts + 1
    }
    else if (rts == "VRS") {
		slacks = slackins + slackouts
		// target coefficient\slackins\slackouts\sum of lamda
        frows = 1 + slackins + slackouts + 1
		// target coefficient,theta,dmus,slackins,slackouts,rhs
        fcols = 1 + 1 + dmus + slackins + slackouts + 1
    }
	else if (rts == "IRS") {
        slacks = slackins + slackouts + 1
        // target coefficient\slackins\slackouts\sum of lamda
        frows = 1 + slackins + slackouts + 1
		// target coefficient,theta,dmus,slackins,slackouts,sum of lamda,rhs
        fcols = 1 + 1 + dmus + slackins + slackouts + 1 + 1
    }
    else if (rts == "DRS") {
        slacks = slackins + slackouts + 1
        // target coefficient\slackins\slackouts\sum of lamda
        frows = 1 + slackins + slackouts + 1
		// target coefficient,theta,dmus,slackins,slackouts,sum of lamda,rhs
        fcols = 1 + 1 + dmus + slackins + slackouts + 1 + 1
    }
    else {
        _error(3498, "invalid rts optoin.")
    }

    // make frame matrix for CRS(CCR)
    F = J(frows, fcols, 0)
    F[1, 1] = 1
    replacesubmat(F, 2, 3, sig * DI)
    replacesubmat(F, 2 + slackins, 3, -sig * DO)
    replacesubmat(F, 2, 3 + dmus, sig * I(slacks))

    // adjustment
	if (rts == "VRS") {
        replacesubmat(F, frows, 3, J(1, dmus, 1))
        F[frows,fcols] = 1
    }
	else if (rts == "IRS") {
        replacesubmat(F, frows, 3, J(1, dmus, 1))
        F[frows,2 + dmus + slacks] = -1
        F[frows,fcols] = 1
    }
    else if (rts == "DRS") {
        replacesubmat(F, frows, 3, J(1, dmus, 1))
        F[frows,2 + dmus + slacks] = 1
        F[frows,fcols] = 1
    }

    // return result
    st_matrix(frameMat, F)
}

/**
 * DEA Loop - Data Envelopment Analysis Loop for DMUs
 */
function _dealp ( string scalar frameMat,
                  string scalar dmuIn,
                  string scalar dmuOut,
                  string scalar rts,
                  string scalar ort,
                  real scalar stagestep,
                  real scalar tol1,
                  real scalar tol2,
                  string scalar minsubscript,
                  string scalar efficientVec,
                  string scalar trace )
{
    real matrix F, DI, DO, M, VARS, LPRSLT, DEALPRSLT, ARTIF
    real scalar dmus, slackins, slackouts, slacks, artificials, artificialrow
    real scalar frows, fcols, isin, isminsubscript
    real colvector effvec, slackidx, skipdmu
    string scalar tracename

	struct BoundCond matrix boundF, boundM
	struct LpParamStruct scalar param

    F  = st_matrix(frameMat)
    DI = st_matrix(dmuIn)
    DO = st_matrix(dmuOut)
    if (cols(DI) != cols(DO)) {
        _error(3200, "in and out count of dmu is not match!")
    }
	if (!(rts == "CRS" || rts == "VRS" || rts == "IRS" || rts == "DRS")) {
		_error(3498, "rts must be one of CRS, VRS, IRS, DRS")
	}

    // basic value setting for artificial variabels
    isin = (ort == "IN")
    isminsubscript = (minsubscript != "nominsubscript")
    frows = rows(F); fcols = cols(F)
    dmus = cols(DI) // or cols(DO), because cols(DI) == cols(DO)
    slackins = rows(DI); slackouts = rows(DO)

    tracename = rts + "-" + (isin ? "IN" : "OUT")
    tracename = tracename + "-" + (stagestep == 1 ? "SI" : "SII")

	// -------------------------------------------------------------------------
	// define number of slacks by rts
	if (rts == "CRS" || rts == "VRS") slacks = slackins + slackouts
	else if (rts == "IRS" || rts == "DRS") slacks = slackins + slackouts + 1

	// define number of artificials by rts, ort, stage
	if (rts == "CRS" || rts == "DRS") {
		if (stagestep == 1) {
			if (isin) {
				artificials = slackins+slackouts; artificialrow = 2;
			}
			else artificials = 0
		}
		else {
			artificials = slackouts; artificialrow = 2+slackins;
		}
	}
	else if (rts == "VRS" || rts == "IRS") {
		if (stagestep == 1) {
			if (isin) {
				artificials = slackins+slackouts+1; artificialrow = 2;
			}
			else {
				artificials = 1; artificialrow = frows //== 2+slackins+slackouts
			}
		}
		else {
			artificials = slackouts+1; artificialrow = 2+slackins
		}
	}
	if (artificials > 0) {
		ARTIF = J(1, artificials, 1) \ J(frows-1, artificials, 0)
		replacesubmat(ARTIF, artificialrow, 1, I(artificials))
		F = F[,1..fcols-1], ARTIF, F[,fcols]
		frows = rows(F); fcols = cols(F) // revise frows, fcols
	}
	// -------------------------------------------------------------------------

	// put artificial for slackins on inactive status
	if (stagestep == 2) {
	    effvec = st_matrix(efficientVec)
		skipdmu = (effvec :== .)
		if (isin) {
		    replacesubmat(F, 2, 3, -F[2..1+slackins,3::2+dmus+slackins])
		}
		else {
		    replacesubmat(F, 2+slackins, 3,
                -F[2+slackins..1+slackins+slackouts,3::2+dmus+slacks])
		}
	}
	else skipdmu = J(1, dmus, 0)
	// -------------------------------------------------------------------------
    boundF = J(1, fcols, BoundCond());
	// set the boundary for the efficiency variable(theta, eta):
	// -INFINITE <= efficiency <= INFINITE
	boundF[1,2].val = 0; boundF[1,2].lower = 0; boundF[1,2].upper = .

	// set boundary value for the weight variable(lamda, mu):
	// 0 <= weight <= INFINITE
	for (i=3; i<dmus+3; i++) {
		boundF[1,i].val = 0; boundF[1,i].lower = 0; boundF[1,i].upper = .
	}

	// set boundary value for the non-structural variable(slack, artificial).
	// 0 <= slacks and atrificials <= INFINITE
	for (i=dmus+3; i<fcols; i++) {
		boundF[1,i].val = 0; boundF[1,i].lower = 0; boundF[1,i].upper = .
	}
	// liststruct(boundF); // for debug

	// set the lp's parameters
	param.rts            = rts
	param.isin           = isin
	param.stagestep      = stagestep
	param.dmus           = dmus
	param.slacks         = slacks
	param.artificials    = artificials
	param.tol1           = tol1
	param.tol2           = tol2
	param.isminsubscript = isminsubscript
	param.trace          = trace
	// liststruct(param); // for debug
	// -------------------------------------------------------------------------
    DEALPRSLT = J(0, 1+ dmus + slacks, 0)
    if (isin) {
        for (dmui=1; dmui<=dmus; dmui++) {
			if (skipdmu[dmui]) {
				LPRSLT = J(1, cols(DEALPRSLT), .)
			}
			else {
				M = F; boundM = boundF
				if (stagestep == 1) replacesubmat(M, 2, 2, DI[,dmui])
				else replacesubmat(M, 2, fcols, DI[,dmui]*effvec[dmui])
				replacesubmat(M, 2+slackins, fcols, DO[,dmui])

				// execute LP
				VARS   = lp_phase1(M, boundM, dmui, tracename, param)
				if (VARS[1,1] == .) {
					LPRSLT = J(1, cols(DEALPRSLT), .)
				}
				else {
					LPRSLT = lp_phase2(M, boundM, VARS, dmui, tracename, param);
				}
			}
            DEALPRSLT = DEALPRSLT \ LPRSLT
        }
    }
    else {
        for (dmui=1; dmui<=dmus; dmui++) {
			if (skipdmu[dmui]) {
				LPRSLT = J(1, cols(DEALPRSLT), .)
			}
			else {
				M = F; boundM = boundF
				replacesubmat(M, 2, fcols, DI[,dmui])
				if (stagestep == 1) {
					if (rts == "CRS" || rts == "DRS") M[1,2] = -1
					replacesubmat(M, 2+slackins, 2, DO[,dmui])
				}
				else replacesubmat(M, 2+slackins, fcols, DO[,dmui]*effvec[dmui])

				// execute LP
				if (artificials == 0) { // if artificials == 0 then skip phase 1
					VARS   = (0, 2+dmus..1+dmus+slacks, 1..1+dmus, 0)
					M = M[,1],
						M[,VARS[,2::cols(VARS)-1] :+ 1],
						M[,cols(M)]
					boundM = boundM[,1],
							 boundM[,VARS[,2::cols(VARS)-1] :+ 1],
							 boundM[,cols(M)]
				}
				else {
					VARS = lp_phase1(M, boundM, dmui, tracename, param)
				}

				if (VARS[1,1] == .) {
					LPRSLT = J(1, cols(DEALPRSLT), .)
				}
				else {
					LPRSLT = lp_phase2(M, boundM, VARS, dmui, tracename, param);
				}
			}
            DEALPRSLT = DEALPRSLT \ LPRSLT
        }
    }

    // adjust efficiency
    if (stagestep == 2) {
        replacesubmat(DEALPRSLT, 1, 1, effvec)
    }
    st_matrix("r(dealprslt)", DEALPRSLT)
}

real matrix function lp_phase1 ( real matrix M,
                                 struct BoundCond matrix boundM,
                                 real scalar dmui,
                                 string scalar aTracename,
                                 struct LpParamStruct scalar param )
{
    real matrix T, VARS
    real scalar mrows, mcols, phase
	real vector reorderidx, bfsidx, nonbfsidx
    string scalar tracename, msg
	struct LpResultStruct scalar lpresult

    mrows = rows(M); mcols = cols(M)
    tracename = aTracename + "-PI"

    // 1st: initialize matrix.
    if (param.trace == "trace") {
        displayas("txt")
        printf("\n\n\n----------[PHASE I]----------")
        printf("\n[DMUi=%g]%s: initialize matrix.\n",
            dmui, tracename); M
    }

    // 2nd: classify basic and nonbasic.
	VARS = (0, 1..1+param.dmus+param.slacks, -1..-param.artificials, 0)
	bfsidx = J(1, mrows-1, .); nonbfsidx = J(1, 0, .)
	for (j = 3+param.dmus; j <= mcols-1; j++) {
		T = (M[2::mrows,j] :== 1)
		if (any(T)) {
		    maxindex(T, 1, i, w); bfsidx[i] = j
		}
		else nonbfsidx = nonbfsidx, j
	}
	reorderidx = (1, bfsidx[1,], 2..2+param.dmus, nonbfsidx[1,], mcols)
	VARS = VARS[,reorderidx];
	M = M[,reorderidx]; boundM = boundM[,reorderidx]

    if (param.trace == "trace") {
        displayas("txt")
        printf("\n[DMUi=%g]%s: classify basic and nonbasic.\n",
            dmui, tracename); M; VARS
    }

    // 3rd: slove the linear programming(LP).
	phase = 1
    lpresult = lp(M, boundM, VARS, dmui, phase, tracename, param)

	if(lpresult.rc) VARS[1,1] = .
    return(VARS)
}

real matrix function lp_phase2 ( real matrix M,
                                 struct BoundCond matrix boundM,
                                 real matrix VARS,
                                 real scalar dmui,
                                 string scalar aTracename,
                                 struct LpParamStruct scalar param )
{
    real matrix T, XB, orgVARS, LPRSLT
    real scalar mrows, mcols, realslacks
	real vector slackidx
    string scalar tracename, msg
	struct LpResultStruct scalar lpresult

    orgVARS = VARS
    mrows = rows(M); mcols = cols(M)

    tracename = aTracename + "-PII"

	// modify target function value:
	M[1,] = J(1,mcols,0); M[1,1] = 1
	if (param.stagestep == 1) { // X = theta
		for (j=2; j<mcols; j++) {
			if (VARS[1,j] == 1) M[1,j] = 1 // because of theta index == 1
		}
	}
	else if (param.stagestep == 2) { // X = S1 + S2 + ... + Sn
	    realslacks = (param.rts == "IRS" || param.rts == "DRS") ?
		        param.slacks-1 : param.slacks;
		slackidx = (2+param.dmus..1+param.dmus+realslacks)
		for (j=2; j<mcols; j++) {
			for (i=1; i<=realslacks; i++) {
				if (VARS[1,j] == slackidx[i] && !allof(M[,j], 0)) M[1,j] = 1
			}
		}
	}

    if (param.trace == "trace") {
        displayas("txt")
        printf("\n----------[PHASE II]----------")
        printf("\n[DMUi=%g]%s: initialize matrix.\n",
            dmui, tracename); M
        printf("\n[DMUi=%g]%s: VARS.\n", dmui, tracename); VARS
    }

    phase = 2
    lpresult = lp(M, boundM, VARS, dmui, phase, tracename, param)

    // -------------------------------------------------------------------------
    // phase 2 final.
    // -------------------------------------------------------------------------
	if(lpresult.rc) {
		LPRSLT = J(1, 1+param.dmus+param.slacks, .)
	}
	else {
		// lpresult = theta(1) + dmus + slacks
		LPRSLT = J(1, 1+param.dmus+param.slacks, 0)
		for (j=1; j<=rows(lpresult.XB) ; j++) {
			if (VARS[1,j+1] > 0) LPRSLT[1, VARS[1,j+1]] =lpresult.XB[j, 1]
		}
		if (param.stagestep == 1 && LPRSLT[1, 1] <= 0) {
			LPRSLT[1, 1] = lpresult.xVal
		}
	}

    if (param.trace == "trace") {
        msg = sprintf("[DMUi=%g]%s-FINAL", dmui, tracename);
        printf("\n%s: original VARS.\n", msg); orgVARS
        printf("\n%s: VARS.\n", msg); VARS
        printf("\n%s: XB.\n", msg); lpresult.XB
        printf("\n%s: LPRSLT.\n", msg); LPRSLT
    }

    return(LPRSLT)
}

/**
 * Declare the LP's tableau structure.
 */
struct LpTableauStruct {
	pointer(real matrix) scalar CB, CNj
	pointer(real matrix) scalar  B,  Nj, b
	pointer(real matrix) scalar Bi,CBBi, rawXB, XB
}

/**
 * return 0: sucess
 * return 1: B inverse error
 * return 2: XB has negative value.
 */
real scalar function decompsition(real matrix M,
								struct BoundCond matrix boundM,
								real scalar mrows,
								real scalar mcols,
								real scalar slacks,
								struct LpTableauStruct scalar tbl,
								struct LpParamStruct scalar param )
{
	real matrix CB, CNj
	real matrix  B,  Nj, b
	real matrix Bi,CBBi, rawXB, XB, BiNjXj
	real scalar Njcols, result

	// set the tableau.
	tbl.CB = &CB; tbl.CNj  = &CNj;
	tbl.B  = &B;  tbl.Nj   = &Nj;   tbl.b = &b
	tbl.Bi = &Bi; tbl.CBBi = &CBBi
	tbl.rawXB = &rawXB; tbl.XB = &XB;

	CB  = M[1,2::1+slacks]; 		CNj = M[1,2+slacks::mcols-1]
	B   = M[2..mrows,2::1+slacks];  Nj  = M[2..mrows,2+slacks::mcols-1]
	b   = M[2..mrows,mcols]

	Bi = lusolve(B, I(rows(B)), 1e-14)
	if (any(Bi :== .)) { // B is singular matrix.
		return (result = 1);
		// Bi = svsolve(B, I(rows(B)), 1e-14)
		// if (any(Bi :== .)) return (result = 1);
	}

	CBBi = CB*Bi
	// BFS(basic feasible solution)
	Njcols = cols(Nj); BiNjXj = J(rows(Nj), 1, 0)
	for (j=1; j<=Njcols; j++) {
		BiNjXj = BiNjXj :+ (Bi*Nj[.,1]*(boundM[1,1+slacks+j].val))
	}
	rawXB = Bi*b - BiNjXj
	XB = edittozerotol(rawXB, param.tol2) // BFS(basic feasible solution)
	// BFS(basic feasible solution) must be nonnegative.
	if (any(XB :< 0)) return (result = 2);

	return (result = 0); // sucess
}

struct LpResultStruct function lp ( real matrix M,
									struct BoundCond matrix boundM,
									real matrix VARS,
									real scalar dmui,
									real scalar phase,
									string scalar tracename,
									struct LpParamStruct scalar param )
{
    real matrix B, CB, Bi, b, XB, rawXB, BiNjXj, CBBi, Nj, CNj, Aj, alpha
    real matrix T, TH1, TH2, valT, lowerT, upperT, LVi, V
    real scalar mrows, mcols, enteringVar, leavingVar, calci, maxiter
    real scalar existArtificial, xVal, alphaVal
	real scalar minYn, tcols, tempVal, minVal, maxVal
	struct BoundCond matrix boundT
	struct LpResultStruct scalar lpresult
	struct LpTableauStruct scalar tbl

	real colvector enterings, leavings
	real scalar    enteringi, leavingi

	real scalar slacks, isin, tol1, tol2, isminsubscript
	string scalar trace

	// -------------------------------------------------------------------------
	slacks = rows(M) - 1 // number of basic feasible solution
	isin = param.isin
	tol1 = param.tol1; tol2 = param.tol2
	isminsubscript = param.isminsubscript
	trace = param.trace
    minYn = 0
	if (param.stagestep == 1) minYn = (phase == 1) ? 1 : isin
	else minYn = (phase == 1)
    // -------------------------------------------------------------------------
	mrows = rows(M); mcols = cols(M)
    LVi = J(slacks, 1, .) // leaving variable index matrix.
    // -------------------------------------------------------------------------
if (trace == "trace") {
    displayas("txt"); msg = "initial tableau in the LP."
    printf("\n[DMUi=%g]%s: %s\n", dmui, tracename, msg); M
}
	lpresult.rc = 0; lpresult.rmsg = ""
    existArtificial = (phase == 2 && any(VARS[,2::1+slacks] :< 0));
    maxiter = st_numscalar("c(maxiter)")
    for (calci=1 ; calci<=maxiter ; calci++) { // prevent infinite loop

if (trace == "trace") {
	printf("\n[DMUi=%g]%s-LOOP[%g] Start...\n", dmui, tracename, calci)
}
        B  = M[2..mrows,2::1+slacks];       CB  = M[1,2::1+slacks]
        Nj = M[2..mrows,2+slacks::mcols-1]; CNj = M[1,2+slacks::mcols-1]
		b  = M[2..mrows,mcols]

		Bi = lusolve(B, I(rows(B)), 1e-14)
		if (any(Bi :== .)) { // B is singular matrix.
		    Bi = svsolve(B, I(rows(B)), 1e-14)
			if (any(Bi :== .)) {
				lpresult.rc = 3498;
				lpresult.rmsg = sprintf("%s[DMUi=%g][LOOP=%g]%s",
					"No Solution(BFS's inverse is not exist):",
					dmui, calci, tracename)
				break;

				/*
				display("B:");B
				display("rank(B) : det(B)"); rank(B), det(B)
				_error(3498, "No Solution(BFS's inverse is not exist):"
						+ "[DMUi=" + strofreal(dmui) + "]"
						+ "[LOOP=" + strofreal(calci) + "]"
						+ tracename)
				*/
			}
		}
		CBBi = CB*Bi

		// BFS(basic feasible solution)
		Njcols = cols(Nj); BiNjXj = J(rows(Nj), 1, 0)
		for (j=1; j<=Njcols; j++) {
			BiNjXj = BiNjXj :+ (Bi*Nj[.,1]*(boundM[1,1+slacks+j].val))
		}
		rawXB = Bi*b - BiNjXj
		XB = edittozerotol(rawXB, tol2) // BFS(basic feasible solution)

		if (any(XB :== .)) {
			lpresult.rc = 3498;
			lpresult.rmsg = sprintf("%s[DMUi=%g][LOOP=%g]%s",
				"No Solution(XB contains missing value):",
				dmui, calci, tracename)
			break;
		}

		// BFS(basic feasible solution) must be nonnegative.

		boundT = boundM[1,2+slacks::mcols-1]
		tcols = cols(boundT); valT = J(1,tcols,.)
		for(j=1; j<=tcols; j++) {
			valT[1,j] = boundT[1,j].val
		}
		xVal = edittozerotol((CB*rawXB + CNj*valT'), tol2) // objective funtion value

        T = VARS[1,2::1+rows(XB)] // rows(XB) equal to number of slacks.

if (trace == "trace") {
    printf("\n[DMUi=%g]%s-LOOP[%g]: CBBi * b = %g\n",
        dmui, tracename, calci, xVal);
    display("Entered index(if value is negative, that's artificial):"); T;
    display("XB = Bi*b - BiNjXj:"); rawXB;
}

        // ---------------------------------------------------------------------
        // loop terminated condition.
        // ---------------------------------------------------------------------
        if (phase == 1) {
            // objective function value is zero
			if (xVal == 0 ) {
				// if all artificals are out or remaining artificals are at zero,
				// stop and go phaseII
				T = (T :< 0) // artificial remain or not?
				if (allof(T, 0) || allof(select(XB, T'), 0)) break;

				// If remaining artificails are not zero, No Solution.
				lpresult.rc = 3498
				lpresult.rmsg = sprintf("%s[DMUi=%g][LOOP=%g]%s",
					"No Solution(Remaining artificails are not zero):",
					dmui, calci, tracename)
				break;
			}
		}

        // ---------------------------------------------------------------------
        // Select entering variable.
        // ---------------------------------------------------------------------
        enteringVar = 0; tempVal = 0; minVal = 0; maxVal = 0; boundi = 0
        Njcols = cols(Nj); T = J(1, Njcols, .)
        for (j=1; j<=Njcols; j++) {
            tempVal = CBBi * Nj[,j] - CNj[1,j]
			// if (abs(tempVal) < tol1) continue;

			boundi = 1 + slacks + j
			if (boundM[1,boundi].val == boundM[1,boundi].lower) { // lower bound
				T[1,j] = tempVal
			}
			else { // upper bound
				T[1,j] = -tempVal
			}
        } // end of for

		if (!minYn) { // maximization.
			T = T :/ (T :< 0)
			if (!allof(T, .)) {
				minindex(T, 1, enterings, w); enteringi = 1
				enteringVar = enterings[enteringi]
				evj = 1+slacks+enteringVar
			}
		}
		else { // minimization.
			T = T :/ (T :> 0)
			if (!allof(T, .)) {
				maxindex(T, 1, enterings, w); enteringi = 1
				enteringVar = enterings[enteringi]
				evj = 1+slacks+enteringVar
			}
		}

        // No more candidate for entering variable.
        if (enteringVar == 0) {
			if (trace == "trace") {
                printf("\n[DMUi=%g]%s-LOOP[%g]:", dmui, tracename, calci)
                printf("No more candidate for entering variable.\n:(CB*Bi*Nj)-Cj\n");T
            }
			if (phase == 1) {
				lpresult.rc = 3498;
				lpresult.rmsg = sprintf("%s[DMUi=%g][LOOP=%g]%s",
					"No Solution(No more candidate for entering variable):",
					dmui, calci, tracename)
			}
            break
        }

if (trace == "trace") {
    displayas("txt"); msg = "Select entering variable."
    printf("\n[DMUi=%g]%s-LOOP[%g]: %s(%g:%g)\n:(CB*Bi*Nj)-Cj\n",
        dmui, tracename, calci, msg, enteringVar, T[enteringVar]); T
}

        // ---------------------------------------------------------------------
        // Select leaving variable.
        // ---------------------------------------------------------------------
        leavingVar = 0
        Aj = Nj[,enteringVar]
        alpha = edittozerotol(Bi*Aj, tol1)
        if (existArtificial) {
            T = VARS[1,2::1+slacks]; tcols = cols(T)
            for (j=1; j<=tcols; j++) {
                if (T[1,j] < 0 && alpha[j,1] != 0) {
                    leavingVar = j; lvj = 1+leavingVar
                    break;
                }
            }
        }

		if (leavingVar == 0) {
			boundT = boundM[1,2::1+slacks]
			tcols = cols(boundT)
			lowerT = upperT = J(1,tcols,.)
			for(j=1; j<=tcols; j++) {
				lowerT[1,j] = boundT[1,j].lower
				upperT[1,j] = boundT[1,j].upper
			}

			// XB=Bi*b
			leavingCase = 0
			if (boundM[1,evj].val ==  boundM[1,evj].lower) {
				minVal = .;
				// 1. alpha's positive min value
				TH1 = boundM[1,evj].lower :+ ((rawXB :- lowerT')
											:/ (alpha :* (alpha :> 0)))
				// TH1 = edittozerotol(TH1, tol1)
				if (any(TH1 :< minVal)) {
					minindex(TH1, 1, mi, w)
					leavingVar = mi[1]; lvj = 1+leavingVar
					if (phase == 1 && w[1,2] >= 2 && VARS[1,lvj] > 0) {
						// if phase 1 and same min ratio test result,
						// artificial variable must leave first.
						tcols = w[1,2]
						for (j=2; j<=tcols; j++) {
							if (VARS[1,mi[j]+1] < 0) {
								leavingVar = mi[j]; lvj = 1+leavingVar
								break;
							}
						}
					}
					minVal = TH1[leavingVar,1];
					leavingCase = 1
				}
				// 2. alpha's negative min value
				TH2 = boundM[1,evj].lower :+ ((rawXB :- upperT')
											:/ (alpha :* (alpha :< 0)))
				// TH2 = edittozerotol(TH2, tol1)
				if (any(TH2 :< minVal)) {
					minindex(TH2, 1, mi, w)
					leavingVar = mi[1]; lvj = 1+leavingVar
					minVal = TH1[leavingVar,1]; leavingCase = 2
				}
				// 3. get the enteringVar's upper value
				if (boundM[1,evj].upper < minVal) {
					minVal = boundM[1,evj].upper; leavingCase = 3
				}

if (trace == "trace") {
    displayas("txt"); msg = "Select leaving variable.[MinVal]"
    printf("\n[DMUi=%g]%s-LOOP[%g]: %s(%g:%g)\n",
        dmui, tracename, calci, msg, leavingVar, minVal)
    display("XB | alpha:(XB=Bi*b, alpha=Bi*Aj):");rawXB,alpha
    display("[MinVal]enteringVar's upper | theta1 | theta2")
	printf("\n[boundM[1,%g].upper:%g][leavingCase:%g]\n",
			evj, boundM[1,evj].upper, leavingCase); TH1,TH2
}

				if (leavingCase == 1) {
					boundM[1,lvj].val = boundM[1,lvj].lower
				}
				else if (leavingCase == 2) {
					boundM[1,lvj].val = boundM[1,lvj].upper
				}
				else { // if (leavingCase == 3)
					boundM[1,evj].val = boundM[1,evj].upper; continue;
				}
			}
			else { // if (boundM[1,evj].val ==  boundM[1,evj].upper)
				maxVal = 0;
				// 1. alpha's positive min value
				TH1 = boundM[1,evj].upper :+ ((rawXB :- upperT')
											:/ (alpha :* (alpha :> 0)))
				// TH1 = edittozerotol(TH1, tol1)
				if (any(TH1 :> maxVal)) {
					maxindex(TH1, 1, mi, w)
					leavingVar = mi[1]; lvj = 1+leavingVar
					maxVal = TH1[leavingVar,1]; leavingCase = 1
				}
				// 2. alpha's negative min value
				TH2 = boundM[1,evj].upper :+ ((rawXB :- lowerT')
											:/ (alpha :* (alpha :< 0)))
				// TH2 = edittozerotol(TH2, tol1)
				if (any(TH2 :> maxVal)) {
					maxindex(TH2, 1, mi, w)
					leavingVar = mi[1]; lvj = 1+leavingVar
					maxVal = TH1[leavingVar,1]; leavingCase = 2
				}
				// 3. get the enteringVar's lower value
				if (boundM[1,evj].lower > maxVal) {
					maxVal = boundM[1,evj].lower; leavingCase = 3
				}

if (trace == "trace") {
    displayas("txt"); msg = "Select leaving variable.[MaxVal]"
    printf("\n[DMUi=%g]%s-LOOP[%g]: %s(%g:%g)\n",
        dmui, tracename, calci, msg, leavingVar, maxVal)
    display("XB | alpha:(XB=Bi*b, alpha=Bi*Aj):");XB,alpha
    display("[MaxVal]enteringVar's lower | theta1 | theta2")
	printf("\n[boundM[1,%g].lower:%g][leavingCase:%g]\n",
			evj, boundM[1,evj].lower, leavingCase); TH1,TH2
}

			    if (leavingCase == 1) {
					boundM[1,lvj].val = boundM[1,lvj].upper
				}
				else if (leavingCase == 2) {
					boundM[1,lvj].val = boundM[1,lvj].lower
				}
				else  { // if (leavingCase == 3)
					boundM[1,evj].val = boundM[1,evj].lower; continue;
				}
			}
        }

        // If no leaving variable exits
        if (leavingVar == 0) {
            if (trace == "trace")
                display("Break: No more candidate for leaving variable.")

			lpresult.rc = 3498;
			lpresult.rmsg = sprintf("%s[DMUi=%g][LOOP=%g]%s",
				"No Solution(No more candidate for leaving variable):",
				dmui, calci, tracename)

            break
        }
        // theta is not leaving at phaseII.
        if (phase == 2 && VARS[,lvj] == 1) {
            if (trace == "trace") display("Break: theta(¥è) is not leaving.")
            break
        }

        // ---------------------------------------------------------------------
        // reply calculatation result.
        // ---------------------------------------------------------------------
		LVi[leavingVar,1] = VARS[,evj]
		_swapcols(M, lvj, evj)
		_swapcols(boundM, lvj, evj)
		_swapcols(VARS, lvj, evj)

        // Clear artificial variable
		if (VARS[,evj] < 0) {
			T = J(1, cols(VARS), 1); T[1,evj] = 0

			VARS   = select(VARS, T)
			M      = select(M, T)
			boundM = select(boundM, T)
			mcols  = cols(M)

			// if exist artificial(phase II)
			if (existArtificial) {
				existArtificial = any(VARS[,2::1+slacks] :< 0)
			}
		}

if (trace == "trace") {
    printf("\n[DMUi=%g]%s-LOOP[%g]: updated tableau.[%g(%g) <--> %g(%g)]\n",
        dmui, tracename, calci, lvj, leavingVar, evj, enteringVar); M
    display("LVi: Entered variable's VARS index value."); LVi
    display("VARS: Variable's index."); VARS
}

    } //end of main for

    // return lpresult
	if (calci > maxiter) {
		lpresult.rc = 3498;
		lpresult.rmsg = sprintf("%s[DMUi=%g][LOOP=%g]%s",
			"No Solution(LOOP greater than maxiter):",
			dmui, calci, tracename)
	}
	if(lpresult.rc) display(lpresult.rmsg)
	lpresult.xVal = xVal
	lpresult.XB = XB
    return(lpresult)
}
/* A[.,lvj] <--> A[.,evj] */
function _swapcols( transmorphic matrix A,
					real scalar lvj,
					real scalar evj )
{
		transmorphic colvector  v

		v = A[., lvj]
		A[., lvj] = A[., evj]
		A[., evj] = v
}


function replacesubmat ( transmorphic matrix M,
                         real scalar row,
                         real scalar col,
                         transmorphic matrix T )
{
    M[|row,col\row + rows(T) - 1, col + cols(T) - 1|] = T
}

function _setup_dearslt_names(string scalar dearsltmat,
                              string scalar dmuinmat,
                              string scalar dmuoutmat )
{
    string matrix DMU_CS     // dmu in matrix column stripes
    string matrix DEARSLT_CS // dea result matrix column stripes
    string matrix DEARSLT_RS // dea result matrix row stripes
    real matrix M
    real scalar mcols, cnt

    M = st_matrix(dearsltmat)
    mcols = cols(M)

    DMU_CS = st_matrixcolstripe(dmuinmat)
    for (i = 1; i <= rows(DMU_CS); i++) {
        DMU_CS[i, 1] = "ref"
    }

    DEARSLT_CS = ("","rank"\"","theta")\DMU_CS\ // column join
        st_matrixrowstripe(dmuinmat)\st_matrixrowstripe(dmuoutmat)
    if (mcols - rows(DEARSLT_CS) > 0) {
        cnt = 0
        for (i = rows(DEARSLT_CS)+1 ; i <= mcols ; i++) {
            DEARSLT_CS = DEARSLT_CS \ ("slack", "slack_" + strofreal(++cnt))
        }
    }

    DEARSLT_RS = st_matrixcolstripe(dmuinmat)

    // name the row and column of dea result matrix
    st_matrixrowstripe(dearsltmat, DEARSLT_RS)
    st_matrixcolstripe(dearsltmat, DEARSLT_CS)
}

/**
 * deamat - dmucount x ( 1(theta) + dmu count + slcak(in, out) count)
 */
function _dmurank( string scalar deamat,
                   real scalar dmuincount,
                   real scalar dmuoutcount,
                   real scalar minrank,
                   real scalar tol )
{
    real matrix M
    real rowvector v, vv, retvec, slcaksum
    real scalar m, mm, row

    M = st_matrix(deamat)
    v = round(M[,1], tol)
    if (minrank) minindex(v, rows(v), i, w)
    else maxindex(v, rows(v), i, w)

    retvec = J(rows(v), 1, .)
    if (allof(w[,2], 1)) {
        retvec[i[1::rows(v)]] = (1::rows(v))
    }
    else {
        // rank correction for ties
        slcaksum = rowsum(M[|1,cols(M) - (dmuincount + dmuoutcount - 1)\.,.|])
        for (m = 1; m <= rows(w); m++) {
            if (w[m,2] >= 2) {
                vv = i[w[m,1]::(w[m,1] + w[m,2] - 1)]
                if (minrank) maxindex(slcaksum[vv], w[m,2], ii, ww)
                else minindex(slcaksum[vv], w[m,2], ii, ww)
                for (mm = 1; mm <= rows(ww); mm++) {
                    for (row = ww[mm,1]; row < ww[mm,1] + ww[mm,2]; row++) {
                        retvec[vv[ii[row]]] = w[m,1] + ww[mm,1] - 1
                    }
                }
            }
            else {
                retvec[i[w[m,1]]] = w[m,1] // row = w[m,1]
            }
        }
    }
    st_matrix("r(rank)", retvec)
}

function maxvecindex( string scalar vecname )
{
    real matrix A

    A = st_matrix(vecname)
    maxindex(A, 1, i, w)

    st_numscalar("r(maxval)", A[i[1]])
    st_numscalar("r(maxindex)", i[1])
    st_matrix("r(maxindexes)", i)
}

function minvecindex( string scalar vecname )
{
    real matrix A

    A = st_matrix(vecname)
    if (sum(A :< .) > 0) {
        minindex(A, 1, i, w)
        st_numscalar("r(minval)", A[i[1]])
        st_numscalar("r(minindex)", i[1])
        st_matrix("r(minindexes)", i)
    }
    // if overall missing value.
    else {
        st_numscalar("r(minval)", .)
        st_numscalar("r(minindex)", 0)
        st_matrix("r(minindexes)", 0)
    }
}

function _roundmat( string scalar matname, real scalar tol )
{
    real matrix A
    A = round(st_matrix(matname), tol)
    st_matrix(matname, A)
}

function _uniqrowmat( string scalar matname, string scalar varname )
{
    st_matrix(matname, sort(uniqrows(st_data(., varname)), 1))
}

function _file_exists( string scalar fn )
{
    st_numscalar("r(fileexists)", fileexists(fn))
}
end

// End of the MATA Definition Area ---------------------------------------------