Disclaimer: I previously posted this on SO, but though it would be more relevant for scicomp. The original post has been deleted.
I have been trying to recreate the functionality provided by MATLABs ordschur function by calling LAPACK directly (MATLAB uses LAPACKS dtrsen function to do its computation). More specifically, the ordschurfunction that orders the eigenvalues according to those on the interior of the unit disc. More information about MATLAB's ordschur function can be found here. More information about LAPACK's dtrsen function can be found here.
However, despite selecting the desired eigenvalues, the output of LAPACKs dtrsen doesn't match that of MATLABs ordschur. Any suggestions for what the issue could be? Changing the values for SELECT manually seems to have no influence on the output.
Below is an example matrix and the associated MATLAB code.
I'm trying to compute the ordered schur factorization of the following sympletic matrix
$Z = \begin{bmatrix}A + BR^{-1}B^T(A^{-1})^TQ & -BR^{-1}B^T(A^{-1})^T \\(A^{-1})^TQ & (A^{-1})^T\end{bmatrix}$
where $A$ and $B$ are matrices for a discrete-time dynamical system, where it ensured that the pair is stabilizable. $Q$ and $R$ are identity matrices. The ordered schur factorization of matrix $Z$ is used to compute the solution to the discrete algebraic riccati equation (see Algebraic Riccati equation wiki)
The following is a sympletic matrix, where the matrices used to generate it fulfil the requirements written above.
A = [1.1047 0.0301;
0.1954 1.0838]
B = [0.2131;
0.2272]
Q = eye(2)
R = 4
Z = [A + B*inv(R)*B'*inv(A)'*Q, -B*inv(R)*B'*inv(A)';
-inv(A)'*Q, inv(A)'];
%Z = [1.1147 0.0394 -0.0100 -0.0094;
% 0.2061 1.0938 -0.0107 -0.0100;
% -0.9097 0.1640 0.9097 -0.1640;
% 0.0252 -0.9272 -0.0252 0.9272]
LAPACK is being called from MATLAB using this interface (lapack.m).
% Computing the Schur factorization using MATLABs schur(A)
[UU,TT] = schur(Z);
P = size(Z,1);
% Settings for LAPACK dtrsen
JOB = 'N';
COMPQ = 'N';
SELECT = abs(ordeig(TT)) < 1; SELECT = SELECT.*ones(size(Z,1),1); % Flagging all eigenvalues within the unit disk (replicating ordshur setting 'udi')
N = P;
T = TT;
LDT = P;
QQ = UU;
LDQ = P;
WR = zeros(P,1);
WI = zeros(P,1);
M = 1 ;
S = 0;
SEP = 0;
WORK = zeros(P,1);
LWORK = P;
IWORK = zeros(P,1);
LIWORK = 1;
INFO = 0;
% Computing ordschur using LAPACK dtrsen
out = lapack('dtrsen',JOB,COMPQ,SELECT,N,T,LDT,QQ,LDQ,WR,WI,M,S,SEP,WORK,LWORK,IWORK,LIWORK,INFO);
TTS = out{5};
UUS = out{7};
% Setting any of following allows for the ordering of the eigenvalues according to:
% 'lhp' left-half plane (real(E)<0)
% 'rhp' right-half plane (real(E)>0)
% 'udi' interior of unit disk (abs(E)<1)
% 'udo' exterior of unit disk (abs(E)>1)
[US,TS] = ordschur(UU,TT,'udi');
if all(all((UUS*TTS*UUS' - Z) < 1e-6))
disp('Correct factorization from FORTRAN');
else
disp('!! Incorrect factorization from FORTRAN');
end
if all(all((US*TS*US' - Z) < 1e-6))
disp('Correct factorization from MATLAB');
else
disp('!! Incorrect factorization from MATLAB');
end
Output from calling LAPACK
TTS =
1.0486 -0.1820 -0.6761 -0.7283
0 1.2304 -0.6545 0.5104
0 0 0.9536 0.1412
0 0 0 0.8128
>> UUS
UUS =
-0.1641 0.1657 -0.9587 -0.1626
-0.3025 0.2208 -0.0669 0.9248
-0.1622 -0.9588 -0.1658 0.1639
0.9248 -0.0665 -0.2211 0.3024
Built-in MATLAB ordschur
TS =
0.9536 -0.1335 0.1414 -0.9838
0 0.8128 -0.8331 0.0047
0 0 1.2304 0.1721
0 0 0 1.0486
>> US
US =
-0.0331 -0.0229 -0.4356 -0.8992
0.0111 -0.0370 -0.8991 0.4361
-0.8994 -0.4352 0.0229 0.0332
0.4357 -0.8993 0.0371 -0.0111
Trying to recompute $Z$ shows that the LAPACK factorization is incorrect, compared to the built-in ordschur
>> UUS*TTS*UUS'
ans =
0.9515 0.1060 -0.0555 0.0077
-0.0650 1.1627 -0.0850 0.0109
-0.6754 -0.4831 0.9916 -0.2879
0.6395 -0.7514 0.1332 0.9396
>> US*TS*US'
ans =
1.1147 0.0394 -0.0100 -0.0094
0.2061 1.0938 -0.0107 -0.0100
-0.9097 0.1640 0.9097 -0.1640
0.0252 -0.9272 -0.0252 0.9272
