Menu
▾
▴
[Octave-cvsupdate] SF.net SVN: octave:[7673] trunk/octave-forge/main/control
|
From: <par...@us...> - 2010-09-07 08:41:04
|
Revision: 7673
http://octave.svn.sourceforge.net/octave/?rev=7673&view=rev
Author: paramaniac
Date: 2010-09-07 08:40:56 +0000 (Tue, 07 Sep 2010)
Log Message:
-----------
control: use SLICOT for ss:minreal
Modified Paths:
--------------
trunk/octave-forge/main/control/DESCRIPTION
trunk/octave-forge/main/control/inst/@lti/minreal.m
trunk/octave-forge/main/control/inst/@ss/__minreal__.m
trunk/octave-forge/main/control/inst/ltimodels.m
trunk/octave-forge/main/control/src/Makefile
Added Paths:
-----------
trunk/octave-forge/main/control/src/AB07MD.f
trunk/octave-forge/main/control/src/TB01PD.f
trunk/octave-forge/main/control/src/TB01UD.f
trunk/octave-forge/main/control/src/TB01XD.f
trunk/octave-forge/main/control/src/sltb01pd.cc
Modified: trunk/octave-forge/main/control/DESCRIPTION
===================================================================
--- trunk/octave-forge/main/control/DESCRIPTION 2010-09-07 06:25:09 UTC (rev 7672)
+++ trunk/octave-forge/main/control/DESCRIPTION 2010-09-07 08:40:56 UTC (rev 7673)
@@ -1,6 +1,6 @@
Name: Control
Version: 2.0.0
-Date: 2010-08-25
+Date: 2010-09-07
Author: Lukas Reichlin <luk...@gm...>
Maintainer: Lukas Reichlin <luk...@gm...>
Title: Control Systems
Modified: trunk/octave-forge/main/control/inst/@lti/minreal.m
===================================================================
--- trunk/octave-forge/main/control/inst/@lti/minreal.m 2010-09-07 06:25:09 UTC (rev 7672)
+++ trunk/octave-forge/main/control/inst/@lti/minreal.m 2010-09-07 08:40:56 UTC (rev 7673)
@@ -1,4 +1,4 @@
-## Copyright (C) 2009 Lukas F. Reichlin
+## Copyright (C) 2009 - 2010 Lukas F. Reichlin
##
## This file is part of LTI Syncope.
##
@@ -23,10 +23,18 @@
## Author: Lukas Reichlin <luk...@gm...>
## Created: October 2009
-## Version: 0.1
+## Version: 0.2
function sys = minreal (sys, tol = "def")
+ if (nargin > 2) # nargin == 0 not possible because minreal is inside @lti
+ print_usage ();
+ endif
+
+ if (isa (tol, "lti") || (tol != "def" && ! isreal (tol) && ! isscalar (tol)))
+ error ("minreal: second argument must be a real scalar");
+ endif
+
sys = __minreal__ (sys, tol);
endfunction
\ No newline at end of file
Modified: trunk/octave-forge/main/control/inst/@ss/__minreal__.m
===================================================================
--- trunk/octave-forge/main/control/inst/@ss/__minreal__.m 2010-09-07 06:25:09 UTC (rev 7672)
+++ trunk/octave-forge/main/control/inst/@ss/__minreal__.m 2010-09-07 08:40:56 UTC (rev 7673)
@@ -1,4 +1,4 @@
-## Copyright (C) 2009 Lukas F. Reichlin
+## Copyright (C) 2009 - 2010 Lukas F. Reichlin
##
## This file is part of LTI Syncope.
##
@@ -17,55 +17,27 @@
## -*- texinfo -*-
## Minimal realization of SS models. The physical meaning of states is lost.
+## Uses SLICOT TB01PD by courtesy of NICONET e.V. <http://www.slicot.org>
-## Algorithm based on sysmin by A. Scottedward Hodel
## Author: Lukas Reichlin <luk...@gm...>
## Created: October 2009
-## Version: 0.1
+## Version: 0.2
function retsys = __minreal__ (sys, tol)
- A = sys.a;
- B = sys.b;
- C = sys.c;
-
- if (! isempty (A))
- if (tol == "def")
- [cflg, Uc] = isctrb (A, B);
- else
- [cflg, Uc] = isctrb (A, B, tol);
- endif
-
- if (! cflg)
- if (! isempty (Uc))
- A = Uc.' * A * Uc;
- B = Uc.' * B;
- C = C * Uc;
- else
- A = B = C = [];
- endif
- endif
+ if (tol == "def")
+ tol = 0;
+ elseif (tol > 1)
+ error ("ss: minreal: require tol <= 1");
endif
- if (! isempty (A))
- if (tol == "def")
- [oflg, Uo] = isobsv (A, C);
- else
- [oflg, Uo] = isobsv (A, C, tol);
- endif
+ [a, b, c, nr] = sltb01pd (sys.a, sys.b, sys.c, tol);
- if (! oflg)
- if (! isempty (Uo))
- A = Uo.' * A * Uo;
- B = Uo.' * B;
- C = C * Uo;
- else
- A = B = C = [];
- endif
- endif
- endif
+ a = a(1:nr, 1:nr);
+ b = b(1:nr, :);
+ c = c(:, 1:nr);
- retsys = ss (A, B, C, sys.d);
+ retsys = ss (a, b, c, sys.d);
retsys.lti = sys.lti; # retain i/o names and tsam
endfunction
Modified: trunk/octave-forge/main/control/inst/ltimodels.m
===================================================================
--- trunk/octave-forge/main/control/inst/ltimodels.m 2010-09-07 06:25:09 UTC (rev 7672)
+++ trunk/octave-forge/main/control/inst/ltimodels.m 2010-09-07 08:40:56 UTC (rev 7673)
@@ -251,13 +251,13 @@
%! 23.0 5.0 7.0 14.0 16.0
%! 4.0 6.0 13.0 20.0 22.0
%! 10.0 12.0 19.0 21.0 3.0
-%! 11.0 18.0 25.0 2.0 9.0];
+%! 11.0 18.0 25.0 2.0 9.0 ];
%!
%! B = [ -1.0 -4.0
%! 4.0 9.0
%! -9.0 -16.0
%! 16.0 25.0
-%! -25.0 -36.0];
+%! -25.0 -36.0 ];
%!
%! tol = 0;
%!
@@ -270,13 +270,50 @@
%! 4.4741 -12.5544 5.3509 5.9403 1.4360
%! 14.4576 7.6855 23.1452 26.3872 -29.9557
%! 0.0000 1.4805 27.4668 22.6564 -0.0072
-%! 0.0000 0.0000 -30.4822 0.6745 18.8680];
+%! 0.0000 0.0000 -30.4822 0.6745 18.8680 ];
%!
%! Bce = [ 31.1199 47.6865
%! 3.2480 0.0000
%! 0.0000 0.0000
%! 0.0000 0.0000
-%! 0.0000 0.0000];
+%! 0.0000 0.0000 ];
%!
%!assert (Ac, Ace, 1e-4);
%!assert (Bc, Bce, 1e-4);
+
+## ss: minreal (SLICOT TB01PD)
+%!shared M, Me
+%! A = [ 1.0 2.0 0.0
+%! 4.0 -1.0 0.0
+%! 0.0 0.0 1.0 ];
+%!
+%! B = [ 1.0
+%! 0.0
+%! 1.0 ];
+%!
+%! C = [ 0.0 1.0 -1.0
+%! 0.0 0.0 1.0 ];
+%!
+%! D = zeros (2, 1);
+%!
+%! sys = ss (A, B, C, D);
+%! sysmin = minreal (sys, 0.0);
+%! [Ar, Br, Cr, Dr] = ssdata (sysmin);
+%! M = [Ar, Br; Cr, Dr];
+%!
+%! Ae = [ 1.0000 -1.4142 1.4142
+%! -2.8284 0.0000 1.0000
+%! 2.8284 1.0000 0.0000 ];
+%!
+%! Be = [-1.0000
+%! 0.7071
+%! 0.7071 ];
+%!
+%! Ce = [ 0.0000 0.0000 -1.4142
+%! 0.0000 0.7071 0.7071 ];
+%!
+%! De = zeros (2, 1);
+%!
+%! Me = [Ae, Be; Ce, De];
+%!
+%!assert (M, Me, 1e-4);
Added: trunk/octave-forge/main/control/src/AB07MD.f
===================================================================
--- trunk/octave-forge/main/control/src/AB07MD.f (rev 0)
+++ trunk/octave-forge/main/control/src/AB07MD.f 2010-09-07 08:40:56 UTC (rev 7673)
@@ -0,0 +1,224 @@
+ SUBROUTINE AB07MD( JOBD, N, M, P, A, LDA, B, LDB, C, LDC, D, LDD,
+ $ 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 find the dual of a given state-space representation.
+C
+C ARGUMENTS
+C
+C Mode Parameters
+C
+C JOBD CHARACTER*1
+C Specifies whether or not a non-zero matrix D appears in
+C the given state space model:
+C = 'D': D is present;
+C = 'Z': D is assumed a zero matrix.
+C
+C Input/Output Parameters
+C
+C N (input) INTEGER
+C The order of the state-space representation. N >= 0.
+C
+C M (input) INTEGER
+C The number of system inputs. M >= 0.
+C
+C P (input) INTEGER
+C The number of system outputs. P >= 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 original state dynamics matrix A.
+C On exit, the leading N-by-N part of this array contains
+C the dual state dynamics matrix A'.
+C
+C LDA INTEGER
+C The leading dimension of array A. LDA >= MAX(1,N).
+C
+C B (input/output) DOUBLE PRECISION array, dimension
+C (LDB,MAX(M,P))
+C On entry, the leading N-by-M part of this array must
+C contain the original input/state matrix B.
+C On exit, the leading N-by-P part of this array contains
+C the dual input/state matrix C'.
+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 P-by-N part of this array must
+C contain the original state/output matrix C.
+C On exit, the leading M-by-N part of this array contains
+C the dual state/output matrix B'.
+C
+C LDC INTEGER
+C The leading dimension of array C.
+C LDC >= MAX(1,M,P) if N > 0.
+C LDC >= 1 if N = 0.
+C
+C D (input/output) DOUBLE PRECISION array, dimension
+C (LDD,MAX(M,P))
+C On entry, if JOBD = 'D', the leading P-by-M part of this
+C array must contain the original direct transmission
+C matrix D.
+C On exit, if JOBD = 'D', the leading M-by-P part of this
+C array contains the dual direct transmission matrix D'.
+C The array D is not referenced if JOBD = 'Z'.
+C
+C LDD INTEGER
+C The leading dimension of array D.
+C LDD >= MAX(1,M,P) if JOBD = 'D'.
+C LDD >= 1 if JOBD = 'Z'.
+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
+C METHOD
+C
+C If the given state-space representation is the M-input/P-output
+C (A,B,C,D), its dual is simply the P-input/M-output (A',C',B',D').
+C
+C REFERENCES
+C
+C None
+C
+C NUMERICAL ASPECTS
+C
+C None
+C
+C CONTRIBUTOR
+C
+C Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, Dec. 1996.
+C Supersedes Release 2.0 routine AB07AD by T.W.C.Williams, Kingston
+C Polytechnic, United Kingdom, March 1982.
+C
+C REVISIONS
+C
+C V. Sima, Research Institute for Informatics, Bucharest, Feb. 2004.
+C
+C KEYWORDS
+C
+C Dual system, state-space model, state-space representation.
+C
+C ******************************************************************
+C
+C .. Scalar Arguments ..
+ CHARACTER JOBD
+ INTEGER INFO, LDA, LDB, LDC, LDD, M, N, P
+C .. Array Arguments ..
+ DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), D(LDD,*)
+C .. Local Scalars ..
+ LOGICAL LJOBD
+ INTEGER J, MINMP, MPLIM
+C .. External functions ..
+ LOGICAL LSAME
+ EXTERNAL LSAME
+C .. External subroutines ..
+ EXTERNAL DCOPY, DSWAP, XERBLA
+C .. Intrinsic Functions ..
+ INTRINSIC MAX, MIN
+C .. Executable Statements ..
+C
+ INFO = 0
+ LJOBD = LSAME( JOBD, 'D' )
+ MPLIM = MAX( M, P )
+ MINMP = MIN( M, P )
+C
+C Test the input scalar arguments.
+C
+ IF( .NOT.LJOBD .AND. .NOT.LSAME( JOBD, 'Z' ) ) THEN
+ INFO = -1
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -2
+ ELSE IF( M.LT.0 ) THEN
+ INFO = -3
+ ELSE IF( P.LT.0 ) THEN
+ INFO = -4
+ ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+ INFO = -6
+ ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
+ INFO = -8
+ ELSE IF( ( N.GT.0 .AND. LDC.LT.MAX( 1, MPLIM ) ) .OR.
+ $ ( N.EQ.0 .AND. LDC.LT.1 ) ) THEN
+ INFO = -10
+ ELSE IF( ( LJOBD .AND. LDD.LT.MAX( 1, MPLIM ) ) .OR.
+ $ ( .NOT.LJOBD .AND. LDD.LT.1 ) ) THEN
+ INFO = -12
+ END IF
+C
+ IF ( INFO.NE.0 ) THEN
+C
+C Error return.
+C
+ CALL XERBLA( 'AB07MD', -INFO )
+ RETURN
+ END IF
+C
+C Quick return if possible.
+C
+ IF ( MAX( N, MINMP ).EQ.0 )
+ $ RETURN
+C
+ IF ( N.GT.0 ) THEN
+C
+C Transpose A, if non-scalar.
+C
+ DO 10 J = 1, N - 1
+ CALL DSWAP( N-J, A(J+1,J), 1, A(J,J+1), LDA )
+ 10 CONTINUE
+C
+C Replace B by C' and C by B'.
+C
+ DO 20 J = 1, MPLIM
+ IF ( J.LE.MINMP ) THEN
+ CALL DSWAP( N, B(1,J), 1, C(J,1), LDC )
+ ELSE IF ( J.GT.P ) THEN
+ CALL DCOPY( N, B(1,J), 1, C(J,1), LDC )
+ ELSE
+ CALL DCOPY( N, C(J,1), LDC, B(1,J), 1 )
+ END IF
+ 20 CONTINUE
+C
+ END IF
+C
+ IF ( LJOBD .AND. MINMP.GT.0 ) THEN
+C
+C Transpose D, if non-scalar.
+C
+ DO 30 J = 1, MPLIM
+ IF ( J.LT.MINMP ) THEN
+ CALL DSWAP( MINMP-J, D(J+1,J), 1, D(J,J+1), LDD )
+ ELSE IF ( J.GT.P ) THEN
+ CALL DCOPY( P, D(1,J), 1, D(J,1), LDD )
+ ELSE IF ( J.GT.M ) THEN
+ CALL DCOPY( M, D(J,1), LDD, D(1,J), 1 )
+ END IF
+ 30 CONTINUE
+C
+ END IF
+C
+ RETURN
+C *** Last line of AB07MD ***
+ END
Modified: trunk/octave-forge/main/control/src/Makefile
===================================================================
--- trunk/octave-forge/main/control/src/Makefile 2010-09-07 06:25:09 UTC (rev 7672)
+++ trunk/octave-forge/main/control/src/Makefile 2010-09-07 08:40:56 UTC (rev 7673)
@@ -1,6 +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
+ slsb04qd.oct slsg03ad.oct slsb02od.oct slab13ad.oct slab01od.oct \
+ sltb01pd.oct
# transmission zeros of state-space models
slab08nd.oct: slab08nd.cc
@@ -111,5 +112,11 @@
mkoctfile slab01od.cc \
AB01OD.f AB01ND.f MB03OY.f MB01PD.f MB01QD.f
+# 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
+
clean:
rm *.o core octave-core *.oct *~
Added: trunk/octave-forge/main/control/src/TB01PD.f
===================================================================
--- trunk/octave-forge/main/control/src/TB01PD.f (rev 0)
+++ trunk/octave-forge/main/control/src/TB01PD.f 2010-09-07 08:40:56 UTC (rev 7673)
@@ -0,0 +1,352 @@
+ SUBROUTINE TB01PD( JOB, EQUIL, N, M, P, A, LDA, B, LDB, C, LDC,
+ $ NR, TOL, 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 find a reduced (controllable, observable, or minimal) state-
+C space representation (Ar,Br,Cr) for any original state-space
+C representation (A,B,C). The matrix Ar is in upper block
+C Hessenberg form.
+C
+C ARGUMENTS
+C
+C Mode Parameters
+C
+C JOB CHARACTER*1
+C Indicates whether the user wishes to remove the
+C uncontrollable and/or unobservable parts as follows:
+C = 'M': Remove both the uncontrollable and unobservable
+C parts to get a minimal state-space representation;
+C = 'C': Remove the uncontrollable part only to get a
+C controllable state-space representation;
+C = 'O': Remove the unobservable part only to get an
+C observable state-space representation.
+C
+C EQUIL CHARACTER*1
+C Specifies whether the user wishes to preliminarily balance
+C the triplet (A,B,C) as follows:
+C = 'S': Perform balancing (scaling);
+C = 'N': Do not perform balancing.
+C
+C Input/Output Parameters
+C
+C N (input) INTEGER
+C The order of the original state-space representation, i.e.
+C the order of the matrix A. N >= 0.
+C
+C M (input) INTEGER
+C The number of system inputs. M >= 0.
+C
+C P (input) INTEGER
+C The number of system outputs. P >= 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 original state dynamics matrix A.
+C On exit, the leading NR-by-NR part of this array contains
+C the upper block Hessenberg state dynamics matrix Ar of a
+C minimal, controllable, or observable realization for the
+C original system, depending on the value of JOB, JOB = 'M',
+C JOB = 'C', or JOB = 'O', respectively.
+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 if JOB = 'C', or (LDB,MAX(M,P)), otherwise.
+C On entry, the leading N-by-M part of this array must
+C contain the original input/state matrix B; if JOB = 'M',
+C or JOB = 'O', the remainder of the leading N-by-MAX(M,P)
+C part is used as internal workspace.
+C On exit, the leading NR-by-M part of this array contains
+C the transformed input/state matrix Br of a minimal,
+C controllable, or observable realization for the original
+C system, depending on the value of JOB, JOB = 'M',
+C JOB = 'C', or JOB = 'O', respectively.
+C If JOB = 'C', only the first IWORK(1) rows of B are
+C nonzero.
+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 P-by-N part of this array must
+C contain the original state/output matrix C; if JOB = 'M',
+C or JOB = 'O', the remainder of the leading MAX(M,P)-by-N
+C part is used as internal workspace.
+C On exit, the leading P-by-NR part of this array contains
+C the transformed state/output matrix Cr of a minimal,
+C controllable, or observable realization for the original
+C system, depending on the value of JOB, JOB = 'M',
+C JOB = 'C', or JOB = 'O', respectively.
+C If JOB = 'M', or JOB = 'O', only the last IWORK(1) columns
+C (in the first NR columns) of C are nonzero.
+C
+C LDC INTEGER
+C The leading dimension of array C.
+C LDC >= MAX(1,M,P) if N > 0.
+C LDC >= 1 if N = 0.
+C
+C NR (output) INTEGER
+C The order of the reduced state-space representation
+C (Ar,Br,Cr) of a minimal, controllable, or observable
+C realization for the original system, depending on
+C JOB = 'M', JOB = 'C', or JOB = 'O'.
+C
+C Tolerances
+C
+C TOL DOUBLE PRECISION
+C The tolerance to be used in rank determination when
+C transforming (A, B, C). If the user sets TOL > 0, then
+C the given value of TOL is used as a lower bound for the
+C reciprocal condition number (see the description of the
+C argument RCOND in the SLICOT routine MB03OD); a
+C (sub)matrix whose estimated condition number is less than
+C 1/TOL is considered to be of full rank. If the user sets
+C TOL <= 0, then an implicitly computed, default tolerance
+C (determined by the SLICOT routine TB01UD) is used instead.
+C
+C Workspace
+C
+C IWORK INTEGER array, dimension (N+MAX(M,P))
+C On exit, if INFO = 0, the first nonzero elements of
+C IWORK(1:N) return the orders of the diagonal blocks of A.
+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 length of the array DWORK.
+C LDWORK >= MAX(1, N + MAX(N, 3*M, 3*P)).
+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
+C METHOD
+C
+C If JOB = 'M', the matrices A and B are operated on by orthogonal
+C similarity transformations (made up of products of Householder
+C transformations) so as to produce an upper block Hessenberg matrix
+C A1 and a matrix B1 with all but its first rank(B) rows zero; this
+C separates out the controllable part of the original system.
+C Applying the same algorithm to the dual of this subsystem,
+C therefore separates out the controllable and observable (i.e.
+C minimal) part of the original system representation, with the
+C final Ar upper block Hessenberg (after using pertransposition).
+C If JOB = 'C', or JOB = 'O', only the corresponding part of the
+C above procedure is applied.
+C
+C REFERENCES
+C
+C [1] Van Dooren, P.
+C The Generalized Eigenstructure Problem in Linear System
+C Theory. (Algorithm 1)
+C IEEE Trans. Auto. Contr., AC-26, pp. 111-129, 1981.
+C
+C NUMERICAL ASPECTS
+C 3
+C The algorithm requires 0(N ) operations and is backward stable.
+C
+C CONTRIBUTOR
+C
+C V. Sima, Katholieke Univ. Leuven, Belgium, Feb. 1998.
+C
+C REVISIONS
+C
+C A. Varga, DLR Oberpfaffenhofen, July 1998.
+C A. Varga, DLR Oberpfaffenhofen, April 28, 1999.
+C V. Sima, Research Institute for Informatics, Bucharest, Mar. 2004.
+C
+C KEYWORDS
+C
+C Hessenberg form, minimal realization, multivariable system,
+C orthogonal transformation, state-space model, state-space
+C representation.
+C
+C ******************************************************************
+C
+C .. Parameters ..
+ INTEGER LDIZ
+ PARAMETER ( LDIZ = 1 )
+ DOUBLE PRECISION ZERO, ONE
+ PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
+C .. Scalar Arguments ..
+ CHARACTER EQUIL, JOB
+ INTEGER INFO, LDA, LDB, LDC, LDWORK, M, N, NR, P
+ DOUBLE PRECISION TOL
+C .. Array Arguments ..
+ INTEGER IWORK(*)
+ DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), DWORK(*)
+C .. Local Scalars ..
+ LOGICAL LEQUIL, LNJOBC, LNJOBO
+ INTEGER I, INDCON, ITAU, IZ, JWORK, KL, MAXMP, NCONT,
+ $ WRKOPT
+ DOUBLE PRECISION MAXRED
+C .. External Functions ..
+ LOGICAL LSAME
+ EXTERNAL LSAME
+C .. External Subroutines ..
+ EXTERNAL AB07MD, TB01ID, TB01UD, TB01XD, XERBLA
+C .. Intrinsic Functions ..
+ INTRINSIC INT, MAX, MIN
+C .. Executable Statements ..
+C
+ INFO = 0
+ MAXMP = MAX( M, P )
+ LNJOBC = .NOT.LSAME( JOB, 'C' )
+ LNJOBO = .NOT.LSAME( JOB, 'O' )
+ LEQUIL = LSAME( EQUIL, 'S' )
+C
+C Test the input scalar arguments.
+C
+ IF( LNJOBC .AND. LNJOBO .AND. .NOT.LSAME( JOB, 'M' ) ) THEN
+ INFO = -1
+ ELSE IF( .NOT.LEQUIL .AND. .NOT.LSAME( EQUIL, 'N' ) ) THEN
+ INFO = -2
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -3
+ ELSE IF( M.LT.0 ) THEN
+ INFO = -4
+ ELSE IF( P.LT.0 ) THEN
+ INFO = -5
+ ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+ INFO = -7
+ ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
+ INFO = -9
+ ELSE IF( LDC.LT.1 .OR. ( N.GT.0 .AND. LDC.LT.MAXMP ) ) THEN
+ INFO = -11
+ ELSE IF( LDWORK.LT.MAX( 1, N + MAX( N, 3*MAXMP ) ) ) THEN
+ INFO = -16
+ END IF
+C
+ IF ( INFO.NE.0 ) THEN
+C
+C Error return.
+C
+ CALL XERBLA( 'TB01PD', -INFO )
+ RETURN
+ END IF
+C
+C Quick return if possible.
+C
+ IF ( N.EQ.0 .OR. ( LNJOBC .AND. MIN( N, P ).EQ.0 ) .OR.
+ $ ( LNJOBO .AND. MIN( N, M ).EQ.0 ) ) THEN
+ NR = 0
+C
+ DO 5 I = 1, N
+ IWORK(I) = 0
+ 5 CONTINUE
+C
+ DWORK(1) = ONE
+ RETURN
+ END IF
+C
+C If required, balance the triplet (A,B,C) (default MAXRED).
+C Workspace: need N.
+C
+C (Note: Comments in the code beginning "Workspace:" describe the
+C minimal amount of real workspace needed at that point in the code,
+C as well as the preferred amount for good performance.)
+C
+ IF ( LEQUIL ) THEN
+ MAXRED = ZERO
+ CALL TB01ID( 'A', N, M, P, MAXRED, A, LDA, B, LDB, C, LDC,
+ $ DWORK, INFO )
+ WRKOPT = N
+ ELSE
+ WRKOPT = 1
+ END IF
+C
+ IZ = 1
+ ITAU = 1
+ JWORK = ITAU + N
+ IF ( LNJOBO ) THEN
+C
+C Separate out controllable subsystem (of order NCONT):
+C A <-- Z'*A*Z, B <-- Z'*B, C <-- C*Z.
+C
+C Workspace: need N + MAX(N, 3*M, P).
+C prefer larger.
+C
+ CALL TB01UD( 'No Z', N, M, P, A, LDA, B, LDB, C, LDC, NCONT,
+ $ INDCON, IWORK, DWORK(IZ), LDIZ, DWORK(ITAU), TOL,
+ $ IWORK(N+1), DWORK(JWORK), LDWORK-JWORK+1, INFO )
+C
+ WRKOPT = INT( DWORK(JWORK) ) + JWORK - 1
+ ELSE
+ NCONT = N
+ END IF
+C
+ IF ( LNJOBC ) THEN
+C
+C Separate out the observable subsystem (of order NR):
+C Form the dual of the subsystem of order NCONT (which is
+C controllable, if JOB = 'M'), leaving rest as it is.
+C
+ CALL AB07MD( 'Z', NCONT, M, P, A, LDA, B, LDB, C, LDC, DWORK,
+ $ 1, INFO )
+C
+C And separate out the controllable part of this dual subsystem.
+C
+C Workspace: need NCONT + MAX(NCONT, 3*P, M).
+C prefer larger.
+C
+ CALL TB01UD( 'No Z', NCONT, P, M, A, LDA, B, LDB, C, LDC, NR,
+ $ INDCON, IWORK, DWORK(IZ), LDIZ, DWORK(ITAU), TOL,
+ $ IWORK(N+1), DWORK(JWORK), LDWORK-JWORK+1, INFO )
+C
+ WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) )+JWORK-1 )
+C
+C Transpose and reorder (to get a block upper Hessenberg
+C matrix A), giving, for JOB = 'M', the controllable and
+C observable (i.e., minimal) part of original system.
+C
+ IF( INDCON.GT.0 ) THEN
+ KL = IWORK(1) - 1
+ IF ( INDCON.GE.2 )
+ $ KL = KL + IWORK(2)
+ ELSE
+ KL = 0
+ END IF
+ CALL TB01XD( 'Zero D', NR, P, M, KL, MAX( 0, NR-1 ), A, LDA,
+ $ B, LDB, C, LDC, DWORK, 1, INFO )
+ ELSE
+ NR = NCONT
+ END IF
+C
+C Annihilate the trailing components of IWORK(1:N).
+C
+ DO 10 I = INDCON + 1, N
+ IWORK(I) = 0
+ 10 CONTINUE
+C
+C Set optimal workspace dimension.
+C
+ DWORK(1) = WRKOPT
+ RETURN
+C *** Last line of TB01PD ***
+ END
Added: trunk/octave-forge/main/control/src/TB01UD.f
===================================================================
--- trunk/octave-forge/main/control/src/TB01UD.f (rev 0)
+++ trunk/octave-forge/main/control/src/TB01UD.f 2010-09-07 08:40:56 UTC (rev 7673)
@@ -0,0 +1,491 @@
+ SUBROUTINE TB01UD( JOBZ, N, M, P, A, LDA, B, LDB, C, LDC, NCONT,
+ $ INDCON, NBLK, Z, LDZ, TAU, TOL, 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 find a controllable realization for the linear time-invariant
+C multi-input system
+C
+C dX/dt = A * X + B * U,
+C Y = C * X,
+C
+C where A, B, and C are N-by-N, N-by-M, and P-by-N matrices,
+C respectively, and A and B are reduced by this routine to
+C orthogonal canonical form using (and optionally accumulating)
+C orthogonal similarity transformations, which are also applied
+C to C. Specifically, the system (A, B, C) is reduced to the
+C triplet (Ac, Bc, Cc), where Ac = Z' * A * Z, Bc = Z' * B,
+C Cc = C * Z, with
+C
+C [ Acont * ] [ Bcont ]
+C Ac = [ ], Bc = [ ],
+C [ 0 Auncont ] [ 0 ]
+C
+C and
+C
+C [ A11 A12 . . . A1,p-1 A1p ] [ B1 ]
+C [ A21 A22 . . . A2,p-1 A2p ] [ 0 ]
+C [ 0 A32 . . . A3,p-1 A3p ] [ 0 ]
+C Acont = [ . . . . . . . ], Bc = [ . ],
+C [ . . . . . . ] [ . ]
+C [ . . . . . ] [ . ]
+C [ 0 0 . . . Ap,p-1 App ] [ 0 ]
+C
+C where the blocks B1, A21, ..., Ap,p-1 have full row ranks and
+C p is the controllability index of the pair. The size of the
+C block Auncont is equal to the dimension of the uncontrollable
+C subspace of the pair (A, B).
+C
+C ARGUMENTS
+C
+C Mode Parameters
+C
+C JOBZ CHARACTER*1
+C Indicates whether the user wishes to accumulate in a
+C matrix Z the orthogonal similarity transformations for
+C reducing the system, as follows:
+C = 'N': Do not form Z and do not store the orthogonal
+C transformations;
+C = 'F': Do not form Z, but store the orthogonal
+C transformations in the factored form;
+C = 'I': Z is initialized to the unit matrix and the
+C orthogonal transformation matrix Z is returned.
+C
+C Input/Output Parameters
+C
+C N (input) INTEGER
+C The order of the original state-space representation,
+C i.e. the order of the matrix A. N >= 0.
+C
+C M (input) INTEGER
+C The number of system inputs, or of columns of B. M >= 0.
+C
+C P (input) INTEGER
+C The number of system outputs, or of rows of C. P >= 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 original state dynamics matrix A.
+C On exit, the leading NCONT-by-NCONT part contains the
+C upper block Hessenberg state dynamics matrix Acont in Ac,
+C given by Z' * A * Z, of a controllable realization for
+C the original system. The elements below the first block-
+C subdiagonal are set to zero. The leading N-by-N part
+C contains the matrix Ac.
+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 N-by-M part of this array must
+C contain the input matrix B.
+C On exit, the leading NCONT-by-M part of this array
+C contains the transformed input matrix Bcont in Bc, given
+C by Z' * B, with all elements but the first block set to
+C zero. The leading N-by-M part contains the matrix Bc.
+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 P-by-N part of this array must
+C contain the output matrix C.
+C On exit, the leading P-by-N part of this array contains
+C the transformed output matrix Cc, given by C * Z.
+C
+C LDC INTEGER
+C The leading dimension of array C. LDC >= MAX(1,P).
+C
+C NCONT (output) INTEGER
+C The order of the controllable state-space representation.
+C
+C INDCON (output) INTEGER
+C The controllability index of the controllable part of the
+C system representation.
+C
+C NBLK (output) INTEGER array, dimension (N)
+C The leading INDCON elements of this array contain the
+C the orders of the diagonal blocks of Acont.
+C
+C Z (output) DOUBLE PRECISION array, dimension (LDZ,N)
+C If JOBZ = 'I', then the leading N-by-N part of this
+C array contains the matrix of accumulated orthogonal
+C similarity transformations which reduces the given system
+C to orthogonal canonical form.
+C If JOBZ = 'F', the elements below the diagonal, with the
+C array TAU, represent the orthogonal transformation matrix
+C as a product of elementary reflectors. The transformation
+C matrix can then be obtained by calling the LAPACK Library
+C routine DORGQR.
+C If JOBZ = 'N', the array Z is not referenced and can be
+C supplied as a dummy array (i.e. set parameter LDZ = 1 and
+C declare this array to be Z(1,1) in the calling program).
+C
+C LDZ INTEGER
+C The leading dimension of array Z. If JOBZ = 'I' or
+C JOBZ = 'F', LDZ >= MAX(1,N); if JOBZ = 'N', LDZ >= 1.
+C
+C TAU (output) DOUBLE PRECISION array, dimension (N)
+C The elements of TAU contain the scalar factors of the
+C elementary reflectors used in the reduction of B and A.
+C
+C Tolerances
+C
+C TOL DOUBLE PRECISION
+C The tolerance to be used in rank determination when
+C transforming (A, B). If the user sets TOL > 0, then
+C the given value of TOL is used as a lower bound for the
+C reciprocal condition number (see the description of the
+C argument RCOND in the SLICOT routine MB03OD); a
+C (sub)matrix whose estimated condition number is less than
+C 1/TOL is considered to be of full rank. If the user sets
+C TOL <= 0, then an implicitly computed, default tolerance,
+C defined by TOLDEF = N*N*EPS, is used instead, where EPS
+C is the machine precision (see LAPACK Library routine
+C DLAMCH).
+C
+C Workspace
+C
+C IWORK INTEGER array, dimension (M)
+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 length of the array DWORK.
+C LDWORK >= MAX(1, N, 3*M, P).
+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
+C METHOD
+C
+C Matrix B is first QR-decomposed and the appropriate orthogonal
+C similarity transformation applied to the matrix A. Leaving the
+C first rank(B) states unchanged, the remaining lower left block
+C of A is then QR-decomposed and the new orthogonal matrix, Q1,
+C is also applied to the right of A to complete the similarity
+C transformation. By continuing in this manner, a completely
+C controllable state-space pair (Acont, Bcont) is found for the
+C given (A, B), where Acont is upper block Hessenberg with each
+C subdiagonal block of full row rank, and Bcont is zero apart from
+C its (independent) first rank(B) rows.
+C All orthogonal transformations determined in this process are also
+C applied to the matrix C, from the right.
+C NOTE that the system controllability indices are easily
+C calculated from the dimensions of the blocks of Acont.
+C
+C REFERENCES
+C
+C [1] Konstantinov, M.M., Petkov, P.Hr. and Christov, N.D.
+C Orthogonal Invariants and Canonical Forms for Linear
+C Controllable Systems.
+C Proc. 8th IFAC World Congress, Kyoto, 1, pp. 49-54, 1981.
+C
+C [2] Paige, C.C.
+C Properties of numerical algorithms related to computing
+C controllablity.
+C IEEE Trans. Auto. Contr., AC-26, pp. 130-138, 1981.
+C
+C [3] Petkov, P.Hr., Konstantinov, M.M., Gu, D.W. and
+C Postlethwaite, I.
+C Optimal Pole Assignment Design of Linear Multi-Input Systems.
+C Leicester University, Report 99-11, May 1996.
+C
+C NUMERICAL ASPECTS
+C 3
+C The algorithm requires 0(N ) operations and is backward stable.
+C
+C FURTHER COMMENTS
+C
+C If the system matrices A and B are badly scaled, it would be
+C useful to scale them with SLICOT routine TB01ID, before calling
+C the routine.
+C
+C CONTRIBUTOR
+C
+C V. Sima, Katholieke Univ. Leuven, Belgium, Feb. 1998.
+C
+C REVISIONS
+C
+C V. Sima, Katholieke Univ. Leuven, Belgium, May 1999, Nov. 2003.
+C A. Varga, DLR Oberpfaffenhofen, March 2002, Nov. 2003.
+C
+C KEYWORDS
+C
+C Controllability, minimal realization, orthogonal canonical form,
+C orthogonal transformation.
+C
+C ******************************************************************
+C
+C .. Parameters ..
+ DOUBLE PRECISION ZERO, ONE
+ PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
+C .. Scalar Arguments ..
+ CHARACTER JOBZ
+ INTEGER INDCON, INFO, LDA, LDB, LDC, LDWORK, LDZ, M, N,
+ $ NCONT, P
+ DOUBLE PRECISION TOL
+C .. Array Arguments ..
+ DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), DWORK(*), TAU(*),
+ $ Z(LDZ,*)
+ INTEGER IWORK(*), NBLK(*)
+C .. Local Scalars ..
+ LOGICAL LJOBF, LJOBI, LJOBZ
+ INTEGER IQR, ITAU, J, MCRT, NBL, NCRT, NI, NJ, RANK,
+ $ WRKOPT
+ DOUBLE PRECISION ANORM, BNORM, FNRM, TOLDEF
+C .. Local Arrays ..
+ DOUBLE PRECISION SVAL(3)
+C .. External Functions ..
+ LOGICAL LSAME
+ DOUBLE PRECISION DLAMCH, DLANGE, DLAPY2
+ EXTERNAL DLAMCH, DLANGE, DLAPY2, LSAME
+C .. External Subroutines ..
+ EXTERNAL DCOPY, DLACPY, DLAPMT, DLASET, DORGQR, DORMQR,
+ $ MB01PD, MB03OY, XERBLA
+C .. Intrinsic Functions ..
+ INTRINSIC DBLE, INT, MAX, MIN
+C ..
+C .. Executable Statements ..
+C
+ INFO = 0
+ LJOBF = LSAME( JOBZ, 'F' )
+ LJOBI = LSAME( JOBZ, 'I' )
+ LJOBZ = LJOBF.OR.LJOBI
+C
+C Test the input scalar arguments.
+C
+ IF( .NOT.LJOBZ .AND. .NOT.LSAME( JOBZ, 'N' ) ) THEN
+ INFO = -1
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -2
+ ELSE IF( M.LT.0 ) THEN
+ INFO = -3
+ ELSE IF( P.LT.0 ) THEN
+ INFO = -4
+ ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+ INFO = -6
+ ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
+ INFO = -8
+ ELSE IF( LDC.LT.MAX( 1, P ) ) THEN
+ INFO = -10
+ ELSE IF( .NOT.LJOBZ .AND. LDZ.LT.1 .OR.
+ $ LJOBZ .AND. LDZ.LT.MAX( 1, N ) ) THEN
+ INFO = -15
+ ELSE IF( LDWORK.LT.MAX( 1, N, 3*M, P ) ) THEN
+ INFO = -20
+ END IF
+C
+ IF ( INFO.NE.0 ) THEN
+C
+C Error return.
+C
+ CALL XERBLA( 'TB01UD', -INFO )
+ RETURN
+ END IF
+C
+ NCONT = 0
+ INDCON = 0
+C
+C Calculate the absolute norms of A and B (used for scaling).
+C
+ ANORM = DLANGE( 'M', N, N, A, LDA, DWORK )
+ BNORM = DLANGE( 'M', N, M, B, LDB, DWORK )
+C
+C Quick return if possible.
+C
+ IF ( MIN( N, M ).EQ.0 .OR. BNORM.EQ.ZERO ) THEN
+ IF( N.GT.0 ) THEN
+ IF ( LJOBI ) THEN
+ CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
+ ELSE IF ( LJOBF ) THEN
+ CALL DLASET( 'Full', N, N, ZERO, ZERO, Z, LDZ )
+ CALL DLASET( 'Full', N, 1, ZERO, ZERO, TAU, N )
+ END IF
+ END IF
+ DWORK(1) = ONE
+ RETURN
+ END IF
+C
+C Scale (if needed) the matrices A and B.
+C
+ CALL MB01PD( 'S', 'G', N, N, 0, 0, ANORM, 0, NBLK, A, LDA, INFO )
+ CALL MB01PD( 'S', 'G', N, M, 0, 0, BNORM, 0, NBLK, B, LDB, INFO )
+C
+C Compute the Frobenius norm of [ B A ] (used for rank estimation).
+C
+ FNRM = DLAPY2( DLANGE( 'F', N, M, B, LDB, DWORK ),
+ $ DLANGE( 'F', N, N, A, LDA, DWORK ) )
+C
+ TOLDEF = TOL
+ IF ( TOLDEF.LE.ZERO ) THEN
+C
+C Use the default tolerance in controllability determination.
+C
+ TOLDEF = DBLE( N*N )*DLAMCH( 'EPSILON' )
+ END IF
+C
+ IF ( FNRM.LT.TOLDEF )
+ $ FNRM = ONE
+C
+ WRKOPT = 1
+ NI = 0
+ ITAU = 1
+ NCRT = N
+ MCRT = M
+ IQR = 1
+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
+ 10 CONTINUE
+C
+C Rank-revealing QR decomposition with column pivoting.
+C The calculation is performed in NCRT rows of B starting from
+C the row IQR (initialized to 1 and then set to rank(B)+1).
+C Workspace: 3*MCRT.
+C
+ CALL MB03OY( NCRT, MCRT, B(IQR,1), LDB, TOLDEF, FNRM, RANK,
+ $ SVAL, IWORK, TAU(ITAU), DWORK, INFO )
+C
+ IF ( RANK.NE.0 ) THEN
+ NJ = NI
+ NI = NCONT
+ NCONT = NCONT + RANK
+ INDCON = INDCON + 1
+ NBLK(INDCON) = RANK
+C
+C Premultiply and postmultiply the appropriate block row
+C and block column of A by Q' and Q, respectively.
+C Workspace: need NCRT;
+C prefer NCRT*NB.
+C
+ CALL DORMQR( 'Left', 'Transpose', NCRT, NCRT, RANK,
+ $ B(IQR,1), LDB, TAU(ITAU), A(NI+1,NI+1), LDA,
+ $ DWORK, LDWORK, INFO )
+ WRKOPT = MAX( WRKOPT, INT( DWORK(1) ) )
+C
+C Workspace: need N;
+C prefer N*NB.
+C
+ CALL DORMQR( 'Right', 'No transpose', N, NCRT, RANK,
+ $ B(IQR,1), LDB, TAU(ITAU), A(1,NI+1), LDA,
+ $ DWORK, LDWORK, INFO )
+ WRKOPT = MAX( WRKOPT, INT( DWORK(1) ) )
+C
+C Postmultiply the appropriate block column of C by Q.
+C Workspace: need P;
+C prefer P*NB.
+C
+ CALL DORMQR( 'Right', 'No transpose', P, NCRT, RANK,
+ $ B(IQR,1), LDB, TAU(ITAU), C(1,NI+1), LDC,
+ $ DWORK, LDWORK, INFO )
+ WRKOPT = MAX( WRKOPT, INT( DWORK(1) ) )
+C
+C If required, save transformations.
+C
+ IF ( LJOBZ.AND.NCRT.GT.1 ) THEN
+ CALL DLACPY( 'L', NCRT-1, MIN( RANK, NCRT-1 ),
+ $ B(IQR+1,1), LDB, Z(NI+2,ITAU), LDZ )
+ END IF
+C
+C Zero the subdiagonal elements of the current matrix.
+C
+ IF ( RANK.GT.1 )
+ $ CALL DLASET( 'L', RANK-1, RANK-1, ZERO, ZERO, B(IQR+1,1),
+ $ LDB )
+C
+C Backward permutation of the columns of B or A.
+C
+ IF ( INDCON.EQ.1 ) THEN
+ CALL DLAPMT( .FALSE., RANK, M, B(IQR,1), LDB, IWORK )
+ IQR = RANK + 1
+ ELSE
+ DO 20 J = 1, MCRT
+ CALL DCOPY( RANK, B(IQR,J), 1, A(NI+1,NJ+IWORK(J)),
+ $ 1 )
+ 20 CONTINUE
+ END IF
+C
+ ITAU = ITAU + RANK
+ IF ( RANK.NE.NCRT ) THEN
+ MCRT = RANK
+ NCRT = NCRT - RANK
+ CALL DLACPY( 'G', NCRT, MCRT, A(NCONT+1,NI+1), LDA,
+ $ B(IQR,1), LDB )
+ CALL DLASET( 'G', NCRT, MCRT, ZERO, ZERO,
+ $ A(NCONT+1,NI+1), LDA )
+ GO TO 10
+ END IF
+ END IF
+C
+C If required, accumulate transformations.
+C Workspace: need N; prefer N*NB.
+C
+ IF ( LJOBI ) THEN
+ CALL DORGQR( N, N, ITAU-1, Z, LDZ, TAU, DWORK,
+ $ LDWORK, INFO )
+ WRKOPT = MAX( WRKOPT, INT( DWORK(1) ) )
+ END IF
+C
+C Annihilate the trailing blocks of B.
+C
+ IF( IQR.LE.N )
+ $ CALL DLASET( 'G', N-IQR+1, M, ZERO, ZERO, B(IQR,1), LDB )
+C
+C Annihilate the trailing elements of TAU, if JOBZ = 'F'.
+C
+ IF ( LJOBF ) THEN
+ DO 30 J = ITAU, N
+ TAU(J) = ZERO
+ 30 CONTINUE
+ END IF
+C
+C Undo scaling of A and B.
+C
+ IF ( INDCON.LT.N ) THEN
+ NBL = INDCON + 1
+ NBLK(NBL) = N - NCONT
+ ELSE
+ NBL = 0
+ END IF
+ CALL MB01PD( 'U', 'H', N, N, 0, 0, ANORM, NBL, NBLK, A, LDA,
+ $ INFO )
+ CALL MB01PD( 'U', 'G', NBLK(1), M, 0, 0, BNORM, 0, NBLK, B, LDB,
+ $ INFO )
+C
+C Set optimal workspace dimension.
+C
+ DWORK(1) = WRKOPT
+ RETURN
+C *** Last line of TB01UD ***
+ END
Added: trunk/octave-forge/main/control/src/TB01XD.f
===================================================================
--- trunk/octave-forge/main/control/src/TB01XD.f (rev 0)
+++ trunk/octave-forge/main/control/src/TB01XD.f 2010-09-07 08:40:56 UTC (rev 7673)
@@ -0,0 +1,284 @@
+ SUBROUTINE TB01XD( JOBD, N, M, P, KL, KU, A, LDA, B, LDB, C, LDC,
+ $ D, LDD, 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 apply a special transformation to a system given as a triple
+C (A,B,C),
+C
+C A <-- P * A' * P, B <-- P * C', C <-- B' * P,
+C
+C where P is a matrix with 1 on the secondary diagonal, and with 0
+C in the other entries. Matrix A can be specified as a band matrix.
+C Optionally, matrix D of the system can be transposed. This
+C transformation is actually a special similarity transformation of
+C the dual system.
+C
+C ARGUMENTS
+C
+C Mode Parameters
+C
+C JOBD CHARACTER*1
+C Specifies whether or not a non-zero matrix D appears in
+C the given state space model:
+C = 'D': D is present;
+C = 'Z': D is assumed a zero matrix.
+C
+C Input/Output Parameters
+C
+C N (input) INTEGER
+C The order of the matrix A, the number of rows of matrix B
+C and the number of columns of matrix C.
+C N represents the dimension of the state vector. N >= 0.
+C
+C M (input) INTEGER.
+C The number of columns of matrix B.
+C M represents the dimension of input vector. M >= 0.
+C
+C P (input) INTEGER.
+C The number of rows of matrix C.
+C P represents the dimension of output vector. P >= 0.
+C
+C KL (input) INTEGER
+C The number of subdiagonals of A to be transformed.
+C MAX( 0, N-1 ) >= KL >= 0.
+C
+C KU (input) INTEGER
+C The number of superdiagonals of A to be transformed.
+C MAX( 0, N-1 ) >= KU >= 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 system state matrix A.
+C On exit, the leading N-by-N part of this array contains
+C the transformed (pertransposed) matrix P*A'*P.
+C
+C LDA INTEGER
+C The leading dimension of the array A. LDA >= MAX(1,N).
+C
+C B (input/output) DOUBLE PRECISION array, dimension
+C (LDB,MAX(M,P))
+C On entry, the leading N-by-M part of this array must
+C contain the original input/state matrix B.
+C On exit, the leading N-by-P part of this array contains
+C the dual input/state matrix P*C'.
+C
+C LDB INTEGER
+C The leading dimension of the array B.
+C LDB >= MAX(1,N) if M > 0 or P > 0.
+C LDB >= 1 if M = 0 and P = 0.
+C
+C C (input/output) DOUBLE PRECISION array, dimension (LDC,N)
+C On entry, the leading P-by-N part of this array must
+C contain the original state/output matrix C.
+C On exit, the leading M-by-N part of this array contains
+C the dual state/output matrix B'*P.
+C
+C LDC INTEGER
+C The leading dimension of array C.
+C LDC >= MAX(1,M,P) if N > 0.
+C LDC >= 1 if N = 0.
+C
+C D (input/output) DOUBLE PRECISION array, dimension
+C (LDD,MAX(M,P))
+C On entry, if JOBD = 'D', the leading P-by-M part of this
+C array must contain the original direct transmission
+C matrix D.
+C On exit, if JOBD = 'D', the leading M-by-P part of this
+C array contains the transposed direct transmission matrix
+C D'. The array D is not referenced if JOBD = 'Z'.
+C
+C LDD INTEGER
+C The leading dimension of array D.
+C LDD >= MAX(1,M,P) if JOBD = 'D'.
+C LDD >= 1 if JOBD = 'Z'.
+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
+C METHOD
+C
+C The rows and/or columns of the matrices of the triplet (A,B,C)
+C and, optionally, of the matrix D are swapped in a special way.
+C
+C NUMERICAL ASPECTS
+C
+C None.
+C
+C CONTRIBUTOR
+C
+C V. Sima, Katholieke Univ. Leuven, Belgium, Feb. 1998.
+C Partly based on routine DMPTR (A. Varga, German Aerospace
+C Research Establishment, DLR, Aug. 1992).
+C
+C
+C REVISIONS
+C
+C 07-31-1998, 04-25-1999, A. Varga.
+C 03-16-2004, V. Sima.
+C
+C KEYWORDS
+C
+C Matrix algebra, matrix operations, similarity transformation.
+C
+C *********************************************************************
+C
+C ..
+C .. Scalar Arguments ..
+ CHARACTER JOBD
+ INTEGER INFO, KL, KU, LDA, LDB, LDC, LDD, M, N, P
+C ..
+C .. Array Arguments ..
+ DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ),
+ $ D( LDD, * )
+C ..
+C .. Local Scalars ..
+ LOGICAL LJOBD
+ INTEGER J, J1, LDA1, MAXMP, MINMP, NM1
+C ..
+C .. External functions ..
+ LOGICAL LSAME
+ EXTERNAL LSAME
+C ..
+C .. External Subroutines ..
+ EXTERNAL DCOPY, DSWAP, XERBLA
+C ..
+C .. Intrinsic Functions ..
+ INTRINSIC MAX, MIN
+C ..
+C .. Executable Statements ..
+C
+C Test the scalar input arguments.
+C
+ INFO = 0
+ LJOBD = LSAME( JOBD, 'D' )
+ MAXMP = MAX( M, P )
+ MINMP = MIN( M, P )
+ NM1 = N - 1
+C
+ IF( .NOT.LJOBD .AND. .NOT.LSAME( JOBD, 'Z' ) ) THEN
+ INFO = -1
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -2
+ ELSE IF( M.LT.0 ) THEN
+ INFO = -3
+ ELSE IF( P.LT.0 ) THEN
+ INFO = -4
+ ELSE IF( KL.LT.0 .OR. KL.GT.MAX( 0, NM1 ) ) THEN
+ INFO = -5
+ ELSE IF( KU.LT.0 .OR. KU.GT.MAX( 0, NM1 ) ) THEN
+ INFO = -6
+ ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+ INFO = -8
+ ELSE IF( ( MAXMP.GT.0 .AND. LDB.LT.MAX( 1, N ) ) .OR.
+ $ ( MINMP.EQ.0 .AND. LDB.LT.1 ) ) THEN
+ INFO = -10
+ ELSE IF( LDC.LT.1 .OR. ( N.GT.0 .AND. LDC.LT.MAXMP ) ) THEN
+ INFO = -12
+ ELSE IF( LDD.LT.1 .OR. ( LJOBD .AND. LDD.LT.MAXMP ) ) THEN
+ INFO = -14
+ END IF
+C
+ IF( INFO.NE.0 ) THEN
+C
+C Error return.
+C
+ CALL XERBLA( 'TB01XD', -INFO )
+ RETURN
+ END IF
+C
+C Quick return if possible.
+C
+ IF ( LJOBD ) THEN
+C
+C Replace D by D', if non-scalar.
+C
+ DO 5 J = 1, MAXMP
+ IF ( J.LT.MINMP ) THEN
+ CALL DSWAP( MINMP-J, D(J+1,J), 1, D(J,J+1), LDD )
+ ELSE IF ( J.GT.P ) THEN
+ CALL DCOPY( P, D(1,J), 1, D(J,1), LDD )
+ ELSE IF ( J.GT.M ) THEN
+ CALL DCOPY( M, D(J,1), LDD, D(1,J), 1 )
+ END IF
+ 5 CONTINUE
+C
+ END IF
+C
+ IF( N.EQ.0 )
+ $ RETURN
+C
+C Replace matrix A by P*A'*P.
+C
+ IF ( KL.EQ.NM1 .AND. KU.EQ.NM1 ) THEN
+C
+C Full matrix A.
+C
+ DO 10 J = 1, NM1
+ CALL DSWAP( N-J, A( 1, J ), 1, A( N-J+1, J+1 ), -LDA )
+ 10 CONTINUE
+C
+ ELSE
+C
+C Band matrix A.
+C
+ LDA1 = LDA + 1
+C
+C Pertranspose the KL subdiagonals.
+C
+ DO 20 J = 1, MIN( KL, N-2 )
+ J1 = ( N - J )/2
+ CALL DSWAP( J1, A(J+1,1), LDA1, A(N-J1+1,N-J1+1-J), -LDA1 )
+ 20 CONTINUE
+C
+C Pertranspose the KU superdiagonals.
+C
+ DO 30 J = 1, MIN( KU, N-2 )
+ J1 = ( N - J )/2
+ CALL DSWAP( J1, A(1,J+1), LDA1, A(N-J1+1-J,N-J1+1), -LDA1 )
+ 30 CONTINUE
+C
+C Pertranspose the diagonal.
+C
+ J1 = N/2
+ CALL DSWAP( J1, A(1,1), LDA1, A(N-J1+1,N-J1+1), -LDA1 )
+C
+ END IF
+C
+C Replace matrix B by P*C' and matrix C by B'*P.
+C
+ DO 40 J = 1, MAXMP
+ IF ( J.LE.MINMP ) THEN
+ CALL DSWAP( N, B(1,J), 1, C(J,1), -LDC )
+ ELSE IF ( J.GT.P ) THEN
+ CALL DCOPY( N, B(1,J), 1, C(J,1), -LDC )
+ ELSE
+ CALL DCOPY( N, C(J,1), -LDC, B(1,J), 1 )
+ END IF
+ 40 CONTINUE
+C
+ RETURN
+C *** Last line of TB01XD ***
+ END
Added: trunk/octave-forge/main/control/src/sltb01pd.cc
===================================================================
--- trunk/octave-forge/main/control/src/sltb01pd.cc (rev 0)
+++ trunk/octave-forge/main/control/src/sltb01pd.cc 2010-09-07 08:40:56 UTC (rev 7673)
@@ -0,0 +1,123 @@
+/*
+
+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/>.
+
+Minimal realization of state-space models.
+Uses SLICOT TB01PD by courtesy of NICONET e.V.
+<http://www.slicot.org>
+
+Author: Lukas Reichlin <luk...@gm...>
+Created: September 2010
+Version: 0.1
+
+*/
+
+#include <octave/oct.h>
+#include <f77-fcn.h>
+#include "common.cc"
+
+extern "C"
+{
+ int F77_FUNC (tb01pd, TB01PD)
+ (char& JOB, char& EQUIL,
+ int& N, int& M, int& P,
+ double* A, int& LDA,
+ double* B, int& LDB,
+ double* C, int& LDC,
+ int& NR,
+ double& TOL,
+ int* IWORK,
+ double* DWORK, int& LDWORK,
+ int& INFO);
+}
+
+DEFUN_DLD (sltb01pd, args, nargout, "Slicot TB01PD Release 5.0")
+{
+ int nargin = args.length ();
+ octave_value_list retval;
+
+ if (nargin != 4)
+ {
+ print_usage ();
+ }
+ else
+ {
+ // arguments in
+ char job = 'M';
+ char equil = 'N';
+
+ Matrix a = args(0).matrix_value ();
+ Matrix b = args(1).matrix_value ();
+ Matrix c = args(2).matrix_value ();
+ double tol = args(3).double_value ();
+
+ int n = a.rows (); // n: number of states
+ int m = b.columns (); // m: number of inputs
+ int p = c.rows (); // p: number of outputs
+
+ int lda = max (1, n);
+ int ldb = max (1, n);
+ int ldc;
+
+ if (n == 0)
+ ldc = 1;
+ else
+ ldc = max (1, m, p);
+
+ // arguments out
+ int nr;
+
+ // workspace
+ int liwork = n + max (m, p);
+ int ldwork = max (1, n + max (n, 3*m, 3*p));
+
+ OCTAVE_LOCAL_BUFFER (int, iwork, liwork);
+ OCTAVE_LOCAL_BUFFER (double, dwork, ldwork);
+
+ // error indicators
+ int info;
+
+
+ // SLICOT routine TB01PD
+ F77_XFCN (tb01pd, TB01PD,
+ (job, equil,
+ n, m, p,
+ a.fortran_vec (), lda,
+ b.fortran_vec (), ldb,
+ c.fortran_vec (), ldc,
+ nr,
+ tol,
+ iwork,
+ dwork, ldwork,
+ info));
+
+ if (f77_exception_encountered)
+ error ("ss: minreal: sltb01pd: exception in SLICOT subroutine TB01PD");
+
+ if (info != 0)
+ error ("ss: minreal: sltb01pd: TB01PD returned info = %d", info);
+
+ // return values
+ retval(0) = a;
+ retval(1) = b;
+ retval(2) = c;
+ retval(3) = octave_value (nr);
+ }
+
+ return retval;
+}
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
