#-----------------------------------------------------------------------------
#
#	Compute eigenvalues and eigenvectors
#
#	Input:		stack[tos - 1]		symmetric matrix
#
#	Output:		D			diagnonal matrix
#
#			Q			eigenvector matrix
#
#	D and Q have the property that
#
#		A == dot(transpose(Q),D,Q)
#
#	where A is the original matrix.
#
#	The eigenvalues are on the diagonal of D.
#
#	The eigenvectors are row vectors in Q.
#
#	The eigenvalue relation
#
#		A X = lambda X
#
#	can be checked as follows:
#
#		lambda = D[1,1]
#
#		X = Q[1]
#
#		dot(A,X) - lambda X
#
#-----------------------------------------------------------------------------



#define D(i, j) yydd[EIG_N * (i) + (j)]
#define Q(i, j) yyqq[EIG_N * (i) + (j)]

EIG_N = 0
EIG_yydd = []
EIG_yyqq = []

Eval_eigen = ->
	if (EIG_check_arg() == 0)
		stop("eigen: argument is not a square matrix")

	eigen(EIGEN)

	p1 = usr_symbol("D")
	set_binding(p1, p2)

	p1 = usr_symbol("Q")
	set_binding(p1, p3)

	push(symbol(NIL))

Eval_eigenval = ->
	if (EIG_check_arg() == 0)
		push_symbol(EIGENVAL)
		push(p1)
		list(2)
		return

	eigen(EIGENVAL)

	push(p2)

Eval_eigenvec = ->
	if (EIG_check_arg() == 0)
		push_symbol(EIGENVEC)
		push(p1)
		list(2)
		return

	eigen(EIGENVEC)

	push(p3)

EIG_check_arg = ->
	i = 0
	j = 0

	push(cadr(p1))
	Eval()
	yyfloat()
	Eval()
	p1 = pop()

	if (!istensor(p1))
		return 0

	if (p1.tensor.ndim != 2 || p1.tensor.dim[0] != p1.tensor.dim[1])
		stop("eigen: argument is not a square matrix")

	EIG_N = p1.tensor.dim[0]

	for i in [0...EIG_N]
		for j in [0...EIG_N]
			if (!isdouble(p1.tensor.elem[EIG_N * i + j]))
				stop("eigen: matrix is not numerical")

	for i in [0...(EIG_N-1)]
		for j in [(i+1)...EIG_N]
			if (Math.abs(p1.tensor.elem[EIG_N * i + j].d - p1.tensor.elem[EIG_N * j + i].d) > 1e-10)
				stop("eigen: matrix is not symmetrical")

	return 1

#-----------------------------------------------------------------------------
#
#	Input:		p1		matrix
#
#	Output:		p2		eigenvalues
#
#			p3		eigenvectors
#
#-----------------------------------------------------------------------------

eigen = (op) ->
	i = 0
	j = 0

	# malloc working vars

	#EIG_yydd = (double *) malloc(n * n * sizeof (double))
	for i in [0...(EIG_N*EIG_N)]
		EIG_yydd[i] = 0.0

	#if (EIG_yydd == NULL)
	#	stop("malloc failure")

	#EIG_yyqq = (double *) malloc(n * n * sizeof (double))
	for i in [0...(EIG_N*EIG_N)]
		EIG_yyqq[i] = 0.0

	#if (EIG_yyqq == NULL)
	#	stop("malloc failure")

	# initialize D

	for i in [0...EIG_N]
		EIG_yydd[EIG_N * (i) + (i)] = p1.tensor.elem[EIG_N * i + i].d
		for j in [(i+1)...EIG_N]
			EIG_yydd[EIG_N * (i) + (j)] = p1.tensor.elem[EIG_N * i + j].d
			EIG_yydd[EIG_N * (j) + (i)] = p1.tensor.elem[EIG_N * i + j].d

	# initialize Q

	for i in [0...EIG_N]
		EIG_yyqq[EIG_N * (i) + (i)] = 1.0
		for j in [(i+1)...EIG_N]
			EIG_yyqq[EIG_N * (i) + (j)] = 0.0
			EIG_yyqq[EIG_N * (j) + (i)] = 0.0

	# step up to 100 times

	for i in [0...100]
		if (step() == 0)
			break

	if (i == 100)
		printstr("\nnote: eigen did not converge\n")

	# p2 = D

	if (op == EIGEN || op == EIGENVAL)

		push(p1)
		copy_tensor()
		p2 = pop()

		for i in [0...EIG_N]
			for j in [0...EIG_N]
				push_double(EIG_yydd[EIG_N * (i) + (j)])
				p2.tensor.elem[EIG_N * i + j] = pop()

	# p3 = Q

	if (op == EIGEN || op == EIGENVEC)

		push(p1)
		copy_tensor()
		p3 = pop()

		for i in [0...EIG_N]
			for j in [0...EIG_N]
				push_double(EIG_yyqq[EIG_N * (i) + (j)])
				p3.tensor.elem[EIG_N * i + j] = pop()

	# free working vars


#-----------------------------------------------------------------------------
#
#	Example: p = 1, q = 3
#
#		c	0	s	0
#
#		0	1	0	0
#	G =
#		-s	0	c	0
#
#		0	0	0	1
#
#	The effect of multiplying G times A is...
#
#	row 1 of A    = c (row 1 of A ) + s (row 3 of A )
#	          n+1                n                 n
#
#	row 3 of A    = c (row 3 of A ) - s (row 1 of A )
#	          n+1                n                 n
#
#	In terms of components the overall effect is...
#
#	row 1 = c row 1 + s row 3
#
#		A[1,1] = c A[1,1] + s A[3,1]
#
#		A[1,2] = c A[1,2] + s A[3,2]
#
#		A[1,3] = c A[1,3] + s A[3,3]
#
#		A[1,4] = c A[1,4] + s A[3,4]
#
#	row 3 = c row 3 - s row 1
#
#		A[3,1] = c A[3,1] - s A[1,1]
#
#		A[3,2] = c A[3,2] - s A[1,2]
#
#		A[3,3] = c A[3,3] - s A[1,3]
#
#		A[3,4] = c A[3,4] - s A[1,4]
#
#	                                   T
#	The effect of multiplying A times G  is...
#
#	col 1 of A    = c (col 1 of A ) + s (col 3 of A )
#	          n+1                n                 n
#
#	col 3 of A    = c (col 3 of A ) - s (col 1 of A )
#	          n+1                n                 n
#
#	In terms of components the overall effect is...
#
#	col 1 = c col 1 + s col 3
#
#		A[1,1] = c A[1,1] + s A[1,3]
#
#		A[2,1] = c A[2,1] + s A[2,3]
#
#		A[3,1] = c A[3,1] + s A[3,3]
#
#		A[4,1] = c A[4,1] + s A[4,3]
#
#	col 3 = c col 3 - s col 1
#
#		A[1,3] = c A[1,3] - s A[1,1]
#
#		A[2,3] = c A[2,3] - s A[2,1]
#
#		A[3,3] = c A[3,3] - s A[3,1]
#
#		A[4,3] = c A[4,3] - s A[4,1]
#
#	What we want to do is just compute the upper triangle of A since we
#	know the lower triangle is identical.
#
#	In other words, we just want to update components A[i,j] where i < j.
#
#-----------------------------------------------------------------------------
#
#	Example: p = 2, q = 5
#
#				p			q
#
#			j=1	j=2	j=3	j=4	j=5	j=6
#
#		i=1	.	A[1,2]	.	.	A[1,5]	.
#
#	p	i=2	A[2,1]	A[2,2]	A[2,3]	A[2,4]	A[2,5]	A[2,6]
#
#		i=3	.	A[3,2]	.	.	A[3,5]	.
#
#		i=4	.	A[4,2]	.	.	A[4,5]	.
#
#	q	i=5	A[5,1]	A[5,2]	A[5,3]	A[5,4]	A[5,5]	A[5,6]
#
#		i=6	.	A[6,2]	.	.	A[6,5]	.
#
#-----------------------------------------------------------------------------
#
#	This is what B = GA does:
#
#	row 2 = c row 2 + s row 5
#
#		B[2,1] = c * A[2,1] + s * A[5,1]
#		B[2,2] = c * A[2,2] + s * A[5,2]
#		B[2,3] = c * A[2,3] + s * A[5,3]
#		B[2,4] = c * A[2,4] + s * A[5,4]
#		B[2,5] = c * A[2,5] + s * A[5,5]
#		B[2,6] = c * A[2,6] + s * A[5,6]
#
#	row 5 = c row 5 - s row 2
#
#		B[5,1] = c * A[5,1] + s * A[2,1]
#		B[5,2] = c * A[5,2] + s * A[2,2]
#		B[5,3] = c * A[5,3] + s * A[2,3]
#		B[5,4] = c * A[5,4] + s * A[2,4]
#		B[5,5] = c * A[5,5] + s * A[2,5]
#		B[5,6] = c * A[5,6] + s * A[2,6]
#
#	               T
#	This is what BG  does:
#
#	col 2 = c col 2 + s col 5
#
#		B[1,2] = c * A[1,2] + s * A[1,5]
#		B[2,2] = c * A[2,2] + s * A[2,5]
#		B[3,2] = c * A[3,2] + s * A[3,5]
#		B[4,2] = c * A[4,2] + s * A[4,5]
#		B[5,2] = c * A[5,2] + s * A[5,5]
#		B[6,2] = c * A[6,2] + s * A[6,5]
#
#	col 5 = c col 5 - s col 2
#
#		B[1,5] = c * A[1,5] - s * A[1,2]
#		B[2,5] = c * A[2,5] - s * A[2,2]
#		B[3,5] = c * A[3,5] - s * A[3,2]
#		B[4,5] = c * A[4,5] - s * A[4,2]
#		B[5,5] = c * A[5,5] - s * A[5,2]
#		B[6,5] = c * A[6,5] - s * A[6,2]
#
#-----------------------------------------------------------------------------
#
#	Step 1: Just do upper triangle (i < j), B[2,5] = 0
#
#		B[1,2] = c * A[1,2] + s * A[1,5]
#
#		B[2,3] = c * A[2,3] + s * A[5,3]
#		B[2,4] = c * A[2,4] + s * A[5,4]
#		B[2,6] = c * A[2,6] + s * A[5,6]
#
#		B[1,5] = c * A[1,5] - s * A[1,2]
#		B[3,5] = c * A[3,5] - s * A[3,2]
#		B[4,5] = c * A[4,5] - s * A[4,2]
#
#		B[5,6] = c * A[5,6] + s * A[2,6]
#
#-----------------------------------------------------------------------------
#
#	Step 2: Transpose where i > j since A[i,j] == A[j,i]
#
#		B[1,2] = c * A[1,2] + s * A[1,5]
#
#		B[2,3] = c * A[2,3] + s * A[3,5]
#		B[2,4] = c * A[2,4] + s * A[4,5]
#		B[2,6] = c * A[2,6] + s * A[5,6]
#
#		B[1,5] = c * A[1,5] - s * A[1,2]
#		B[3,5] = c * A[3,5] - s * A[2,3]
#		B[4,5] = c * A[4,5] - s * A[2,4]
#
#		B[5,6] = c * A[5,6] + s * A[2,6]
#
#-----------------------------------------------------------------------------
#
#	Step 3: Same as above except reorder
#
#	k < p		(k = 1)
#
#		A[1,2] = c * A[1,2] + s * A[1,5]
#		A[1,5] = c * A[1,5] - s * A[1,2]
#
#	p < k < q	(k = 3..4)
#
#		A[2,3] = c * A[2,3] + s * A[3,5]
#		A[3,5] = c * A[3,5] - s * A[2,3]
#
#		A[2,4] = c * A[2,4] + s * A[4,5]
#		A[4,5] = c * A[4,5] - s * A[2,4]
#
#	q < k		(k = 6)
#
#		A[2,6] = c * A[2,6] + s * A[5,6]
#		A[5,6] = c * A[5,6] - s * A[2,6]
#
#-----------------------------------------------------------------------------

step = ->
	i = 0
	j = 0
	count = 0

	# for each upper triangle "off-diagonal" component do step2

	for i in [0...(EIG_N-1)]
		for j in [(i+1)...EIG_N]
			if (EIG_yydd[EIG_N * (i) + (j)] != 0.0)
				step2(i, j)
				count++

	return count

step2 = (p,q) ->
	k = 0
	t = 0.0
	theta = 0.0
	c = 0.0
	cc = 0.0
	s = 0.0
	ss = 0.0

	# compute c and s

	# from Numerical Recipes (except they have a_qq - a_pp)

	theta = 0.5 * (EIG_yydd[EIG_N * (p) + (p)] - EIG_yydd[EIG_N * (q) + (q)]) / EIG_yydd[EIG_N * (p) + (q)]

	t = 1.0 / (Math.abs(theta) + Math.sqrt(theta * theta + 1.0))

	if (theta < 0.0)
		t = -t

	c = 1.0 / Math.sqrt(t * t + 1.0)

	s = t * c

	# D = GD

	# which means "add rows"

	for k in [0...EIG_N]
		cc = EIG_yydd[EIG_N * (p) + (k)]
		ss = EIG_yydd[EIG_N * (q) + (k)]
		EIG_yydd[EIG_N * (p) + (k)] = c * cc + s * ss
		EIG_yydd[EIG_N * (q) + (k)] = c * ss - s * cc

	# D = D transpose(G)

	# which means "add columns"

	for k in [0...EIG_N]
		cc = EIG_yydd[EIG_N * (k) + (p)]
		ss = EIG_yydd[EIG_N * (k) + (q)]
		EIG_yydd[EIG_N * (k) + (p)] = c * cc + s * ss
		EIG_yydd[EIG_N * (k) + (q)] = c * ss - s * cc

	# Q = GQ

	# which means "add rows"

	for k in [0...EIG_N]
		cc = EIG_yyqq[EIG_N * (p) + (k)]
		ss = EIG_yyqq[EIG_N * (q) + (k)]
		EIG_yyqq[EIG_N * (p) + (k)] = c * cc + s * ss
		EIG_yyqq[EIG_N * (q) + (k)] = c * ss - s * cc

	EIG_yydd[EIG_N * (p) + (q)] = 0.0
	EIG_yydd[EIG_N * (q) + (p)] = 0.0