Menu
▾
▴
[Octave-cvsupdate] SF.net SVN: octave:[7713] trunk/octave-forge/main/control
|
From: <par...@us...> - 2010-09-14 10:59:45
|
Revision: 7713
http://octave.svn.sourceforge.net/octave/?rev=7713&view=rev
Author: paramaniac
Date: 2010-09-14 10:59:36 +0000 (Tue, 14 Sep 2010)
Log Message:
-----------
control: add lyapchol and dlyapchol solvers
Modified Paths:
--------------
trunk/octave-forge/main/control/INDEX
trunk/octave-forge/main/control/inst/dlyap.m
trunk/octave-forge/main/control/inst/lyap.m
trunk/octave-forge/main/control/inst/test_control.m
trunk/octave-forge/main/control/src/Makefile
Added Paths:
-----------
trunk/octave-forge/main/control/inst/dlyapchol.m
trunk/octave-forge/main/control/inst/lyapchol.m
trunk/octave-forge/main/control/src/SB03OD.f
trunk/octave-forge/main/control/src/SG03BD.f
trunk/octave-forge/main/control/src/SG03BU.f
trunk/octave-forge/main/control/src/SG03BV.f
trunk/octave-forge/main/control/src/SG03BW.f
trunk/octave-forge/main/control/src/SG03BX.f
trunk/octave-forge/main/control/src/SG03BY.f
trunk/octave-forge/main/control/src/slsb03od.cc
trunk/octave-forge/main/control/src/slsg03bd.cc
Modified: trunk/octave-forge/main/control/INDEX
===================================================================
--- trunk/octave-forge/main/control/INDEX 2010-09-14 09:11:45 UTC (rev 7712)
+++ trunk/octave-forge/main/control/INDEX 2010-09-14 10:59:36 UTC (rev 7713)
@@ -92,7 +92,9 @@
care
dare
dlyap
+ dlyapchol
lyap
+ lyapchol
Octave-only
isctrb
isobsv
Modified: trunk/octave-forge/main/control/inst/dlyap.m
===================================================================
--- trunk/octave-forge/main/control/inst/dlyap.m 2010-09-14 09:11:45 UTC (rev 7712)
+++ trunk/octave-forge/main/control/inst/dlyap.m 2010-09-14 10:59:36 UTC (rev 7713)
@@ -32,6 +32,8 @@
## AXA' - EXE' + B = 0 (Generalized Lyapunov Equation)
## @end group
## @end example
+##
+## @seealso{dlyapchol, lyap, lyapchol}
## @end deftypefn
## Author: Lukas Reichlin <luk...@gm...>
Added: trunk/octave-forge/main/control/inst/dlyapchol.m
===================================================================
--- trunk/octave-forge/main/control/inst/dlyapchol.m (rev 0)
+++ trunk/octave-forge/main/control/inst/dlyapchol.m 2010-09-14 10:59:36 UTC (rev 7713)
@@ -0,0 +1,95 @@
+## Copyright (C) 2010 Lukas F. Reichlin
+##
+## This file is part of LTI Syncope.
+##
+## LTI Syncope is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## LTI Syncope is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn{Function File} {@var{u} =} dlyapchol (@var{a}, @var{b})
+## @deftypefnx{Function File} {@var{u} =} dlyapchol (@var{a}, @var{b}, @var{e})
+## Compute Cholesky factor of discrete-time Lyapunov equations.
+## Uses SLICOT SB03OD and SG03BD by courtesy of NICONET e.V.
+## <http://www.slicot.org>
+##
+## @example
+## @group
+## A U' U A' - X + B B' = 0 (Lyapunov Equation)
+##
+## A U' U A' - E U' U E' + B B' = 0 (Generalized Lyapunov Equation)
+## @end group
+## @end example
+##
+## @seealso{dlyap, lyap, lyapchol}
+## @end deftypefn
+
+## Author: Lukas Reichlin <luk...@gm...>
+## Created: January 2010
+## Version: 0.2
+
+function [u, scale] = dlyapchol (a, b, e)
+
+ switch (nargin)
+ case 2
+
+ if (! isreal (a) || isempty (a) || ! issquare (a))
+ error ("dlyapchol: a must be real and square");
+ endif
+
+ if (! isreal (b) || isempty (b))
+ error ("dlyapchol: b must be real and square")
+ endif
+
+ if (rows (a) != rows (b))
+ error ("dlyapchol: a and b must have the same number of rows");
+ endif
+
+ [u, scale] = slsb03od (a.', b.', true);
+
+ ## NOTE: TRANS = 'T' not suitable because we need U' U, not U U'
+
+ case 3
+
+ if (! isreal (a) || isempty (a) || ! issquare (a))
+ error ("dlyapchol: a must be real and square");
+ endif
+
+ if (! isreal (e) || isempty (e) || ! issquare (e))
+ error ("dlyapchol: e must be real and square");
+ endif
+
+ if (! isreal (b) || isempty (b))
+ error ("dlyapchol: b must be real");
+ endif
+
+ if (rows (b) != rows (a) || rows (e) != rows (a))
+ error ("dlyapchol: a, b, e not conformal");
+ endif
+
+ [u, scale] = slsg03bd (a.', e.', b.', true);
+
+ ## NOTE: TRANS = 'T' not suitable because we need U' U, not U U'
+
+ otherwise
+ print_usage ();
+
+ endswitch
+
+ if (scale < 1)
+ warning ("dlyapchol: solution scaled by %g to prevent overflow", scale);
+ endif
+
+endfunction
+
+
+## TODO: add tests
Modified: trunk/octave-forge/main/control/inst/lyap.m
===================================================================
--- trunk/octave-forge/main/control/inst/lyap.m 2010-09-14 09:11:45 UTC (rev 7712)
+++ trunk/octave-forge/main/control/inst/lyap.m 2010-09-14 10:59:36 UTC (rev 7713)
@@ -32,6 +32,8 @@
## AXE' + EXA' + B = 0 (Generalized Lyapunov Equation)
## @end group
## @end example
+##
+## @seealso{lyapchol, dlyap, dlyapchol}
## @end deftypefn
## Author: Lukas Reichlin <luk...@gm...>
Added: trunk/octave-forge/main/control/inst/lyapchol.m
===================================================================
--- trunk/octave-forge/main/control/inst/lyapchol.m (rev 0)
+++ trunk/octave-forge/main/control/inst/lyapchol.m 2010-09-14 10:59:36 UTC (rev 7713)
@@ -0,0 +1,143 @@
+## Copyright (C) 2010 Lukas F. Reichlin
+##
+## This file is part of LTI Syncope.
+##
+## LTI Syncope is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## LTI Syncope is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn{Function File} {@var{u} =} lyapchol (@var{a}, @var{b})
+## @deftypefnx{Function File} {@var{u} =} lyapchol (@var{a}, @var{b}, @var{e})
+## Compute Cholesky factor of continuous-time Lyapunov equations.
+## Uses SLICOT SB03OD and SG03BD by courtesy of NICONET e.V.
+## <http://www.slicot.org>
+##
+## @example
+## @group
+## A U' U + U' U A' + B B' = 0 (Lyapunov Equation)
+##
+## A U' U E' + E U' U A' + B B' = 0 (Generalized Lyapunov Equation)
+## @end group
+## @end example
+##
+## @seealso{lyap, dlyap, dlyapchol}
+## @end deftypefn
+
+## Author: Lukas Reichlin <luk...@gm...>
+## Created: January 2010
+## Version: 0.2
+
+function [u, scale] = lyapchol (a, b, e)
+
+ switch (nargin)
+ case 2
+
+ if (! isreal (a) || isempty (a) || ! issquare (a))
+ error ("lyapchol: a must be real and square");
+ endif
+
+ if (! isreal (b) || isempty (b))
+ error ("lyapchol: b must be real and square")
+ endif
+
+ if (rows (a) != rows (b))
+ error ("lyapchol: a and b must have the same number of rows");
+ endif
+
+ [u, scale] = slsb03od (a.', b.', false);
+
+ ## NOTE: TRANS = 'T' not suitable because we need U' U, not U U'
+
+ case 3
+
+ if (! isreal (a) || isempty (a) || ! issquare (a))
+ error ("lyapchol: a must be real and square");
+ endif
+
+ if (! isreal (e) || isempty (e) || ! issquare (e))
+ error ("lyapchol: e must be real and square");
+ endif
+
+ if (! isreal (b) || isempty (b))
+ error ("lyapchol: b must be real");
+ endif
+
+ if (rows (b) != rows (a) || rows (e) != rows (a))
+ error ("lyapchol: a, b, e not conformal");
+ endif
+
+ [u, scale] = slsg03bd (a.', e.', b.', false);
+
+ ## NOTE: TRANS = 'T' not suitable because we need U' U, not U U'
+
+ otherwise
+ print_usage ();
+
+ endswitch
+
+ if (scale < 1)
+ warning ("lyapchol: solution scaled by %g to prevent overflow", scale);
+ endif
+
+endfunction
+
+
+%!shared U, U_exp, X, X_exp
+%!
+%! A = [ -1.0 37.0 -12.0 -12.0
+%! -1.0 -10.0 0.0 4.0
+%! 2.0 -4.0 7.0 -6.0
+%! 2.0 2.0 7.0 -9.0 ].';
+%!
+%! B = [ 1.0 2.5 1.0 3.5
+%! 0.0 1.0 0.0 1.0
+%! -1.0 -2.5 -1.0 -1.5
+%! 1.0 2.5 4.0 -5.5
+%! -1.0 -2.5 -4.0 3.5 ].';
+%!
+%! U = lyapchol (A, B);
+%!
+%! X = U.' * U; # use lyap at home!
+%!
+%! U_exp = [ 1.0000 0.0000 0.0000 0.0000
+%! 3.0000 1.0000 0.0000 0.0000
+%! 2.0000 -1.0000 1.0000 0.0000
+%! -1.0000 1.0000 -2.0000 1.0000 ].';
+%!
+%! X_exp = [ 1.0000 3.0000 2.0000 -1.0000
+%! 3.0000 10.0000 5.0000 -2.0000
+%! 2.0000 5.0000 6.0000 -5.0000
+%! -1.0000 -2.0000 -5.0000 7.0000 ];
+%!
+%!assert (U, U_exp, 1e-4);
+%!assert (X, X_exp, 1e-4);
+
+%!shared U, U_exp, X, X_exp
+%!
+%! A = [ -1.0 3.0 -4.0
+%! 0.0 5.0 -2.0
+%! -4.0 4.0 1.0 ].';
+%!
+%! E = [ 2.0 1.0 3.0
+%! 2.0 0.0 1.0
+%! 4.0 5.0 1.0 ].';
+%!
+%! B = [ 2.0 -1.0 7.0 ].';
+%!
+%! U = lyapchol (A, B, E);
+%!
+%! U_exp = [ 1.6003 -0.4418 -0.1523
+%! 0.0000 0.6795 -0.2499
+%! 0.0000 0.0000 0.2041 ];
+%!
+%!assert (U, U_exp, 1e-4);
Modified: trunk/octave-forge/main/control/inst/test_control.m
===================================================================
--- trunk/octave-forge/main/control/inst/test_control.m 2010-09-14 09:11:45 UTC (rev 7712)
+++ trunk/octave-forge/main/control/inst/test_control.m 2010-09-14 10:59:36 UTC (rev 7713)
@@ -41,6 +41,8 @@
test dlyap
test gram
test covar
+test lyapchol
+## test dlyapchol # TODO: add tests
## various oct-files
test place
Modified: trunk/octave-forge/main/control/src/Makefile
===================================================================
--- trunk/octave-forge/main/control/src/Makefile 2010-09-14 09:11:45 UTC (rev 7712)
+++ trunk/octave-forge/main/control/src/Makefile 2010-09-14 10:59:36 UTC (rev 7713)
@@ -1,7 +1,7 @@
all: slab08nd.oct slab13dd.oct slsb10hd.oct slsb10ed.oct slab13bd.oct \
slsb01bd.oct slsb10fd.oct slsb10dd.oct slsb03md.oct slsb04md.oct \
slsb04qd.oct slsg03ad.oct slsb02od.oct slab13ad.oct slab01od.oct \
- sltb01pd.oct
+ sltb01pd.oct slsb03od.oct slsg03bd.oct
# transmission zeros of state-space models
slab08nd.oct: slab08nd.cc
@@ -107,16 +107,29 @@
SB03OR.f SB03OY.f SB04PX.f MB04NY.f MB04OY.f \
SB03OV.f
-# Staircase form using orthogonal transformations
+# staircase form using orthogonal transformations
slab01od.oct: slab01od.cc
mkoctfile slab01od.cc \
AB01OD.f AB01ND.f MB03OY.f MB01PD.f MB01QD.f
-# Minimal realization of state-space models
+# minimal realization of state-space models
sltb01pd.oct: sltb01pd.cc
mkoctfile sltb01pd.cc \
TB01PD.f TB01XD.f TB01ID.f AB07MD.f TB01UD.f \
MB03OY.f MB01PD.f MB01QD.f
+# Cholesky factor of Lyapunov equations
+slsb03od.oct: slsb03od.cc
+ mkoctfile slsb03od.cc \
+ SB03OD.f select.f SB03OU.f SB03OT.f MB04ND.f \
+ MB04OD.f SB03OR.f SB03OY.f SB04PX.f MB04NY.f \
+ MB04OY.f SB03OV.f
+
+# Cholesky factor of generalized Lyapunov equations
+slsg03bd.oct: slsg03bd.cc
+ mkoctfile slsg03bd.cc \
+ SG03BD.f SG03BV.f SG03BU.f SG03BW.f SG03BX.f \
+ SG03BY.f MB02UU.f MB02UV.f
+
clean:
rm *.o core octave-core *.oct *~
Added: trunk/octave-forge/main/control/src/SB03OD.f
===================================================================
--- trunk/octave-forge/main/control/src/SB03OD.f (rev 0)
+++ trunk/octave-forge/main/control/src/SB03OD.f 2010-09-14 10:59:36 UTC (rev 7713)
@@ -0,0 +1,662 @@
+ SUBROUTINE SB03OD( DICO, FACT, TRANS, N, M, A, LDA, Q, LDQ, B,
+ $ LDB, SCALE, WR, WI, DWORK, LDWORK, INFO )
+C
+C SLICOT RELEASE 5.0.
+C
+C Copyright (c) 2002-2009 NICONET e.V.
+C
+C This program is free software: you can redistribute it and/or
+C modify it under the terms of the GNU General Public License as
+C published by the Free Software Foundation, either version 2 of
+C the License, or (at your option) any later version.
+C
+C This program is distributed in the hope that it will be useful,
+C but WITHOUT ANY WARRANTY; without even the implied warranty of
+C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+C GNU General Public License for more details.
+C
+C You should have received a copy of the GNU General Public License
+C along with this program. If not, see
+C <http://www.gnu.org/licenses/>.
+C
+C PURPOSE
+C
+C To solve for X = op(U)'*op(U) either the stable non-negative
+C definite continuous-time Lyapunov equation
+C 2
+C op(A)'*X + X*op(A) = -scale *op(B)'*op(B) (1)
+C
+C or the convergent non-negative definite discrete-time Lyapunov
+C equation
+C 2
+C op(A)'*X*op(A) - X = -scale *op(B)'*op(B) (2)
+C
+C where op(K) = K or K' (i.e., the transpose of the matrix K), A is
+C an N-by-N matrix, op(B) is an M-by-N matrix, U is an upper
+C triangular matrix containing the Cholesky factor of the solution
+C matrix X, X = op(U)'*op(U), and scale is an output scale factor,
+C set less than or equal to 1 to avoid overflow in X. If matrix B
+C has full rank then the solution matrix X will be positive-definite
+C and hence the Cholesky factor U will be nonsingular, but if B is
+C rank deficient then X may be only positive semi-definite and U
+C will be singular.
+C
+C In the case of equation (1) the matrix A must be stable (that
+C is, all the eigenvalues of A must have negative real parts),
+C and for equation (2) the matrix A must be convergent (that is,
+C all the eigenvalues of A must lie inside the unit circle).
+C
+C ARGUMENTS
+C
+C Mode Parameters
+C
+C DICO CHARACTER*1
+C Specifies the type of Lyapunov equation to be solved as
+C follows:
+C = 'C': Equation (1), continuous-time case;
+C = 'D': Equation (2), discrete-time case.
+C
+C FACT CHARACTER*1
+C Specifies whether or not the real Schur factorization
+C of the matrix A is supplied on entry, as follows:
+C = 'F': On entry, A and Q contain the factors from the
+C real Schur factorization of the matrix A;
+C = 'N': The Schur factorization of A will be computed
+C and the factors will be stored in A and Q.
+C
+C TRANS CHARACTER*1
+C Specifies the form of op(K) to be used, as follows:
+C = 'N': op(K) = K (No transpose);
+C = 'T': op(K) = K**T (Transpose).
+C
+C Input/Output Parameters
+C
+C N (input) INTEGER
+C The order of the matrix A and the number of columns in
+C matrix op(B). N >= 0.
+C
+C M (input) INTEGER
+C The number of rows in matrix op(B). M >= 0.
+C
+C A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
+C On entry, the leading N-by-N part of this array must
+C contain the matrix A. If FACT = 'F', then A contains
+C an upper quasi-triangular matrix S in Schur canonical
+C form; the elements below the upper Hessenberg part of the
+C array A are not referenced.
+C On exit, the leading N-by-N upper Hessenberg part of this
+C array contains the upper quasi-triangular matrix S in
+C Schur canonical form from the Shur factorization of A.
+C The contents of array A is not modified if FACT = 'F'.
+C
+C LDA INTEGER
+C The leading dimension of array A. LDA >= MAX(1,N).
+C
+C Q (input or output) DOUBLE PRECISION array, dimension
+C (LDQ,N)
+C On entry, if FACT = 'F', then the leading N-by-N part of
+C this array must contain the orthogonal matrix Q of the
+C Schur factorization of A.
+C Otherwise, Q need not be set on entry.
+C On exit, the leading N-by-N part of this array contains
+C the orthogonal matrix Q of the Schur factorization of A.
+C The contents of array Q is not modified if FACT = 'F'.
+C
+C LDQ INTEGER
+C The leading dimension of array Q. LDQ >= MAX(1,N).
+C
+C B (input/output) DOUBLE PRECISION array, dimension (LDB,N)
+C if TRANS = 'N', and dimension (LDB,max(M,N)), if
+C TRANS = 'T'.
+C On entry, if TRANS = 'N', the leading M-by-N part of this
+C array must contain the coefficient matrix B of the
+C equation.
+C On entry, if TRANS = 'T', the leading N-by-M part of this
+C array must contain the coefficient matrix B of the
+C equation.
+C On exit, the leading N-by-N part of this array contains
+C the upper triangular Cholesky factor U of the solution
+C matrix X of the problem, X = op(U)'*op(U).
+C If M = 0 and N > 0, then U is set to zero.
+C
+C LDB INTEGER
+C The leading dimension of array B.
+C LDB >= MAX(1,N,M), if TRANS = 'N';
+C LDB >= MAX(1,N), if TRANS = 'T'.
+C
+C SCALE (output) DOUBLE PRECISION
+C The scale factor, scale, set less than or equal to 1 to
+C prevent the solution overflowing.
+C
+C WR (output) DOUBLE PRECISION array, dimension (N)
+C WI (output) DOUBLE PRECISION array, dimension (N)
+C If FACT = 'N', and INFO >= 0 and INFO <= 2, WR and WI
+C contain the real and imaginary parts, respectively, of
+C the eigenvalues of A.
+C If FACT = 'F', WR and WI are not referenced.
+C
+C Workspace
+C
+C DWORK DOUBLE PRECISION array, dimension (LDWORK)
+C On exit, if INFO = 0, or INFO = 1, DWORK(1) returns the
+C optimal value of LDWORK.
+C
+C LDWORK INTEGER
+C The length of the array DWORK.
+C If M > 0, LDWORK >= MAX(1,4*N + MIN(M,N));
+C If M = 0, LDWORK >= 1.
+C For optimum performance LDWORK should sometimes be larger.
+C
+C Error Indicator
+C
+C INFO INTEGER
+C = 0: successful exit;
+C < 0: if INFO = -i, the i-th argument had an illegal
+C value;
+C = 1: if the Lyapunov equation is (nearly) singular
+C (warning indicator);
+C if DICO = 'C' this means that while the matrix A
+C (or the factor S) has computed eigenvalues with
+C negative real parts, it is only just stable in the
+C sense that small perturbations in A can make one or
+C more of the eigenvalues have a non-negative real
+C part;
+C if DICO = 'D' this means that while the matrix A
+C (or the factor S) has computed eigenvalues inside
+C the unit circle, it is nevertheless only just
+C convergent, in the sense that small perturbations
+C in A can make one or more of the eigenvalues lie
+C outside the unit circle;
+C perturbed values were used to solve the equation;
+C = 2: if FACT = 'N' and DICO = 'C', but the matrix A is
+C not stable (that is, one or more of the eigenvalues
+C of A has a non-negative real part), or DICO = 'D',
+C but the matrix A is not convergent (that is, one or
+C more of the eigenvalues of A lies outside the unit
+C circle); however, A will still have been factored
+C and the eigenvalues of A returned in WR and WI.
+C = 3: if FACT = 'F' and DICO = 'C', but the Schur factor S
+C supplied in the array A is not stable (that is, one
+C or more of the eigenvalues of S has a non-negative
+C real part), or DICO = 'D', but the Schur factor S
+C supplied in the array A is not convergent (that is,
+C one or more of the eigenvalues of S lies outside the
+C unit circle);
+C = 4: if FACT = 'F' and the Schur factor S supplied in
+C the array A has two or more consecutive non-zero
+C elements on the first sub-diagonal, so that there is
+C a block larger than 2-by-2 on the diagonal;
+C = 5: if FACT = 'F' and the Schur factor S supplied in
+C the array A has a 2-by-2 diagonal block with real
+C eigenvalues instead of a complex conjugate pair;
+C = 6: if FACT = 'N' and the LAPACK Library routine DGEES
+C has failed to converge. This failure is not likely
+C to occur. The matrix B will be unaltered but A will
+C be destroyed.
+C
+C METHOD
+C
+C The method used by the routine is based on the Bartels and Stewart
+C method [1], except that it finds the upper triangular matrix U
+C directly without first finding X and without the need to form the
+C normal matrix op(B)'*op(B).
+C
+C The Schur factorization of a square matrix A is given by
+C
+C A = QSQ',
+C
+C where Q is orthogonal and S is an N-by-N block upper triangular
+C matrix with 1-by-1 and 2-by-2 blocks on its diagonal (which
+C correspond to the eigenvalues of A). If A has already been
+C factored prior to calling the routine however, then the factors
+C Q and S may be supplied and the initial factorization omitted.
+C
+C If TRANS = 'N', the matrix B is factored as (QR factorization)
+C _ _ _ _ _
+C B = P ( R ), M >= N, B = P ( R Z ), M < N,
+C ( 0 )
+C _ _
+C where P is an M-by-M orthogonal matrix and R is a square upper
+C _ _ _ _ _
+C triangular matrix. Then, the matrix B = RQ, or B = ( R Z )Q (if
+C M < N) is factored as
+C _ _
+C B = P ( R ), M >= N, B = P ( R Z ), M < N.
+C
+C If TRANS = 'T', the matrix B is factored as (RQ factorization)
+C _
+C _ _ ( Z ) _
+C B = ( 0 R ) P, M >= N, B = ( _ ) P, M < N,
+C ( R )
+C _ _
+C where P is an M-by-M orthogonal matrix and R is a square upper
+C _ _ _ _ _
+C triangular matrix. Then, the matrix B = Q'R, or B = Q'( Z' R' )'
+C (if M < N) is factored as
+C _ _
+C B = ( R ) P, M >= N, B = ( Z ) P, M < N.
+C ( R )
+C
+C These factorizations are utilised to either transform the
+C continuous-time Lyapunov equation to the canonical form
+C 2
+C op(S)'*op(V)'*op(V) + op(V)'*op(V)*op(S) = -scale *op(F)'*op(F),
+C
+C or the discrete-time Lyapunov equation to the canonical form
+C 2
+C op(S)'*op(V)'*op(V)*op(S) - op(V)'*op(V) = -scale *op(F)'*op(F),
+C
+C where V and F are upper triangular, and
+C
+C F = R, M >= N, F = ( R Z ), M < N, if TRANS = 'N';
+C ( 0 0 )
+C
+C F = R, M >= N, F = ( 0 Z ), M < N, if TRANS = 'T'.
+C ( 0 R )
+C
+C The transformed equation is then solved for V, from which U is
+C obtained via the QR factorization of V*Q', if TRANS = 'N', or
+C via the RQ factorization of Q*V, if TRANS = 'T'.
+C
+C REFERENCES
+C
+C [1] Bartels, R.H. and Stewart, G.W.
+C Solution of the matrix equation A'X + XB = C.
+C Comm. A.C.M., 15, pp. 820-826, 1972.
+C
+C [2] Hammarling, S.J.
+C Numerical solution of the stable, non-negative definite
+C Lyapunov equation.
+C IMA J. Num. Anal., 2, pp. 303-325, 1982.
+C
+C NUMERICAL ASPECTS
+C 3
+C The algorithm requires 0(N ) operations and is backward stable.
+C
+C FURTHER COMMENTS
+C
+C The Lyapunov equation may be very ill-conditioned. In particular,
+C if A is only just stable (or convergent) then the Lyapunov
+C equation will be ill-conditioned. A symptom of ill-conditioning
+C is "large" elements in U relative to those of A and B, or a
+C "small" value for scale. A condition estimate can be computed
+C using SLICOT Library routine SB03MD.
+C
+C SB03OD routine can be also used for solving "unstable" Lyapunov
+C equations, i.e., when matrix A has all eigenvalues with positive
+C real parts, if DICO = 'C', or with moduli greater than one,
+C if DICO = 'D'. Specifically, one may solve for X = op(U)'*op(U)
+C either the continuous-time Lyapunov equation
+C 2
+C op(A)'*X + X*op(A) = scale *op(B)'*op(B), (3)
+C
+C or the discrete-time Lyapunov equation
+C 2
+C op(A)'*X*op(A) - X = scale *op(B)'*op(B), (4)
+C
+C provided, for equation (3), the given matrix A is replaced by -A,
+C or, for equation (4), the given matrices A and B are replaced by
+C inv(A) and B*inv(A), if TRANS = 'N' (or inv(A)*B, if TRANS = 'T'),
+C respectively. Although the inversion generally can rise numerical
+C problems, in case of equation (4) it is expected that the matrix A
+C is enough well-conditioned, having only eigenvalues with moduli
+C greater than 1. However, if A is ill-conditioned, it could be
+C preferable to use the more general SLICOT Lyapunov solver SB03MD.
+C
+C CONTRIBUTOR
+C
+C Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, Aug. 1997.
+C Supersedes Release 2.0 routine SB03CD by Sven Hammarling,
+C NAG Ltd, United Kingdom.
+C
+C REVISIONS
+C
+C Dec. 1997, April 1998, May 1998, May 1999, Oct. 2001 (V. Sima).
+C March 2002 (A. Varga).
+C
+C KEYWORDS
+C
+C Lyapunov equation, orthogonal transformation, real Schur form,
+C Sylvester equation.
+C
+C ******************************************************************
+C
+C .. Parameters ..
+ DOUBLE PRECISION ZERO, ONE
+ PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
+C .. Scalar Arguments ..
+ CHARACTER DICO, FACT, TRANS
+ INTEGER INFO, LDA, LDB, LDQ, LDWORK, M, N
+ DOUBLE PRECISION SCALE
+C .. Array Arguments ..
+ DOUBLE PRECISION A(LDA,*), B(LDB,*), DWORK(*), Q(LDQ,*), WI(*),
+ $ WR(*)
+C .. Local Scalars ..
+ LOGICAL CONT, LTRANS, NOFACT
+ INTEGER I, IFAIL, INFORM, ITAU, J, JWORK, K, L, MINMN,
+ $ NE, SDIM, WRKOPT
+ DOUBLE PRECISION EMAX, TEMP
+C .. Local Arrays ..
+ LOGICAL BWORK(1)
+C .. External Functions ..
+ LOGICAL LSAME, SELECT
+ DOUBLE PRECISION DLAPY2
+ EXTERNAL DLAPY2, LSAME, SELECT
+C .. External Subroutines ..
+ EXTERNAL DCOPY, DGEES, DGEMM, DGEMV, DGEQRF, DGERQF,
+ $ DLACPY, DLASET, DTRMM, SB03OU, XERBLA
+C .. Intrinsic Functions ..
+ INTRINSIC INT, MAX, MIN
+C .. Executable Statements ..
+C
+C Test the input scalar arguments.
+C
+ CONT = LSAME( DICO, 'C' )
+ NOFACT = LSAME( FACT, 'N' )
+ LTRANS = LSAME( TRANS, 'T' )
+ MINMN = MIN( M, N )
+C
+ INFO = 0
+ IF( .NOT.CONT .AND. .NOT.LSAME( DICO, 'D' ) ) THEN
+ INFO = -1
+ ELSE IF( .NOT.NOFACT .AND. .NOT.LSAME( FACT, 'F' ) ) THEN
+ INFO = -2
+ ELSE IF( .NOT.LTRANS .AND. .NOT.LSAME( TRANS, 'N' ) ) THEN
+ INFO = -3
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -4
+ ELSE IF( M.LT.0 ) THEN
+ INFO = -5
+ ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+ INFO = -7
+ ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
+ INFO = -9
+ ELSE IF( ( LDB.LT.MAX( 1, N ) ) .OR.
+ $ ( LDB.LT.MAX( 1, N, M ) .AND. .NOT.LTRANS ) ) THEN
+ INFO = -11
+ ELSE IF( LDWORK.LT.1 .OR. ( M.GT.0 .AND. LDWORK.LT.4*N + MINMN ) )
+ $ THEN
+ INFO = -16
+ END IF
+C
+ IF ( INFO.NE.0 ) THEN
+C
+C Error return.
+C
+ CALL XERBLA( 'SB03OD', -INFO )
+ RETURN
+ END IF
+C
+C Quick return if possible.
+C
+ IF( MINMN.EQ.0 ) THEN
+ IF( M.EQ.0 )
+ $ CALL DLASET( 'Full', N, N, ZERO, ZERO, B, LDB )
+ SCALE = ONE
+ DWORK(1) = ONE
+ RETURN
+ END IF
+C
+C Start the solution.
+C
+C (Note: Comments in the code beginning "Workspace:" describe the
+C minimal amount of real workspace needed at that point in the
+C code, as well as the preferred amount for good performance.
+C NB refers to the optimal block size for the immediately
+C following subroutine, as returned by ILAENV.)
+C
+ IF ( NOFACT ) THEN
+C
+C Find the Schur factorization of A, A = Q*S*Q'.
+C Workspace: need 3*N;
+C prefer larger.
+C
+ CALL DGEES( 'Vectors', 'Not ordered', SELECT, N, A, LDA, SDIM,
+ $ WR, WI, Q, LDQ, DWORK, LDWORK, BWORK, INFORM )
+ IF ( INFORM.NE.0 ) THEN
+ INFO = 6
+ RETURN
+ END IF
+ WRKOPT = DWORK(1)
+C
+C Check the eigenvalues for stability.
+C
+ IF ( CONT ) THEN
+ EMAX = WR(1)
+C
+ DO 20 J = 2, N
+ IF ( WR(J).GT.EMAX )
+ $ EMAX = WR(J)
+ 20 CONTINUE
+C
+ ELSE
+ EMAX = DLAPY2( WR(1), WI(1) )
+C
+ DO 40 J = 2, N
+ TEMP = DLAPY2( WR(J), WI(J) )
+ IF ( TEMP.GT.EMAX )
+ $ EMAX = TEMP
+ 40 CONTINUE
+C
+ END IF
+C
+ IF ( ( CONT ) .AND. ( EMAX.GE.ZERO ) .OR.
+ $ ( .NOT.CONT ) .AND. ( EMAX.GE.ONE ) ) THEN
+ INFO = 2
+ RETURN
+ END IF
+ ELSE
+ WRKOPT = 0
+ END IF
+C
+C Perform the QR or RQ factorization of B,
+C _ _ _ _ _
+C B = P ( R ), or B = P ( R Z ), if TRANS = 'N', or
+C ( 0 )
+C _
+C _ _ ( Z ) _
+C B = ( 0 R ) P, or B = ( _ ) P, if TRANS = 'T'.
+C ( R )
+C Workspace: need MIN(M,N) + N;
+C prefer MIN(M,N) + N*NB.
+C
+ ITAU = 1
+ JWORK = ITAU + MINMN
+ IF ( LTRANS ) THEN
+ CALL DGERQF( N, M, B, LDB, DWORK(ITAU), DWORK(JWORK),
+ $ LDWORK-JWORK+1, IFAIL )
+ WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) ) + JWORK - 1, MINMN*N)
+ JWORK = ITAU
+C
+C Form in B
+C _ _ _ _ _ _
+C B := Q'R, m >= n, B := Q'*( Z' R' )', m < n, with B an
+C n-by-min(m,n) matrix.
+C Use a BLAS 3 operation if enough workspace, and BLAS 2,
+C _
+C otherwise: B is formed column by column.
+C
+ IF ( LDWORK.GE.JWORK+MINMN*N-1 ) THEN
+ K = JWORK
+C
+ DO 60 I = 1, MINMN
+ CALL DCOPY( N, Q(N-MINMN+I,1), LDQ, DWORK(K), 1 )
+ K = K + N
+ 60 CONTINUE
+C
+ CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Non-unit',
+ $ N, MINMN, ONE, B(N-MINMN+1,M-MINMN+1), LDB,
+ $ DWORK(JWORK), N )
+ IF ( M.LT.N )
+ $ CALL DGEMM( 'Transpose', 'No transpose', N, M, N-M,
+ $ ONE, Q, LDQ, B, LDB, ONE, DWORK(JWORK), N )
+ CALL DLACPY( 'Full', N, MINMN, DWORK(JWORK), N, B, LDB )
+ ELSE
+ NE = N - MINMN
+C
+ DO 80 J = 1, MINMN
+ NE = NE + 1
+ CALL DCOPY( NE, B(1,M-MINMN+J), 1, DWORK(JWORK), 1 )
+ CALL DGEMV( 'Transpose', NE, N, ONE, Q, LDQ,
+ $ DWORK(JWORK), 1, ZERO, B(1,J), 1 )
+ 80 CONTINUE
+C
+ END IF
+ ELSE
+ CALL DGEQRF( M, N, B, LDB, DWORK(ITAU), DWORK(JWORK),
+ $ LDWORK-JWORK+1, IFAIL )
+ WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) ) + JWORK - 1, MINMN*N)
+ JWORK = ITAU
+C
+C Form in B
+C _ _ _ _ _ _
+C B := RQ, m >= n, B := ( R Z )*Q, m < n, with B an
+C min(m,n)-by-n matrix.
+C Use a BLAS 3 operation if enough workspace, and BLAS 2,
+C _
+C otherwise: B is formed row by row.
+C
+ IF ( LDWORK.GE.JWORK+MINMN*N-1 ) THEN
+ CALL DLACPY( 'Full', MINMN, N, Q, LDQ, DWORK(JWORK), MINMN )
+ CALL DTRMM( 'Left', 'Upper', 'No transpose', 'Non-unit',
+ $ MINMN, N, ONE, B, LDB, DWORK(JWORK), MINMN )
+ IF ( M.LT.N )
+ $ CALL DGEMM( 'No transpose', 'No transpose', M, N, N-M,
+ $ ONE, B(1,M+1), LDB, Q(M+1,1), LDQ, ONE,
+ $ DWORK(JWORK), MINMN )
+ CALL DLACPY( 'Full', MINMN, N, DWORK(JWORK), MINMN, B, LDB )
+ ELSE
+ NE = MINMN + MAX( 0, N-M )
+C
+ DO 100 J = 1, MINMN
+ CALL DCOPY( NE, B(J,J), LDB, DWORK(JWORK), 1 )
+ CALL DGEMV( 'Transpose', NE, N, ONE, Q(J,1), LDQ,
+ $ DWORK(JWORK), 1, ZERO, B(J,1), LDB )
+ NE = NE - 1
+ 100 CONTINUE
+C
+ END IF
+ END IF
+ JWORK = ITAU + MINMN
+C
+C Solve for U the transformed Lyapunov equation
+C 2 _ _
+C op(S)'*op(U)'*op(U) + op(U)'*op(U)*op(S) = -scale *op(B)'*op(B),
+C
+C or
+C 2 _ _
+C op(S)'*op(U)'*op(U)*op(S) - op(U)'*op(U) = -scale *op(B)'*op(B)
+C
+C Workspace: need MIN(M,N) + 4*N;
+C prefer larger.
+C
+ CALL SB03OU( .NOT.CONT, LTRANS, N, MINMN, A, LDA, B, LDB,
+ $ DWORK(ITAU), B, LDB, SCALE, DWORK(JWORK),
+ $ LDWORK-JWORK+1, INFO )
+ IF ( INFO.GT.1 ) THEN
+ INFO = INFO + 1
+ RETURN
+ END IF
+ WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) ) + JWORK - 1 )
+ JWORK = ITAU
+C
+C Form U := U*Q' or U := Q*U in the array B.
+C Use a BLAS 3 operation if enough workspace, and BLAS 2, otherwise.
+C Workspace: need N;
+C prefer N*N;
+C
+ IF ( LDWORK.GE.JWORK+N*N-1 ) THEN
+ IF ( LTRANS ) THEN
+ CALL DLACPY( 'Full', N, N, Q, LDQ, DWORK(JWORK), N )
+ CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Non-unit', N,
+ $ N, ONE, B, LDB, DWORK(JWORK), N )
+ ELSE
+ K = JWORK
+C
+ DO 120 I = 1, N
+ CALL DCOPY( N, Q(1,I), 1, DWORK(K), N )
+ K = K + 1
+ 120 CONTINUE
+C
+ CALL DTRMM( 'Left', 'Upper', 'No transpose', 'Non-unit', N,
+ $ N, ONE, B, LDB, DWORK(JWORK), N )
+ END IF
+ CALL DLACPY( 'Full', N, N, DWORK(JWORK), N, B, LDB )
+ WRKOPT = MAX( WRKOPT, JWORK + N*N - 1 )
+ ELSE
+ IF ( LTRANS ) THEN
+C
+C U is formed column by column ( U := Q*U ).
+C
+ DO 140 I = 1, N
+ CALL DCOPY( I, B(1,I), 1, DWORK(JWORK), 1 )
+ CALL DGEMV( 'No transpose', N, I, ONE, Q, LDQ,
+ $ DWORK(JWORK), 1, ZERO, B(1,I), 1 )
+ 140 CONTINUE
+ ELSE
+C
+C U is formed row by row ( U' := Q*U' ).
+C
+ DO 160 I = 1, N
+ CALL DCOPY( N-I+1, B(I,I), LDB, DWORK(JWORK), 1 )
+ CALL DGEMV( 'No transpose', N, N-I+1, ONE, Q(1,I), LDQ,
+ $ DWORK(JWORK), 1, ZERO, B(I,1), LDB )
+ 160 CONTINUE
+ END IF
+ END IF
+C
+C Lastly find the QR or RQ factorization of U, overwriting on B,
+C to give the required Cholesky factor.
+C Workspace: need 2*N;
+C prefer N + N*NB;
+C
+ JWORK = ITAU + N
+ IF ( LTRANS ) THEN
+ CALL DGERQF( N, N, B, LDB, DWORK(ITAU), DWORK(JWORK),
+ $ LDWORK-JWORK+1, IFAIL )
+ ELSE
+ CALL DGEQRF( N, N, B, LDB, DWORK(ITAU), DWORK(JWORK),
+ $ LDWORK-JWORK+1, IFAIL )
+ END IF
+ WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) ) + JWORK - 1 )
+C
+C Make the diagonal elements of U non-negative.
+C
+ IF ( LTRANS ) THEN
+C
+ DO 200 J = 1, N
+ IF ( B(J,J).LT.ZERO ) THEN
+C
+ DO 180 I = 1, J
+ B(I,J) = -B(I,J)
+ 180 CONTINUE
+C
+ END IF
+ 200 CONTINUE
+C
+ ELSE
+ K = JWORK
+C
+ DO 240 J = 1, N
+ DWORK(K) = B(J,J)
+ L = JWORK
+C
+ DO 220 I = 1, J
+ IF ( DWORK(L).LT.ZERO ) B(I,J) = -B(I,J)
+ L = L + 1
+ 220 CONTINUE
+C
+ K = K + 1
+ 240 CONTINUE
+ END IF
+C
+ IF( N.GT.1 )
+ $ CALL DLASET( 'Lower', N-1, N-1, ZERO, ZERO, B(2,1), LDB )
+C
+C Set the optimal workspace.
+C
+ DWORK(1) = WRKOPT
+C
+ RETURN
+C *** Last line of SB03OD ***
+ END
Added: trunk/octave-forge/main/control/src/SG03BD.f
===================================================================
--- trunk/octave-forge/main/control/src/SG03BD.f (rev 0)
+++ trunk/octave-forge/main/control/src/SG03BD.f 2010-09-14 10:59:36 UTC (rev 7713)
@@ -0,0 +1,814 @@
+ SUBROUTINE SG03BD( DICO, FACT, TRANS, N, M, A, LDA, E, LDE, Q,
+ $ LDQ, Z, LDZ, B, LDB, SCALE, ALPHAR, ALPHAI,
+ $ BETA, DWORK, LDWORK, INFO )
+C
+C SLICOT RELEASE 5.0.
+C
+C Copyright (c) 2002-2009 NICONET e.V.
+C
+C This program is free software: you can redistribute it and/or
+C modify it under the terms of the GNU General Public License as
+C published by the Free Software Foundation, either version 2 of
+C the License, or (at your option) any later version.
+C
+C This program is distributed in the hope that it will be useful,
+C but WITHOUT ANY WARRANTY; without even the implied warranty of
+C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+C GNU General Public License for more details.
+C
+C You should have received a copy of the GNU General Public License
+C along with this program. If not, see
+C <http://www.gnu.org/licenses/>.
+C
+C PURPOSE
+C
+C To compute the Cholesky factor U of the matrix X,
+C
+C T
+C X = op(U) * op(U),
+C
+C which is the solution of either the generalized
+C c-stable continuous-time Lyapunov equation
+C
+C T T
+C op(A) * X * op(E) + op(E) * X * op(A)
+C
+C 2 T
+C = - SCALE * op(B) * op(B), (1)
+C
+C or the generalized d-stable discrete-time Lyapunov equation
+C
+C T T
+C op(A) * X * op(A) - op(E) * X * op(E)
+C
+C 2 T
+C = - SCALE * op(B) * op(B), (2)
+C
+C without first finding X and without the need to form the matrix
+C op(B)**T * op(B).
+C
+C op(K) is either K or K**T for K = A, B, E, U. A and E are N-by-N
+C matrices, op(B) is an M-by-N matrix. The resulting matrix U is an
+C N-by-N upper triangular matrix with non-negative entries on its
+C main diagonal. SCALE is an output scale factor set to avoid
+C overflow in U.
+C
+C In the continuous-time case (1) the pencil A - lambda * E must be
+C c-stable (that is, all eigenvalues must have negative real parts).
+C In the discrete-time case (2) the pencil A - lambda * E must be
+C d-stable (that is, the moduli of all eigenvalues must be smaller
+C than one).
+C
+C ARGUMENTS
+C
+C Mode Parameters
+C
+C DICO CHARACTER*1
+C Specifies which type of the equation is considered:
+C = 'C': Continuous-time equation (1);
+C = 'D': Discrete-time equation (2).
+C
+C FACT CHARACTER*1
+C Specifies whether the generalized real Schur
+C factorization of the pencil A - lambda * E is supplied
+C on entry or not:
+C = 'N': Factorization is not supplied;
+C = 'F': Factorization is supplied.
+C
+C TRANS CHARACTER*1
+C Specifies whether the transposed equation is to be solved
+C or not:
+C = 'N': op(A) = A, op(E) = E;
+C = 'T': op(A) = A**T, op(E) = E**T.
+C
+C Input/Output Parameters
+C
+C N (input) INTEGER
+C The order of the matrix A. N >= 0.
+C
+C M (input) INTEGER
+C The number of rows in the matrix op(B). M >= 0.
+C
+C A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
+C On entry, if FACT = 'F', then the leading N-by-N upper
+C Hessenberg part of this array must contain the
+C generalized Schur factor A_s of the matrix A (see
+C definition (3) in section METHOD). A_s must be an upper
+C quasitriangular matrix. The elements below the upper
+C Hessenberg part of the array A are not referenced.
+C If FACT = 'N', then the leading N-by-N part of this
+C array must contain the matrix A.
+C On exit, the leading N-by-N part of this array contains
+C the generalized Schur factor A_s of the matrix A. (A_s is
+C an upper quasitriangular matrix.)
+C
+C LDA INTEGER
+C The leading dimension of the array A. LDA >= MAX(1,N).
+C
+C E (input/output) DOUBLE PRECISION array, dimension (LDE,N)
+C On entry, if FACT = 'F', then the leading N-by-N upper
+C triangular part of this array must contain the
+C generalized Schur factor E_s of the matrix E (see
+C definition (4) in section METHOD). The elements below the
+C upper triangular part of the array E are not referenced.
+C If FACT = 'N', then the leading N-by-N part of this
+C array must contain the coefficient matrix E of the
+C equation.
+C On exit, the leading N-by-N part of this array contains
+C the generalized Schur factor E_s of the matrix E. (E_s is
+C an upper triangular matrix.)
+C
+C LDE INTEGER
+C The leading dimension of the array E. LDE >= MAX(1,N).
+C
+C Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
+C On entry, if FACT = 'F', then the leading N-by-N part of
+C this array must contain the orthogonal matrix Q from
+C the generalized Schur factorization (see definitions (3)
+C and (4) in section METHOD).
+C If FACT = 'N', Q need not be set on entry.
+C On exit, the leading N-by-N part of this array contains
+C the orthogonal matrix Q from the generalized Schur
+C factorization.
+C
+C LDQ INTEGER
+C The leading dimension of the array Q. LDQ >= MAX(1,N).
+C
+C Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
+C On entry, if FACT = 'F', then the leading N-by-N part of
+C this array must contain the orthogonal matrix Z from
+C the generalized Schur factorization (see definitions (3)
+C and (4) in section METHOD).
+C If FACT = 'N', Z need not be set on entry.
+C On exit, the leading N-by-N part of this array contains
+C the orthogonal matrix Z from the generalized Schur
+C factorization.
+C
+C LDZ INTEGER
+C The leading dimension of the array Z. LDZ >= MAX(1,N).
+C
+C B (input/output) DOUBLE PRECISION array, dimension (LDB,N1)
+C On entry, if TRANS = 'T', the leading N-by-M part of this
+C array must contain the matrix B and N1 >= MAX(M,N).
+C If TRANS = 'N', the leading M-by-N part of this array
+C must contain the matrix B and N1 >= N.
+C On exit, the leading N-by-N part of this array contains
+C the Cholesky factor U of the solution matrix X of the
+C problem, X = op(U)**T * op(U).
+C If M = 0 and N > 0, then U is set to zero.
+C
+C LDB INTEGER
+C The leading dimension of the array B.
+C If TRANS = 'T', LDB >= MAX(1,N).
+C If TRANS = 'N', LDB >= MAX(1,M,N).
+C
+C SCALE (output) DOUBLE PRECISION
+C The scale factor set to avoid overflow in U.
+C 0 < SCALE <= 1.
+C
+C ALPHAR (output) DOUBLE PRECISION array, dimension (N)
+C ALPHAI (output) DOUBLE PRECISION array, dimension (N)
+C BETA (output) DOUBLE PRECISION array, dimension (N)
+C If INFO = 0, 3, 5, 6, or 7, then
+C (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, are the
+C eigenvalues of the matrix pencil A - lambda * E.
+C
+C Workspace
+C
+C DWORK DOUBLE PRECISION array, dimension (LDWORK)
+C On exit, if INFO = 0, DWORK(1) returns the optimal value
+C of LDWORK.
+C
+C LDWORK INTEGER
+C The dimension of the array DWORK.
+C LDWORK >= MAX(1,4*N,6*N-6), if FACT = 'N';
+C LDWORK >= MAX(1,2*N,6*N-6), if FACT = 'F'.
+C For good performance, LDWORK should be larger.
+C
+C Error indicator
+C
+C INFO INTEGER
+C = 0: successful exit;
+C < 0: if INFO = -i, the i-th argument had an illegal
+C value;
+C = 1: the pencil A - lambda * E is (nearly) singular;
+C perturbed values were used to solve the equation
+C (but the reduced (quasi)triangular matrices A and E
+C are unchanged);
+C = 2: FACT = 'F' and the matrix contained in the upper
+C Hessenberg part of the array A is not in upper
+C quasitriangular form;
+C = 3: FACT = 'F' and there is a 2-by-2 block on the main
+C diagonal of the pencil A_s - lambda * E_s whose
+C eigenvalues are not conjugate complex;
+C = 4: FACT = 'N' and the pencil A - lambda * E cannot be
+C reduced to generalized Schur form: LAPACK routine
+C DGEGS has failed to converge;
+C = 5: DICO = 'C' and the pencil A - lambda * E is not
+C c-stable;
+C = 6: DICO = 'D' and the pencil A - lambda * E is not
+C d-stable;
+C = 7: the LAPACK routine DSYEVX utilized to factorize M3
+C failed to converge in the discrete-time case (see
+C section METHOD for SLICOT Library routine SG03BU).
+C This error is unlikely to occur.
+C
+C METHOD
+C
+C An extension [2] of Hammarling's method [1] to generalized
+C Lyapunov equations is utilized to solve (1) or (2).
+C
+C First the pencil A - lambda * E is reduced to real generalized
+C Schur form A_s - lambda * E_s by means of orthogonal
+C transformations (QZ-algorithm):
+C
+C A_s = Q**T * A * Z (upper quasitriangular) (3)
+C
+C E_s = Q**T * E * Z (upper triangular). (4)
+C
+C If the pencil A - lambda * E has already been factorized prior to
+C calling the routine however, then the factors A_s, E_s, Q and Z
+C may be supplied and the initial factorization omitted.
+C
+C Depending on the parameters TRANS and M the N-by-N upper
+C triangular matrix B_s is defined as follows. In any case Q_B is
+C an M-by-M orthogonal matrix, which need not be accumulated.
+C
+C 1. If TRANS = 'N' and M < N, B_s is the upper triangular matrix
+C from the QR-factorization
+C
+C ( Q_B O ) ( B * Z )
+C ( ) * B_s = ( ),
+C ( O I ) ( O )
+C
+C where the O's are zero matrices of proper size and I is the
+C identity matrix of order N-M.
+C
+C 2. If TRANS = 'N' and M >= N, B_s is the upper triangular matrix
+C from the (rectangular) QR-factorization
+C
+C ( B_s )
+C Q_B * ( ) = B * Z,
+C ( O )
+C
+C where O is the (M-N)-by-N zero matrix.
+C
+C 3. If TRANS = 'T' and M < N, B_s is the upper triangular matrix
+C from the RQ-factorization
+C
+C ( Q_B O )
+C (B_s O ) * ( ) = ( Q**T * B O ).
+C ( O I )
+C
+C 4. If TRANS = 'T' and M >= N, B_s is the upper triangular matrix
+C from the (rectangular) RQ-factorization
+C
+C ( B_s O ) * Q_B = Q**T * B,
+C
+C where O is the N-by-(M-N) zero matrix.
+C
+C Assuming SCALE = 1, the transformation of A, E and B described
+C above leads to the reduced continuous-time equation
+C
+C T T
+C op(A_s) op(U_s) op(U_s) op(E_s)
+C
+C T T
+C + op(E_s) op(U_s) op(U_s) op(A_s)
+C
+C T
+C = - op(B_s) op(B_s) (5)
+C
+C or to the reduced discrete-time equation
+C
+C T T
+C op(A_s) op(U_s) op(U_s) op(A_s)
+C
+C T T
+C - op(E_s) op(U_s) op(U_s) op(E_s)
+C
+C T
+C = - op(B_s) op(B_s). (6)
+C
+C For brevity we restrict ourself to equation (5) and the case
+C TRANS = 'N'. The other three cases can be treated in a similar
+C fashion.
+C
+C We use the following partitioning for the matrices A_s, E_s, B_s
+C and U_s
+C
+C ( A11 A12 ) ( E11 E12 )
+C A_s = ( ), E_s = ( ),
+C ( 0 A22 ) ( 0 E22 )
+C
+C ( B11 B12 ) ( U11 U12 )
+C B_s = ( ), U_s = ( ). (7)
+C ( 0 B22 ) ( 0 U22 )
+C
+C The size of the (1,1)-blocks is 1-by-1 (iff A_s(2,1) = 0.0) or
+C 2-by-2.
+C
+C We compute U11 and U12**T in three steps.
+C
+C Step I:
+C
+C From (5) and (7) we get the 1-by-1 or 2-by-2 equation
+C
+C T T T T
+C A11 * U11 * U11 * E11 + E11 * U11 * U11 * A11
+C
+C T
+C = - B11 * B11.
+C
+C For brevity, details are omitted here. See [2]. The technique
+C for computing U11 is similar to those applied to standard
+C Lyapunov equations in Hammarling's algorithm ([1], section 6).
+C
+C Furthermore, the auxiliary matrices M1 and M2 defined as
+C follows
+C
+C -1 -1
+C M1 = U11 * A11 * E11 * U11
+C
+C -1 -1
+C M2 = B11 * E11 * U11
+C
+C are computed in a numerically reliable way.
+C
+C Step II:
+C
+C The generalized Sylvester equation
+C
+C T T T T
+C A22 * U12 + E22 * U12 * M1 =
+C
+C T T T T T
+C - B12 * M2 - A12 * U11 - E12 * U11 * M1
+C
+C is solved for U12**T.
+C
+C Step III:
+C
+C It can be shown that
+C
+C T T T T
+C A22 * U22 * U22 * E22 + E22 * U22 * U22 * A22 =
+C
+C T T
+C - B22 * B22 - y * y (8)
+C
+C holds, where y is defined as
+C
+C T T T T T T
+C y = B12 - ( E12 * U11 + E22 * U12 ) * M2 .
+C
+C If B22_tilde is the square triangular matrix arising from the
+C (rectangular) QR-factorization
+C
+C ( B22_tilde ) ( B22 )
+C Q_B_tilde * ( ) = ( ),
+C ( O ) ( y**T )
+C
+C where Q_B_tilde is an orthogonal matrix of order N, then
+C
+C T T T
+C - B22 * B22 - y * y = - B22_tilde * B22_tilde.
+C
+C Replacing the right hand side in (8) by the term
+C - B22_tilde**T * B22_tilde leads to a reduced generalized
+C Lyapunov equation of lower dimension compared to (5).
+C
+C The recursive application of the steps I to III yields the
+C solution U_s of the equation (5).
+C
+C It remains to compute the solution matrix U of the original
+C problem (1) or (2) from the matrix U_s. To this end we transform
+C the solution back (with respect to the transformation that led
+C from (1) to (5) (from (2) to (6)) and apply the QR-factorization
+C (RQ-factorization). The upper triangular solution matrix U is
+C obtained by
+C
+C Q_U * U = U_s * Q**T (if TRANS = 'N')
+C
+C or
+C
+C U * Q_U = Z * U_s (if TRANS = 'T')
+C
+C where Q_U is an N-by-N orthogonal matrix. Again, the orthogonal
+C matrix Q_U need not be accumulated.
+C
+C REFERENCES
+C
+C [1] Hammarling, S.J.
+C Numerical solution of the stable, non-negative definite
+C Lyapunov equation.
+C IMA J. Num. Anal., 2, pp. 303-323, 1982.
+C
+C [2] Penzl, T.
+C Numerical solution of generalized Lyapunov equations.
+C Advances in Comp. Math., vol. 8, pp. 33-48, 1998.
+C
+C NUMERICAL ASPECTS
+C
+C The number of flops required by the routine is given by the
+C following table. Note that we count a single floating point
+C arithmetic operation as one flop.
+C
+C | FACT = 'F' FACT = 'N'
+C ---------+--------------------------------------------------
+C M <= N | (13*N**3+6*M*N**2 (211*N**3+6*M*N**2
+C | +6*M**2*N-2*M**3)/3 +6*M**2*N-2*M**3)/3
+C |
+C M > N | (11*N**3+12*M*N**2)/3 (209*N**3+12*M*N**2)/3
+C
+C FURTHER COMMENTS
+C
+C The Lyapunov equation may be very ill-conditioned. In particular,
+C if DICO = 'D' and the pencil A - lambda * E has a pair of almost
+C reciprocal eigenvalues, or DICO = 'C' and the pencil has an almost
+C degenerate pair of eigenvalues, then the Lyapunov equation will be
+C ill-conditioned. Perturbed values were used to solve the equation.
+C A condition estimate can be obtained from the routine SG03AD.
+C When setting the error indicator INFO, the routine does not test
+C for near instability in the equation but only for exact
+C instability.
+C
+C CONTRIBUTOR
+C
+C T. Penzl, Technical University Chemnitz, Germany, Aug. 1998.
+C
+C REVISIONS
+C
+C Sep. 1998 (V. Sima).
+C May 1999 (V. Sima).
+C March 2002 (A. Varga).
+C Feb. 2004 (V. Sima).
+C
+C KEYWORDS
+C
+C Lyapunov equation
+C
+C ******************************************************************
+C
+C .. Parameters ..
+ DOUBLE PRECISION MONE, ONE, TWO, ZERO
+ PARAMETER ( MONE = -1.0D+0, ONE = 1.0D+0, TWO = 2.0D+0,
+ $ ZERO = 0.0D+0 )
+C .. Scalar Arguments ..
+ DOUBLE PRECISION SCALE
+ INTEGER INFO, LDA, LDB, LDE, LDQ, LDWORK, LDZ, M, N
+ CHARACTER DICO, FACT, TRANS
+C .. Array Arguments ..
+ DOUBLE PRECISION A(LDA,*), ALPHAI(*), ALPHAR(*), B(LDB,*),
+ $ BETA(*), DWORK(*), E(LDE,*), Q(LDQ,*), Z(LDZ,*)
+C .. Local Scalars ..
+ DOUBLE PRECISION S1, S2, SAFMIN, WI, WR1, WR2
+ INTEGER I, INFO1, MINMN, MINWRK, OPTWRK
+ LOGICAL ISDISC, ISFACT, ISTRAN
+C .. Local Arrays ..
+ DOUBLE PRECISION E1(2,2)
+C .. External Functions ..
+ DOUBLE PRECISION DLAMCH, DLAPY2
+ LOGICAL LSAME
+ EXTERNAL DLAMCH, DLAPY2, LSAME
+C .. External Subroutines ..
+ EXTERNAL DCOPY, DGEGS, DGEMM, DGEMV, DGEQRF, DGERQF,
+ $ DLACPY, DLAG2, DLASET, DSCAL, DTRMM, SG03BU,
+ $ SG03BV, XERBLA
+C .. Intrinsic Functions ..
+ INTRINSIC ABS, DBLE, INT, MAX, MIN, SIGN
+C .. Executable Statements ..
+C
+C Decode input parameters.
+C
+ ISDISC = LSAME( DICO, 'D' )
+ ISFACT = LSAME( FACT, 'F' )
+ ISTRAN = LSAME( TRANS, 'T' )
+C
+C Compute minimal workspace.
+C
+ IF (ISFACT ) THEN
+ MINWRK = MAX( 1, 2*N, 6*N-6 )
+ ELSE
+ MINWRK = MAX( 1, 4*N, 6*N-6 )
+ END IF
+C
+C Check the scalar input parameters.
+C
+ IF ( .NOT.( ISDISC .OR. LSAME( DICO, 'C' ) ) ) THEN
+ INFO = -1
+ ELSEIF ( .NOT.( ISFACT .OR. LSAME( FACT, 'N' ) ) ) THEN
+ INFO = -2
+ ELSEIF ( .NOT.( ISTRAN .OR. LSAME( TRANS, 'N' ) ) ) THEN
+ INFO = -3
+ ELSEIF ( N .LT. 0 ) THEN
+ INFO = -4
+ ELSEIF ( M .LT. 0 ) THEN
+ INFO = -5
+ ELSEIF ( LDA .LT. MAX( 1, N ) ) THEN
+ INFO = -7
+ ELSEIF ( LDE .LT. MAX( 1, N ) ) THEN
+ INFO = -9
+ ELSEIF ( LDQ .LT. MAX( 1, N ) ) THEN
+ INFO = -11
+ ELSEIF ( LDZ .LT. MAX( 1, N ) ) THEN
+ INFO = -13
+ ELSEIF ( ( ISTRAN .AND. ( LDB .LT. MAX( 1, N ) ) ) .OR.
+ $ ( .NOT.ISTRAN .AND. ( LDB .LT. MAX( 1, M, N ) ) ) ) THEN
+ INFO = -15
+ ELSEIF ( LDWORK .LT. MINWRK ) THEN
+ INFO = -21
+ ELSE
+ INFO = 0
+ END IF
+ IF ( INFO .NE. 0 ) THEN
+ CALL XERBLA( 'SG03BD', -INFO )
+ RETURN
+ END IF
+C
+ SCALE = ONE
+C
+C Quick return if possible.
+C
+ MINMN = MIN( M, N )
+ IF ( MINMN .EQ. 0 ) THEN
+ IF ( N.GT.0 )
+ $ CALL DLASET( 'Full', N, N, ZERO, ZERO, B, LDB )
+ DWORK(1) = ONE
+ RETURN
+ ENDIF
+C
+ IF ( ISFACT ) THEN
+C
+C Make sure the upper Hessenberg part of A is quasitriangular.
+C
+ DO 20 I = 1, N-2
+ IF ( A(I+1,I).NE.ZERO .AND. A(I+2,I+1).NE.ZERO ) THEN
+ INFO = 2
+ RETURN
+ END IF
+ 20 CONTINUE
+ END IF
+C
+ IF ( .NOT.ISFACT ) THEN
+C
+C Reduce the pencil A - lambda * E to generalized Schur form.
+C
+C A := Q**T * A * Z (upper quasitriangular)
+C E := Q**T * E * Z (upper triangular)
+C
+C ( Workspace: >= MAX(1,4*N) )
+C
+ CALL DGEGS( 'Vectors', 'Vectors', N, A, LDA, E, LDE, ALPHAR,
+ $ ALPHAI, BETA, Q, LDQ, Z, LDZ, DWORK, LDWORK,
+ $ INFO1 )
+ IF...
[truncated message content] |
