TQLI.hp42s

LBL "TQLI"
@ Calculates eigenvalues and eigenvectors of a symmetric tridiagonal matrix.
@ Optionally (FLAG 02) results are sorted in ascending eigenvalue order.
@ FLAGS 03 and 04 serve for compatibility with the LATTICE program.
@
@ flags:
@  FLAG 01: omit calculation of eigenvectors if SET
@  FLAG 02: sort eigenvalues in ascending order and eigenvectors accordingly
@  FLAG 03: expect last element of VE to be zero if SET (instead of first one)
@  FLAG 04: create an identity MZ matrix if SET
@
@ FLAGS 03 and 04 will be cleared by TQLI.
@
@ input:
@  VD: 1xN vector of diagonal elements
@  VE: 1xN vector of off-diagonal elements (see also FLAG 03 for zero padding)
@  MZ: NxN matrix as prepared by tred2 or unity matrix 
@      (only required input if FLAG 04 not set)
@ output:
@  VD: 1xN vector of eigenvalues
@  MZ: NxN matrix of eigenvectors (if FLAG 01 is not set)
@
@ internally used:
@ REG 01: dimension of matrix (n)
@ REG 02: iterator of outermost loop (l)
@ REG 03: counter for number of iterations (iter)
@ REG 04: inner loop iterator (m)
@ REG 05: inner loop iterator (i)
@ REG 06: inner loop iterator (k)
@ REG 07: calculation helper (s)
@ REG 08: calculation helper (c)
@ REG 09: calculation helper (r)
@ REG 10: calculation helper (p)
@ REG 11: calculation helper (g)
@ REG 12: calculation helper (f)
@ REG 13: calculation helper (b)
RCL "VD"
DIM?
STO 01		@ n = dim(VD)
FC?C 03		@ FLAG 03 not set?
	XEQ 30		@ shift zero element of VE
FS?C 04		@ FLAG 04 set?
	XEQ 32		@ initialize MZ matrix
RCL 01          @ n
1
RCL 01          @ n
STOIJ
0
STOEL           @ e[n] = 0
RCL 01		@ n
1ᴇ3
÷
1
+
STO 02		@ l=1; l<=n; l++
LBL 11
	0
	STO 03		@ iter = 0
	LBL 12		@ do {
		RCL 01		@ n
		1
		-		@ n-1
		1ᴇ3
		÷
		RCL 02		@ l iterator
		IP		@ l
		+
		STO 04		@ m=l; m<=n-1; m++
		LBL 13
		INDEX "VD"
		1
		RCL 04		@ m iterator
		STOIJ
		RCLEL		@ d[m]
		ABS
		J+
		RCLEL		@ d[m+1]
		ABS
		+		@ fabs(d[m]) + fabs(d[m+1])
		1ᴇ-10		@ EPS = 1E-10
		INDEX "VE"
		1
		RCL 04		@ m iterator
		STOIJ
		R↓
		R↓
		RCLEL		@ e[m]
		ABS
		X<=Y?		@ if (fabs(e[m]) <= EPS*dd)
			GTO 14		@ break
		ISG 04		@ while (m <= n-1)
			GTO 13
		LBL 14
			RCL 04		@ m iterator
			IP		@ m
			RCL 02		@ l iterator
			IP		@ l
			X=Y?		@ if (m == l)
				GTO 21		@ exit do ... while loop
			1
			STO+ 03		@ iter++
			RCL 03		@ iter
			30
			X=Y?		@ if (iter == 30)
				GTO 91		@ nerror
			INDEX "VD"
			1
			1
			RCL+ 02	        @ += l
			STOIJ
			RCLEL		@ d[l+1]
			J-
			RCLEL		@ d[l]
			-
			2
			÷		@ (d[l+1]-d[l]) / 2.0
			RCLIJ
			INDEX "VE"
			STOIJ
			R↓
			R↓
			RCLEL		@ e[l]
			÷
			STO 11		@ g = (d[l+1]-d[l]) / (2.0*e[l])
			X^2
			1
			+
			SQRT		@ sqrt(g^2 + 1)
			STO 09		@ r = pythag(g,1.0);
			INDEX "VD"
			1
			RCL 04		@ m
			STOIJ
			RCLEL		@ d[m]
			STO 11		@ g = d[m]
			1
			RCL 02		@ l
			STOIJ
			RCLEL		@ d[l]
			STO- 11         @ g = d[m]-d[l]
			RCLIJ
			INDEX "VE"
			STOIJ
			RCLEL		@ e[l]
			RCL 09		@ r
			ABS		@ abs(r)
			RCL 11		@ g
			SIGN		@ sign(g)
			×		@ abs(r) * sign(g)
			RCL 11		@ g
			+		@ g+SIGN(r,g)
			÷		@ e[l] / (g+SIGN(r,g))
			STO+ 11		@ g = d[m]-d[l] + e[l] / (g+SIGN(r,g))
			1
			STO 07		@ s = 1.0
			STO 08		@ c = 1.0
			0
			STO 10		@ p = 0.0
			RCL 02		@ l iterator
			IP		@ l
			1
			-		@ l-1
			1ᴇ3
			÷
			-1
			RCL+ 04		@ m iterator -1
			IP		@ m-1
			+
			STO 05		@ i=m-1; i>l-1; i--
			LBL 15
				INDEX "VE"
				1
				RCL 05		@ i
				STOIJ
				RCLEL		@ e[i]
				RCL× 07		@ *= s
				STO 12		@ f = s*e[i]
				RCLEL		@ e[i]
				RCL× 08		@ *= c
				STO 13		@ b = c*e[i]
				RCL 12		@ f
				X^2
				RCL 11		@ g
				X^2
				+
				SQRT		@ sqrt(f^2 + g^2)
				STO 09		@ r = pythag(f,g)
				J+
				STOEL		@ e[i+1] = r
				X≠0?		@ r != 0
					GTO 16
				RCLIJ
				INDEX "VD"
				STOIJ
				RCLEL		@ d[i+1]
				RCL- 10		@ -= p
				STOEL		@ d[i+1] -= p
				INDEX "VE"
				1
				RCL 04		@ m
				STOIJ
				0
				STOEL		@ e[m] = 0.0
				GTO 19		@ break
				LBL 16
				RCL 12		@ f
				RCL÷ 09		@ /= r
				STO 07		@ s = f/r
				RCL 11		@ g
				RCL÷ 09		@ /= r
				STO 08		@ c = g/r
				INDEX "VD"
				1
				RCL 05		@ i
				1
				+
				STOIJ
				RCLEL		@ d[i+1]
				RCL- 10		@ -= p
				STO 11		@ g = d[i+1]-p
				J-
				RCLEL		@ d[i]
				RCL- 11		@ -= g
				RCL× 07		@ (d[i]-g)*s
				2
				RCL× 08		@ *= c
				RCL× 13		@ 2.0*c*b
				+
				STO 09		@ r = (d[i]-g)*s + 2.0*c*b
				RCL× 07		@ *= s
				STO 10		@ p = s*r
				RCL+ 11		@ += g
				J+
				STOEL		@ d[i+1] = g+p
				RCL 08		@ c
				RCL× 09		@ *= r
				RCL- 13		@ -= b
				STO 11		@ g = c*r-b
				FS? 01		@ skip eigenvector calculation?
					GTO 18
				RCL 01		@ n
				1ᴇ3
				÷
				1
				+
				STO 06		@ k=1; k<=n; k++
				INDEX "MZ"
				LBL 17
					RCL 06		@ k
					1
					RCL+ 05		@ i+1
					STOIJ
					RCLEL		@ z[k][i+1]
					STO 12		@ f = z[k][i+1]
					J-
					RCLEL		@ z[k][i]
					RCL× 07		@ *= s
					RCL 08		@ c
					RCL× 12		@ *= f
					+		@ s*z[k][i] + c*f
					J+
					STOEL		@ z[k][i+1] = s*z[k][i] + c*f
					J-
					RCLEL		@ z[k][i]
					RCL× 08		@ *= c
					RCL 07		@ s
					RCL× 12		@ *= f
					-
					STOEL		@ z[k][i] = c*z[k][i] - s*f
					ISG 06		@ for (k=1; k<=n; k++)
						GTO 17
				LBL 18
				DSE 05		@ for (i=m-1; i>=l; i--)
					GTO 15
			LBL 19
			RCL 09		@ r
			X≠0?		@ if (r != 0)
				GTO 20
			RCL 02		@ l iterator
			IP		@ l
			RCL 05		@ i iterator
			IP		@ i
			X≥Y?		@ if (i >= l)
				GTO 14		@ continue
			LBL 20
			INDEX "VD"
			1
			RCL 02		@ l
			STOIJ
			RCLEL		@ d[l]
			RCL- 10		@ -= p
			STOEL		@ d[l] -= p
			RCLIJ
			INDEX "VE"
			STOIJ
			RCL 11		@ g
			STOEL		@ e[l] = g
			1
			RCL 04		@ m
			STOIJ		@ e[m]
			0
			STOEL		@ e[m] = 0.0
			GTO 12		@ next do ... while loop iteration
	LBL 21
	ISG 02		@ for (l=1; l<=n; l++)
		GTO 11
	FS? 02		@ FLAG 02 set?
		XEQ 34		@ sort eigenvalues and -vectors
	GTO 99
LBL 30		@ shift zero element of VE from beginning to end
	INDEX "VE"
	1
	2
	STOIJ
	RCL 01		@ n
	1ᴇ3
	÷
	2
	+
	STO 05          @ i=2; i<=n; i++
	LBL 31          @ do {
		RCLEL           @ e[i]
		J-
		STOEL           @ e[i-1] = e[i]
		J+
		J+
		ISG 05          @ } while (i<=n)
			GTO 31
	RTN
LBL 32		@ initialize MZ identity matrix
	RCL 01		@ n
	ENTER
	NEWMAT
	STO "MZ"	@ z = n x n matrix
	INDEX "MZ"
	RCL 01		@ n
	1ᴇ3
	÷
	1
	+
	STO 02		@ l=1; l<=n; l++
	LBL 33		@ do {
		RCL 02
		ENTER
		STOIJ
		1
		STOEL		@ z[l][l] = 1
		ISG 02		@ } while (l<=n)
			GTO 33
	RTN
LBL 34		@ sort eigenvalues and -vectors (based on nrecipes eigsrt)
	RCL "MZ"
	TRANS
	STO "MZ"	@ need eigenvectors as row vectors
	RCL 01		@ n
	1
	-
	1ᴇ3
	÷
	1
	+
	STO 05		@ i=1; i<n; i++
	LBL 35		@ do {
		RCL 05		@ i
		STO 06		@ k = i
		INDEX "VD"
		1
		X<>Y
		STOIJ
		RCLEL		@ d[k]
		STO 10		@ p = d[k=i]
		RCL 01		@ n
		1ᴇ3
		÷
		1
		RCL 05		@ i
		IP
		+
		+
		STO 04		@ m=i+1; m<=n; m++ (use m instead of j)
		LBL 36		@ do {
			1
			RCL 04		@ m
			STOIJ
			RCL 10		@ p
			RCLEL		@ d[m]
			X>Y?		@ if !(d[m] <= p)
				GTO 37		@ skip if case step
			1
			RCL 04		@ m
			STO 06		@ k = m
			STOIJ
			RCLEL		@ d[k=m]
			STO 10		@ p = d[k=m]
			LBL 37
			ISG 04		@ } while (m<=n)
				GTO 36
		RCL 05		@ i iterator
		IP		@ i
		RCL 06		@ k iterator
		IP		@ k
		X=Y?		@ if !(k != i)
			GTO 38		@ no swapping necessary
		1
		RCL 05		@ i
		STOIJ
		RCLEL		@ d[i]
		1
		RCL 06		@ k
		STOIJ
		R↓
		R↓
		STOEL		@ d[k] = d[i]
		1
		RCL 05		@ i
		STOIJ
		RCL 10		@ p
		STOEL		@ d[i] = 0
		FS? 01		@ FLAG 01 set?
			GTO 38		@ skip eigenvector sorting
		RCL 05		@ i
		RCL 06		@ k
		INDEX "MZ"
		R<>R		@ swap rows i and k of MZ
		LBL 38
		ISG 05		@ } while (i<n)
			GTO 35
	RCL "MZ"
	TRANS
	STO "MZ"	@ restore eigenvectors back as column vectors
	RTN
LBL 91		@ nrerror("Too many iterations in tqli")
	"Err: Too "
	├"many iter"
	AVIEW
LBL 99
	RTN
END