Skip to content

Commit f9dc340

Browse files
VasileSima4VasileSima4
committed
Populated the directory
1 parent b7945d4 commit f9dc340

332 files changed

Lines changed: 59961 additions & 0 deletions

File tree

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

‎Contents.m‎

Lines changed: 431 additions & 0 deletions
Large diffs are not rendered by default.

‎HaeigZ.m‎

Lines changed: 147 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,147 @@
1+
% HAEIGZ.F - Gateway function for computing a balancing transformation
2+
% or the eigenvalues of a complex Hamiltonian matrix, using
3+
% SLICOT Library routines MB04DZ and MB03XZ.
4+
% The gateway accepts real or complex input matrices.
5+
%
6+
% Matlab call:
7+
% [W(,Ao,Go(,U1,U2(,l,scal)))] = HaeigZ(A,QG(,job,jobu,balanc))
8+
% [W(,Se,De(,Ue1,Ue2(,l,scal)))] = HaeigZ(A,QG(,job,jobu,balanc))
9+
% [l,scal(,Ab,QGb)] = HaeigZ(A,QG,job,balanc)
10+
%
11+
% [l,scal] = HaeigZ(A,QG,-1,balanc)
12+
% [l,scal,Ab,QGb] = HaeigZ(A,QG,-1,balanc)
13+
% [W] = HaeigZ(A,QG)
14+
% [W,Se,De] = HaeigZ(A,QG)
15+
% [W,Se,De,Ue1,Ue2] = HaeigZ(A,QG,0,1)
16+
% [W,Se,De,Ue1,Ue2,l,scal] = HaeigZ(A,QG,0,1,balanc)
17+
% [W,Ao] = HaeigZ(A,QG,1)
18+
% [W,Ao,U1,U2] = HaeigZ(A,QG,1,1)
19+
% [W,Ao,Go] = HaeigZ(A,QG,2)
20+
% [W,Ao,Go,U1,U2] = HaeigZ(A,QG,2,1)
21+
% [W,Ao,Go,U1,U2,l,scal] = HaeigZ(A,QG,2,1,balanc)
22+
%
23+
% HAEIGZ computes the eigenvalues of a complex n-by-n Hamiltonian matrix H,
24+
% with
25+
%
26+
% [ A G ] H H
27+
% H = [ H ], G = G , Q = Q , (1)
28+
% [ Q -A ]
29+
%
30+
% where A, G and Q are complex n-by-n matrices.
31+
%
32+
% Due to the structure of H, if lambda is an eigenvalue, then
33+
% -conjugate(lambda) is also an eigenvalue. This does not mean that
34+
% purely imaginary eigenvalues are necessarily multiple. The function
35+
% computes the eigenvalues of H using an embedding to a skew-
36+
% Hamiltonian matrix He,
37+
%
38+
% [ Ae Ge ] T T
39+
% He = [ T ], Ge = -Ge , Qe = -Qe , (2)
40+
% [ Qe Ae ]
41+
%
42+
% where Ae, Ge, and Qe are real 2*n-by-2*n matrices. Then, an
43+
% orthogonal symplectic matrix Ue is used to reduce He to the
44+
% structured real Schur form,
45+
%
46+
% T [ Se De ] T
47+
% Ue He Ue = [ T ], De = -De , (3)
48+
% [ 0 Se ]
49+
%
50+
% where Ue is a 4n-by-4n real symplectic matrix, and Se is upper
51+
% quasi-triangular (real Schur form).
52+
%
53+
% Optionally, if job > 0, the matrix i*He is further transformed to
54+
% the structured complex Schur form
55+
%
56+
% H [ Ao Go ] H
57+
% U (i*He) U = [ H ], Go = Go , (4)
58+
% [ 0 -Ao ]
59+
%
60+
% where U is a 4n-by-4n unitary symplectic matrix, and Ao is upper
61+
% triangular (Schur form). Optionally, if jobu = 1, the unitary
62+
% symplectic transformation matrix
63+
%
64+
% ( U1 U2 )
65+
% U = ( )
66+
% ( -U2 U1 )
67+
%
68+
% is computed.
69+
%
70+
% If job = -1, an accurate similarity transformation T such that
71+
% Hb = T\H*T has, as nearly as possible, approximately equal row and
72+
% column norms. T is a permutation of a diagonal matrix and
73+
% symplectic. T is stored in an n-vector scal as described in MB04DZ.
74+
%
75+
% Description of input parameters:
76+
% A - the n-by-n matrix A.
77+
% QG - an n-by-(n+1) matrix containing the triangles of the
78+
% Hermitian matrices G and Q, as follows:
79+
% the leading n-by-n lower triangular part contains the
80+
% lower triangle of the matrix Q, and the n-by-n upper
81+
% triangular part of the submatrix in the columns 2 to n+1
82+
% contains the upper triangle of the matrix G of H in (1).
83+
% So, if i >= j, then Q(i,j) = conj(Q(j,i)) is stored in
84+
% QG(i,j) and G(j,i) = conj(G(i,j)) is stored in QG(j,i+1).
85+
% QG is an empty matrix if n = 0.
86+
% job - (optional) scalar indicating the computation to be
87+
% performed, as follows:
88+
% = -1 : compute a balancing transformation only;
89+
% = 0 : compute the eigenvalues only (default);
90+
% = 1 : compute the eigenvalues and the matrix Ao in (4);
91+
% = 2 : compute the eigenvalues and the matrices Ao and Go
92+
% in (4).
93+
% jobu - (optional) if job > 0, scalar indicating whether the
94+
% unitary transformation matrix U is returned, as follows:
95+
% = 0 : U is not required (default);
96+
% = 1 : on exit, U contains the unitary transformation
97+
% matrix.
98+
% balanc - determines whether H should be permuted (balanc = 1),
99+
% scaled (balanc = 2), or permuted and scaled (balanc = 3)
100+
% prior to eigenvalue computations. Otherwise balanc = 0
101+
% (default). This parameter is optional if job >= 0, but
102+
% compulsory if job < 0.
103+
%
104+
% Description of output parameters:
105+
% l - if job = -1 or balanc > 0, an integer determined when H was
106+
% balanced. The balanced A(I,J) = 0 if I > J and J = 1:l-1.
107+
% The balanced Q(I,J) = 0 if J = 1:l-1 or I = 1:l-1.
108+
% scal - if job = -1 or balanc > 0, an n-vector containing details
109+
% of the permutation and/or scaling factors applied when
110+
% balancing. See MB04DZ, for details.
111+
% Ab - if job = -1, the matrix A of the balanced Hamiltonian Hb.
112+
% The lower triangular part of the first l-1 columns of A is
113+
% zero.
114+
% QGb - if job = -1, the matrices Q and G of the balanced
115+
% Hamiltonian Hb, stored compactly. The lower triangular and
116+
% diagonal part of the first l-1 columns of QGb is zero.
117+
% W - the 2n-vector of the eigenvalues of the matrix H.
118+
% Se - if job = 0, the computed 2n-by-2n upper real Schur
119+
% submatrix Se in (3).
120+
% De - if job = 0, the computed 2n-by-2n skew-symmetric submatrix
121+
% De in (3).
122+
% Ue1 - if job = 0, the computed 2n-by-2n (1,1) block of the matrix
123+
% Ue in (3).
124+
% Ue2 - if job = 0, the computed 2n-by-2n (2,1) block of the matrix
125+
% Ue in (3).
126+
% Ao - if job > 0, the computed 2n-by-2n upper triangular
127+
% submatrix Ao in (4).
128+
% Go - if job = 2, the computed 2n-by-2n Hermitian matrix Go
129+
% in (4).
130+
% U1 - if jobu = 1, a 2n-by-2n matrix containing the computed
131+
% matrix U1.
132+
% U2 - if jobu = 1, an 2n-by-2n matrix containing the computed
133+
% matrix U2.
134+
%
135+
% Contributors:
136+
% V. Sima, Research Institute for Informatics, Bucharest, Nov. 2011.
137+
%
138+
139+
% RELEASE 2.0 of SLICOT Basic Systems and Control Toolbox.
140+
% Based on SLICOT RELEASE 5.7, Copyright (c) 2002-2020 NICONET e.V.
141+
%
142+
% Contributor:
143+
% V. Sima, Research Institute for Informatics, Bucharest, July 2012.
144+
%
145+
% Revisions:
146+
% V. Sima, Nov. 2012.
147+
%

‎HaeigZ.mexw64‎

240 KB
Binary file not shown.

‎Hameig.m‎

Lines changed: 72 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,72 @@
1+
function W = Hameig( A, QG, scale )
2+
%HAMEIG computes the eigenvalues of a Hamiltonian matrix given in a
3+
% compressed form.
4+
%
5+
% W = HAMEIG(A,QG,SCALE) computes the eigenvalues of a Hamiltonian
6+
% matrix
7+
%
8+
% ( A G )
9+
% H = ( T ), with A, Q, and G n-by-n. (1)
10+
% ( Q -A )
11+
%
12+
% The symmetric matrices Q and G are stored in the n-by-(n+1) array QG,
13+
% as follows: the leading n-by-n lower triangular part contains the
14+
% lower triangle of the matrix Q, and the n-by-n upper triangular part
15+
% of the submatrix in the columns 2 to n+1 contains the upper triangle
16+
% of the matrix G of H in (1). So, if i >= j, then Q(i,j) = Q(j,i)
17+
% is stored in QG(i,j) and G(i,j) = G(j,i) is stored in QG(j,i+1).
18+
%
19+
% SCALE is a scalar specifying whether or not balancing operations
20+
% should be performed when computing the eigenvalues, as follows:
21+
% SCALE = 0 : do not use balancing;
22+
% SCALE = 1 : do scaling in order to equilibrate the rows and columns.
23+
% Default: SCALE = 1.
24+
%
25+
% W is an 2*n-vector containing the eigenvalues of H.
26+
%
27+
% Comments
28+
% Eigenvalues are stored in W as follows: The first n entries of W are
29+
% the computed eigenvalues of H with non-negative real part. They are
30+
% stored in decreasing order of magnitude of the real parts, i.e.,
31+
% real( W(I) ) >= real( W(I+1) ). (In particular, an eigenvalue closest
32+
% to the imaginary axis is W(N).) In addition, eigenvalues with zero real
33+
% part are sorted in decreasing order of magnitude of imaginary parts.
34+
% Note that non-real eigenvalues with non-zero real part appear in
35+
% complex conjugate pairs, but eigenvalues with zero real part (in the
36+
% first n entries of W) do not, in general, appear in complex conjugate
37+
% pairs. The remaining n eigenvalues are the negatives of the first
38+
% n eigenvalues.
39+
%
40+
% See also HAMILEIG
41+
%
42+
43+
% RELEASE 2.0 of SLICOT Basic Systems and Control Toolbox.
44+
% Based on SLICOT RELEASE 5.7, Copyright (c) 2002-2020 NICONET e.V.
45+
%
46+
% V. Sima, Nov. 2002.
47+
% Revised: Mar. 2009.
48+
%
49+
50+
ni = nargin;
51+
%
52+
if ni < 2 || nargout < 1,
53+
error( [ 'Usage: W = Hameig(A,QG)', sprintf('\n'),...
54+
' W = Hameig(A,QG,0)' ] )
55+
end
56+
%
57+
[n, m] = size( A ); [n1, m1] = size( QG );
58+
if ~( m == n && n1 == n && m1 == n + min( n, 1 ) ),
59+
error( 'Arrays A and/or QG have wrong or incompatible sizes' )
60+
end
61+
%
62+
if ni < 3, scale = 1; end
63+
%
64+
job = 0;
65+
[ Wr, Wi ] = Hamileig( A, QG, job, scale );
66+
W = Wr + 1i*Wi; W(n+1:2*n) = -W;
67+
for ii = n + 1 : 2*n,
68+
if ~( Wr(ii-n) == 0 ), W(ii) = conj( W(ii) ); end
69+
end
70+
W = reshape( W, 2*n, min( n, 1 ) );
71+
%
72+
% end Hameig

‎Hamileig.m‎

Lines changed: 118 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,118 @@
1+
% HAMILEIG.F - MEX-function to compute the eigenvalues of a
2+
% Hamiltonian matrix using the square-reduced approach
3+
% implemented in SLICOT routines MB03SD and MB04ZD.
4+
%
5+
% [Ao,QGo(,U)] = Hamileig(A,QG,job(,compu,S)), job = -1;
6+
% [WR,WI(,Ao,QGo,U)] = Hamileig(A,QG(,job,jobscl,compu,S)), job = 0,
7+
% or job = 1.
8+
%
9+
% [Ao,QGo(,U)] = Hamileig(A,QG,-1(,compu,S))
10+
% [WR,WI(,Ao,QGo,U)] = Hamileig(A,QG(,0, jobscl,compu,S))
11+
% [WR,WI] = Hamileig(A,QG, 1(,jobscl))
12+
%
13+
% HAMILEIG transforms a Hamiltonian matrix
14+
%
15+
% ( A G )
16+
% H = ( T ) (1)
17+
% ( Q -A )
18+
%
19+
% (with G and Q symmetric matrices) into a square-reduced Hamiltonian
20+
% matrix
21+
%
22+
% ( A' G' )
23+
% H' = ( T ), (2)
24+
% ( Q' -A' )
25+
% T
26+
% by an orthogonal symplectic similarity transformation H' = U H U,
27+
% where
28+
%
29+
% ( U1 U2 )
30+
% U = ( ), (3)
31+
% ( -U2 U1 )
32+
%
33+
% and computes the eigenvalues of H. Therefore, H' is such that
34+
%
35+
% 2 ( A'' G'' )
36+
% H' = ( T ), (4)
37+
% ( 0 A'' )
38+
%
39+
% with A'' upper Hessenberg and G'' skew symmetric. The square roots
40+
% of the eigenvalues of A'' = A'*A' + G'*Q' are the eigenvalues of H.
41+
%
42+
% Description of input parameters:
43+
% A - the n-by-n matrix A.
44+
% QG - an n-by-(n+1) matrix containing the triangles of the
45+
% symmetric matrices Q and G, as follows:
46+
% the leading n-by-n lower triangular part contains the lower
47+
% triangle of the matrix Q, and the n-by-n upper triangular
48+
% part of the submatrix in the columns 2 to n+1 contains the
49+
% upper triangle of the matrix G of H in (1).
50+
% So, if i >= j, then Q(i,j) = Q(j,i) is stored in QG(i,j)
51+
% and G(i,j) = G(j,i) is stored in QG(j,i+1).
52+
% QG is an empty matrix if n = 0.
53+
% job - (optional) scalar indicating the computation to be
54+
% performed, as follows:
55+
% =-1 : compute the square-reduced matrix H' in (2);
56+
% = 0 : compute the eigenvalues of H (default);
57+
% = 1 : compute the eigenvalues of H, assuming that the
58+
% given Hamiltonian matrix is already in the reduced
59+
% form (2).
60+
% jobscl - (optional) if job >= 0, scalar specifying whether or not
61+
% balancing operations should be performed when computing the
62+
% eigenvalues of A'', as follows:
63+
% = 0 : do not use balancing;
64+
% = 1 : do scaling in order to equilibrate the rows
65+
% and columns of A'' (default).
66+
% If job = -1, jobscl is not used.
67+
% compu - (optional) if job <= 0, scalar indicating whether the
68+
% orthogonal symplectic similarity transformation matrix U is
69+
% returned or accumulated into an orthogonal symplectic
70+
% matrix, or if the transformation matrix is not required,
71+
% as follows:
72+
% = 0 : U is not required (default);
73+
% = 1 : on entry, U need not be set;
74+
% on exit, U contains the orthogonal symplectic
75+
% matrix U;
76+
% = 2 : the orthogonal symplectic similarity transformations
77+
% are accumulated into U;
78+
% on input, U must contain an orthogonal symplectic
79+
% matrix S;
80+
% on exit, U contains S*U.
81+
% If job = 1, compu is not used.
82+
% S - (optional) if job <= 0 and compu = 2, an n-by-2*n matrix
83+
% containing the first n rows of the given orthogonal
84+
% symplectic matrix S.
85+
%
86+
% Description of output parameters:
87+
% WR,WI - if job >= 0, the n-vectors of real parts and imaginary
88+
% parts, respectively, of the n computed eigenvalues of H'
89+
% with non-negative real part. The remaining n eigenvalues
90+
% are the negatives of these eigenvalues.
91+
% Eigenvalues are stored in WR and WI in decreasing order of
92+
% magnitude of the real parts, i.e., WR(I) >= WR(I+1).
93+
% (In particular, an eigenvalue closest to the imaginary
94+
% axis is WR(N)+WI(N)i.)
95+
% In addition, eigenvalues with zero real part are sorted in
96+
% decreasing order of magnitude of imaginary parts. Note
97+
% that non-real eigenvalues with non-zero real part appear
98+
% in complex conjugate pairs, but eigenvalues with zero real
99+
% part do not, in general, appear in complex conjugate
100+
% pairs.
101+
% Ao - if job <= 0, the computed n-by-n submatrix A' of H' in (2).
102+
% QGo - if job <= 0, the computed n-by-(n+1) matrix containing the
103+
% triangles of the symmetric matrices Q' and G' of H' in (2),
104+
% stored in the same way as the initial matrices Q and G.
105+
% U - if job <= 0 and compu > 0, an n-by-2*n matrix containing
106+
% the first n rows of the computed orthogonal symplectic
107+
% matrix U.
108+
%
109+
110+
% RELEASE 2.0 of SLICOT Basic Systems and Control Toolbox.
111+
% Based on SLICOT RELEASE 5.7, Copyright (c) 2002-2020 NICONET e.V.
112+
%
113+
% Contributor:
114+
% V. Sima, Research Institute for Informatics, Bucharest, Nov. 2002.
115+
%
116+
% Revisions:
117+
% -
118+
%

‎Hamileig.mexw64‎

27 KB
Binary file not shown.

‎Hesscond.m‎

Lines changed: 41 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,41 @@
1+
function rcnd = Hesscond( H, normh )
2+
%HESSCOND Computes an estimate of the reciprocal of the condition
3+
% number of an upper Hessenberg matrix H.
4+
%
5+
% RCND = HESSCOND(H) computes an estimate of the reciprocal
6+
% of the condition number of H in the 1-norm.
7+
%
8+
% RCND = HESSCOND(H,NORMH) has an additional input argument:
9+
%
10+
% NORMH specifies whether the 1-norm or the infinity-norm
11+
% reciprocal condition number is required.
12+
% NORMH = 1: 1-norm;
13+
% NORMH = 2: infinity-norm.
14+
% Default: NORMH = 1.
15+
%
16+
% See also HESSOL, HESSL
17+
%
18+
19+
% RELEASE 2.0 of SLICOT Basic Systems and Control Toolbox.
20+
% Based on SLICOT RELEASE 5.7, Copyright (c) 2002-2020 NICONET e.V.
21+
%
22+
% V. Sima, Dec. 2002.
23+
% Revised: Oct. 2004, Mar. 2009.
24+
%
25+
26+
ni = nargin;
27+
%
28+
if ni < 1,
29+
error('Usage: RCND = HESSCOND(H,NORMH)')
30+
elseif ni == 1,
31+
normh = 1;
32+
end
33+
%
34+
if isreal( H ) == 1,
35+
mtype = 0;
36+
else
37+
mtype = 1;
38+
end
39+
rcnd = Hessol(2,H,mtype,normh);
40+
%
41+
% end Hesscond

0 commit comments

Comments
 (0)