From: <par...@us...> - 2010-01-10 16:19:07
|
Revision: 6730 http://octave.svn.sourceforge.net/octave/?rev=6730&view=rev Author: paramaniac Date: 2010-01-10 16:18:59 +0000 (Sun, 10 Jan 2010) Log Message: ----------- control-oo: add support for Sylvester equations Modified Paths: -------------- trunk/octave-forge/extra/control-oo/inst/dlyap.m trunk/octave-forge/extra/control-oo/inst/lyap.m trunk/octave-forge/extra/control-oo/src/Makefile Added Paths: ----------- trunk/octave-forge/extra/control-oo/src/SB04MD.f trunk/octave-forge/extra/control-oo/src/SB04MR.f trunk/octave-forge/extra/control-oo/src/SB04MU.f trunk/octave-forge/extra/control-oo/src/SB04MW.f trunk/octave-forge/extra/control-oo/src/SB04MY.f trunk/octave-forge/extra/control-oo/src/SB04QD.f trunk/octave-forge/extra/control-oo/src/SB04QR.f trunk/octave-forge/extra/control-oo/src/SB04QU.f trunk/octave-forge/extra/control-oo/src/SB04QY.f trunk/octave-forge/extra/control-oo/src/slsb04md.cc trunk/octave-forge/extra/control-oo/src/slsb04qd.cc Modified: trunk/octave-forge/extra/control-oo/inst/dlyap.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/dlyap.m 2010-01-10 15:16:53 UTC (rev 6729) +++ trunk/octave-forge/extra/control-oo/inst/dlyap.m 2010-01-10 16:18:59 UTC (rev 6730) @@ -17,10 +17,15 @@ ## -*- texinfo -*- ## @deftypefn{Function File} {@var{x} =} dlyap (@var{a}, @var{b}) -## +## @deftypefnx{Function File} {@var{x} =} dlyap (@var{a}, @var{b}, @var{c}) +## Solve discrete-time Lyapunov or Sylvester equations. +## Uses SLICOT SB03MD and SB04QD by courtesy of NICONET e.V. +## <http://www.slicot.org> ## @example ## @group +## AXA' - X + B = 0 (Lyapunov Equation) ## +## AXB' - X + C = 0 (Sylvester Equation) ## @end group ## @end example ## @end deftypefn @@ -29,29 +34,51 @@ ## Created: January 2010 ## Version: 0.1 -function x = dlyap (a, b) +function x = dlyap (a, b, c) - if (nargin != 2) - print_usage (); - endif + if (nargin == 2) - na = issquare (a); - nb = issquare (b); + na = issquare (a); + nb = issquare (b); - if (! na) - error ("lyap: a must be square"); - endif + if (! na) + error ("lyap: a must be square"); + endif - if (! nb) - error ("lyap: b must be square") - endif + if (! nb) + error ("lyap: b must be square") + endif - if (na != nb) - error ("lyap: a and b must be of identical size"); - endif + if (na != nb) + error ("lyap: a and b must be of identical size"); + endif - x = slsb03md (a, -b, true); # AXA' - X = -B + x = slsb03md (a, -b, true); # AXA' - X = -B + + elseif (nargin == 3) + + n = issquare (a); + m = issquare (b); + [crows, ccols] = size (c); + + if (! n) + error ("dlyap: a must be square"); + endif + + if (! m) + error ("dlyap: b must be square"); + endif + + if (crows != n || ccols != m) + error ("dlyap: c must be a (%dx%d) matrix", n, m); + endif + + x = slsb04qd (-a, b, c); # AXB' - X = -C + else + print_usage (); + endif + endfunction @@ -71,4 +98,26 @@ %! 1.0000 3.0000 0.0000 %! 1.0000 0.0000 4.0000]; %! +%!assert (X, X_exp, 1e-4); + +## Sylvester +%!shared X, X_exp +%! A = [1.0 2.0 3.0 +%! 6.0 7.0 8.0 +%! 9.0 2.0 3.0]; +%! +%! B = [7.0 2.0 3.0 +%! 2.0 1.0 2.0 +%! 3.0 4.0 1.0]; +%! +%! C = [271.0 135.0 147.0 +%! 923.0 494.0 482.0 +%! 578.0 383.0 287.0]; +%! +%! X = dlyap (-A, B, C); +%! +%! X_exp = [2.0000 3.0000 6.0000 +%! 4.0000 7.0000 1.0000 +%! 5.0000 3.0000 2.0000]; +%! %!assert (X, X_exp, 1e-4); \ No newline at end of file Modified: trunk/octave-forge/extra/control-oo/inst/lyap.m =================================================================== --- trunk/octave-forge/extra/control-oo/inst/lyap.m 2010-01-10 15:16:53 UTC (rev 6729) +++ trunk/octave-forge/extra/control-oo/inst/lyap.m 2010-01-10 16:18:59 UTC (rev 6730) @@ -16,11 +16,16 @@ ## along with this program. If not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn{Function File} {@var{x} =} lyap (@var{a}, @var{q}) -## +## @deftypefn{Function File} {@var{x} =} lyap (@var{a}, @var{b}) +## @deftypefnx{Function File} {@var{x} =} lyap (@var{a}, @var{b}, @var{c}) +## Solve continuous-time Lyapunov or Sylvester equations. +## Uses SLICOT SB03MD and SB04MD by courtesy of NICONET e.V. +## <http://www.slicot.org> ## @example ## @group +## AX + XA' + B = 0 (Lyapunov Equation) ## +## AX + XB + C = 0 (Sylvester Equation) ## @end group ## @end example ## @end deftypefn @@ -29,36 +34,80 @@ ## Created: January 2010 ## Version: 0.1 -function x = lyap (a, q) +function x = lyap (a, b, c) - if (nargin != 2) - print_usage (); - endif + if (nargin == 2) - na = issquare (a); - nq = issquare (q); + na = issquare (a); + nb = issquare (b); - if (! na) - error ("lyap: a must be square"); - endif + if (! na) + error ("lyap: a must be square"); + endif - if (! nq) - error ("lyap: q must be square") - endif + if (! nb) + error ("lyap: b must be square") + endif - if (na != nq) - error ("lyap: a and q must be of identical size"); + if (na != nb) + error ("lyap: a and b must be of identical size"); + endif + + x = slsb03md (a, -b, false); # AX + XA' = -B + + elseif (nargin == 3) + + n = issquare (a); + m = issquare (b); + [crows, ccols] = size (c); + + if (! n) + error ("lyap: a must be square"); + endif + + if (! m) + error ("lyap: b must be square"); + endif + + if (crows != n || ccols != m) + error ("lyap: c must be a (%dx%d) matrix", n, m); + endif + + x = slsb04md (a, b, -c); # AX + XB = -C + + else + print_usage (); endif - - x = slsb03md (a, -q, false); endfunction +## Lyapunov %!shared X, X_exp %! A = [1, 2; -3, -4]; %! Q = [3, 1; 1, 1]; %! X = lyap (A, Q); %! X_exp = [ 6.1667, -3.8333; %! -3.8333, 3.0000]; -%!assert (X, X_exp, 1e-4); \ No newline at end of file +%!assert (X, X_exp, 1e-4); + +## Sylvester +%!shared X, X_exp +%! A = [2.0 1.0 3.0 +%! 0.0 2.0 1.0 +%! 6.0 1.0 2.0]; +%! +%! B = [2.0 1.0 +%! 1.0 6.0]; +%! +%! C = [2.0 1.0 +%! 1.0 4.0 +%! 0.0 5.0]; +%! +%! X = lyap (A, B, -C); +%! +%! X_exp = [-2.7685 0.5498 +%! -1.0531 0.6865 +%! 4.5257 -0.4389]; +%! +%!assert (X, X_exp, 1e-4); Modified: trunk/octave-forge/extra/control-oo/src/Makefile =================================================================== --- trunk/octave-forge/extra/control-oo/src/Makefile 2010-01-10 15:16:53 UTC (rev 6729) +++ trunk/octave-forge/extra/control-oo/src/Makefile 2010-01-10 16:18:59 UTC (rev 6730) @@ -1,5 +1,6 @@ all: slab08nd.oct slab13dd.oct slsb10hd.oct slsb10ed.oct slab13bd.oct \ - slsb01bd.oct slsb10fd.oct slsb10dd.oct slsb03md.oct + slsb01bd.oct slsb10fd.oct slsb10dd.oct slsb03md.oct slsb04md.oct \ + slsb04qd.oct # transmission zeros of state-space models slab08nd.oct: slab08nd.cc @@ -74,5 +75,15 @@ SB03MD.f select.f SB03MX.f SB03MY.f MB01RD.f \ SB03MV.f SB03MW.f SB04PX.f +# Sylvester equations - continuous-time +slsb04md.oct: slsb04md.cc + mkoctfile slsb04md.cc \ + SB04MD.f SB04MU.f SB04MY.f SB04MR.f SB04MW.f + +# Sylvester equations - discrete-time +slsb04qd.oct: slsb04qd.cc + mkoctfile slsb04qd.cc \ + SB04QD.f SB04QU.f SB04QY.f SB04MW.f SB04QR.f + clean: rm *.o core octave-core *.oct *~ Added: trunk/octave-forge/extra/control-oo/src/SB04MD.f =================================================================== --- trunk/octave-forge/extra/control-oo/src/SB04MD.f (rev 0) +++ trunk/octave-forge/extra/control-oo/src/SB04MD.f 2010-01-10 16:18:59 UTC (rev 6730) @@ -0,0 +1,347 @@ + SUBROUTINE SB04MD( N, M, A, LDA, B, LDB, C, LDC, Z, LDZ, IWORK, + $ 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 the continuous-time Sylvester equation +C +C AX + XB = C +C +C where A, B, C and X are general N-by-N, M-by-M, N-by-M and +C N-by-M matrices respectively. +C +C ARGUMENTS +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 order of the matrix 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 coefficient matrix A of the equation. +C On exit, the leading N-by-N upper Hessenberg part of this +C array contains the matrix H, and the remainder of the +C leading N-by-N part, together with the elements 2,3,...,N +C of array DWORK, contain the orthogonal transformation +C matrix U (stored in factored form). +C +C LDA INTEGER +C The leading dimension of array A. LDA >= MAX(1,N). +C +C B (input/output) DOUBLE PRECISION array, dimension (LDB,M) +C On entry, the leading M-by-M part of this array must +C contain the coefficient matrix B of the equation. +C On exit, the leading M-by-M part of this array contains +C the quasi-triangular Schur factor S of the matrix B'. +C +C LDB INTEGER +C The leading dimension of array B. LDB >= MAX(1,M). +C +C C (input/output) DOUBLE PRECISION array, dimension (LDC,M) +C On entry, the leading N-by-M part of this array must +C contain the coefficient matrix C of the equation. +C On exit, the leading N-by-M part of this array contains +C the solution matrix X of the problem. +C +C LDC INTEGER +C The leading dimension of array C. LDC >= MAX(1,N). +C +C Z (output) DOUBLE PRECISION array, dimension (LDZ,M) +C The leading M-by-M part of this array contains the +C orthogonal matrix Z used to transform B' to real upper +C Schur form. +C +C LDZ INTEGER +C The leading dimension of array Z. LDZ >= MAX(1,M). +C +C Workspace +C +C IWORK INTEGER array, dimension (4*N) +C +C DWORK DOUBLE PRECISION array, dimension (LDWORK) +C On exit, if INFO = 0, DWORK(1) returns the optimal value +C of LDWORK, and DWORK(2), DWORK(3),..., DWORK(N) contain +C the scalar factors of the elementary reflectors used to +C reduce A to upper Hessenberg form, as returned by LAPACK +C Library routine DGEHRD. +C +C LDWORK INTEGER +C The length of the array DWORK. +C LDWORK = MAX(1, 2*N*N + 8*N, 5*M, N + M). +C For optimum 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 > 0: if INFO = i, 1 <= i <= M, the QR algorithm failed to +C compute all the eigenvalues (see LAPACK Library +C routine DGEES); +C > M: if a singular matrix was encountered whilst solving +C for the (INFO-M)-th column of matrix X. +C +C METHOD +C +C The matrix A is transformed to upper Hessenberg form H = U'AU by +C the orthogonal transformation matrix U; matrix B' is transformed +C to real upper Schur form S = Z'B'Z using the orthogonal +C transformation matrix Z. The matrix C is also multiplied by the +C transformations, F = U'CZ, and the solution matrix Y of the +C transformed system +C +C HY + YS' = F +C +C is computed by back substitution. Finally, the matrix Y is then +C multiplied by the orthogonal transformation matrices, X = UYZ', in +C order to obtain the solution matrix X to the original problem. +C +C REFERENCES +C +C [1] Golub, G.H., Nash, S. and Van Loan, C.F. +C A Hessenberg-Schur method for the problem AX + XB = C. +C IEEE Trans. Auto. Contr., AC-24, pp. 909-913, 1979. +C +C NUMERICAL ASPECTS +C 3 3 2 2 +C The algorithm requires about (5/3) N + 10 M + 5 N M + 2.5 M N +C operations and is backward stable. +C +C CONTRIBUTORS +C +C Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, Aug. 1997. +C Supersedes Release 2.0 routine SB04AD by G. Golub, S. Nash, and +C C. Van Loan, Stanford University, California, United States of +C America, January 1982. +C +C REVISIONS +C +C V. Sima, Katholieke Univ. Leuven, Belgium, June 2000, Aug. 2000. +C +C KEYWORDS +C +C Hessenberg form, 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 .. + INTEGER INFO, LDA, LDB, LDC, LDWORK, LDZ, M, N +C .. Array Arguments .. + INTEGER IWORK(*) + DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), DWORK(*), Z(LDZ,*) +C .. Local Scalars .. + INTEGER I, IEIG, IFAIL, IHI, ILO, IND, ITAU, JWORK, + $ SDIM, WRKOPT +C .. Local Scalars .. + LOGICAL SELECT +C .. Local Arrays .. + LOGICAL BWORK(1) +C .. External Subroutines .. + EXTERNAL DCOPY, DGEES, DGEHRD, DGEMM, DGEMV, DLACPY, + $ DORMHR, DSWAP, SB04MU, SB04MY, XERBLA +C .. Intrinsic Functions .. + INTRINSIC INT, MAX +C .. Executable Statements .. +C + INFO = 0 +C +C Test the input scalar arguments. +C + IF( N.LT.0 ) THEN + INFO = -1 + ELSE IF( M.LT.0 ) THEN + INFO = -2 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -4 + ELSE IF( LDB.LT.MAX( 1, M ) ) THEN + INFO = -6 + ELSE IF( LDC.LT.MAX( 1, N ) ) THEN + INFO = -8 + ELSE IF( LDZ.LT.MAX( 1, M ) ) THEN + INFO = -10 + ELSE IF( LDWORK.LT.MAX( 1, 2*N*N + 8*N, 5*M, N + M ) ) THEN + INFO = -13 + END IF +C + IF ( INFO.NE.0 ) THEN +C +C Error return. +C + CALL XERBLA( 'SB04MD', -INFO ) + RETURN + END IF +C +C Quick return if possible. +C + IF ( N.EQ.0 .OR. M.EQ.0 ) THEN + DWORK(1) = ONE + RETURN + END IF +C + ILO = 1 + IHI = N + WRKOPT = 1 +C +C Step 1 : Reduce A to upper Hessenberg and B' to quasi-upper +C triangular. That is, H = U' * A * U (store U in factored +C form) and S = Z' * B' * Z (save Z). +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 + DO 20 I = 2, M + CALL DSWAP( I-1, B(1,I), 1, B(I,1), LDB ) + 20 CONTINUE +C +C Workspace: need 5*M; +C prefer larger. +C + IEIG = M + 1 + JWORK = IEIG + M + CALL DGEES( 'Vectors', 'Not ordered', SELECT, M, B, LDB, + $ SDIM, DWORK, DWORK(IEIG), Z, LDZ, DWORK(JWORK), + $ LDWORK-JWORK+1, BWORK, INFO ) + IF ( INFO.NE.0 ) + $ RETURN + WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) )+JWORK-1 ) +C +C Workspace: need 2*N; +C prefer N + N*NB. +C + ITAU = 2 + JWORK = ITAU + N - 1 + CALL DGEHRD( N, ILO, IHI, A, LDA, DWORK(ITAU), DWORK(JWORK), + $ LDWORK-JWORK+1, IFAIL ) + WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) )+JWORK-1 ) +C +C Step 2 : Form F = ( U' * C ) * Z. Use BLAS 3, if enough space. +C +C Workspace: need N + M; +C prefer N + M*NB. +C + CALL DORMHR( 'Left', 'Transpose', N, M, ILO, IHI, A, LDA, + $ DWORK(ITAU), C, LDC, DWORK(JWORK), LDWORK-JWORK+1, + $ IFAIL ) + WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) )+JWORK-1 ) +C + IF ( LDWORK.GE.JWORK - 1 + N*M ) THEN + CALL DGEMM( 'No transpose', 'No transpose', N, M, M, ONE, C, + $ LDC, Z, LDZ, ZERO, DWORK(JWORK), N ) + CALL DLACPY( 'Full', N, M, DWORK(JWORK), N, C, LDC ) + WRKOPT = MAX( WRKOPT, JWORK - 1 + N*M ) + ELSE +C + DO 40 I = 1, N + CALL DGEMV( 'Transpose', M, M, ONE, Z, LDZ, C(I,1), LDC, + $ ZERO, DWORK(JWORK), 1 ) + CALL DCOPY( M, DWORK(JWORK), 1, C(I,1), LDC ) + 40 CONTINUE +C + END IF +C + IND = M + 60 CONTINUE + IF ( IND.GT.1 ) THEN +C +C Step 3 : Solve H * Y + Y * S' = F for Y. +C + IF ( B(IND,IND-1).EQ.ZERO ) THEN +C +C Solve a special linear algebraic system of order N. +C Workspace: N*(N+1)/2 + 3*N. +C + CALL SB04MY( M, N, IND, A, LDA, B, LDB, C, LDC, + $ DWORK(JWORK), IWORK, INFO ) +C + IF ( INFO.NE.0 ) THEN + INFO = INFO + M + RETURN + END IF + WRKOPT = MAX( WRKOPT, JWORK + N*( N + 1 )/2 + 2*N - 1 ) + IND = IND - 1 + ELSE +C +C Solve a special linear algebraic system of order 2*N. +C Workspace: 2*N*N + 8*N; +C + CALL SB04MU( M, N, IND, A, LDA, B, LDB, C, LDC, + $ DWORK(JWORK), IWORK, INFO ) +C + IF ( INFO.NE.0 ) THEN + INFO = INFO + M + RETURN + END IF + WRKOPT = MAX( WRKOPT, JWORK + 2*N*N + 7*N - 1 ) + IND = IND - 2 + END IF + GO TO 60 + ELSE IF ( IND.EQ.1 ) THEN +C +C Solve a special linear algebraic system of order N. +C Workspace: N*(N+1)/2 + 3*N; +C + CALL SB04MY( M, N, IND, A, LDA, B, LDB, C, LDC, + $ DWORK(JWORK), IWORK, INFO ) + IF ( INFO.NE.0 ) THEN + INFO = INFO + M + RETURN + END IF + WRKOPT = MAX( WRKOPT, JWORK + N*( N + 1 )/2 + 2*N - 1 ) + END IF +C +C Step 4 : Form C = ( U * Y ) * Z'. Use BLAS 3, if enough space. +C +C Workspace: need N + M; +C prefer N + M*NB. +C + CALL DORMHR( 'Left', 'No transpose', N, M, ILO, IHI, A, LDA, + $ DWORK(ITAU), C, LDC, DWORK(JWORK), LDWORK-JWORK+1, + $ IFAIL ) + WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) )+JWORK-1 ) +C + IF ( LDWORK.GE.JWORK - 1 + N*M ) THEN + CALL DGEMM( 'No transpose', 'Transpose', N, M, M, ONE, C, LDC, + $ Z, LDZ, ZERO, DWORK(JWORK), N ) + CALL DLACPY( 'Full', N, M, DWORK(JWORK), N, C, LDC ) + ELSE +C + DO 80 I = 1, N + CALL DGEMV( 'No transpose', M, M, ONE, Z, LDZ, C(I,1), LDC, + $ ZERO, DWORK(JWORK), 1 ) + CALL DCOPY( M, DWORK(JWORK), 1, C(I,1), LDC ) + 80 CONTINUE + END IF +C + RETURN +C *** Last line of SB04MD *** + END Added: trunk/octave-forge/extra/control-oo/src/SB04MR.f =================================================================== --- trunk/octave-forge/extra/control-oo/src/SB04MR.f (rev 0) +++ trunk/octave-forge/extra/control-oo/src/SB04MR.f 2010-01-10 16:18:59 UTC (rev 6730) @@ -0,0 +1,222 @@ + SUBROUTINE SB04MR( M, D, IPR, 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 a linear algebraic system of order M whose coefficient +C matrix has zeros below the second subdiagonal. The matrix is +C stored compactly, row-wise. +C +C ARGUMENTS +C +C Input/Output Parameters +C +C M (input) INTEGER +C The order of the system. M >= 0. +C Note that parameter M should have twice the value in the +C original problem (see SLICOT Library routine SB04MU). +C +C D (input/output) DOUBLE PRECISION array, dimension +C (M*(M+1)/2+3*M) +C On entry, the first M*(M+1)/2 + 2*M elements of this array +C must contain the coefficient matrix, stored compactly, +C row-wise, and the next M elements must contain the right +C hand side of the linear system, as set by SLICOT Library +C routine SB04MU. +C On exit, the content of this array is updated, the last M +C elements containing the solution with components +C interchanged (see IPR). +C +C IPR (output) INTEGER array, dimension (2*M) +C The leading M elements contain information about the +C row interchanges performed for solving the system. +C Specifically, the i-th component of the solution is +C specified by IPR(i). +C +C Error Indicator +C +C INFO INTEGER +C = 0: successful exit; +C = 1: if a singular matrix was encountered. +C +C METHOD +C +C Gaussian elimination with partial pivoting is used. The rows of +C the matrix are not actually permuted, only their indices are +C interchanged in array IPR. +C +C REFERENCES +C +C [1] Golub, G.H., Nash, S. and Van Loan, C.F. +C A Hessenberg-Schur method for the problem AX + XB = C. +C IEEE Trans. Auto. Contr., AC-24, pp. 909-913, 1979. +C +C NUMERICAL ASPECTS +C +C None. +C +C CONTRIBUTORS +C +C Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, Sep. 1997. +C Supersedes Release 2.0 routine SB04AR by G. Golub, S. Nash, and +C C. Van Loan, Stanford University, California, United States of +C America, January 1982. +C +C REVISIONS +C +C - +C +C KEYWORDS +C +C Hessenberg form, orthogonal transformation, real Schur form, +C Sylvester equation. +C +C ****************************************************************** +C + DOUBLE PRECISION ZERO + PARAMETER ( ZERO = 0.0D0 ) +C .. Scalar Arguments .. + INTEGER INFO, M +C .. Array Arguments .. + INTEGER IPR(*) + DOUBLE PRECISION D(*) +C .. Local Scalars .. + INTEGER I, I1, I2, IPRM, IPRM1, J, K, L, M1, MPI, MPI1, + $ MPI2 + DOUBLE PRECISION D1, D2, D3, DMAX +C .. External Subroutines .. + EXTERNAL DAXPY +C .. Intrinsic Functions .. + INTRINSIC ABS +C .. Executable Statements .. +C + INFO = 0 + I2 = ( M*( M + 5 ) )/2 + MPI = M + IPRM = I2 + M1 = M + I1 = 1 +C + DO 20 I = 1, M + MPI = MPI + 1 + IPRM = IPRM + 1 + IPR(MPI) = I1 + IPR(I) = IPRM + I1 = I1 + M1 + IF ( I.GE.3 ) M1 = M1 - 1 + 20 CONTINUE +C + M1 = M - 1 + MPI1 = M + 1 +C +C Reduce to upper triangular form. +C + DO 80 I = 1, M1 + MPI = MPI1 + MPI1 = MPI1 + 1 + IPRM = IPR(MPI) + D1 = D(IPRM) + I1 = 2 + IF ( I.EQ.M1 ) I1 = 1 + MPI2 = MPI + I1 + L = 0 + DMAX = ABS( D1 ) +C + DO 40 J = MPI1, MPI2 + D2 = D(IPR(J)) + D3 = ABS( D2 ) + IF ( D3.GT.DMAX ) THEN + DMAX = D3 + D1 = D2 + L = J - MPI + END IF + 40 CONTINUE +C +C Check singularity. +C + IF ( DMAX.EQ.ZERO ) THEN + INFO = 1 + RETURN + END IF +C + IF ( L.GT.0 ) THEN +C +C Permute the row indices. +C + K = IPRM + J = MPI + L + IPRM = IPR(J) + IPR(J) = K + IPR(MPI) = IPRM + K = IPR(I) + I2 = I + L + IPR(I) = IPR(I2) + IPR(I2) = K + END IF + IPRM = IPRM + 1 +C +C Annihilate the subdiagonal elements of the matrix. +C + I2 = I + D3 = D(IPR(I)) +C + DO 60 J = MPI1, MPI2 + I2 = I2 + 1 + IPRM1 = IPR(J) + DMAX = -D(IPRM1)/D1 + D(IPR(I2)) = D(IPR(I2)) + DMAX*D3 + CALL DAXPY( M-I, DMAX, D(IPRM), 1, D(IPRM1+1), 1 ) + 60 CONTINUE +C + IPR(MPI1) = IPR(MPI1) + 1 + IF ( I.NE.M1 ) IPR(MPI2) = IPR(MPI2) + 1 + 80 CONTINUE +C + MPI = M + M + IPRM = IPR(MPI) +C +C Check singularity. +C + IF ( D(IPRM).EQ.ZERO ) THEN + INFO = 1 + RETURN + END IF +C +C Back substitution. +C + D(IPR(M)) = D(IPR(M))/D(IPRM) +C + DO 120 I = M1, 1, -1 + MPI = MPI - 1 + IPRM = IPR(MPI) + IPRM1 = IPRM + DMAX = ZERO +C + DO 100 K = I+1, M + IPRM1 = IPRM1 + 1 + DMAX = DMAX + D(IPR(K))*D(IPRM1) + 100 CONTINUE +C + D(IPR(I)) = ( D(IPR(I)) - DMAX )/D(IPRM) + 120 CONTINUE +C + RETURN +C *** Last line of SB04MR *** + END Added: trunk/octave-forge/extra/control-oo/src/SB04MU.f =================================================================== --- trunk/octave-forge/extra/control-oo/src/SB04MU.f (rev 0) +++ trunk/octave-forge/extra/control-oo/src/SB04MU.f 2010-01-10 16:18:59 UTC (rev 6730) @@ -0,0 +1,190 @@ + SUBROUTINE SB04MU( N, M, IND, A, LDA, B, LDB, C, LDC, D, IPR, + $ 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 construct and solve a linear algebraic system of order 2*M +C whose coefficient matrix has zeros below the second subdiagonal. +C Such systems appear when solving continuous-time Sylvester +C equations using the Hessenberg-Schur method. +C +C ARGUMENTS +C +C Input/Output Parameters +C +C N (input) INTEGER +C The order of the matrix B. N >= 0. +C +C M (input) INTEGER +C The order of the matrix A. M >= 0. +C +C IND (input) INTEGER +C IND and IND - 1 specify the indices of the columns in C +C to be computed. IND > 1. +C +C A (input) DOUBLE PRECISION array, dimension (LDA,M) +C The leading M-by-M part of this array must contain an +C upper Hessenberg matrix. +C +C LDA INTEGER +C The leading dimension of array A. LDA >= MAX(1,M). +C +C B (input) DOUBLE PRECISION array, dimension (LDB,N) +C The leading N-by-N part of this array must contain a +C matrix in real Schur form. +C +C LDB INTEGER +C The leading dimension of array B. LDB >= MAX(1,N). +C +C C (input/output) DOUBLE PRECISION array, dimension (LDC,N) +C On entry, the leading M-by-N part of this array must +C contain the coefficient matrix C of the equation. +C On exit, the leading M-by-N part of this array contains +C the matrix C with columns IND-1 and IND updated. +C +C LDC INTEGER +C The leading dimension of array C. LDC >= MAX(1,M). +C +C Workspace +C +C D DOUBLE PRECISION array, dimension (2*M*M+7*M) +C +C IPR INTEGER array, dimension (4*M) +C +C Error Indicator +C +C INFO INTEGER +C = 0: successful exit; +C > 0: if INFO = IND, a singular matrix was encountered. +C +C METHOD +C +C A special linear algebraic system of order 2*M, whose coefficient +C matrix has zeros below the second subdiagonal is constructed and +C solved. The coefficient matrix is stored compactly, row-wise. +C +C REFERENCES +C +C [1] Golub, G.H., Nash, S. and Van Loan, C.F. +C A Hessenberg-Schur method for the problem AX + XB = C. +C IEEE Trans. Auto. Contr., AC-24, pp. 909-913, 1979. +C +C NUMERICAL ASPECTS +C +C None. +C +C CONTRIBUTORS +C +C Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, Sep. 1997. +C Supersedes Release 2.0 routine SB04AU by G. Golub, S. Nash, and +C C. Van Loan, Stanford University, California, United States of +C America, January 1982. +C +C REVISIONS +C +C - +C +C KEYWORDS +C +C Hessenberg form, orthogonal transformation, real Schur form, +C Sylvester equation. +C +C ****************************************************************** +C + DOUBLE PRECISION ZERO + PARAMETER ( ZERO = 0.0D0 ) +C .. Scalar Arguments .. + INTEGER INFO, IND, LDA, LDB, LDC, M, N +C .. Array Arguments .. + INTEGER IPR(*) + DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), D(*) +C .. Local Scalars .. + INTEGER I, I2, IND1, J, K, K1, K2, M2 + DOUBLE PRECISION TEMP +C .. External Subroutines .. + EXTERNAL DAXPY, SB04MR +C .. Intrinsic Functions .. + INTRINSIC MAX, MIN +C .. Executable Statements .. +C + IND1 = IND - 1 +C + DO 20 I = IND + 1, N + CALL DAXPY( M, -B(IND1,I), C(1,I), 1, C(1,IND1), 1 ) + CALL DAXPY( M, -B(IND,I), C(1,I), 1, C(1,IND), 1 ) + 20 CONTINUE +C +C Construct the linear algebraic system of order 2*M. +C + K1 = -1 + M2 = 2*M + I2 = M*(M2 + 5) + K = M2 +C + DO 60 I = 1, M +C + DO 40 J = MAX( 1, I - 1 ), M + K1 = K1 + 2 + K2 = K1 + K + TEMP = A(I,J) + IF ( I.NE.J ) THEN + D(K1) = TEMP + D(K1+1) = ZERO + IF ( J.GT.I ) D(K2) = ZERO + D(K2+1) = TEMP + ELSE + D(K1) = TEMP + B(IND1,IND1) + D(K1+1) = B(IND1,IND) + D(K2) = B(IND,IND1) + D(K2+1) = TEMP + B(IND,IND) + END IF + 40 CONTINUE +C + K1 = K2 + K = K - MIN( 2, I ) +C +C Store the right hand side. +C + I2 = I2 + 2 + D(I2) = C(I,IND) + D(I2-1) = C(I,IND1) + 60 CONTINUE +C +C Solve the linear algebraic system and store the solution in C. +C + CALL SB04MR( M2, D, IPR, INFO ) +C + IF ( INFO.NE.0 ) THEN + INFO = IND + ELSE + I2 = 0 +C + DO 80 I = 1, M + I2 = I2 + 2 + C(I,IND1) = D(IPR(I2-1)) + C(I,IND) = D(IPR(I2)) + 80 CONTINUE +C + END IF +C + RETURN +C *** Last line of SB04MU *** + END Added: trunk/octave-forge/extra/control-oo/src/SB04MW.f =================================================================== --- trunk/octave-forge/extra/control-oo/src/SB04MW.f (rev 0) +++ trunk/octave-forge/extra/control-oo/src/SB04MW.f 2010-01-10 16:18:59 UTC (rev 6730) @@ -0,0 +1,194 @@ + SUBROUTINE SB04MW( M, D, IPR, 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 a linear algebraic system of order M whose coefficient +C matrix is in upper Hessenberg form, stored compactly, row-wise. +C +C ARGUMENTS +C +C Input/Output Parameters +C +C M (input) INTEGER +C The order of the system. M >= 0. +C +C D (input/output) DOUBLE PRECISION array, dimension +C (M*(M+1)/2+2*M) +C On entry, the first M*(M+1)/2 + M elements of this array +C must contain an upper Hessenberg matrix, stored compactly, +C row-wise, and the next M elements must contain the right +C hand side of the linear system, as set by SLICOT Library +C routine SB04MY. +C On exit, the content of this array is updated, the last M +C elements containing the solution with components +C interchanged (see IPR). +C +C IPR (output) INTEGER array, dimension (2*M) +C The leading M elements contain information about the +C row interchanges performed for solving the system. +C Specifically, the i-th component of the solution is +C specified by IPR(i). +C +C Error Indicator +C +C INFO INTEGER +C = 0: successful exit; +C = 1: if a singular matrix was encountered. +C +C METHOD +C +C Gaussian elimination with partial pivoting is used. The rows of +C the matrix are not actually permuted, only their indices are +C interchanged in array IPR. +C +C REFERENCES +C +C [1] Golub, G.H., Nash, S. and Van Loan, C.F. +C A Hessenberg-Schur method for the problem AX + XB = C. +C IEEE Trans. Auto. Contr., AC-24, pp. 909-913, 1979. +C +C NUMERICAL ASPECTS +C +C None. +C +C CONTRIBUTORS +C +C Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, Sep. 1997. +C Supersedes Release 2.0 routine SB04AW by G. Golub, S. Nash, and +C C. Van Loan, Stanford University, California, United States of +C America, January 1982. +C +C REVISIONS +C +C - +C +C KEYWORDS +C +C Hessenberg form, orthogonal transformation, real Schur form, +C Sylvester equation. +C +C ****************************************************************** +C + DOUBLE PRECISION ZERO + PARAMETER ( ZERO = 0.0D0 ) +C .. Scalar Arguments .. + INTEGER INFO, M +C .. Array Arguments .. + INTEGER IPR(*) + DOUBLE PRECISION D(*) +C .. Local Scalars .. + INTEGER I, I1, IPRM, IPRM1, K, M1, M2, MPI + DOUBLE PRECISION D1, D2, MULT +C .. External Subroutines .. + EXTERNAL DAXPY +C .. Intrinsic Functions .. + INTRINSIC ABS +C .. Executable Statements .. +C + INFO = 0 + M1 = ( M*( M + 3 ) )/2 + M2 = M + M + MPI = M + IPRM = M1 + M1 = M + I1 = 1 +C + DO 20 I = 1, M + MPI = MPI + 1 + IPRM = IPRM + 1 + IPR(MPI) = I1 + IPR(I) = IPRM + I1 = I1 + M1 + IF ( I.GT.1 ) M1 = M1 - 1 + 20 CONTINUE +C + M1 = M - 1 + MPI = M +C +C Reduce to upper triangular form. +C + DO 40 I = 1, M1 + I1 = I + 1 + MPI = MPI + 1 + IPRM = IPR(MPI) + IPRM1 = IPR(MPI+1) + D1 = D(IPRM) + D2 = D(IPRM1) + IF ( ABS( D1 ).LE.ABS( D2 ) ) THEN +C +C Permute the row indices. +C + K = IPRM + IPR(MPI) = IPRM1 + IPRM = IPRM1 + IPRM1 = K + K = IPR(I) + IPR(I) = IPR(I1) + IPR(I1) = K + D1 = D2 + END IF +C +C Check singularity. +C + IF ( D1.EQ.ZERO ) THEN + INFO = 1 + RETURN + END IF +C + MULT = -D(IPRM1)/D1 + IPRM1 = IPRM1 + 1 + IPR(MPI+1) = IPRM1 +C +C Annihilate the subdiagonal elements of the matrix. +C + D(IPR(I1)) = D(IPR(I1)) + MULT*D(IPR(I)) + CALL DAXPY( M-I, MULT, D(IPRM+1), 1, D(IPRM1), 1 ) + 40 CONTINUE +C +C Check singularity. +C + IF ( D(IPR(M2)).EQ.ZERO ) THEN + INFO = 1 + RETURN + END IF +C +C Back substitution. +C + D(IPR(M)) = D(IPR(M))/D(IPR(M2)) + MPI = M2 +C + DO 80 I = M1, 1, -1 + MPI = MPI - 1 + IPRM = IPR(MPI) + IPRM1 = IPRM + MULT = ZERO +C + DO 60 I1 = I + 1, M + IPRM1 = IPRM1 + 1 + MULT = MULT + D(IPR(I1))*D(IPRM1) + 60 CONTINUE +C + D(IPR(I)) = ( D(IPR(I)) - MULT )/D(IPRM) + 80 CONTINUE +C + RETURN +C *** Last line of SB04MW *** + END Added: trunk/octave-forge/extra/control-oo/src/SB04MY.f =================================================================== --- trunk/octave-forge/extra/control-oo/src/SB04MY.f (rev 0) +++ trunk/octave-forge/extra/control-oo/src/SB04MY.f 2010-01-10 16:18:59 UTC (rev 6730) @@ -0,0 +1,168 @@ + SUBROUTINE SB04MY( N, M, IND, A, LDA, B, LDB, C, LDC, D, IPR, + $ 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 construct and solve a linear algebraic system of order M whose +C coefficient matrix is in upper Hessenberg form. Such systems +C appear when solving Sylvester equations using the Hessenberg-Schur +C method. +C +C ARGUMENTS +C +C Input/Output Parameters +C +C N (input) INTEGER +C The order of the matrix B. N >= 0. +C +C M (input) INTEGER +C The order of the matrix A. M >= 0. +C +C IND (input) INTEGER +C The index of the column in C to be computed. IND >= 1. +C +C A (input) DOUBLE PRECISION array, dimension (LDA,M) +C The leading M-by-M part of this array must contain an +C upper Hessenberg matrix. +C +C LDA INTEGER +C The leading dimension of array A. LDA >= MAX(1,M). +C +C B (input) DOUBLE PRECISION array, dimension (LDB,N) +C The leading N-by-N part of this array must contain a +C matrix in real Schur form. +C +C LDB INTEGER +C The leading dimension of array B. LDB >= MAX(1,N). +C +C C (input/output) DOUBLE PRECISION array, dimension (LDC,N) +C On entry, the leading M-by-N part of this array must +C contain the coefficient matrix C of the equation. +C On exit, the leading M-by-N part of this array contains +C the matrix C with column IND updated. +C +C LDC INTEGER +C The leading dimension of array C. LDC >= MAX(1,M). +C +C Workspace +C +C D DOUBLE PRECISION array, dimension (M*(M+1)/2+2*M) +C +C IPR INTEGER array, dimension (2*M) +C +C Error Indicator +C +C INFO INTEGER +C = 0: successful exit; +C > 0: if INFO = IND, a singular matrix was encountered. +C +C METHOD +C +C A special linear algebraic system of order M, with coefficient +C matrix in upper Hessenberg form is constructed and solved. The +C coefficient matrix is stored compactly, row-wise. +C +C REFERENCES +C +C [1] Golub, G.H., Nash, S. and Van Loan, C.F. +C A Hessenberg-Schur method for the problem AX + XB = C. +C IEEE Trans. Auto. Contr., AC-24, pp. 909-913, 1979. +C +C NUMERICAL ASPECTS +C +C None. +C +C CONTRIBUTORS +C +C Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, Sep. 1997. +C Supersedes Release 2.0 routine SB04AY by G. Golub, S. Nash, and +C C. Van Loan, Stanford University, California, United States of +C America, January 1982. +C +C REVISIONS +C +C - +C +C KEYWORDS +C +C Hessenberg form, orthogonal transformation, real Schur form, +C Sylvester equation. +C +C ****************************************************************** +C +C .. Scalar Arguments .. + INTEGER INFO, IND, LDA, LDB, LDC, M, N +C .. Array Arguments .. + INTEGER IPR(*) + DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), D(*) +C .. Local Scalars .. + INTEGER I, I2, J, K, K1, K2, M1 +C .. External Subroutines .. + EXTERNAL DAXPY, DCOPY, SB04MW +C .. Intrinsic Functions .. + INTRINSIC MAX +C .. Executable Statements .. +C + DO 20 I = IND + 1, N + CALL DAXPY( M, -B(IND,I), C(1,I), 1, C(1,IND), 1 ) + 20 CONTINUE +C + M1 = M + 1 + I2 = ( M*M1 )/2 + M1 + K2 = 1 + K = M +C +C Construct the linear algebraic system of order M. +C + DO 40 I = 1, M + J = M1 - K + CALL DCOPY ( K, A(I,J), LDA, D(K2), 1 ) + K1 = K2 + K2 = K2 + K + IF ( I.GT.1 ) THEN + K1 = K1 + 1 + K = K - 1 + END IF + D(K1) = D(K1) + B(IND,IND) +C +C Store the right hand side. +C + D(I2) = C(I,IND) + I2 = I2 + 1 + 40 CONTINUE +C +C Solve the linear algebraic system and store the solution in C. +C + CALL SB04MW( M, D, IPR, INFO ) +C + IF ( INFO.NE.0 ) THEN + INFO = IND + ELSE +C + DO 60 I = 1, M + C(I,IND) = D(IPR(I)) + 60 CONTINUE +C + END IF +C + RETURN +C *** Last line of SB04MY *** + END Added: trunk/octave-forge/extra/control-oo/src/SB04QD.f =================================================================== --- trunk/octave-forge/extra/control-oo/src/SB04QD.f (rev 0) +++ trunk/octave-forge/extra/control-oo/src/SB04QD.f 2010-01-10 16:18:59 UTC (rev 6730) @@ -0,0 +1,376 @@ + SUBROUTINE SB04QD( N, M, A, LDA, B, LDB, C, LDC, Z, LDZ, IWORK, + $ 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 the discrete-time Sylvester equation +C +C X + AXB = C, +C +C where A, B, C and X are general N-by-N, M-by-M, N-by-M and +C N-by-M matrices respectively. A Hessenberg-Schur method, which +C reduces A to upper Hessenberg form, H = U'AU, and B' to real +C Schur form, S = Z'B'Z (with U, Z orthogonal matrices), is used. +C +C ARGUMENTS +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 order of the matrix 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 coefficient matrix A of the equation. +C On exit, the leading N-by-N upper Hessenberg part of this +C array contains the matrix H, and the remainder of the +C leading N-by-N part, together with the elements 2,3,...,N +C of array DWORK, contain the orthogonal transformation +C matrix U (stored in factored form). +C +C LDA INTEGER +C The leading dimension of array A. LDA >= MAX(1,N). +C +C B (input/output) DOUBLE PRECISION array, dimension (LDB,M) +C On entry, the leading M-by-M part of this array must +C contain the coefficient matrix B of the equation. +C On exit, the leading M-by-M part of this array contains +C the quasi-triangular Schur factor S of the matrix B'. +C +C LDB INTEGER +C The leading dimension of array B. LDB >= MAX(1,M). +C +C C (input/output) DOUBLE PRECISION array, dimension (LDC,M) +C On entry, the leading N-by-M part of this array must +C contain the coefficient matrix C of the equation. +C On exit, the leading N-by-M part of this array contains +C the solution matrix X of the problem. +C +C LDC INTEGER +C The leading dimension of array C. LDC >= MAX(1,N). +C +C Z (output) DOUBLE PRECISION array, dimension (LDZ,M) +C The leading M-by-M part of this array contains the +C orthogonal matrix Z used to transform B' to real upper +C Schur form. +C +C LDZ INTEGER +C The leading dimension of array Z. LDZ >= MAX(1,M). +C +C Workspace +C +C IWORK INTEGER array, dimension (4*N) +C +C DWORK DOUBLE PRECISION array, dimension (LDWORK) +C On exit, if INFO = 0, DWORK(1) returns the optimal value +C of LDWORK, and DWORK(2), DWORK(3),..., DWORK(N) contain +C the scalar factors of the elementary reflectors used to +C reduce A to upper Hessenberg form, as returned by LAPACK +C Library routine DGEHRD. +C +C LDWORK INTEGER +C The length of the array DWORK. +C LDWORK = MAX(1, 2*N*N + 9*N, 5*M, N + M). +C For optimum 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 > 0: if INFO = i, 1 <= i <= M, the QR algorithm failed to +C compute all the eigenvalues of B (see LAPACK Library +C routine DGEES); +C > M: if a singular matrix was encountered whilst solving +C for the (INFO-M)-th column of matrix X. +C +C METHOD +C +C The matrix A is transformed to upper Hessenberg form H = U'AU by +C the orthogonal transformation matrix U; matrix B' is transformed +C to real upper Schur form S = Z'B'Z using the orthogonal +C transformation matrix Z. The matrix C is also multiplied by the +C transformations, F = U'CZ, and the solution matrix Y of the +C transformed system +C +C Y + HYS' = F +C +C is computed by back substitution. Finally, the matrix Y is then +C multiplied by the orthogonal transformation matrices, X = UYZ', in +C order to obtain the solution matrix X to the original problem. +C +C REFERENCES +C +C [1] Golub, G.H., Nash, S. and Van Loan, C.F. +C A Hessenberg-Schur method for the problem AX + XB = C. +C IEEE Trans. Auto. Contr., AC-24, pp. 909-913, 1979. +C +C [2] Sima, V. +C Algorithms for Linear-quadratic Optimization. +C Marcel Dekker, Inc., New York, 1996. +C +C NUMERICAL ASPECTS +C 3 3 2 2 +C The algorithm requires about (5/3) N + 10 M + 5 N M + 2.5 M N +C operations and is backward stable. +C +C CONTRIBUTORS +C +C D. Sima, University of Bucharest, May 2000, Aug. 2000. +C +C REVISIONS +C +C V. Sima, Research Institute for Informatics, Bucharest, May 2000. +C +C KEYWORDS +C +C Hessenberg form, orthogonal transformation, real Schur form, +C Sylvester equation. +C +C ****************************************************************** +C +C .. Parameters .. + DOUBLE PRECISION ONE, ZERO + PARAMETER ( ONE = 1.0D0, ZERO = 0.0D0 ) +C .. Scalar Arguments .. + INTEGER INFO, LDA, LDB, LDC, LDWORK, LDZ, M, N +C .. Array Arguments .. + INTEGER IWORK(*) + DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), DWORK(*), Z(LDZ,*) +C .. Local Scalars .. + INTEGER BL, CHUNK, I, IEIG, IFAIL, IHI, ILO, IND, ITAU, + $ JWORK, SDIM, WRKOPT +C .. Local Scalars .. + LOGICAL BLAS3, BLOCK +C .. Local Arrays .. + LOGICAL BWORK(1) +C .. External Functions .. + LOGICAL SELECT +C .. External Subroutines .. + EXTERNAL DCOPY, DGEES, DGEHRD, DGEMM, DGEMV, DLACPY, + $ DORMHR, DSWAP, SB04QU, SB04QY, XERBLA +C .. Intrinsic Functions .. + INTRINSIC INT, MAX, MIN +C .. Executable Statements .. +C + INFO = 0 +C +C Test the input scalar arguments. +C + IF( N.LT.0 ) THEN + INFO = -1 + ELSE IF( M.LT.0 ) THEN + INFO = -2 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -4 + ELSE IF( LDB.LT.MAX( 1, M ) ) THEN + INFO = -6 + ELSE IF( LDC.LT.MAX( 1, N ) ) THEN + INFO = -8 + ELSE IF( LDZ.LT.MAX( 1, M ) ) THEN + INFO = -10 + ELSE IF( LDWORK.LT.MAX( 1, 2*N*N + 9*N, 5*M, N + M ) ) THEN + INFO = -13 + END IF +C + IF ( INFO.NE.0 ) THEN +C +C Error return. +C + CALL XERBLA( 'SB04QD', -INFO ) + RETURN + END IF +C +C Quick return if possible. +C + IF ( N.EQ.0 .OR. M.EQ.0 ) THEN + DWORK(1) = ONE + RETURN + END IF +C + ILO = 1 + IHI = N + WRKOPT = 2*N*N + 9*N +C +C Step 1 : Reduce A to upper Hessenberg and B' to quasi-upper +C triangular. That is, H = U' * A * U (store U in factored +C form) and S = Z' * B' * Z (save Z). +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 + DO 20 I = 2, M + CALL DSWAP( I-1, B(1,I), 1, B(I,1), LDB ) + 20 CONTINUE +C +C Workspace: need 5*M; +C prefer larger. +C + IEIG = M + 1 + JWORK = IEIG + M + CALL DGEES( 'Vectors', 'Not ordered', SELECT, M, B, LDB, + $ SDIM, DWORK, DWORK(IEIG), Z, LDZ, DWORK(JWORK), + $ LDWORK-JWORK+1, BWORK, INFO ) + IF ( INFO.NE.0 ) + $ RETURN + WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) )+JWORK-1 ) +C +C Workspace: need 2*N; +C prefer N + N*NB. +C + ITAU = 2 + JWORK = ITAU + N - 1 + CALL DGEHRD( N, ILO, IHI, A, LDA, DWORK(ITAU), DWORK(JWORK), + $ LDWORK-JWORK+1, IFAIL ) + WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) )+JWORK-1 ) +C +C Step 2 : Form F = ( U' * C ) * Z. Use BLAS 3, if enough space. +C +C Workspace: need N + M; +C prefer N + M*NB. +C + CALL DORMHR( 'Left', 'Transpose', N, M, ILO, IHI, A, LDA, + $ DWORK(ITAU), C, LDC, DWORK(JWORK), LDWORK-JWORK+1, + $ IFAIL ) + WRKOPT = MAX( WRKOPT, MAX( INT( DWORK(JWORK) ), N*M )+JWORK-1 ) +C + CHUNK = ( LDWORK - JWORK + 1 ) / M + BLOCK = MIN( CHUNK, N ).GT.1 + BLAS3 = CHUNK.GE.N .AND. BLOCK +C + IF ( BLAS3 ) THEN + CALL DGEMM( 'No transpose', 'No transpose', N, M, M, ONE, C, + $ LDC, Z, LDZ, ZERO, DWORK(JWORK), N ) + CALL DLACPY( 'Full', N, M, DWORK(JWORK), N, C, LDC ) +C + ELSE IF ( BLOCK ) THEN +C +C Use as many rows of C as possible. +C + DO 40 I = 1, N, CHUNK + BL = MIN( N-I+1, CHUNK ) + CALL DGEMM( 'NoTranspose', 'NoTranspose', BL, M, M, ONE, + $ C(I,1), LDC, Z, LDZ, ZERO, DWORK(JWORK), BL ) + CALL DLACPY( 'Full', BL, M, DWORK(JWORK), BL, C(I,1), LDC ) + 40 CONTINUE +C + ELSE +C + DO 60 I = 1, N + CALL DGEMV( 'Transpose', M, M, ONE, Z, LDZ, C(I,1), LDC, + $ ZERO, DWORK(JWORK), 1 ) + CALL DCOPY( M, DWORK(JWORK), 1, C(I,1), LDC ) + 60 CONTINUE +C + END IF +C +C Step 3 : Solve Y + H * Y * S' = F for Y. +C + IND = M + 80 CONTINUE +C + IF ( IND.GT.1 ) THEN + IF ( B(IND,IND-1).EQ.ZERO ) THEN +C +C Solve a special linear algebraic system of order N. +C Workspace: N*(N+1)/2 + 3*N. +C + CALL SB04QY( M, N, IND, A, LDA, B, LDB, C, LDC, + $ DWORK(JWORK), IWORK, INFO ) +C + IF ( INFO.NE.0 ) THEN + INFO = INFO + M + RETURN + END IF + IND = IND - 1 + ELSE +C +C Solve a special linear algebraic system of order 2*N. +C Workspace: 2*N*N + 9*N; +C + CALL SB04QU( M, N, IND, A, LDA, B, LDB, C, LDC, + $ DWORK(JWORK), IWORK, INFO ) +C + IF ( INFO.NE.0 ) THEN + INFO = INFO + M + RETURN + END IF + IND = IND - 2 + END IF + GO TO 80 + ELSE IF ( IND.EQ.1 ) THEN +C +C Solve a special linear algebraic system of order N. +C Workspace: N*(N+1)/2 + 3*N; +C + CALL SB04QY( M, N, IND, A, LDA, B, LDB, C, LDC, + $ DWORK(JWORK), IWORK, INFO ) + IF ( INFO.NE.0 ) THEN + INFO = INFO + M + RETURN + END IF + END IF +C +C Step 4 : Form C = ( U * Y ) * Z'. Use BLAS 3, if enough space. +C +C Workspace: need N + M; +C prefer N + M*NB. +C + CALL DORMHR( 'Left', 'No transpose', N, M, ILO, IHI, A, LDA, + $ DWORK(ITAU), C, LDC, DWORK(JWORK), LDWORK-JWORK+1, + $ IFAIL ) + WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) )+JWORK-1 ) +C + IF ( BLAS3 ) THEN + CALL DGEMM( 'No transpose', 'Transpose', N, M, M, ONE, C, LDC, + $ Z, LDZ, ZERO, DWORK(JWORK), N ) + CALL DLACPY( 'Full', N, M, DWORK(JWORK), N, C, LDC ) +C + ELSE IF ( BLOCK ) THEN +C +C Use as many rows of C as possible. +C + DO 100 I = 1, N, CHUNK + BL = MIN( N-I+1, CHUNK ) + CALL DGEMM( 'NoTranspose', 'Transpose', BL, M, M, ONE, + $ C(I,1), LDC, Z, LDZ, ZERO, DWORK(JWORK), BL ) + CALL DLACPY( 'Full', BL, M, DWORK(JWORK), BL, C(I,1), LDC ) + 100 CONTINUE +C + ELSE +C + DO 120 I = 1, N + CALL DGEMV( 'No transpose', M, M, ONE, Z, LDZ, C(I,1), LDC, + $ ZERO, DWORK(JWORK), 1 ) + CALL DCOPY( M, DWORK(JWORK), 1, C(I,1), LDC ) + 120 CONTINUE + END IF +C + RETURN +C *** Last line of SB04QD *** + END Added: trunk/octave-forge/extra/control-oo/src/SB04QR.f =================================================================== --- trunk/octave-forge/extra/control-oo/src/SB04QR.f (rev 0) +++ trunk/octave-forge/extra/control-oo/src/SB04QR.f 2010-01-10 16:18:59 UTC (rev 6730) @@ -0,0 +1,224 @@ + SUBROUTINE SB04QR( M, D, IPR, 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 a linear algebraic system of order M whose coefficient +C matrix has zeros below the third subdiagonal and zero elements on +C the third subdiagonal with even column indices. The matrix is +C stored compactly, row-wise. +C +C ARGUMENTS +C +C Input/Output Parameters +C +C M (input) INTEGER +C The order of the system. M >= 0, M even. +C Note that parameter M should have twice the value in the +C original problem (see SLICOT Library routine SB04QU). +C +C D (input/output) DOUBLE PRECISION array, dimension +C (M*M/2+4*M) +C On entry, the first M*M/2 + 3*M elements of this array +C must contain the coefficient matrix, stored compactly, +C row-wise, and the next M elements must contain the right +C hand side of the linear system, as set by SLICOT Library +C routine SB04QU. +C On exit, the content of this array is updated, the last M +C elements containing the solution with components +C interchanged (see IPR). +C +C IPR (output) INTEGER array, dimension (2*M) +C The leading M elements contain information about the +C row interchanges performed for solving the system. +C Specifically, the i-th component of the solution is +C specified by IPR(i). +C +C Error Indicator +C +C INFO INTEGER +C = 0: successful exit; +C = 1: if a singular matrix was encountered. +C +C METHOD +C +C Gaussian elimination with partial pivoting is used. The rows of +C the matrix are not actually permuted, only their indices are +C interchanged in array IPR. +C +C REFERENCES +C +C [1] Golub, G.H., Nash, S. and Van Loan, C.F. +C A Hessenberg-Schur metho... [truncated message content] |