3
$\begingroup$

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
$\endgroup$
2
  • 1
    $\begingroup$ Please, check if scaling is the source of the difference. Eigenvalues/eigenvectors are not unique. $\endgroup$
    – Anton Menshov
    Aug 7, 2023 at 19:34
  • $\begingroup$ @AntonMenshov thanks for the reply. I'm unsure where this scaling should take place, when ordering the eigenvalues, it seems that the only requirement is that the ones inside the disk is to be selected. Nothing about scale. $\endgroup$
    – two_Thomas
    Aug 8, 2023 at 11:52

2 Answers 2

3
$\begingroup$

First of all, please include a definition of Z in your minimal example so that we can reproduce your output.

The bug clearly seems to be in the Fortran code. Have you verified that $Z \approx UTU^*$ for both outputs?

A possible issue is that you are passing a vector of Matlab's bools as SELECT, while Fortran expects a logical value. Fortran's logical values do not necessarily map to 0 and 1 as in C or Matlab; rather, the mapping is compiler-dependent. How did you compile and link the Lapack interface? With which compiler?

I am not sure this will solve the problem; I am just suggesting some ideas.

$\endgroup$
1
  • 1
    $\begingroup$ Thanks for the reply. I have added additional context in the post about the matrix $Z$. You're correct, I can see now that the FORTRAN/LAPACK output for some matrices does not yield a correct factorisation (see code above and the outputs). I have tried several ways of setting the SELECTvalues, ranging from integers to true/false and also MATLABs logical values. The LAPACK library was compiled using mingw. Essentially I'm using this library to call the LAPACK functions. Using the instructions i compiled it to a mex file. $\endgroup$
    – two_Thomas
    Aug 8, 2023 at 11:57
2
$\begingroup$

You should use COMPQ = 'V'; instead of COMPQ = 'N';

accroding to the official documentation of LAPACK dtrsen,

COMPQ = 'V' : update the matrix Q of Schur vectors;(i.e. the 7th input argument in dtrsen function)

COMPQ = 'N' : do not update Q.(i.e. the 7th input argument in dtrsen )

The updated MATLAB script is:

clear all;close all;clc;

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]


% Computing the Schur factorization using MATLABs schur(A)
[UU,TT] = schur(Z);
P = size(Z,1);

% Settings for LAPACK dtrsen
JOB = 'N';
COMPQ = 'V';
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
$\endgroup$
3
  • $\begingroup$ Thanks! This fixed the factorization issue. However, the resulting factorization from dtrsen still doesn't result in the same matrices as MATLAB's ordschur. At least, it doesn't seem to be selecting all eigenvalues within the unit disc. Any suggestions? $\endgroup$
    – two_Thomas
    Aug 28, 2023 at 6:24
  • $\begingroup$ @two_Thomas If I set SELECT = [false;false;false;false;false;true;true;]; the result seems true. I don't know how to explain it. The only thing I am certain of is SELECT does matter to get correct result. $\endgroup$
    – 138 Aspen
    Aug 28, 2023 at 9:06
  • $\begingroup$ I'm really unsure how this variable then works. Especially its dimension and which eigenvalues need to be selected. Do you have any intuition regarding how you selected these values? $\endgroup$
    – two_Thomas
    Aug 28, 2023 at 18:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.