BALANCE2 Balancing generalized eigenvalue problem

Balancing generalized eigenvalue problem using Ward's algorithm.

BALANCE2 Balancing generalized eigenvalue problem
[T1 T2] = BALANCE(A,B) computes matrices T1 and T2 that balance
matrices A and B for the generalized eigenvalue problem by T1*A*T2 and
T1*B*T2 according to [War1981].

[T1 T2] = BALANCE(A,B,rho) uses rho as computer radix. The default is
rho = 2.

[T1 T2 ILO IHI] = BALANCE(A,B) returns the indices ILO and IHI.


Example 1, see [War1981, pp.148f.]:

A = [-2.0e+1 -1.0e+4 -2.0e+0 -1.0e+6 -1.0e+1 -2.0e+5
6.0e-3 4.0e+0 6.0e-4 2.0e+2 3.0e-3 3.0e+1
-2.0e-1 -3.0e+2 -4.0e-2 -1.0e+4 0.0e+0 3.0e+3
6.0e-5 4.0e-2 9.0e-6 9.0e+0 3.0e-5 5.0e-1
6.0e-2 5.0e+1 8.0e-3 -4.0e+3 8.0e-2 0.0e+0
0.0e+0 1.0e+3 7.0e-1 -2.0e+5 1.3e+1 -6.0e+4 ];

B = [-2.0e+1 -1.0e+4 2.0e+0 -2.0e+6 1.0e+1 -1.0e+5
5.0e-3 3.0e+0 -2.0e-4 4.0e+2 -1.0e-3 3.0e+1
0.0e+0 -1.0e+2 -8.0e-2 2.0e+4 -4.0e-1 0.0e+0
5.0e-5 3.0e-2 2.0e-6 4.0e+0 2.0e-5 1.0e-1
4.0e-2 3.0e+1 -1.0e-3 3.0e+3 -1.0e-2 6.0e+2
-1.0e+0 0.0e+0 4.0e-1 -1.0e+5 4.0e+0 2.0e+4 ];

[T1 T2] = balance2(A,B);
lambda1 = sort(eig(A,B)); lambda2 = sort(eig(T1*A*T2,T1*B*T2));
fprintf(1,'%-16s %-16s\n','unbalanced','balanced');
for i = 1:size(A,1),
fprintf(1,'%+16.15f %+16.15f\n',lambda1(i),lambda2(i));
end


Example 2, see [War1981, pp.148f.]:

A = full(spdiags(10.^[1:2:13; 0:2:12; -1:2:11]',[-1:1],7,7));
B = eye(7);

[T1 T2] = balance2(A,B);
lambda1 = eig(A,B); lambda2 = eig(T1*A*T2,T1*B*T2);
fprintf(1,'%-16s %-16s\n','unbalanced','balanced');
for i = 1:size(A,1),
fprintf(1,'%+16.15e %+16.15e\n',lambda1(i),lambda2(i));
end


References:
[War1981] Ward, R.C.: "Balancing the generalized eigenvalue problem",
SIAM Journal on Scientific and Statistical Computing, Vol.
2, No. 2, pp.141-152, 1981.

See also balance and eig.