			% P = transition matrix such that A * P = P * D
			% A = diagonal matrix of eigenvalues
			function [P, D, count] = eig_cj(A)
			[Ni, Nj] = size(A);
			assert(Ni == Nj, "The input must be diagonal");
			N = Ni;
			P = eye(N);
			D = A; % copy matrix A

			count = 0;

			while off(D, N) > (1e-9 * N^2)
			for p = 1:(N-1)
			for q = (p+1):N
			[c, s] = from(p, q, D);
			[D, P] = transition(D, P, p, q, c, s, N);
			end
			end
			count = count + 1;
			end
			end

			% off norm of A, A is a square matrix, N = size of A
			function out = off(A, N)
			out = 0;
			for i = 1:N
			for j = 1:N
			if i ~= j
			out = out + A(i,j)^2;
			end
			end
			end

			out = sqrt(out);
			end

			% cos and sin deduction from p, q and a square matrix A
			function [c, s] = from(p, q, A)
			if abs(A(p,q)) > eps
			tau = (A(q,q) - A(p,p)) / (2 * A(p,q));
			if tau >= 0
			tau = -tau + sqrt(1 + tau^2);
			else
			tau = -tau - sqrt(1 + tau^2);
			end
			c = 1.0 / sqrt(1 + tau^2);
			s = tau * c;
			else
			c = 1.0;
			s = 0.0;
			end
			end

			% apply a rotation to A: J(p,q)' A J(p,q)
			% hope it takes A and P by reference
			function [A, P] = transition(A, P, p, q, c, s, N)
			% left multiplication by J'
			for j = 1:N
			t = A(p,j);
			A(p,j) = c * t - s * A(q,j);
			A(q,j) = s * t + c * A(q,j);
			end

			% right multiplication by J
			for i = 1:N
			t = A(i,p);
			A(i,p) = c * t - s * A(i,q);
			A(i,q) = s * t + c * A(i,q);

			tp = P(i,p);
			P(i,p) = c * tp - s * P(i,q);
			P(i,q) = s * tp + c * P(i,q);
			end
			end

compphys/exo3/eig_cj.m
33bba04159d6master

eig_cj.m
View Options