*! geneigen  1.0.2    28Dec2002  C F Baum
*  makes use of J.Weesie's qsortidx
* 1.0.1 2C04 Correct line 218
* 1.0.2 2C28 Correct first loop elmhes
*
program define geneigen, rclass
version 7.0
capt local n = rowsof(`1')
if _rc > 0 {
	di in r "usage: geneigen matrixname"
	error 198
	}
local nc = colsof(`1')
if `n' == 0 | `n' ~= `nc' {
	di in r "error: geneigen must be called with a square matrix"
	error 198
	}
tempname mat1 mat2 mat3 wr wi modulus res prod
mat `mat1' = `1'
_balanc `mat1'
mat `mat2'= r(_balanc)
_elmhes `mat2'
mat `mat3' = r(_elmhes)
_hqr `mat3'
mat `wr' = r(lambdar)
mat `wi' = r(lambdai)
mat `modulus' = J(`n',1,0)
forv i=1/`n' {
	mat `modulus'[`i',1] = sqrt(`wr'[`i',1]^2+`wi'[`i',1]^2)
	local m = `modulus'[`i',1]
	local mod "`mod' `m'"
	}
qsortidx `mod'
local ord = "`s(order)'"
matrix `res' = J(`n',4,0)
scalar `prod' = 1
forv i=1/`n' {
	local elt: word `i' of `ord'
	mat `res'[`i',1]=`modulus'[`elt',1]
	mat `res'[`i',2]=`wr'[`elt',1]
	mat `res'[`i',3]=`wi'[`elt',1]
	mat `res'[`i',4]=`res'[`i',1]
	if `i'>1 { mat `res'[`i',4]=`res'[`i',4]+`res'[`i'-1,4] }
	scalar `prod' = `prod' *  `modulus'[`elt',1]
	}
mat coln `res' = Modulus Real Imag Cum
return scalar n = `n'
return scalar det = `prod'
return matrix eigen `res'
end

* NumRecip balanc
program define _balanc,rclass
version 7.0
tempname a r c g f s
mat `a' = `1'
local n = rowsof(`a')
scalar radix = 2
scalar sqrdx = radix*radix
local last = 0
while `last' == 0 {
	local last = 1
	forv i = 1/`n' {
		scalar `r' = 0
		scalar `c' = 0
		forv j = 1/`n' {
			if `j' ~= `i' {
				scalar `c' = `c' + abs(`a'[`j',`i'])
				scalar `r' = `r' + abs(`a'[`i',`j'])
				}
			}
		if `c' ~= 0 & `r' ~= 0 {
			scalar `g' = `r' / radix
			scalar `f' = 1
			scalar `s' = `c' + `r'
			while `c' < `g' {
				scalar `f' = `f' * radix
				scalar `c' = `c' * sqrdx
				}
			scalar `g' = `r' * radix
			while `c' > `g' {
				scalar `f' = `f' / radix
				scalar `c' = `c' / sqrdx
			}
		if (`c'+`r')/`f' < 0.95 * `s' {
			local last = 0
			scalar `g' = 1.0/`f'
			forv j=1/`n' {
				mat `a'[`i',`j'] = `a'[`i',`j'] * `g'
				mat `a'[`j',`i'] = `a'[`j',`i'] * `f'
				}
			}
		}
	}
}
return matrix _balanc  `a'
end

* NumRecip elmhes
program define _elmhes,rclass
version 7.0
tempname a x y t
mat `a' = `1'
local n = rowsof(`a')
local en1 = `n'-1
forv m = 2/`en1' {
	scalar `x' = 0.0
	local i = `m'
	forv j = `m'/`n' {
		if abs(`a'[`j',`m'-1]) > abs(`x') {
			scalar `x' = `a'[`j',`m'-1]
			local i = `j'
			}
		}
	if `i' ~= `m' {
		local mm1 = `m'-1
		forv j = `mm1'/`n' {
			scalar `t' = `a'[`i',`j']
			mat `a'[`i',`j'] = `a'[`m',`j']
			mat `a'[`m',`j'] = `t'
			}
		forv j = 1/`n' {
			scalar `t' = `a'[`j',`i']
			mat `a'[`j',`i'] = `a'[`j',`m']
			mat `a'[`j',`m'] = `t'
			}
		}
	if `x' ~= 0 {
		local m1 = `m'+1
		forv i = `m1'/`n' {
			scalar `y' = `a'[`i',`m'-1]
			if `y' ~= 0.0 {
				scalar `y' = `y' / `x'
				mat `a'[`i',`m'-1] = `y'
				forv j = `m'/`n' {
					mat `a'[`i',`j'] = `a'[`i',`j'] - `y'*`a'[`m',`j']
					}
				forv j = 1/`n' {
					mat `a'[`j',`m'] = `a'[`j',`m'] + `y'*`a'[`j',`i']
					}
				}
			}
		}
	}
return matrix _elmhes `a'
end

* NumRecip hqr
program define _hqr,rclass
version 7.0
tempname a wr wi
mat `a' = `1'
local n = rowsof(`a')
mat `wr' = J(`n',1,0)
mat `wi' = J(`n',1,0)
tempname z y w x v u t s r q p anorm
scalar `anorm' = abs(`a'[1,1])
forv i = 2/`n' {
	local im1 = `i'-1
	forv j = `im1'/`n' {
		scalar `anorm' = `anorm' + abs(`a'[1,`j'])
		}
	}
local nn = `n'
scalar `t' = 0.0
while `nn' >= 1 {
	local its = 0
	local ell = 1
*	di in r "** ell, nn `ell' `nn'"
	while `ell' <= `nn'  {
*	while `ell' < `nn' - 1 {
		forv l = `nn'(-1)2 {
			scalar `s' = abs(`a'[`l'-1,`l'-1]) + abs(`a'[`l',`l'])
			if `s' == 0.0 { scalar `s' = `anorm' }
*		di in r "s " `s' " `l' " `a'[`l',`l'-1]
			if abs(`a'[`l',`l'-1]) + `s' == `s' {
					local ell = `l'
					continue, break
					}
			}
		scalar `x' = `a'[`nn',`nn']
*	di in r "ell, nn, x `ell' `nn' " `x'
		if `ell' == `nn' {
			mat `wr'[`nn',1] = `x' + `t'
			mat `wi'[`nn',1] = 0.0
			local nn = `nn' - 1
*	di in r "ell, nn now `ell' `nn'"
			}
		else {
			scalar `y' = `a'[`nn'-1,`nn'-1]
			scalar `w' = `a'[`nn',`nn'-1] * `a'[`nn'-1,`nn']
			if `ell' == `nn'-1 {
				scalar `p' = 0.5 * (`y'-`x')
				scalar `q' = `p'*`p'+`w'
				scalar `z' = sqrt(abs(`q'))
				scalar `x' = `x' + `t'
				if `q' >= 0.0 {
					scalar `z' = `p' + `z'*sign(`p')
					mat `wr'[`nn'-1,1] = `x' + `z'
					mat `wr'[`nn',1] = `x' + `z'
					if `z' ~= 0 { mat `wr'[`nn',1] = `x' - `w' / `z' }
					mat `wi'[`nn'-1,1] = 0.0
					mat `wi'[`nn',1] = 0.0
					}
				else {
					mat `wr'[`nn'-1,1] = `x' + `p'
					mat `wr'[`nn',1] = `x' + `p'
					mat `wi'[`nn'-1,1] = -`z'
					mat `wi'[`nn',1] = `z'
					}
				local nn = `nn' - 2
			}
			else {
				if `its' == 30 {
					di in r "Too many iterations"
					error 198
					}
				if `its' == 10 | `its' == 20 {
					scalar `t' = `t' + `x'
* corr: forv 1/nn
					forv i = 1/`nn' { mat `a'[`i',`i'] = `a'[`i',`i'] - `x' }
					scalar `s' = abs(`a'[`nn',`nn'-1]) + abs(`a'[`nn'-1,`nn'-2])
					scalar `y' = 0.75 * `s'
					scalar `x' = 0.75 * `s'
					scalar `w' = -0.4375 * `s' * `s'
					}
				local its = `its' + 1
*				di in r "its is `its' nn is `nn'"
				local nn2 = `nn' - 2
*		di in r "`nn2'|`ell'"
				forv m = `nn2'(-1)`ell' {
					local emm = `m'
*		di in r "`m'|`emm'"
					scalar `z' = `a'[`m',`m']
					scalar `r' = `x' - `z'
					scalar `s' = `y' - `z'
					scalar `p' = (`r'*`s'-`w')/`a'[`m'+1,`m']+`a'[`m',`m'+1]
					scalar `q' = `a'[`m'+1,`m'+1]-`z'-`r'-`s'
					scalar `r' = `a'[`m'+2,`m'+1]
					scalar `s' = abs(`p')+abs(`q')+abs(`r')
					scalar `p' = `p' / `s'
					scalar `q' = `q' / `s'
					scalar `r' = `r' / `s'
					if `m' == `ell' { continue, break }
					scalar `u' = abs(`a'[`m',`m'-1])*(abs(`q')+abs(`r'))
					scalar `v' = abs(`p')*(abs(`a'[`m'-1,`m'-1])+abs(`z')+abs(`a'[`m'+1,`m'+1]))
					if `u' + `v' == `v' { continue, break }
					}
*		di in r "`m'|`emm'"
				local m2 = `emm' + 2
				forv i = `m2'/ `nn' {
					mat `a'[`i',`i'-2] = 0.0
					if `i' ~= `m2' { mat `a'[`i',`i'-3] = 0.0 }
					}
				local nn1 = `nn' - 1
*	di in r "`emm' `nn1' "
				forv k = `emm'/`nn1' {
					if `k' ~= `emm' {
						scalar `p' = `a'[`k',`k'-1]
						scalar `q' = `a'[`k'+1,`k'-1]
						scalar `r' = 0.0
						if `k' ~= `nn1' { scalar `r' = `a'[`k'+2,`k'-1] }
						scalar `x' = abs(`p') + abs(`q') + abs(`r')
						if `x' ~= 0.0 {
							scalar `p' = `p' / `x'
							scalar `q' = `q' / `x'
							scalar `r' = `r' / `x'
							}
						}
					scalar `s' = sqrt(`p'*`p'+`q'*`q'+`r'*`r')*sign(`p')
					if `s' ~= 0.0 {
						if `k' ~= `emm' {
							mat `a'[`k',`k'-1] = -`s'*`x'
							}
						else {
							if `ell' ~= `emm' { mat `a'[`k',`k'-1] = -`a'[`k',`k'-1] }
							}
						scalar `p' = `p' + `s'
						scalar `x' = `p' / `s'
						scalar `y' = `q' / `s'
						scalar `z' = `r' / `s'
						scalar `q' = `q' / `p'
						scalar `r' = `r' / `p'
						forv j = `k'/`nn' {
							scalar `p' = `a'[`k',`j'] + `q' * `a'[`k'+1,`j']
							if `k' ~= `nn'-1 {
								scalar `p' = `p' + `r' * `a'[`k'+2,`j']
								mat `a'[`k'+2,`j'] = `a'[`k'+2,`j'] - `p'*`z'
								}
							mat `a'[`k'+1,`j'] = `a'[`k'+1,`j'] - `p'*`y'
							mat `a'[`k',`j'] = `a'[`k',`j'] - `p'*`x'
							}
						local mmin = min(`nn',`k'+3)
						forv i = `ell'/`mmin' {
							scalar `p' = `x'*`a'[`i',`k']+`y'*`a'[`i',`k'+1]
							if `k' ~= `nn'-1 {
								scalar `p' = `p' + `z'*`a'[`i',`k'+2]
								mat `a'[`i',`k'+2] = `a'[`i',`k'+2] - `p'*`r'
								}
							mat `a'[`i',`k'+1] = `a'[`i',`k'+1] - `p'*`q'
							mat `a'[`i',`k'] = `a'[`i',`k'] - `p'
							}  /*  18 */
					}	/* end if s ~= 0 */
				}  	/*  19 */
				continue  /* equivalent to goto 2 */
			}  		/* light blue */
		}  		/* dark blue */
	}  		/* end of while ell < nn-1  */
}  		/* end of while nn >= 1  */
* mat list `wr'
* mat list `wi'
return matrix lambdar `wr'
return matrix lambdai `wi'
end