From: <par...@us...> - 2011-08-22 18:25:46
|
Revision: 8474 http://octave.svn.sourceforge.net/octave/?rev=8474&view=rev Author: paramaniac Date: 2011-08-22 18:25:37 +0000 (Mon, 22 Aug 2011) Log Message: ----------- control: add code for Slicot TD04AD tf2ss conversion Added Paths: ----------- trunk/octave-forge/main/control/devel/tf2ss/ trunk/octave-forge/main/control/devel/tf2ss/AB07MD.f trunk/octave-forge/main/control/devel/tf2ss/MB01PD.f trunk/octave-forge/main/control/devel/tf2ss/MB01QD.f trunk/octave-forge/main/control/devel/tf2ss/MB03OY.f trunk/octave-forge/main/control/devel/tf2ss/TB01ID.f trunk/octave-forge/main/control/devel/tf2ss/TB01PD.f trunk/octave-forge/main/control/devel/tf2ss/TB01UD.f trunk/octave-forge/main/control/devel/tf2ss/TB01XD.f trunk/octave-forge/main/control/devel/tf2ss/TD03AY.f trunk/octave-forge/main/control/devel/tf2ss/TD04AD.dat trunk/octave-forge/main/control/devel/tf2ss/TD04AD.f trunk/octave-forge/main/control/devel/tf2ss/TD04AD.html trunk/octave-forge/main/control/devel/tf2ss/TD04AD.res trunk/octave-forge/main/control/devel/tf2ss/makefile_tf2ss.m Added: trunk/octave-forge/main/control/devel/tf2ss/AB07MD.f =================================================================== --- trunk/octave-forge/main/control/devel/tf2ss/AB07MD.f (rev 0) +++ trunk/octave-forge/main/control/devel/tf2ss/AB07MD.f 2011-08-22 18:25:37 UTC (rev 8474) @@ -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-2010 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 Property changes on: trunk/octave-forge/main/control/devel/tf2ss/AB07MD.f ___________________________________________________________________ Added: svn:executable + * Added: trunk/octave-forge/main/control/devel/tf2ss/MB01PD.f =================================================================== --- trunk/octave-forge/main/control/devel/tf2ss/MB01PD.f (rev 0) +++ trunk/octave-forge/main/control/devel/tf2ss/MB01PD.f 2011-08-22 18:25:37 UTC (rev 8474) @@ -0,0 +1,271 @@ + SUBROUTINE MB01PD( SCUN, TYPE, M, N, KL, KU, ANRM, NBL, NROWS, A, + $ LDA, INFO ) +C +C SLICOT RELEASE 5.0. +C +C Copyright (c) 2002-2010 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 scale a matrix or undo scaling. Scaling is performed, if +C necessary, so that the matrix norm will be in a safe range of +C representable numbers. +C +C ARGUMENTS +C +C Mode Parameters +C +C SCUN CHARACTER*1 +C SCUN indicates the operation to be performed. +C = 'S': scale the matrix. +C = 'U': undo scaling of the matrix. +C +C TYPE CHARACTER*1 +C TYPE indicates the storage type of the input matrix. +C = 'G': A is a full matrix. +C = 'L': A is a (block) lower triangular matrix. +C = 'U': A is an (block) upper triangular matrix. +C = 'H': A is an (block) upper Hessenberg matrix. +C = 'B': A is a symmetric band matrix with lower bandwidth +C KL and upper bandwidth KU and with the only the +C lower half stored. +C = 'Q': A is a symmetric band matrix with lower bandwidth +C KL and upper bandwidth KU and with the only the +C upper half stored. +C = 'Z': A is a band matrix with lower bandwidth KL and +C upper bandwidth KU. +C +C Input/Output Parameters +C +C M (input) INTEGER +C The number of rows of the matrix A. M >= 0. +C +C N (input) INTEGER +C The number of columns of the matrix A. N >= 0. +C +C KL (input) INTEGER +C The lower bandwidth of A. Referenced only if TYPE = 'B', +C 'Q' or 'Z'. +C +C KU (input) INTEGER +C The upper bandwidth of A. Referenced only if TYPE = 'B', +C 'Q' or 'Z'. +C +C ANRM (input) DOUBLE PRECISION +C The norm of the initial matrix A. ANRM >= 0. +C When ANRM = 0 then an immediate return is effected. +C ANRM should be preserved between the call of the routine +C with SCUN = 'S' and the corresponding one with SCUN = 'U'. +C +C NBL (input) INTEGER +C The number of diagonal blocks of the matrix A, if it has a +C block structure. To specify that matrix A has no block +C structure, set NBL = 0. NBL >= 0. +C +C NROWS (input) INTEGER array, dimension max(1,NBL) +C NROWS(i) contains the number of rows and columns of the +C i-th diagonal block of matrix A. The sum of the values +C NROWS(i), for i = 1: NBL, should be equal to min(M,N). +C The elements of the array NROWS are not referenced if +C NBL = 0. +C +C A (input/output) DOUBLE PRECISION array, dimension (LDA,N) +C On entry, the leading M by N part of this array must +C contain the matrix to be scaled/unscaled. +C On exit, the leading M by N part of A will contain +C the modified matrix. +C The storage mode of A is specified by TYPE. +C +C LDA (input) INTEGER +C The leading dimension of the array A. LDA >= max(1,M). +C +C Error Indicator +C +C INFO (output) INTEGER +C = 0: successful exit +C < 0: if INFO = -i, the i-th argument had an illegal +C value. +C +C METHOD +C +C Denote by ANRM the norm of the matrix, and by SMLNUM and BIGNUM, +C two positive numbers near the smallest and largest safely +C representable numbers, respectively. The matrix is scaled, if +C needed, such that the norm of the result is in the range +C [SMLNUM, BIGNUM]. The scaling factor is represented as a ratio +C of two numbers, one of them being ANRM, and the other one either +C SMLNUM or BIGNUM, depending on ANRM being less than SMLNUM or +C larger than BIGNUM, respectively. For undoing the scaling, the +C norm is again compared with SMLNUM or BIGNUM, and the reciprocal +C of the previous scaling factor is used. +C +C CONTRIBUTOR +C +C V. Sima, Katholieke Univ. Leuven, Belgium, Nov. 1996. +C +C REVISIONS +C +C Oct. 2001, V. Sima, Research Institute for Informatics, Bucharest. +C +C ****************************************************************** +C +C .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) +C .. Scalar Arguments .. + CHARACTER SCUN, TYPE + INTEGER INFO, KL, KU, LDA, M, MN, N, NBL + DOUBLE PRECISION ANRM +C .. Array Arguments .. + INTEGER NROWS ( * ) + DOUBLE PRECISION A( LDA, * ) +C .. Local Scalars .. + LOGICAL FIRST, LSCALE + INTEGER I, ISUM, ITYPE + DOUBLE PRECISION BIGNUM, SMLNUM +C .. External Functions .. + LOGICAL LSAME + DOUBLE PRECISION DLAMCH + EXTERNAL DLAMCH, LSAME +C .. +C .. External Subroutines .. + EXTERNAL DLABAD, MB01QD, XERBLA +C .. Intrinsic Functions .. + INTRINSIC MAX, MIN +C .. Save statement .. + SAVE BIGNUM, FIRST, SMLNUM +C .. Data statements .. + DATA FIRST/.TRUE./ +C .. +C .. Executable Statements .. +C +C Test the input scalar arguments. +C + INFO = 0 + LSCALE = LSAME( SCUN, 'S' ) + IF( LSAME( TYPE, 'G' ) ) THEN + ITYPE = 0 + ELSE IF( LSAME( TYPE, 'L' ) ) THEN + ITYPE = 1 + ELSE IF( LSAME( TYPE, 'U' ) ) THEN + ITYPE = 2 + ELSE IF( LSAME( TYPE, 'H' ) ) THEN + ITYPE = 3 + ELSE IF( LSAME( TYPE, 'B' ) ) THEN + ITYPE = 4 + ELSE IF( LSAME( TYPE, 'Q' ) ) THEN + ITYPE = 5 + ELSE IF( LSAME( TYPE, 'Z' ) ) THEN + ITYPE = 6 + ELSE + ITYPE = -1 + END IF +C + MN = MIN( M, N ) +C + ISUM = 0 + IF( NBL.GT.0 ) THEN + DO 10 I = 1, NBL + ISUM = ISUM + NROWS(I) + 10 CONTINUE + END IF +C + IF( .NOT.LSCALE .AND. .NOT.LSAME( SCUN, 'U' ) ) THEN + INFO = -1 + ELSE IF( ITYPE.EQ.-1 ) THEN + INFO = -2 + ELSE IF( M.LT.0 ) THEN + INFO = -3 + ELSE IF( N.LT.0 .OR. + $ ( ( ITYPE.EQ.4 .OR. ITYPE.EQ.5 ) .AND. N.NE.M ) ) THEN + INFO = -4 + ELSE IF( ANRM.LT.ZERO ) THEN + INFO = -7 + ELSE IF( NBL.LT.0 ) THEN + INFO = -8 + ELSE IF( NBL.GT.0 .AND. ISUM.NE.MN ) THEN + INFO = -9 + ELSE IF( ITYPE.LE.3 .AND. LDA.LT.MAX( 1, M ) ) THEN + INFO = -11 + ELSE IF( ITYPE.GE.4 ) THEN + IF( KL.LT.0 .OR. KL.GT.MAX( M-1, 0 ) ) THEN + INFO = -5 + ELSE IF( KU.LT.0 .OR. KU.GT.MAX( N-1, 0 ) .OR. + $ ( ( ITYPE.EQ.4 .OR. ITYPE.EQ.5 ) .AND. KL.NE.KU ) ) + $ THEN + INFO = -6 + ELSE IF( ( ITYPE.EQ.4 .AND. LDA.LT.KL+1 ) .OR. + $ ( ITYPE.EQ.5 .AND. LDA.LT.KU+1 ) .OR. + $ ( ITYPE.EQ.6 .AND. LDA.LT.2*KL+KU+1 ) ) THEN + INFO = -11 + END IF + END IF +C + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'MB01PD', -INFO ) + RETURN + END IF +C +C Quick return if possible. +C + IF( MN.EQ.0 .OR. ANRM.EQ.ZERO ) + $ RETURN +C + IF ( FIRST ) THEN +C +C Get machine parameters. +C + SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' ) + BIGNUM = ONE / SMLNUM + CALL DLABAD( SMLNUM, BIGNUM ) + FIRST = .FALSE. + END IF +C + IF ( LSCALE ) THEN +C +C Scale A, if its norm is outside range [SMLNUM,BIGNUM]. +C + IF( ANRM.LT.SMLNUM ) THEN +C +C Scale matrix norm up to SMLNUM. +C + CALL MB01QD( TYPE, M, N, KL, KU, ANRM, SMLNUM, NBL, NROWS, + $ A, LDA, INFO ) + ELSE IF( ANRM.GT.BIGNUM ) THEN +C +C Scale matrix norm down to BIGNUM. +C + CALL MB01QD( TYPE, M, N, KL, KU, ANRM, BIGNUM, NBL, NROWS, + $ A, LDA, INFO ) + END IF +C + ELSE +C +C Undo scaling. +C + IF( ANRM.LT.SMLNUM ) THEN + CALL MB01QD( TYPE, M, N, KL, KU, SMLNUM, ANRM, NBL, NROWS, + $ A, LDA, INFO ) + ELSE IF( ANRM.GT.BIGNUM ) THEN + CALL MB01QD( TYPE, M, N, KL, KU, BIGNUM, ANRM, NBL, NROWS, + $ A, LDA, INFO ) + END IF + END IF +C + RETURN +C *** Last line of MB01PD *** + END Property changes on: trunk/octave-forge/main/control/devel/tf2ss/MB01PD.f ___________________________________________________________________ Added: svn:executable + * Added: trunk/octave-forge/main/control/devel/tf2ss/MB01QD.f =================================================================== --- trunk/octave-forge/main/control/devel/tf2ss/MB01QD.f (rev 0) +++ trunk/octave-forge/main/control/devel/tf2ss/MB01QD.f 2011-08-22 18:25:37 UTC (rev 8474) @@ -0,0 +1,334 @@ + SUBROUTINE MB01QD( TYPE, M, N, KL, KU, CFROM, CTO, NBL, NROWS, A, + $ LDA, INFO ) +C +C SLICOT RELEASE 5.0. +C +C Copyright (c) 2002-2010 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 multiply the M by N real matrix A by the real scalar CTO/CFROM. +C This is done without over/underflow as long as the final result +C CTO*A(I,J)/CFROM does not over/underflow. TYPE specifies that +C A may be full, (block) upper triangular, (block) lower triangular, +C (block) upper Hessenberg, or banded. +C +C ARGUMENTS +C +C Mode Parameters +C +C TYPE CHARACTER*1 +C TYPE indices the storage type of the input matrix. +C = 'G': A is a full matrix. +C = 'L': A is a (block) lower triangular matrix. +C = 'U': A is a (block) upper triangular matrix. +C = 'H': A is a (block) upper Hessenberg matrix. +C = 'B': A is a symmetric band matrix with lower bandwidth +C KL and upper bandwidth KU and with the only the +C lower half stored. +C = 'Q': A is a symmetric band matrix with lower bandwidth +C KL and upper bandwidth KU and with the only the +C upper half stored. +C = 'Z': A is a band matrix with lower bandwidth KL and +C upper bandwidth KU. +C +C Input/Output Parameters +C +C M (input) INTEGER +C The number of rows of the matrix A. M >= 0. +C +C N (input) INTEGER +C The number of columns of the matrix A. N >= 0. +C +C KL (input) INTEGER +C The lower bandwidth of A. Referenced only if TYPE = 'B', +C 'Q' or 'Z'. +C +C KU (input) INTEGER +C The upper bandwidth of A. Referenced only if TYPE = 'B', +C 'Q' or 'Z'. +C +C CFROM (input) DOUBLE PRECISION +C CTO (input) DOUBLE PRECISION +C The matrix A is multiplied by CTO/CFROM. A(I,J) is +C computed without over/underflow if the final result +C CTO*A(I,J)/CFROM can be represented without over/ +C underflow. CFROM must be nonzero. +C +C NBL (input) INTEGER +C The number of diagonal blocks of the matrix A, if it has a +C block structure. To specify that matrix A has no block +C structure, set NBL = 0. NBL >= 0. +C +C NROWS (input) INTEGER array, dimension max(1,NBL) +C NROWS(i) contains the number of rows and columns of the +C i-th diagonal block of matrix A. The sum of the values +C NROWS(i), for i = 1: NBL, should be equal to min(M,N). +C The array NROWS is not referenced if NBL = 0. +C +C A (input/output) DOUBLE PRECISION array, dimension (LDA,N) +C The matrix to be multiplied by CTO/CFROM. See TYPE for +C the storage type. +C +C LDA (input) INTEGER +C The leading dimension of the array A. LDA >= max(1,M). +C +C Error Indicator +C +C INFO INTEGER +C Not used in this implementation. +C +C METHOD +C +C Matrix A is multiplied by the real scalar CTO/CFROM, taking into +C account the specified storage mode of the matrix. +C MB01QD is a version of the LAPACK routine DLASCL, modified for +C dealing with block triangular, or block Hessenberg matrices. +C For efficiency, no tests of the input scalar parameters are +C performed. +C +C CONTRIBUTOR +C +C V. Sima, Katholieke Univ. Leuven, Belgium, Nov. 1996. +C +C ****************************************************************** +C +C .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) +C .. +C .. Scalar Arguments .. + CHARACTER TYPE + INTEGER INFO, KL, KU, LDA, M, N, NBL + DOUBLE PRECISION CFROM, CTO +C .. +C .. Array Arguments .. + INTEGER NROWS ( * ) + DOUBLE PRECISION A( LDA, * ) +C .. +C .. Local Scalars .. + LOGICAL DONE, NOBLC + INTEGER I, IFIN, ITYPE, J, JFIN, JINI, K, K1, K2, K3, + $ K4 + DOUBLE PRECISION BIGNUM, CFROM1, CFROMC, CTO1, CTOC, MUL, SMLNUM +C .. +C .. External Functions .. + LOGICAL LSAME + DOUBLE PRECISION DLAMCH + EXTERNAL LSAME, DLAMCH +C .. +C .. Intrinsic Functions .. + INTRINSIC ABS, MAX, MIN +C .. +C .. Executable Statements .. +C + IF( LSAME( TYPE, 'G' ) ) THEN + ITYPE = 0 + ELSE IF( LSAME( TYPE, 'L' ) ) THEN + ITYPE = 1 + ELSE IF( LSAME( TYPE, 'U' ) ) THEN + ITYPE = 2 + ELSE IF( LSAME( TYPE, 'H' ) ) THEN + ITYPE = 3 + ELSE IF( LSAME( TYPE, 'B' ) ) THEN + ITYPE = 4 + ELSE IF( LSAME( TYPE, 'Q' ) ) THEN + ITYPE = 5 + ELSE + ITYPE = 6 + END IF +C +C Quick return if possible. +C + IF( MIN( M, N ).EQ.0 ) + $ RETURN +C +C Get machine parameters. +C + SMLNUM = DLAMCH( 'S' ) + BIGNUM = ONE / SMLNUM +C + CFROMC = CFROM + CTOC = CTO +C + 10 CONTINUE + CFROM1 = CFROMC*SMLNUM + CTO1 = CTOC / BIGNUM + IF( ABS( CFROM1 ).GT.ABS( CTOC ) .AND. CTOC.NE.ZERO ) THEN + MUL = SMLNUM + DONE = .FALSE. + CFROMC = CFROM1 + ELSE IF( ABS( CTO1 ).GT.ABS( CFROMC ) ) THEN + MUL = BIGNUM + DONE = .FALSE. + CTOC = CTO1 + ELSE + MUL = CTOC / CFROMC + DONE = .TRUE. + END IF +C + NOBLC = NBL.EQ.0 +C + IF( ITYPE.EQ.0 ) THEN +C +C Full matrix +C + DO 30 J = 1, N + DO 20 I = 1, M + A( I, J ) = A( I, J )*MUL + 20 CONTINUE + 30 CONTINUE +C + ELSE IF( ITYPE.EQ.1 ) THEN +C + IF ( NOBLC ) THEN +C +C Lower triangular matrix +C + DO 50 J = 1, N + DO 40 I = J, M + A( I, J ) = A( I, J )*MUL + 40 CONTINUE + 50 CONTINUE +C + ELSE +C +C Block lower triangular matrix +C + JFIN = 0 + DO 80 K = 1, NBL + JINI = JFIN + 1 + JFIN = JFIN + NROWS( K ) + DO 70 J = JINI, JFIN + DO 60 I = JINI, M + A( I, J ) = A( I, J )*MUL + 60 CONTINUE + 70 CONTINUE + 80 CONTINUE + END IF +C + ELSE IF( ITYPE.EQ.2 ) THEN +C + IF ( NOBLC ) THEN +C +C Upper triangular matrix +C + DO 100 J = 1, N + DO 90 I = 1, MIN( J, M ) + A( I, J ) = A( I, J )*MUL + 90 CONTINUE + 100 CONTINUE +C + ELSE +C +C Block upper triangular matrix +C + JFIN = 0 + DO 130 K = 1, NBL + JINI = JFIN + 1 + JFIN = JFIN + NROWS( K ) + IF ( K.EQ.NBL ) JFIN = N + DO 120 J = JINI, JFIN + DO 110 I = 1, MIN( JFIN, M ) + A( I, J ) = A( I, J )*MUL + 110 CONTINUE + 120 CONTINUE + 130 CONTINUE + END IF +C + ELSE IF( ITYPE.EQ.3 ) THEN +C + IF ( NOBLC ) THEN +C +C Upper Hessenberg matrix +C + DO 150 J = 1, N + DO 140 I = 1, MIN( J+1, M ) + A( I, J ) = A( I, J )*MUL + 140 CONTINUE + 150 CONTINUE +C + ELSE +C +C Block upper Hessenberg matrix +C + JFIN = 0 + DO 180 K = 1, NBL + JINI = JFIN + 1 + JFIN = JFIN + NROWS( K ) +C + IF ( K.EQ.NBL ) THEN + JFIN = N + IFIN = N + ELSE + IFIN = JFIN + NROWS( K+1 ) + END IF +C + DO 170 J = JINI, JFIN + DO 160 I = 1, MIN( IFIN, M ) + A( I, J ) = A( I, J )*MUL + 160 CONTINUE + 170 CONTINUE + 180 CONTINUE + END IF +C + ELSE IF( ITYPE.EQ.4 ) THEN +C +C Lower half of a symmetric band matrix +C + K3 = KL + 1 + K4 = N + 1 + DO 200 J = 1, N + DO 190 I = 1, MIN( K3, K4-J ) + A( I, J ) = A( I, J )*MUL + 190 CONTINUE + 200 CONTINUE +C + ELSE IF( ITYPE.EQ.5 ) THEN +C +C Upper half of a symmetric band matrix +C + K1 = KU + 2 + K3 = KU + 1 + DO 220 J = 1, N + DO 210 I = MAX( K1-J, 1 ), K3 + A( I, J ) = A( I, J )*MUL + 210 CONTINUE + 220 CONTINUE +C + ELSE IF( ITYPE.EQ.6 ) THEN +C +C Band matrix +C + K1 = KL + KU + 2 + K2 = KL + 1 + K3 = 2*KL + KU + 1 + K4 = KL + KU + 1 + M + DO 240 J = 1, N + DO 230 I = MAX( K1-J, K2 ), MIN( K3, K4-J ) + A( I, J ) = A( I, J )*MUL + 230 CONTINUE + 240 CONTINUE +C + END IF +C + IF( .NOT.DONE ) + $ GO TO 10 +C + RETURN +C *** Last line of MB01QD *** + END Property changes on: trunk/octave-forge/main/control/devel/tf2ss/MB01QD.f ___________________________________________________________________ Added: svn:executable + * Added: trunk/octave-forge/main/control/devel/tf2ss/MB03OY.f =================================================================== --- trunk/octave-forge/main/control/devel/tf2ss/MB03OY.f (rev 0) +++ trunk/octave-forge/main/control/devel/tf2ss/MB03OY.f 2011-08-22 18:25:37 UTC (rev 8474) @@ -0,0 +1,388 @@ + SUBROUTINE MB03OY( M, N, A, LDA, RCOND, SVLMAX, RANK, SVAL, JPVT, + $ TAU, DWORK, INFO ) +C +C SLICOT RELEASE 5.0. +C +C Copyright (c) 2002-2010 NICONET e.V. +C +C This program is free software: you can redistribute it and/or +C modify it under the terms of the GNU General Public License as +C published by the Free Software Foundation, either version 2 of +C the License, or (at your option) any later version. +C +C This program is distributed in the hope that it will be useful, +C but WITHOUT ANY WARRANTY; without even the implied warranty of +C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +C GNU General Public License for more details. +C +C You should have received a copy of the GNU General Public License +C along with this program. If not, see +C <http://www.gnu.org/licenses/>. +C +C PURPOSE +C +C To compute a rank-revealing QR factorization of a real general +C M-by-N matrix A, which may be rank-deficient, and estimate its +C effective rank using incremental condition estimation. +C +C The routine uses a truncated QR factorization with column pivoting +C [ R11 R12 ] +C A * P = Q * R, where R = [ ], +C [ 0 R22 ] +C with R11 defined as the largest leading upper triangular submatrix +C whose estimated condition number is less than 1/RCOND. The order +C of R11, RANK, is the effective rank of A. Condition estimation is +C performed during the QR factorization process. Matrix R22 is full +C (but of small norm), or empty. +C +C MB03OY does not perform any scaling of the matrix A. +C +C ARGUMENTS +C +C Input/Output Parameters +C +C M (input) INTEGER +C The number of rows of the matrix A. M >= 0. +C +C N (input) INTEGER +C The number of columns of the matrix A. N >= 0. +C +C A (input/output) DOUBLE PRECISION array, dimension +C ( LDA, N ) +C On entry, the leading M-by-N part of this array must +C contain the given matrix A. +C On exit, the leading RANK-by-RANK upper triangular part +C of A contains the triangular factor R11, and the elements +C below the diagonal in the first RANK columns, with the +C array TAU, represent the orthogonal matrix Q as a product +C of RANK elementary reflectors. +C The remaining N-RANK columns contain the result of the +C QR factorization process used. +C +C LDA INTEGER +C The leading dimension of the array A. LDA >= max(1,M). +C +C RCOND (input) DOUBLE PRECISION +C RCOND is used to determine the effective rank of A, which +C is defined as the order of the largest leading triangular +C submatrix R11 in the QR factorization with pivoting of A, +C whose estimated condition number is less than 1/RCOND. +C 0 <= RCOND <= 1. +C NOTE that when SVLMAX > 0, the estimated rank could be +C less than that defined above (see SVLMAX). +C +C SVLMAX (input) DOUBLE PRECISION +C If A is a submatrix of another matrix B, and the rank +C decision should be related to that matrix, then SVLMAX +C should be an estimate of the largest singular value of B +C (for instance, the Frobenius norm of B). If this is not +C the case, the input value SVLMAX = 0 should work. +C SVLMAX >= 0. +C +C RANK (output) INTEGER +C The effective (estimated) rank of A, i.e., the order of +C the submatrix R11. +C +C SVAL (output) DOUBLE PRECISION array, dimension ( 3 ) +C The estimates of some of the singular values of the +C triangular factor R: +C SVAL(1): largest singular value of R(1:RANK,1:RANK); +C SVAL(2): smallest singular value of R(1:RANK,1:RANK); +C SVAL(3): smallest singular value of R(1:RANK+1,1:RANK+1), +C if RANK < MIN( M, N ), or of R(1:RANK,1:RANK), +C otherwise. +C If the triangular factorization is a rank-revealing one +C (which will be the case if the leading columns were well- +C conditioned), then SVAL(1) will also be an estimate for +C the largest singular value of A, and SVAL(2) and SVAL(3) +C will be estimates for the RANK-th and (RANK+1)-st singular +C values of A, respectively. +C By examining these values, one can confirm that the rank +C is well defined with respect to the chosen value of RCOND. +C The ratio SVAL(1)/SVAL(2) is an estimate of the condition +C number of R(1:RANK,1:RANK). +C +C JPVT (output) INTEGER array, dimension ( N ) +C If JPVT(i) = k, then the i-th column of A*P was the k-th +C column of A. +C +C TAU (output) DOUBLE PRECISION array, dimension ( MIN( M, N ) ) +C The leading RANK elements of TAU contain the scalar +C factors of the elementary reflectors. +C +C Workspace +C +C DWORK DOUBLE PRECISION array, dimension ( 3*N-1 ) +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 routine computes a truncated QR factorization with column +C pivoting of A, A * P = Q * R, with R defined above, and, +C during this process, finds the largest leading submatrix whose +C estimated condition number is less than 1/RCOND, taking the +C possible positive value of SVLMAX into account. This is performed +C using the LAPACK incremental condition estimation scheme and a +C slightly modified rank decision test. The factorization process +C stops when RANK has been determined. +C +C The matrix Q is represented as a product of elementary reflectors +C +C Q = H(1) H(2) . . . H(k), where k = rank <= min(m,n). +C +C Each H(i) has the form +C +C H = I - tau * v * v' +C +C where tau is a real scalar, and v is a real vector with +C v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in +C A(i+1:m,i), and tau in TAU(i). +C +C The matrix P is represented in jpvt as follows: If +C jpvt(j) = i +C then the jth column of P is the ith canonical unit vector. +C +C REFERENCES +C +C [1] Bischof, C.H. and P. Tang. +C Generalizing Incremental Condition Estimation. +C LAPACK Working Notes 32, Mathematics and Computer Science +C Division, Argonne National Laboratory, UT, CS-91-132, +C May 1991. +C +C [2] Bischof, C.H. and P. Tang. +C Robust Incremental Condition Estimation. +C LAPACK Working Notes 33, Mathematics and Computer Science +C Division, Argonne National Laboratory, UT, CS-91-133, +C May 1991. +C +C NUMERICAL ASPECTS +C +C The algorithm is backward stable. +C +C CONTRIBUTOR +C +C V. Sima, Katholieke Univ. Leuven, Belgium, Feb. 1998. +C +C REVISIONS +C +C V. Sima, Research Institute for Informatics, Bucharest, Jan. 2009. +C V. Sima, Jan. 2010, following Bujanovic and Drmac's suggestion. +C +C KEYWORDS +C +C Eigenvalue problem, matrix operations, orthogonal transformation, +C singular values. +C +C ****************************************************************** +C +C .. Parameters .. + INTEGER IMAX, IMIN + PARAMETER ( IMAX = 1, IMIN = 2 ) + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) +C .. Scalar Arguments .. + INTEGER INFO, LDA, M, N, RANK + DOUBLE PRECISION RCOND, SVLMAX +C .. Array Arguments .. + INTEGER JPVT( * ) + DOUBLE PRECISION A( LDA, * ), DWORK( * ), SVAL( 3 ), TAU( * ) +C .. +C .. Local Scalars .. + INTEGER I, ISMAX, ISMIN, ITEMP, J, MN, PVT + DOUBLE PRECISION AII, C1, C2, S1, S2, SMAX, SMAXPR, SMIN, + $ SMINPR, TEMP, TEMP2, TOLZ +C .. +C .. External Functions .. + INTEGER IDAMAX + DOUBLE PRECISION DLAMCH, DNRM2 + EXTERNAL DLAMCH, DNRM2, IDAMAX +C .. External Subroutines .. + EXTERNAL DLAIC1, DLARF, DLARFG, DSCAL, DSWAP, XERBLA +C .. Intrinsic Functions .. + INTRINSIC ABS, MAX, MIN, SQRT +C .. +C .. Executable Statements .. +C +C Test the input scalar arguments. +C + INFO = 0 + IF( M.LT.0 ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( LDA.LT.MAX( 1, M ) ) THEN + INFO = -4 + ELSE IF( RCOND.LT.ZERO .OR. RCOND.GT.ONE ) THEN + INFO = -5 + ELSE IF( SVLMAX.LT.ZERO ) THEN + INFO = -6 + END IF +C + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'MB03OY', -INFO ) + RETURN + END IF +C +C Quick return if possible. +C + MN = MIN( M, N ) + IF( MN.EQ.0 ) THEN + RANK = 0 + SVAL( 1 ) = ZERO + SVAL( 2 ) = ZERO + SVAL( 3 ) = ZERO + RETURN + END IF +C + TOLZ = SQRT( DLAMCH( 'Epsilon' ) ) + ISMIN = 1 + ISMAX = ISMIN + N +C +C Initialize partial column norms and pivoting vector. The first n +C elements of DWORK store the exact column norms. The already used +C leading part is then overwritten by the condition estimator. +C + DO 10 I = 1, N + DWORK( I ) = DNRM2( M, A( 1, I ), 1 ) + DWORK( N+I ) = DWORK( I ) + JPVT( I ) = I + 10 CONTINUE +C +C Compute factorization and determine RANK using incremental +C condition estimation. +C + RANK = 0 +C + 20 CONTINUE + IF( RANK.LT.MN ) THEN + I = RANK + 1 +C +C Determine ith pivot column and swap if necessary. +C + PVT = ( I-1 ) + IDAMAX( N-I+1, DWORK( I ), 1 ) +C + IF( PVT.NE.I ) THEN + CALL DSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 ) + ITEMP = JPVT( PVT ) + JPVT( PVT ) = JPVT( I ) + JPVT( I ) = ITEMP + DWORK( PVT ) = DWORK( I ) + DWORK( N+PVT ) = DWORK( N+I ) + END IF +C +C Save A(I,I) and generate elementary reflector H(i). +C + IF( I.LT.M ) THEN + AII = A( I, I ) + CALL DLARFG( M-I+1, A( I, I ), A( I+1, I ), 1, TAU( I ) ) + ELSE + TAU( M ) = ZERO + END IF +C + IF( RANK.EQ.0 ) THEN +C +C Initialize; exit if matrix is zero (RANK = 0). +C + SMAX = ABS( A( 1, 1 ) ) + IF ( SMAX.EQ.ZERO ) THEN + SVAL( 1 ) = ZERO + SVAL( 2 ) = ZERO + SVAL( 3 ) = ZERO + RETURN + END IF + SMIN = SMAX + SMAXPR = SMAX + SMINPR = SMIN + C1 = ONE + C2 = ONE + ELSE +C +C One step of incremental condition estimation. +C + CALL DLAIC1( IMIN, RANK, DWORK( ISMIN ), SMIN, A( 1, I ), + $ A( I, I ), SMINPR, S1, C1 ) + CALL DLAIC1( IMAX, RANK, DWORK( ISMAX ), SMAX, A( 1, I ), + $ A( I, I ), SMAXPR, S2, C2 ) + END IF +C + IF( SVLMAX*RCOND.LE.SMAXPR ) THEN + IF( SVLMAX*RCOND.LE.SMINPR ) THEN + IF( SMAXPR*RCOND.LE.SMINPR ) THEN +C +C Continue factorization, as rank is at least RANK. +C + IF( I.LT.N ) THEN +C +C Apply H(i) to A(i:m,i+1:n) from the left. +C + AII = A( I, I ) + A( I, I ) = ONE + CALL DLARF( 'Left', M-I+1, N-I, A( I, I ), 1, + $ TAU( I ), A( I, I+1 ), LDA, + $ DWORK( 2*N+1 ) ) + A( I, I ) = AII + END IF +C +C Update partial column norms. +C + DO 30 J = I + 1, N + IF( DWORK( J ).NE.ZERO ) THEN + TEMP = ABS( A( I, J ) ) / DWORK( J ) + TEMP = MAX( ( ONE + TEMP )*( ONE - TEMP ), ZERO) + TEMP2 = TEMP*( DWORK( J ) / DWORK( N+J ) )**2 + IF( TEMP2.LE.TOLZ ) THEN + IF( M-I.GT.0 ) THEN + DWORK( J ) = DNRM2( M-I, A( I+1, J ), 1 ) + DWORK( N+J ) = DWORK( J ) + ELSE + DWORK( J ) = ZERO + DWORK( N+J ) = ZERO + END IF + ELSE + DWORK( J ) = DWORK( J )*SQRT( TEMP ) + END IF + END IF + 30 CONTINUE +C + DO 40 I = 1, RANK + DWORK( ISMIN+I-1 ) = S1*DWORK( ISMIN+I-1 ) + DWORK( ISMAX+I-1 ) = S2*DWORK( ISMAX+I-1 ) + 40 CONTINUE +C + DWORK( ISMIN+RANK ) = C1 + DWORK( ISMAX+RANK ) = C2 + SMIN = SMINPR + SMAX = SMAXPR + RANK = RANK + 1 + GO TO 20 + END IF + END IF + END IF + END IF +C +C Restore the changed part of the (RANK+1)-th column and set SVAL. +C + IF ( RANK.LT.N ) THEN + IF ( I.LT.M ) THEN + CALL DSCAL( M-I, -A( I, I )*TAU( I ), A( I+1, I ), 1 ) + A( I, I ) = AII + END IF + END IF + IF ( RANK.EQ.0 ) THEN + SMIN = ZERO + SMINPR = ZERO + END IF + SVAL( 1 ) = SMAX + SVAL( 2 ) = SMIN + SVAL( 3 ) = SMINPR +C + RETURN +C *** Last line of MB03OY *** + END Property changes on: trunk/octave-forge/main/control/devel/tf2ss/MB03OY.f ___________________________________________________________________ Added: svn:executable + * Added: trunk/octave-forge/main/control/devel/tf2ss/TB01ID.f =================================================================== --- trunk/octave-forge/main/control/devel/tf2ss/TB01ID.f (rev 0) +++ trunk/octave-forge/main/control/devel/tf2ss/TB01ID.f 2011-08-22 18:25:37 UTC (rev 8474) @@ -0,0 +1,402 @@ + SUBROUTINE TB01ID( JOB, N, M, P, MAXRED, A, LDA, B, LDB, C, LDC, + $ SCALE, INFO ) +C +C SLICOT RELEASE 5.0. +C +C Copyright (c) 2002-2010 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 reduce the 1-norm of a system matrix +C +C S = ( A B ) +C ( C 0 ) +C +C corresponding to the triple (A,B,C), by balancing. This involves +C a diagonal similarity transformation inv(D)*A*D applied +C iteratively to A to make the rows and columns of +C -1 +C diag(D,I) * S * diag(D,I) +C +C as close in norm as possible. +C +C The balancing can be performed optionally on the following +C particular system matrices +C +C S = A, S = ( A B ) or S = ( A ) +C ( C ) +C +C ARGUMENTS +C +C Mode Parameters +C +C JOB CHARACTER*1 +C Indicates which matrices are involved in balancing, as +C follows: +C = 'A': All matrices are involved in balancing; +C = 'B': B and A matrices are involved in balancing; +C = 'C': C and A matrices are involved in balancing; +C = 'N': B and C matrices are not involved in balancing. +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 MAXRED (input/output) DOUBLE PRECISION +C On entry, the maximum allowed reduction in the 1-norm of +C S (in an iteration) if zero rows or columns are +C encountered. +C If MAXRED > 0.0, MAXRED must be larger than one (to enable +C the norm reduction). +C If MAXRED <= 0.0, then the value 10.0 for MAXRED is +C used. +C On exit, if the 1-norm of the given matrix S is non-zero, +C the ratio between the 1-norm of the given matrix and the +C 1-norm of the balanced matrix. +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 balanced matrix inv(D)*A*D. +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 (LDB,M) +C On entry, if M > 0, the leading N-by-M part of this array +C must contain the system input matrix B. +C On exit, if M > 0, the leading N-by-M part of this array +C contains the balanced matrix inv(D)*B. +C The array B is not referenced if M = 0. +C +C LDB INTEGER +C The leading dimension of the array B. +C LDB >= MAX(1,N) if M > 0. +C LDB >= 1 if M = 0. +C +C C (input/output) DOUBLE PRECISION array, dimension (LDC,N) +C On entry, if P > 0, the leading P-by-N part of this array +C must contain the system output matrix C. +C On exit, if P > 0, the leading P-by-N part of this array +C contains the balanced matrix C*D. +C The array C is not referenced if P = 0. +C +C LDC INTEGER +C The leading dimension of the array C. LDC >= MAX(1,P). +C +C SCALE (output) DOUBLE PRECISION array, dimension (N) +C The scaling factors applied to S. If D(j) is the scaling +C factor applied to row and column j, then SCALE(j) = D(j), +C for j = 1,...,N. +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 Balancing consists of applying a diagonal similarity +C transformation +C -1 +C diag(D,I) * S * diag(D,I) +C +C to make the 1-norms of each row of the first N rows of S and its +C corresponding column nearly equal. +C +C Information about the diagonal matrix D is returned in the vector +C SCALE. +C +C REFERENCES +C +C [1] Anderson, E., Bai, Z., Bischof, C., Demmel, J., Dongarra, J., +C Du Croz, J., Greenbaum, A., Hammarling, S., McKenney, A., +C Ostrouchov, S., and Sorensen, D. +C LAPACK Users' Guide: Second Edition. +C SIAM, Philadelphia, 1995. +C +C NUMERICAL ASPECTS +C +C None. +C +C CONTRIBUTOR +C +C Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, Jan. 1998. +C This subroutine is based on LAPACK routine DGEBAL, and routine +C BALABC (A. Varga, German Aerospace Research Establishment, DLR). +C +C REVISIONS +C +C - +C +C KEYWORDS +C +C Balancing, eigenvalue, matrix algebra, matrix operations, +C similarity transformation. +C +C ********************************************************************* +C +C .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) + DOUBLE PRECISION SCLFAC + PARAMETER ( SCLFAC = 1.0D+1 ) + DOUBLE PRECISION FACTOR, MAXR + PARAMETER ( FACTOR = 0.95D+0, MAXR = 10.0D+0 ) +C .. +C .. Scalar Arguments .. + CHARACTER JOB + INTEGER INFO, LDA, LDB, LDC, M, N, P + DOUBLE PRECISION MAXRED +C .. +C .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ), + $ SCALE( * ) +C .. +C .. Local Scalars .. + LOGICAL NOCONV, WITHB, WITHC + INTEGER I, ICA, IRA, J + DOUBLE PRECISION CA, CO, F, G, MAXNRM, RA, RO, S, SFMAX1, + $ SFMAX2, SFMIN1, SFMIN2, SNORM, SRED +C .. +C .. External Functions .. + LOGICAL LSAME + INTEGER IDAMAX + DOUBLE PRECISION DASUM, DLAMCH + EXTERNAL DASUM, DLAMCH, IDAMAX, LSAME +C .. +C .. External Subroutines .. + EXTERNAL DSCAL, XERBLA +C .. +C .. Intrinsic Functions .. + INTRINSIC ABS, MAX, MIN +C .. +C .. Executable Statements .. +C +C Test the scalar input arguments. +C + INFO = 0 + WITHB = LSAME( JOB, 'A' ) .OR. LSAME( JOB, 'B' ) + WITHC = LSAME( JOB, 'A' ) .OR. LSAME( JOB, 'C' ) +C + IF( .NOT.WITHB .AND. .NOT.WITHC .AND. .NOT.LSAME( JOB, '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( MAXRED.GT.ZERO .AND. MAXRED.LT.ONE ) THEN + INFO = -5 + ELSE IF( LDA.LT.MAX( 1, N ) ) THEN + INFO = -7 + ELSE IF( ( M.GT.0 .AND. LDB.LT.MAX( 1, N ) ) .OR. + $ ( M.EQ.0 .AND. LDB.LT.1 ) ) THEN + INFO = -9 + ELSE IF( LDC.LT.MAX( 1, P ) ) THEN + INFO = -11 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'TB01ID', -INFO ) + RETURN + END IF +C + IF( N.EQ.0 ) + $ RETURN +C +C Compute the 1-norm of the required part of matrix S and exit if +C it is zero. +C + SNORM = ZERO +C + DO 10 J = 1, N + SCALE( J ) = ONE + CO = DASUM( N, A( 1, J ), 1 ) + IF( WITHC .AND. P.GT.0 ) + $ CO = CO + DASUM( P, C( 1, J ), 1 ) + SNORM = MAX( SNORM, CO ) + 10 CONTINUE +C + IF( WITHB ) THEN +C + DO 20 J = 1, M + SNORM = MAX( SNORM, DASUM( N, B( 1, J ), 1 ) ) + 20 CONTINUE +C + END IF +C + IF( SNORM.EQ.ZERO ) + $ RETURN +C +C Set some machine parameters and the maximum reduction in the +C 1-norm of S if zero rows or columns are encountered. +C + SFMIN1 = DLAMCH( 'S' ) / DLAMCH( 'P' ) + SFMAX1 = ONE / SFMIN1 + SFMIN2 = SFMIN1*SCLFAC + SFMAX2 = ONE / SFMIN2 +C + SRED = MAXRED + IF( SRED.LE.ZERO ) SRED = MAXR +C + MAXNRM = MAX( SNORM/SRED, SFMIN1 ) +C +C Balance the matrix. +C +C Iterative loop for norm reduction. +C + 30 CONTINUE + NOCONV = .FALSE. +C + DO 90 I = 1, N + CO = ZERO + RO = ZERO +C + DO 40 J = 1, N + IF( J.EQ.I ) + $ GO TO 40 + CO = CO + ABS( A( J, I ) ) + RO = RO + ABS( A( I, J ) ) + 40 CONTINUE +C + ICA = IDAMAX( N, A( 1, I ), 1 ) + CA = ABS( A( ICA, I ) ) + IRA = IDAMAX( N, A( I, 1 ), LDA ) + RA = ABS( A( I, IRA ) ) +C + IF( WITHC .AND. P.GT.0 ) THEN + CO = CO + DASUM( P, C( 1, I ), 1 ) + ICA = IDAMAX( P, C( 1, I ), 1 ) + CA = MAX( CA, ABS( C( ICA, I ) ) ) + END IF +C + IF( WITHB .AND. M.GT.0 ) THEN + RO = RO + DASUM( M, B( I, 1 ), LDB ) + IRA = IDAMAX( M, B( I, 1 ), LDB ) + RA = MAX( RA, ABS( B( I, IRA ) ) ) + END IF +C +C Special case of zero CO and/or RO. +C + IF( CO.EQ.ZERO .AND. RO.EQ.ZERO ) + $ GO TO 90 + IF( CO.EQ.ZERO ) THEN + IF( RO.LE.MAXNRM ) + $ GO TO 90 + CO = MAXNRM + END IF + IF( RO.EQ.ZERO ) THEN + IF( CO.LE.MAXNRM ) + $ GO TO 90 + RO = MAXNRM + END IF +C +C Guard against zero CO or RO due to underflow. +C + G = RO / SCLFAC + F = ONE + S = CO + RO + 50 CONTINUE + IF( CO.GE.G .OR. MAX( F, CO, CA ).GE.SFMAX2 .OR. + $ MIN( RO, G, RA ).LE.SFMIN2 )GO TO 60 + F = F*SCLFAC + CO = CO*SCLFAC + CA = CA*SCLFAC + G = G / SCLFAC + RO = RO / SCLFAC + RA = RA / SCLFAC + GO TO 50 +C + 60 CONTINUE + G = CO / SCLFAC + 70 CONTINUE + IF( G.LT.RO .OR. MAX( RO, RA ).GE.SFMAX2 .OR. + $ MIN( F, CO, G, CA ).LE.SFMIN2 )GO TO 80 + F = F / SCLFAC + CO = CO / SCLFAC + CA = CA / SCLFAC + G = G / SCLFAC + RO = RO*SCLFAC + RA = RA*SCLFAC + GO TO 70 +C +C Now balance. +C + 80 CONTINUE + IF( ( CO+RO ).GE.FACTOR*S ) + $ GO TO 90 + IF( F.LT.ONE .AND. SCALE( I ).LT.ONE ) THEN + IF( F*SCALE( I ).LE.SFMIN1 ) + $ GO TO 90 + END IF + IF( F.GT.ONE .AND. SCALE( I ).GT.ONE ) THEN + IF( SCALE( I ).GE.SFMAX1 / F ) + $ GO TO 90 + END IF + G = ONE / F + SCALE( I ) = SCALE( I )*F + NOCONV = .TRUE. +C + CALL DSCAL( N, G, A( I, 1 ), LDA ) + CALL DSCAL( N, F, A( 1, I ), 1 ) + IF( M.GT.0 ) CALL DSCAL( M, G, B( I, 1 ), LDB ) + IF( P.GT.0 ) CALL DSCAL( P, F, C( 1, I ), 1 ) +C + 90 CONTINUE +C + IF( NOCONV ) + $ GO TO 30 +C +C Set the norm reduction parameter. +C + MAXRED = SNORM + SNORM = ZERO +C + DO 100 J = 1, N + CO = DASUM( N, A( 1, J ), 1 ) + IF( WITHC .AND. P.GT.0 ) + $ CO = CO + DASUM( P, C( 1, J ), 1 ) + SNORM = MAX( SNORM, CO ) + 100 CONTINUE +C + IF( WITHB ) THEN +C + DO 110 J = 1, M + SNORM = MAX( SNORM, DASUM( N, B( 1, J ), 1 ) ) + 110 CONTINUE +C + END IF + MAXRED = MAXRED/SNORM + RETURN +C *** Last line of TB01ID *** + END Property changes on: trunk/octave-forge/main/control/devel/tf2ss/TB01ID.f ___________________________________________________________________ Added: svn:executable + * Added: trunk/octave-forge/main/control/devel/tf2ss/TB01PD.f =================================================================== --- trunk/octave-forge/main/control/devel/tf2ss/TB01PD.f (rev 0) +++ trunk/octave-forge/main/control/devel/tf2ss/TB01PD.f 2011-08-22 18:25:37 UTC (rev 8474) @@ -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-2010 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,... [truncated message content] |
From: <par...@us...> - 2011-09-03 16:37:50
|
Revision: 8487 http://octave.svn.sourceforge.net/octave/?rev=8487&view=rev Author: paramaniac Date: 2011-09-03 16:37:43 +0000 (Sat, 03 Sep 2011) Log Message: ----------- control: delete *.o for developer makefiles Modified Paths: -------------- trunk/octave-forge/main/control/devel/makefile_all.m trunk/octave-forge/main/control/devel/makefile_chol.m trunk/octave-forge/main/control/devel/makefile_h2syn.m trunk/octave-forge/main/control/devel/makefile_hankel.m trunk/octave-forge/main/control/devel/makefile_helpers.m trunk/octave-forge/main/control/devel/makefile_hinfsyn.m trunk/octave-forge/main/control/devel/makefile_lqr.m trunk/octave-forge/main/control/devel/makefile_lyap.m trunk/octave-forge/main/control/devel/makefile_minreal.m trunk/octave-forge/main/control/devel/makefile_ncfsyn.m trunk/octave-forge/main/control/devel/makefile_norm.m trunk/octave-forge/main/control/devel/makefile_place.m trunk/octave-forge/main/control/devel/makefile_scale.m trunk/octave-forge/main/control/devel/makefile_ss2tf.m trunk/octave-forge/main/control/devel/makefile_staircase.m trunk/octave-forge/main/control/devel/makefile_zero.m trunk/octave-forge/main/control/devel/slred/slconred/makefile_slconred.m trunk/octave-forge/main/control/devel/slred/slmodred/makefile_slmodred.m trunk/octave-forge/main/control/devel/tf2ss/makefile_tf2ss.m Modified: trunk/octave-forge/main/control/devel/makefile_all.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_all.m 2011-09-02 21:30:40 UTC (rev 8486) +++ trunk/octave-forge/main/control/devel/makefile_all.m 2011-09-03 16:37:43 UTC (rev 8487) @@ -27,4 +27,5 @@ makefile_staircase makefile_zero +system ("rm *.o"); cd (homedir); Modified: trunk/octave-forge/main/control/devel/makefile_chol.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_chol.m 2011-09-02 21:30:40 UTC (rev 8486) +++ trunk/octave-forge/main/control/devel/makefile_chol.m 2011-09-03 16:37:43 UTC (rev 8487) @@ -21,4 +21,5 @@ SG03BD.f SG03BV.f SG03BU.f SG03BW.f SG03BX.f \ SG03BY.f MB02UU.f MB02UV.f +system ("rm *.o"); cd (homedir); Modified: trunk/octave-forge/main/control/devel/makefile_h2syn.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_h2syn.m 2011-09-02 21:30:40 UTC (rev 8486) +++ trunk/octave-forge/main/control/devel/makefile_h2syn.m 2011-09-03 16:37:43 UTC (rev 8487) @@ -26,4 +26,5 @@ MA02ED.f select.f SB03MX.f SB02MR.f SB02MV.f \ MB01UD.f SB03MV.f SB04PX.f +system ("rm *.o"); cd (homedir); Modified: trunk/octave-forge/main/control/devel/makefile_hankel.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_hankel.m 2011-09-02 21:30:40 UTC (rev 8486) +++ trunk/octave-forge/main/control/devel/makefile_hankel.m 2011-09-03 16:37:43 UTC (rev 8487) @@ -22,4 +22,5 @@ SB03OR.f SB03OY.f SB04PX.f MB04NY.f MB04OY.f \ SB03OV.f +system ("rm *.o"); cd (homedir); Modified: trunk/octave-forge/main/control/devel/makefile_helpers.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_helpers.m 2011-09-02 21:30:40 UTC (rev 8486) +++ trunk/octave-forge/main/control/devel/makefile_helpers.m 2011-09-03 16:37:43 UTC (rev 8487) @@ -19,4 +19,5 @@ mkoctfile is_real_square_matrix.cc +system ("rm *.o"); cd (homedir); \ No newline at end of file Modified: trunk/octave-forge/main/control/devel/makefile_hinfsyn.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_hinfsyn.m 2011-09-02 21:30:40 UTC (rev 8486) +++ trunk/octave-forge/main/control/devel/makefile_hinfsyn.m 2011-09-03 16:37:43 UTC (rev 8487) @@ -26,4 +26,5 @@ SB03MX.f SB02MR.f SB02MV.f MB01UD.f SB03MV.f \ SB04PX.f +system ("rm *.o"); cd (homedir); Modified: trunk/octave-forge/main/control/devel/makefile_lqr.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_lqr.m 2011-09-02 21:30:40 UTC (rev 8486) +++ trunk/octave-forge/main/control/devel/makefile_lqr.m 2011-09-03 16:37:43 UTC (rev 8487) @@ -19,4 +19,5 @@ SG02AD.f SB02OU.f SB02OV.f SB02OW.f SB02OY.f \ MB01SD.f MB02VD.f MB02PD.f MA02GD.f +system ("rm *.o"); cd (homedir); Modified: trunk/octave-forge/main/control/devel/makefile_lyap.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_lyap.m 2011-09-02 21:30:40 UTC (rev 8486) +++ trunk/octave-forge/main/control/devel/makefile_lyap.m 2011-09-03 16:37:43 UTC (rev 8487) @@ -30,4 +30,5 @@ SG03AD.f MB01RW.f MB01RD.f SG03AX.f SG03AY.f \ MB02UU.f MB02UV.f +system ("rm *.o"); cd (homedir); Modified: trunk/octave-forge/main/control/devel/makefile_minreal.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_minreal.m 2011-09-02 21:30:40 UTC (rev 8486) +++ trunk/octave-forge/main/control/devel/makefile_minreal.m 2011-09-03 16:37:43 UTC (rev 8487) @@ -18,4 +18,5 @@ mkoctfile sltg01jd.cc \ TG01JD.f TG01AD.f TB01XD.f MA02CD.f TG01HX.f +system ("rm *.o"); cd (homedir); Modified: trunk/octave-forge/main/control/devel/makefile_ncfsyn.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_ncfsyn.m 2011-09-02 21:30:40 UTC (rev 8486) +++ trunk/octave-forge/main/control/devel/makefile_ncfsyn.m 2011-09-03 16:37:43 UTC (rev 8487) @@ -31,4 +31,5 @@ MB02VD.f SB02OY.f SB02OW.f SB02OV.f SB02OU.f \ SB02MR.f MA02GD.f SB02MV.f +system ("rm *.o"); cd (homedir); Modified: trunk/octave-forge/main/control/devel/makefile_norm.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_norm.m 2011-09-02 21:30:40 UTC (rev 8486) +++ trunk/octave-forge/main/control/devel/makefile_norm.m 2011-09-03 16:37:43 UTC (rev 8487) @@ -27,4 +27,5 @@ MB04QF.f MB03YA.f MB03YD.f MB02RZ.f MB04QU.f \ MB02SZ.f MB03YT.f +system ("rm *.o"); cd (homedir); Modified: trunk/octave-forge/main/control/devel/makefile_place.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_place.m 2011-09-02 21:30:40 UTC (rev 8486) +++ trunk/octave-forge/main/control/devel/makefile_place.m 2011-09-03 16:37:43 UTC (rev 8487) @@ -15,4 +15,5 @@ SB01BD.f MB03QD.f MB03QY.f SB01BX.f SB01BY.f \ select.f +system ("rm *.o"); cd (homedir); Modified: trunk/octave-forge/main/control/devel/makefile_scale.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_scale.m 2011-09-02 21:30:40 UTC (rev 8486) +++ trunk/octave-forge/main/control/devel/makefile_scale.m 2011-09-03 16:37:43 UTC (rev 8487) @@ -19,4 +19,5 @@ mkoctfile sltg01ad.cc \ TG01AD.f +system ("rm *.o"); cd (homedir); Modified: trunk/octave-forge/main/control/devel/makefile_ss2tf.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_ss2tf.m 2011-09-02 21:30:40 UTC (rev 8486) +++ trunk/octave-forge/main/control/devel/makefile_ss2tf.m 2011-09-03 16:37:43 UTC (rev 8487) @@ -16,4 +16,5 @@ TB04BX.f MA02AD.f MB02RD.f MB01PD.f MB02SD.f \ MB01QD.f +system ("rm *.o"); cd (homedir); Modified: trunk/octave-forge/main/control/devel/makefile_staircase.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_staircase.m 2011-09-02 21:30:40 UTC (rev 8486) +++ trunk/octave-forge/main/control/devel/makefile_staircase.m 2011-09-03 16:37:43 UTC (rev 8487) @@ -30,4 +30,5 @@ TG01ID.f TB01XD.f MA02CD.f AB07MD.f TG01HX.f \ MA02BD.f +system ("rm *.o"); cd (homedir); \ No newline at end of file Modified: trunk/octave-forge/main/control/devel/makefile_zero.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_zero.m 2011-09-02 21:30:40 UTC (rev 8486) +++ trunk/octave-forge/main/control/devel/makefile_zero.m 2011-09-03 16:37:43 UTC (rev 8487) @@ -24,4 +24,5 @@ mkoctfile sltg04bx.cc \ TG04BX.f MB02RD.f MB02SD.f +system ("rm *.o"); cd (homedir); Modified: trunk/octave-forge/main/control/devel/slred/slconred/makefile_slconred.m =================================================================== --- trunk/octave-forge/main/control/devel/slred/slconred/makefile_slconred.m 2011-09-02 21:30:40 UTC (rev 8486) +++ trunk/octave-forge/main/control/devel/slred/slconred/makefile_slconred.m 2011-09-03 16:37:43 UTC (rev 8487) @@ -14,4 +14,6 @@ mkoctfile SB16CD.f SB16CY.f AB09IX.f SB03OD.f MB02UD.f \ AB09DD.f MA02AD.f MB03UD.f select.f SB03OU.f \ MB01SD.f SB03OT.f MB04ND.f MB04OD.f SB03OR.f \ - SB03OY.f SB04PX.f MB04NY.f MB04OY.f SB03OV.f \ No newline at end of file + SB03OY.f SB04PX.f MB04NY.f MB04OY.f SB03OV.f + +system ("rm *.o"); Modified: trunk/octave-forge/main/control/devel/slred/slmodred/makefile_slmodred.m =================================================================== --- trunk/octave-forge/main/control/devel/slred/slmodred/makefile_slmodred.m 2011-09-02 21:30:40 UTC (rev 8486) +++ trunk/octave-forge/main/control/devel/slred/slmodred/makefile_slmodred.m 2011-09-03 16:37:43 UTC (rev 8487) @@ -24,3 +24,5 @@ SB02MR.f SB02MS.f MB03QD.f SB02MU.f SB02MV.f \ SB02MW.f MB04ND.f MB04OD.f MB03QY.f SB03OR.f \ SB03OY.f SB04PX.f MB04NY.f MB04OY.f SB03OV.f + +system ("rm *.o"); Modified: trunk/octave-forge/main/control/devel/tf2ss/makefile_tf2ss.m =================================================================== --- trunk/octave-forge/main/control/devel/tf2ss/makefile_tf2ss.m 2011-09-02 21:30:40 UTC (rev 8486) +++ trunk/octave-forge/main/control/devel/tf2ss/makefile_tf2ss.m 2011-09-03 16:37:43 UTC (rev 8487) @@ -2,3 +2,5 @@ sltd04ad.cc \ TD04AD.f TD03AY.f TB01PD.f TB01XD.f AB07MD.f \ TB01UD.f TB01ID.f MB01PD.f MB03OY.f MB01QD.f + +system ("rm *.o"); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-09-21 19:44:55
|
Revision: 8578 http://octave.svn.sourceforge.net/octave/?rev=8578&view=rev Author: paramaniac Date: 2011-09-21 19:44:48 +0000 (Wed, 21 Sep 2011) Log Message: ----------- control: add zpk (sort of) Added Paths: ----------- trunk/octave-forge/main/control/devel/zpk.m trunk/octave-forge/main/control/devel/zpktest.m Added: trunk/octave-forge/main/control/devel/zpk.m =================================================================== --- trunk/octave-forge/main/control/devel/zpk.m (rev 0) +++ trunk/octave-forge/main/control/devel/zpk.m 2011-09-21 19:44:48 UTC (rev 8578) @@ -0,0 +1,102 @@ +## Copyright (C) 2011 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 LTI Syncope. If not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {@var{s} =} zpk (@var{"s"}) +## @deftypefnx {Function File} {@var{z} =} zpk (@var{"z"}, @var{tsam}) +## @deftypefnx {Function File} {@var{sys} =} zpk (@var{sys}) +## @deftypefnx {Function File} {@var{sys} =} zpk (@var{k}) +## @deftypefnx {Function File} {@var{sys} =} zpk (@var{z}, @var{p}, @var{k}, @dots{}) +## @deftypefnx {Function File} {@var{sys} =} zpk @var{z}, @var{p}, @var{k}, @var{tsam}, @dots{}) +## @deftypefnx {Function File} {@var{sys} =} zpk (@var{z}, @var{p}, @var{k}, @var{tsam}, @dots{}) +## Create transfer function model from zero-pole-gain data. +## +## @strong{Inputs} +## @table @var +## @item sys +## LTI model to be converted to transfer function. +## @item num +## Numerator or cell of numerators. Each numerator must be a row vector +## containing the coefficients of the polynomial in descending powers of +## the transfer function variable. +## @item den +## Denominator or cell of denominators. Each denominator must be a row vector +## containing the coefficients of the polynomial in descending powers of +## the transfer function variable. +## @item tsam +## Sampling time in seconds. If @var{tsam} is not specified, a continuous-time +## model is assumed. +## @item @dots{} +## Optional pairs of properties and values. +## Type @command{set (tf)} for more information. +## @end table +## +## @strong{Outputs} +## @table @var +## @item sys +## Transfer function model. +## @end table +## +## @strong{Example} +## @example +## @group +## @end group +## @end example +## +## @seealso{tf, ss, dss, frd} +## @end deftypefn + +## Author: Lukas Reichlin <luk...@gm...> +## Created: September 2011 +## Version: 0.1 + +function sys = zpk (z = {}, p = {}, k = [], varargin) + + switch (nargin) + case 0 + sys = tf (); + return; + + case 1 + if (isa (z, "lti") || is_real_matrix (z) || ischar (z)) + sys = tf (z); + return; + else + print_usage (); + endif + + case 2 + if (ischar (z) && issample (p, -1)) + sys = tf (z, p); + else + print_usage (); + endif + + otherwise + if (! iscell (z)) + z = {z}; + endif + if (! iscell (p)) + p = {p}; + endif + num = cellfun (@(zer, gain) real (gain * poly (zer)), z, num2cell (k), "uniformoutput", false); + den = cellfun (@(pol) real (poly (pol)), p, "uniformoutput", false); + sys = tf (num, den, varargin{:}); + endswitch + +endfunction + Added: trunk/octave-forge/main/control/devel/zpktest.m =================================================================== --- trunk/octave-forge/main/control/devel/zpktest.m (rev 0) +++ trunk/octave-forge/main/control/devel/zpktest.m 2011-09-21 19:44:48 UTC (rev 8578) @@ -0,0 +1,20 @@ +zpk (0, [], 1) + +zpk () + +zpk ([1,2;3,4]) + +zpk ("s") + +zpk ("z", 0.3) + +zpk ([1,2,3], [4,5,6], 7) + +zpk ([1,2,3], [4,5,6], 7, 0.1) + + +z = {[1, 5, 7], [1]; [1, 7], [1, 5, 5]}; +p = {[1, 5, 6], [1, 2]; [1, 8, 6], [1, 3, 2]}; +k = [1, 2; 3, 4]; + +zpk (z, p, k) This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-10-04 20:41:34
|
Revision: 8674 http://octave.svn.sourceforge.net/octave/?rev=8674&view=rev Author: paramaniac Date: 2011-10-04 20:41:28 +0000 (Tue, 04 Oct 2011) Log Message: ----------- control: add draft code Added Paths: ----------- trunk/octave-forge/main/control/devel/__dss_bilin__.m trunk/octave-forge/main/control/devel/test_dss_bilin.m Added: trunk/octave-forge/main/control/devel/__dss_bilin__.m =================================================================== --- trunk/octave-forge/main/control/devel/__dss_bilin__.m (rev 0) +++ trunk/octave-forge/main/control/devel/__dss_bilin__.m 2011-10-04 20:41:28 UTC (rev 8674) @@ -0,0 +1,64 @@ +## Copyright (C) 2011 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 LTI Syncope. If not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## 1. Discrete -> continuous +## _ +## E = alpha*E + A +## _ +## A = beta (A - alpha*E) +## _ +## B = sqrt(2*alpha*beta) * B +## _ -1 +## C = sqrt(2*alpha*beta) * C * (alpha*E + A) * E +## _ -1 +## D = D - C * (alpha*E + A) * B +## +## +## 2. Continuous -> discrete +## _ +## E = beta*E - A +## _ +## A = alpha (beta*E + A) +## _ +## B = sqrt(2*alpha*beta) * B +## _ -1 +## C = sqrt(2*alpha*beta) * C * (beta*E - A) * E +## _ -1 +## D = D + C * (beta*E - A) * B + +## Author: Lukas Reichlin <luk...@gm...> +## Created: October 2011 +## Version: 0.1 + +function [Ar, Br, Cr, Dr, Er] = __dss_bilin__ (A, B, C, D, E, beta, discrete) + + if (discrete) + Er = E + A; + Ar = beta * (A - E); + Br = sqrt (2*beta) * B; + Cr = sqrt (2*beta) * C / (E + A) * E; + Dr = D - C / (E + A) * B; + else + Er = beta*E - A; + Ar = beta*E + A; + Br = sqrt (2*beta) * B; + Cr = sqrt (2*beta) * C / (beta*E - A) * E; + Dr = D + C / (beta*E - A) * B; + endif + +endfunction Added: trunk/octave-forge/main/control/devel/test_dss_bilin.m =================================================================== --- trunk/octave-forge/main/control/devel/test_dss_bilin.m (rev 0) +++ trunk/octave-forge/main/control/devel/test_dss_bilin.m 2011-10-04 20:41:28 UTC (rev 8674) @@ -0,0 +1,5 @@ +sys = inv (Boeing707); + +d = c2d (sys, 0.2, "bilin") + +c = d2c (d, "bilin") This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-10-07 08:04:42
|
Revision: 8707 http://octave.svn.sourceforge.net/octave/?rev=8707&view=rev Author: paramaniac Date: 2011-10-07 08:04:32 +0000 (Fri, 07 Oct 2011) Log Message: ----------- control: remove cruft in devel folder Added Paths: ----------- trunk/octave-forge/main/control/devel/test_tf2ss.m trunk/octave-forge/main/control/devel/tf2ss_draft.m Removed Paths: ------------- trunk/octave-forge/main/control/devel/tf2ss/ Copied: trunk/octave-forge/main/control/devel/test_tf2ss.m (from rev 8706, trunk/octave-forge/main/control/devel/tf2ss/test_tf2ss.m) =================================================================== --- trunk/octave-forge/main/control/devel/test_tf2ss.m (rev 0) +++ trunk/octave-forge/main/control/devel/test_tf2ss.m 2011-10-07 08:04:32 UTC (rev 8707) @@ -0,0 +1,19 @@ +index = [ 3 3 ]; + +dcoeff = [ 1.0 6.0 11.0 6.0 + 1.0 6.0 11.0 6.0 ]; + +ucoeff = zeros (2, 2, 4); + +u11 = [ 1.0 6.0 12.0 7.0 ]; +u12 = [ 0.0 1.0 4.0 3.0 ]; +u21 = [ 0.0 0.0 1.0 1.0 ]; +u22 = [ 1.0 8.0 20.0 15.0 ]; + +ucoeff(1,1,:) = u11; +ucoeff(1,2,:) = u12; +ucoeff(2,1,:) = u21; +ucoeff(2,2,:) = u22; + +[a, b, c, d] = sltd04ad (ucoeff, dcoeff, index, 0) + Copied: trunk/octave-forge/main/control/devel/tf2ss_draft.m (from rev 8706, trunk/octave-forge/main/control/devel/tf2ss/tf2ss.m) =================================================================== --- trunk/octave-forge/main/control/devel/tf2ss_draft.m (rev 0) +++ trunk/octave-forge/main/control/devel/tf2ss_draft.m 2011-10-07 08:04:32 UTC (rev 8707) @@ -0,0 +1,76 @@ +function retsys = tf2ss_draft () + +num = {[1, 5, 7], [1]; [1, 7], [1, 5, 5]}; +den = {[1, 5, 6], [1, 2]; [1, 8, 6], [1, 3, 2]}; +sys = tf (num, den); + +%{ +sys = tf (1, [1, 0]) +sys = tf (1, [1, 1]) + +sys = tf (1, conv ([1, 1, 1], [1, 4, 6, 4, 1])) + +sys = tf () +sys = tf ("s") +%} + +sys = tf (WestlandLynx); +tol = sqrt (eps) + + [p, m] = size (sys); + [num, den] = tfdata (sys); + + numc = cell (p, m); + denc = cell (p, 1); + + ## multiply all denominators in a row and + ## update each numerator accordingly + for i = 1 : p + denc(i) = __conv__ (den{i,:}); + for j = 1 : m + idx = setdiff (1:m, j); + numc(i,j) = __conv__ (num{i,j}, den{i,idx}); + endfor + endfor + + len_numc = cellfun (@length, numc); + len_denc = cellfun (@length, denc); + + ## tfpoly ensures that there are no leading zeros + tmp = len_numc > repmat (len_denc, 1, m); + if (any (tmp(:))) + error ("tf: tf2ss: system must be proper"); + endif + + max_len_denc = max (len_denc(:)); + ucoeff = zeros (p, m, max_len_denc); + dcoeff = zeros (p, max_len_denc); + index = len_denc-1; + + for i = 1 : p + len = len_denc(i); + dcoeff(i, 1:len) = denc{i}; + for j = 1 : m + ucoeff(i, j, len-len_numc(i,j)+1 : len) = numc{i,j}; + endfor + endfor +index, prod (index), eps*prod (index) +min (sqrt (eps), eps*prod (index)) + [a, b, c, d] = sltd04ad (ucoeff, dcoeff, index, tol); + + retsys = ss (a, b, c, d); + +endfunction + + +function vec = __conv__ (vec, varargin) + + if (nargin == 1) + return; + else + for k = 1 : nargin-1 + vec = conv (vec, varargin{k}); + endfor + endif + +endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-10-11 15:30:45
|
Revision: 8730 http://octave.svn.sourceforge.net/octave/?rev=8730&view=rev Author: paramaniac Date: 2011-10-11 15:30:34 +0000 (Tue, 11 Oct 2011) Log Message: ----------- control: add some draft code to devel folder Added Paths: ----------- trunk/octave-forge/main/control/devel/@iddata/ trunk/octave-forge/main/control/devel/@iddata/iddata.m trunk/octave-forge/main/control/devel/__iddata_dim__.m Added: trunk/octave-forge/main/control/devel/@iddata/iddata.m =================================================================== --- trunk/octave-forge/main/control/devel/@iddata/iddata.m (rev 0) +++ trunk/octave-forge/main/control/devel/@iddata/iddata.m 2011-10-11 15:30:34 UTC (rev 8730) @@ -0,0 +1,34 @@ +## Author: Lukas Reichlin <luk...@gm...> +## Created: October 2011 +## Version: 0.1 + +function dat = iddata (y = [], u = [], tsam = [], varargin) + + if (nargin == 1 && isa (y, "iddata")) + dat = y; + return; + elseif (nargin < 3) + print_usage (); + endif + + if (! issample (tsam, 1)) + error ("iddata: invalid sampling time"); + endif + + [p, m] = __iddata_dim__ (y, u); + + outname = repmat ({""}, p, 1); + inname = repmat ({""}, m, 1); + + dat = struct ("y", y, "outname", {outname}, "outunit", {outname}, + "u", u, "inname", {inname}, "inunit", {inname}, + "tsam", tsam, "timeunit", {""}, + "name", "", "notes", {{}}, "userdata", []); + + dat = class (dat, "iddata"); + + if (nargin > 3) + dat = set (dat, varargin{:}); + endif + +endfunction \ No newline at end of file Added: trunk/octave-forge/main/control/devel/__iddata_dim__.m =================================================================== --- trunk/octave-forge/main/control/devel/__iddata_dim__.m (rev 0) +++ trunk/octave-forge/main/control/devel/__iddata_dim__.m 2011-10-11 15:30:34 UTC (rev 8730) @@ -0,0 +1,26 @@ +## Author: Lukas Reichlin <luk...@gm...> +## Created: October 2011 +## Version: 0.1 + +function [p, m] = __iddata_dim__ (y, u) + + if (! is_real_matrix (y, u)) + error ("iddata: inputs and outputs must be real"); + endif + + [ly, p] = size (y); + [lu, m] = size (u); + + if (ly != lu) + error ("iddata: matrices ""y"" and ""u"" must have the same number of samples (rows)"); + endif + + if (ly < p) + warning ("iddata: more outputs than samples - matrice ""y"" should probably be transposed"); + endif + + if (lu < m) + warning ("iddata: more inputs than samples - matrice ""u"" should probably be transposed"); + endif + +endfunction \ No newline at end of file This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-10-17 18:58:02
|
Revision: 8766 http://octave.svn.sourceforge.net/octave/?rev=8766&view=rev Author: paramaniac Date: 2011-10-17 18:57:55 +0000 (Mon, 17 Oct 2011) Log Message: ----------- control: add a test case, remove cruft Added Paths: ----------- trunk/octave-forge/main/control/devel/test_tf2dss.m Removed Paths: ------------- trunk/octave-forge/main/control/devel/__sys2ss__.m Deleted: trunk/octave-forge/main/control/devel/__sys2ss__.m =================================================================== --- trunk/octave-forge/main/control/devel/__sys2ss__.m 2011-10-17 18:55:20 UTC (rev 8765) +++ trunk/octave-forge/main/control/devel/__sys2ss__.m 2011-10-17 18:57:55 UTC (rev 8766) @@ -1,159 +0,0 @@ -## Copyright (C) 2009, 2011 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 LTI Syncope. If not, see <http://www.gnu.org/licenses/>. - -## -*- texinfo -*- -## TF to SS conversion. - -## Author: Lukas Reichlin <luk...@gm...> -## Created: October 2009 -## Version: 0.3 - -function [retsys, retlti] = __sys2ss__ (sys) - - ## TODO: determine appropriate tolerance from number of inputs - ## (since we multiply all denominators in a row), index, ... - ## default tolerance from TB01UD is TOLDEF = N*N*EPS - - ## SECRET WISH: a routine which accepts individual denominators for - ## each channel and which supports descriptor systems - - [p, m] = size (sys); - [num, den] = tfdata (sys); - - len_num = cellfun (@length, num); - len_den = cellfun (@length, den); - - ## check for properness - ## tfpoly ensures that there are no leading zeros - tmp = len_num > len_den; - if (any (tmp(:))) - ## separation into strictly proper and polynomial part - [numq, numr] = cellfun (@deconv, num, den, "uniformoutput", false); - numq = cellfun (@__remove_leading_zeros__, numq, "uniformoutput", false); - numr = cellfun (@__remove_leading_zeros__, numr, "uniformoutput", false); - - ## minimal state-space realization for the proper part - [a1, b1, c1] = __proper_tf2ss__ (numr, den, p, m); - e1 = eye (size (a1)); - - ## minimal realization for the polynomial part - [e2, a2, b2, c2] = __polynomial_tf2ss__ (numq, p, m); - - ## assemble irreducible descriptor realization - e = blkdiag (e1, e2); - a = blkdiag (a1, a2); - b = vertcat (b1, b2); - c = horzcat (c1, c2); - retsys = dss (a, b, c, [], e); - else - [a, b, c, d] = __proper_tf2ss__ (num, den, p, m); - retsys = ss (a, b, c, d); - endif - - retlti = sys.lti; # preserve lti properties such as tsam - -endfunction - - -## transfer function to state-space conversion for proper models -function [a, b, c, d] = __proper_tf2ss__ (num, den, p, m) -num, den - ## new cells for the TF of same row denominators - numc = cell (p, m); - denc = cell (p, 1); - - ## multiply all denominators in a row and - ## update each numerator accordingly - for i = 1 : p - denc(i) = __conv__ (den{i,:}); - for j = 1 : m - idx = setdiff (1:m, j); - numc(i,j) = __conv__ (num{i,j}, den{i,idx}); - endfor - endfor - - len_numc = cellfun (@length, numc); - len_denc = cellfun (@length, denc); - - ## check for properness - ## tfpoly ensures that there are no leading zeros - tmp = len_numc > repmat (len_denc, 1, m); - if (any (tmp(:))) - error ("tf: tf2ss: system must be proper"); - endif - - ## create arrays and fill in the data - ## in a way that Slicot TD04AD can use - max_len_denc = max (len_denc(:)); - ucoeff = zeros (p, m, max_len_denc); - dcoeff = zeros (p, max_len_denc); - index = len_denc-1; - - for i = 1 : p - len = len_denc(i); - dcoeff(i, 1:len) = denc{i}; - for j = 1 : m - ucoeff(i, j, len-len_numc(i,j)+1 : len) = numc{i,j}; - endfor - endfor - - [a, b, c, d] = sltd04ad (ucoeff, dcoeff, index, sqrt (eps)); - -endfunction - - -function [e2, a2, b2, c2] = __polynomial_tf2ss__ (numq, p, m) - - len_numq = cellfun (@length, numq); - max_len_numq = max (len_numq(:)); - numq = cellfun (@(x) prepad (x, max_len_numq, 0, 2), numq, "uniformoutput", false); - f = @(y) cellfun (@(x) x(y), numq); - s = num2cell (1 : max_len_numq); - D = cellfun (f, s, "uniformoutput", false) - - e2 = diag (ones (p*max_len_numq, 1), -p); - a2 = eye (p*max_len_numq) - b2 = vertcat (D{:}) - c2 = horzcat (zeros (p, (p-1)*max_len_numq), -eye (p)) - - ## TODO: remove uncontrollable part - -endfunction - -## convolution for more than two arguments -function vec = __conv__ (vec, varargin) - - if (nargin == 1) - return; - else - for k = 1 : nargin-1 - vec = conv (vec, varargin{k}); - endfor - endif - -endfunction - - -function p = __remove_leading_zeros__ (p) - - idx = find (p != 0); - - if (! isempty (idx) && idx(1) > 1) - p = p(idx(1) : end); # p(idx) would remove all zeros - endif - -endfunction \ No newline at end of file Added: trunk/octave-forge/main/control/devel/test_tf2dss.m =================================================================== --- trunk/octave-forge/main/control/devel/test_tf2dss.m (rev 0) +++ trunk/octave-forge/main/control/devel/test_tf2dss.m 2011-10-17 18:57:55 UTC (rev 8766) @@ -0,0 +1,15 @@ +num = {[1, 5, 6], [1, 2], [1, 2, -1, 2, 2]; [1, 8, 6], [1, 3, 2], [3, 4, 0]}; +den = {[1, 5, 7], [1], [1, -1, 1]; [1, 7], [1, 5, 5], [1, 1]}; + +%num = {[1, 5, 6], [1, 2], [1, 2, -1, 2, 2]; [1, 8, 6], [1, 3, 2], [3, 4, 0]}.'; +%den = {[1, 5, 7], [1], [1, -1, 1]; [1, 7], [1, 5, 5], [1, 1]}.'; + +%num = {[1, 5, 6], [1, 2]; [1, 8, 6], [1, 3, 2]}; +%den = {[1, 5, 7], [1]; [1, 7], [1, 5, 5]}; + + +sys = tf (num, den) + +P = ss (sys) + + This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-10-21 17:05:43
|
Revision: 8814 http://octave.svn.sourceforge.net/octave/?rev=8814&view=rev Author: paramaniac Date: 2011-10-21 17:05:35 +0000 (Fri, 21 Oct 2011) Log Message: ----------- control: add draft code Modified Paths: -------------- trunk/octave-forge/main/control/devel/ctrbf.m Added Paths: ----------- trunk/octave-forge/main/control/devel/ctrbf/ trunk/octave-forge/main/control/devel/ctrbf/MB01PD.f trunk/octave-forge/main/control/devel/ctrbf/MB01QD.f trunk/octave-forge/main/control/devel/ctrbf/MB03OY.f trunk/octave-forge/main/control/devel/ctrbf/TB01UD.dat trunk/octave-forge/main/control/devel/ctrbf/TB01UD.f trunk/octave-forge/main/control/devel/ctrbf/TB01UD.res trunk/octave-forge/main/control/devel/ctrbf/common.cc trunk/octave-forge/main/control/devel/ctrbf/makefile_ctrbf.m trunk/octave-forge/main/control/devel/ctrbf/sltb01ud.cc trunk/octave-forge/main/control/devel/ctrbf/test_ctrbf.m Added: trunk/octave-forge/main/control/devel/ctrbf/MB01PD.f =================================================================== --- trunk/octave-forge/main/control/devel/ctrbf/MB01PD.f (rev 0) +++ trunk/octave-forge/main/control/devel/ctrbf/MB01PD.f 2011-10-21 17:05:35 UTC (rev 8814) @@ -0,0 +1,271 @@ + SUBROUTINE MB01PD( SCUN, TYPE, M, N, KL, KU, ANRM, NBL, NROWS, A, + $ LDA, INFO ) +C +C SLICOT RELEASE 5.0. +C +C Copyright (c) 2002-2010 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 scale a matrix or undo scaling. Scaling is performed, if +C necessary, so that the matrix norm will be in a safe range of +C representable numbers. +C +C ARGUMENTS +C +C Mode Parameters +C +C SCUN CHARACTER*1 +C SCUN indicates the operation to be performed. +C = 'S': scale the matrix. +C = 'U': undo scaling of the matrix. +C +C TYPE CHARACTER*1 +C TYPE indicates the storage type of the input matrix. +C = 'G': A is a full matrix. +C = 'L': A is a (block) lower triangular matrix. +C = 'U': A is an (block) upper triangular matrix. +C = 'H': A is an (block) upper Hessenberg matrix. +C = 'B': A is a symmetric band matrix with lower bandwidth +C KL and upper bandwidth KU and with the only the +C lower half stored. +C = 'Q': A is a symmetric band matrix with lower bandwidth +C KL and upper bandwidth KU and with the only the +C upper half stored. +C = 'Z': A is a band matrix with lower bandwidth KL and +C upper bandwidth KU. +C +C Input/Output Parameters +C +C M (input) INTEGER +C The number of rows of the matrix A. M >= 0. +C +C N (input) INTEGER +C The number of columns of the matrix A. N >= 0. +C +C KL (input) INTEGER +C The lower bandwidth of A. Referenced only if TYPE = 'B', +C 'Q' or 'Z'. +C +C KU (input) INTEGER +C The upper bandwidth of A. Referenced only if TYPE = 'B', +C 'Q' or 'Z'. +C +C ANRM (input) DOUBLE PRECISION +C The norm of the initial matrix A. ANRM >= 0. +C When ANRM = 0 then an immediate return is effected. +C ANRM should be preserved between the call of the routine +C with SCUN = 'S' and the corresponding one with SCUN = 'U'. +C +C NBL (input) INTEGER +C The number of diagonal blocks of the matrix A, if it has a +C block structure. To specify that matrix A has no block +C structure, set NBL = 0. NBL >= 0. +C +C NROWS (input) INTEGER array, dimension max(1,NBL) +C NROWS(i) contains the number of rows and columns of the +C i-th diagonal block of matrix A. The sum of the values +C NROWS(i), for i = 1: NBL, should be equal to min(M,N). +C The elements of the array NROWS are not referenced if +C NBL = 0. +C +C A (input/output) DOUBLE PRECISION array, dimension (LDA,N) +C On entry, the leading M by N part of this array must +C contain the matrix to be scaled/unscaled. +C On exit, the leading M by N part of A will contain +C the modified matrix. +C The storage mode of A is specified by TYPE. +C +C LDA (input) INTEGER +C The leading dimension of the array A. LDA >= max(1,M). +C +C Error Indicator +C +C INFO (output) INTEGER +C = 0: successful exit +C < 0: if INFO = -i, the i-th argument had an illegal +C value. +C +C METHOD +C +C Denote by ANRM the norm of the matrix, and by SMLNUM and BIGNUM, +C two positive numbers near the smallest and largest safely +C representable numbers, respectively. The matrix is scaled, if +C needed, such that the norm of the result is in the range +C [SMLNUM, BIGNUM]. The scaling factor is represented as a ratio +C of two numbers, one of them being ANRM, and the other one either +C SMLNUM or BIGNUM, depending on ANRM being less than SMLNUM or +C larger than BIGNUM, respectively. For undoing the scaling, the +C norm is again compared with SMLNUM or BIGNUM, and the reciprocal +C of the previous scaling factor is used. +C +C CONTRIBUTOR +C +C V. Sima, Katholieke Univ. Leuven, Belgium, Nov. 1996. +C +C REVISIONS +C +C Oct. 2001, V. Sima, Research Institute for Informatics, Bucharest. +C +C ****************************************************************** +C +C .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) +C .. Scalar Arguments .. + CHARACTER SCUN, TYPE + INTEGER INFO, KL, KU, LDA, M, MN, N, NBL + DOUBLE PRECISION ANRM +C .. Array Arguments .. + INTEGER NROWS ( * ) + DOUBLE PRECISION A( LDA, * ) +C .. Local Scalars .. + LOGICAL FIRST, LSCALE + INTEGER I, ISUM, ITYPE + DOUBLE PRECISION BIGNUM, SMLNUM +C .. External Functions .. + LOGICAL LSAME + DOUBLE PRECISION DLAMCH + EXTERNAL DLAMCH, LSAME +C .. +C .. External Subroutines .. + EXTERNAL DLABAD, MB01QD, XERBLA +C .. Intrinsic Functions .. + INTRINSIC MAX, MIN +C .. Save statement .. + SAVE BIGNUM, FIRST, SMLNUM +C .. Data statements .. + DATA FIRST/.TRUE./ +C .. +C .. Executable Statements .. +C +C Test the input scalar arguments. +C + INFO = 0 + LSCALE = LSAME( SCUN, 'S' ) + IF( LSAME( TYPE, 'G' ) ) THEN + ITYPE = 0 + ELSE IF( LSAME( TYPE, 'L' ) ) THEN + ITYPE = 1 + ELSE IF( LSAME( TYPE, 'U' ) ) THEN + ITYPE = 2 + ELSE IF( LSAME( TYPE, 'H' ) ) THEN + ITYPE = 3 + ELSE IF( LSAME( TYPE, 'B' ) ) THEN + ITYPE = 4 + ELSE IF( LSAME( TYPE, 'Q' ) ) THEN + ITYPE = 5 + ELSE IF( LSAME( TYPE, 'Z' ) ) THEN + ITYPE = 6 + ELSE + ITYPE = -1 + END IF +C + MN = MIN( M, N ) +C + ISUM = 0 + IF( NBL.GT.0 ) THEN + DO 10 I = 1, NBL + ISUM = ISUM + NROWS(I) + 10 CONTINUE + END IF +C + IF( .NOT.LSCALE .AND. .NOT.LSAME( SCUN, 'U' ) ) THEN + INFO = -1 + ELSE IF( ITYPE.EQ.-1 ) THEN + INFO = -2 + ELSE IF( M.LT.0 ) THEN + INFO = -3 + ELSE IF( N.LT.0 .OR. + $ ( ( ITYPE.EQ.4 .OR. ITYPE.EQ.5 ) .AND. N.NE.M ) ) THEN + INFO = -4 + ELSE IF( ANRM.LT.ZERO ) THEN + INFO = -7 + ELSE IF( NBL.LT.0 ) THEN + INFO = -8 + ELSE IF( NBL.GT.0 .AND. ISUM.NE.MN ) THEN + INFO = -9 + ELSE IF( ITYPE.LE.3 .AND. LDA.LT.MAX( 1, M ) ) THEN + INFO = -11 + ELSE IF( ITYPE.GE.4 ) THEN + IF( KL.LT.0 .OR. KL.GT.MAX( M-1, 0 ) ) THEN + INFO = -5 + ELSE IF( KU.LT.0 .OR. KU.GT.MAX( N-1, 0 ) .OR. + $ ( ( ITYPE.EQ.4 .OR. ITYPE.EQ.5 ) .AND. KL.NE.KU ) ) + $ THEN + INFO = -6 + ELSE IF( ( ITYPE.EQ.4 .AND. LDA.LT.KL+1 ) .OR. + $ ( ITYPE.EQ.5 .AND. LDA.LT.KU+1 ) .OR. + $ ( ITYPE.EQ.6 .AND. LDA.LT.2*KL+KU+1 ) ) THEN + INFO = -11 + END IF + END IF +C + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'MB01PD', -INFO ) + RETURN + END IF +C +C Quick return if possible. +C + IF( MN.EQ.0 .OR. ANRM.EQ.ZERO ) + $ RETURN +C + IF ( FIRST ) THEN +C +C Get machine parameters. +C + SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' ) + BIGNUM = ONE / SMLNUM + CALL DLABAD( SMLNUM, BIGNUM ) + FIRST = .FALSE. + END IF +C + IF ( LSCALE ) THEN +C +C Scale A, if its norm is outside range [SMLNUM,BIGNUM]. +C + IF( ANRM.LT.SMLNUM ) THEN +C +C Scale matrix norm up to SMLNUM. +C + CALL MB01QD( TYPE, M, N, KL, KU, ANRM, SMLNUM, NBL, NROWS, + $ A, LDA, INFO ) + ELSE IF( ANRM.GT.BIGNUM ) THEN +C +C Scale matrix norm down to BIGNUM. +C + CALL MB01QD( TYPE, M, N, KL, KU, ANRM, BIGNUM, NBL, NROWS, + $ A, LDA, INFO ) + END IF +C + ELSE +C +C Undo scaling. +C + IF( ANRM.LT.SMLNUM ) THEN + CALL MB01QD( TYPE, M, N, KL, KU, SMLNUM, ANRM, NBL, NROWS, + $ A, LDA, INFO ) + ELSE IF( ANRM.GT.BIGNUM ) THEN + CALL MB01QD( TYPE, M, N, KL, KU, BIGNUM, ANRM, NBL, NROWS, + $ A, LDA, INFO ) + END IF + END IF +C + RETURN +C *** Last line of MB01PD *** + END Property changes on: trunk/octave-forge/main/control/devel/ctrbf/MB01PD.f ___________________________________________________________________ Added: svn:executable + * Added: trunk/octave-forge/main/control/devel/ctrbf/MB01QD.f =================================================================== --- trunk/octave-forge/main/control/devel/ctrbf/MB01QD.f (rev 0) +++ trunk/octave-forge/main/control/devel/ctrbf/MB01QD.f 2011-10-21 17:05:35 UTC (rev 8814) @@ -0,0 +1,334 @@ + SUBROUTINE MB01QD( TYPE, M, N, KL, KU, CFROM, CTO, NBL, NROWS, A, + $ LDA, INFO ) +C +C SLICOT RELEASE 5.0. +C +C Copyright (c) 2002-2010 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 multiply the M by N real matrix A by the real scalar CTO/CFROM. +C This is done without over/underflow as long as the final result +C CTO*A(I,J)/CFROM does not over/underflow. TYPE specifies that +C A may be full, (block) upper triangular, (block) lower triangular, +C (block) upper Hessenberg, or banded. +C +C ARGUMENTS +C +C Mode Parameters +C +C TYPE CHARACTER*1 +C TYPE indices the storage type of the input matrix. +C = 'G': A is a full matrix. +C = 'L': A is a (block) lower triangular matrix. +C = 'U': A is a (block) upper triangular matrix. +C = 'H': A is a (block) upper Hessenberg matrix. +C = 'B': A is a symmetric band matrix with lower bandwidth +C KL and upper bandwidth KU and with the only the +C lower half stored. +C = 'Q': A is a symmetric band matrix with lower bandwidth +C KL and upper bandwidth KU and with the only the +C upper half stored. +C = 'Z': A is a band matrix with lower bandwidth KL and +C upper bandwidth KU. +C +C Input/Output Parameters +C +C M (input) INTEGER +C The number of rows of the matrix A. M >= 0. +C +C N (input) INTEGER +C The number of columns of the matrix A. N >= 0. +C +C KL (input) INTEGER +C The lower bandwidth of A. Referenced only if TYPE = 'B', +C 'Q' or 'Z'. +C +C KU (input) INTEGER +C The upper bandwidth of A. Referenced only if TYPE = 'B', +C 'Q' or 'Z'. +C +C CFROM (input) DOUBLE PRECISION +C CTO (input) DOUBLE PRECISION +C The matrix A is multiplied by CTO/CFROM. A(I,J) is +C computed without over/underflow if the final result +C CTO*A(I,J)/CFROM can be represented without over/ +C underflow. CFROM must be nonzero. +C +C NBL (input) INTEGER +C The number of diagonal blocks of the matrix A, if it has a +C block structure. To specify that matrix A has no block +C structure, set NBL = 0. NBL >= 0. +C +C NROWS (input) INTEGER array, dimension max(1,NBL) +C NROWS(i) contains the number of rows and columns of the +C i-th diagonal block of matrix A. The sum of the values +C NROWS(i), for i = 1: NBL, should be equal to min(M,N). +C The array NROWS is not referenced if NBL = 0. +C +C A (input/output) DOUBLE PRECISION array, dimension (LDA,N) +C The matrix to be multiplied by CTO/CFROM. See TYPE for +C the storage type. +C +C LDA (input) INTEGER +C The leading dimension of the array A. LDA >= max(1,M). +C +C Error Indicator +C +C INFO INTEGER +C Not used in this implementation. +C +C METHOD +C +C Matrix A is multiplied by the real scalar CTO/CFROM, taking into +C account the specified storage mode of the matrix. +C MB01QD is a version of the LAPACK routine DLASCL, modified for +C dealing with block triangular, or block Hessenberg matrices. +C For efficiency, no tests of the input scalar parameters are +C performed. +C +C CONTRIBUTOR +C +C V. Sima, Katholieke Univ. Leuven, Belgium, Nov. 1996. +C +C ****************************************************************** +C +C .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) +C .. +C .. Scalar Arguments .. + CHARACTER TYPE + INTEGER INFO, KL, KU, LDA, M, N, NBL + DOUBLE PRECISION CFROM, CTO +C .. +C .. Array Arguments .. + INTEGER NROWS ( * ) + DOUBLE PRECISION A( LDA, * ) +C .. +C .. Local Scalars .. + LOGICAL DONE, NOBLC + INTEGER I, IFIN, ITYPE, J, JFIN, JINI, K, K1, K2, K3, + $ K4 + DOUBLE PRECISION BIGNUM, CFROM1, CFROMC, CTO1, CTOC, MUL, SMLNUM +C .. +C .. External Functions .. + LOGICAL LSAME + DOUBLE PRECISION DLAMCH + EXTERNAL LSAME, DLAMCH +C .. +C .. Intrinsic Functions .. + INTRINSIC ABS, MAX, MIN +C .. +C .. Executable Statements .. +C + IF( LSAME( TYPE, 'G' ) ) THEN + ITYPE = 0 + ELSE IF( LSAME( TYPE, 'L' ) ) THEN + ITYPE = 1 + ELSE IF( LSAME( TYPE, 'U' ) ) THEN + ITYPE = 2 + ELSE IF( LSAME( TYPE, 'H' ) ) THEN + ITYPE = 3 + ELSE IF( LSAME( TYPE, 'B' ) ) THEN + ITYPE = 4 + ELSE IF( LSAME( TYPE, 'Q' ) ) THEN + ITYPE = 5 + ELSE + ITYPE = 6 + END IF +C +C Quick return if possible. +C + IF( MIN( M, N ).EQ.0 ) + $ RETURN +C +C Get machine parameters. +C + SMLNUM = DLAMCH( 'S' ) + BIGNUM = ONE / SMLNUM +C + CFROMC = CFROM + CTOC = CTO +C + 10 CONTINUE + CFROM1 = CFROMC*SMLNUM + CTO1 = CTOC / BIGNUM + IF( ABS( CFROM1 ).GT.ABS( CTOC ) .AND. CTOC.NE.ZERO ) THEN + MUL = SMLNUM + DONE = .FALSE. + CFROMC = CFROM1 + ELSE IF( ABS( CTO1 ).GT.ABS( CFROMC ) ) THEN + MUL = BIGNUM + DONE = .FALSE. + CTOC = CTO1 + ELSE + MUL = CTOC / CFROMC + DONE = .TRUE. + END IF +C + NOBLC = NBL.EQ.0 +C + IF( ITYPE.EQ.0 ) THEN +C +C Full matrix +C + DO 30 J = 1, N + DO 20 I = 1, M + A( I, J ) = A( I, J )*MUL + 20 CONTINUE + 30 CONTINUE +C + ELSE IF( ITYPE.EQ.1 ) THEN +C + IF ( NOBLC ) THEN +C +C Lower triangular matrix +C + DO 50 J = 1, N + DO 40 I = J, M + A( I, J ) = A( I, J )*MUL + 40 CONTINUE + 50 CONTINUE +C + ELSE +C +C Block lower triangular matrix +C + JFIN = 0 + DO 80 K = 1, NBL + JINI = JFIN + 1 + JFIN = JFIN + NROWS( K ) + DO 70 J = JINI, JFIN + DO 60 I = JINI, M + A( I, J ) = A( I, J )*MUL + 60 CONTINUE + 70 CONTINUE + 80 CONTINUE + END IF +C + ELSE IF( ITYPE.EQ.2 ) THEN +C + IF ( NOBLC ) THEN +C +C Upper triangular matrix +C + DO 100 J = 1, N + DO 90 I = 1, MIN( J, M ) + A( I, J ) = A( I, J )*MUL + 90 CONTINUE + 100 CONTINUE +C + ELSE +C +C Block upper triangular matrix +C + JFIN = 0 + DO 130 K = 1, NBL + JINI = JFIN + 1 + JFIN = JFIN + NROWS( K ) + IF ( K.EQ.NBL ) JFIN = N + DO 120 J = JINI, JFIN + DO 110 I = 1, MIN( JFIN, M ) + A( I, J ) = A( I, J )*MUL + 110 CONTINUE + 120 CONTINUE + 130 CONTINUE + END IF +C + ELSE IF( ITYPE.EQ.3 ) THEN +C + IF ( NOBLC ) THEN +C +C Upper Hessenberg matrix +C + DO 150 J = 1, N + DO 140 I = 1, MIN( J+1, M ) + A( I, J ) = A( I, J )*MUL + 140 CONTINUE + 150 CONTINUE +C + ELSE +C +C Block upper Hessenberg matrix +C + JFIN = 0 + DO 180 K = 1, NBL + JINI = JFIN + 1 + JFIN = JFIN + NROWS( K ) +C + IF ( K.EQ.NBL ) THEN + JFIN = N + IFIN = N + ELSE + IFIN = JFIN + NROWS( K+1 ) + END IF +C + DO 170 J = JINI, JFIN + DO 160 I = 1, MIN( IFIN, M ) + A( I, J ) = A( I, J )*MUL + 160 CONTINUE + 170 CONTINUE + 180 CONTINUE + END IF +C + ELSE IF( ITYPE.EQ.4 ) THEN +C +C Lower half of a symmetric band matrix +C + K3 = KL + 1 + K4 = N + 1 + DO 200 J = 1, N + DO 190 I = 1, MIN( K3, K4-J ) + A( I, J ) = A( I, J )*MUL + 190 CONTINUE + 200 CONTINUE +C + ELSE IF( ITYPE.EQ.5 ) THEN +C +C Upper half of a symmetric band matrix +C + K1 = KU + 2 + K3 = KU + 1 + DO 220 J = 1, N + DO 210 I = MAX( K1-J, 1 ), K3 + A( I, J ) = A( I, J )*MUL + 210 CONTINUE + 220 CONTINUE +C + ELSE IF( ITYPE.EQ.6 ) THEN +C +C Band matrix +C + K1 = KL + KU + 2 + K2 = KL + 1 + K3 = 2*KL + KU + 1 + K4 = KL + KU + 1 + M + DO 240 J = 1, N + DO 230 I = MAX( K1-J, K2 ), MIN( K3, K4-J ) + A( I, J ) = A( I, J )*MUL + 230 CONTINUE + 240 CONTINUE +C + END IF +C + IF( .NOT.DONE ) + $ GO TO 10 +C + RETURN +C *** Last line of MB01QD *** + END Property changes on: trunk/octave-forge/main/control/devel/ctrbf/MB01QD.f ___________________________________________________________________ Added: svn:executable + * Added: trunk/octave-forge/main/control/devel/ctrbf/MB03OY.f =================================================================== --- trunk/octave-forge/main/control/devel/ctrbf/MB03OY.f (rev 0) +++ trunk/octave-forge/main/control/devel/ctrbf/MB03OY.f 2011-10-21 17:05:35 UTC (rev 8814) @@ -0,0 +1,388 @@ + SUBROUTINE MB03OY( M, N, A, LDA, RCOND, SVLMAX, RANK, SVAL, JPVT, + $ TAU, DWORK, INFO ) +C +C SLICOT RELEASE 5.0. +C +C Copyright (c) 2002-2010 NICONET e.V. +C +C This program is free software: you can redistribute it and/or +C modify it under the terms of the GNU General Public License as +C published by the Free Software Foundation, either version 2 of +C the License, or (at your option) any later version. +C +C This program is distributed in the hope that it will be useful, +C but WITHOUT ANY WARRANTY; without even the implied warranty of +C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +C GNU General Public License for more details. +C +C You should have received a copy of the GNU General Public License +C along with this program. If not, see +C <http://www.gnu.org/licenses/>. +C +C PURPOSE +C +C To compute a rank-revealing QR factorization of a real general +C M-by-N matrix A, which may be rank-deficient, and estimate its +C effective rank using incremental condition estimation. +C +C The routine uses a truncated QR factorization with column pivoting +C [ R11 R12 ] +C A * P = Q * R, where R = [ ], +C [ 0 R22 ] +C with R11 defined as the largest leading upper triangular submatrix +C whose estimated condition number is less than 1/RCOND. The order +C of R11, RANK, is the effective rank of A. Condition estimation is +C performed during the QR factorization process. Matrix R22 is full +C (but of small norm), or empty. +C +C MB03OY does not perform any scaling of the matrix A. +C +C ARGUMENTS +C +C Input/Output Parameters +C +C M (input) INTEGER +C The number of rows of the matrix A. M >= 0. +C +C N (input) INTEGER +C The number of columns of the matrix A. N >= 0. +C +C A (input/output) DOUBLE PRECISION array, dimension +C ( LDA, N ) +C On entry, the leading M-by-N part of this array must +C contain the given matrix A. +C On exit, the leading RANK-by-RANK upper triangular part +C of A contains the triangular factor R11, and the elements +C below the diagonal in the first RANK columns, with the +C array TAU, represent the orthogonal matrix Q as a product +C of RANK elementary reflectors. +C The remaining N-RANK columns contain the result of the +C QR factorization process used. +C +C LDA INTEGER +C The leading dimension of the array A. LDA >= max(1,M). +C +C RCOND (input) DOUBLE PRECISION +C RCOND is used to determine the effective rank of A, which +C is defined as the order of the largest leading triangular +C submatrix R11 in the QR factorization with pivoting of A, +C whose estimated condition number is less than 1/RCOND. +C 0 <= RCOND <= 1. +C NOTE that when SVLMAX > 0, the estimated rank could be +C less than that defined above (see SVLMAX). +C +C SVLMAX (input) DOUBLE PRECISION +C If A is a submatrix of another matrix B, and the rank +C decision should be related to that matrix, then SVLMAX +C should be an estimate of the largest singular value of B +C (for instance, the Frobenius norm of B). If this is not +C the case, the input value SVLMAX = 0 should work. +C SVLMAX >= 0. +C +C RANK (output) INTEGER +C The effective (estimated) rank of A, i.e., the order of +C the submatrix R11. +C +C SVAL (output) DOUBLE PRECISION array, dimension ( 3 ) +C The estimates of some of the singular values of the +C triangular factor R: +C SVAL(1): largest singular value of R(1:RANK,1:RANK); +C SVAL(2): smallest singular value of R(1:RANK,1:RANK); +C SVAL(3): smallest singular value of R(1:RANK+1,1:RANK+1), +C if RANK < MIN( M, N ), or of R(1:RANK,1:RANK), +C otherwise. +C If the triangular factorization is a rank-revealing one +C (which will be the case if the leading columns were well- +C conditioned), then SVAL(1) will also be an estimate for +C the largest singular value of A, and SVAL(2) and SVAL(3) +C will be estimates for the RANK-th and (RANK+1)-st singular +C values of A, respectively. +C By examining these values, one can confirm that the rank +C is well defined with respect to the chosen value of RCOND. +C The ratio SVAL(1)/SVAL(2) is an estimate of the condition +C number of R(1:RANK,1:RANK). +C +C JPVT (output) INTEGER array, dimension ( N ) +C If JPVT(i) = k, then the i-th column of A*P was the k-th +C column of A. +C +C TAU (output) DOUBLE PRECISION array, dimension ( MIN( M, N ) ) +C The leading RANK elements of TAU contain the scalar +C factors of the elementary reflectors. +C +C Workspace +C +C DWORK DOUBLE PRECISION array, dimension ( 3*N-1 ) +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 routine computes a truncated QR factorization with column +C pivoting of A, A * P = Q * R, with R defined above, and, +C during this process, finds the largest leading submatrix whose +C estimated condition number is less than 1/RCOND, taking the +C possible positive value of SVLMAX into account. This is performed +C using the LAPACK incremental condition estimation scheme and a +C slightly modified rank decision test. The factorization process +C stops when RANK has been determined. +C +C The matrix Q is represented as a product of elementary reflectors +C +C Q = H(1) H(2) . . . H(k), where k = rank <= min(m,n). +C +C Each H(i) has the form +C +C H = I - tau * v * v' +C +C where tau is a real scalar, and v is a real vector with +C v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in +C A(i+1:m,i), and tau in TAU(i). +C +C The matrix P is represented in jpvt as follows: If +C jpvt(j) = i +C then the jth column of P is the ith canonical unit vector. +C +C REFERENCES +C +C [1] Bischof, C.H. and P. Tang. +C Generalizing Incremental Condition Estimation. +C LAPACK Working Notes 32, Mathematics and Computer Science +C Division, Argonne National Laboratory, UT, CS-91-132, +C May 1991. +C +C [2] Bischof, C.H. and P. Tang. +C Robust Incremental Condition Estimation. +C LAPACK Working Notes 33, Mathematics and Computer Science +C Division, Argonne National Laboratory, UT, CS-91-133, +C May 1991. +C +C NUMERICAL ASPECTS +C +C The algorithm is backward stable. +C +C CONTRIBUTOR +C +C V. Sima, Katholieke Univ. Leuven, Belgium, Feb. 1998. +C +C REVISIONS +C +C V. Sima, Research Institute for Informatics, Bucharest, Jan. 2009. +C V. Sima, Jan. 2010, following Bujanovic and Drmac's suggestion. +C +C KEYWORDS +C +C Eigenvalue problem, matrix operations, orthogonal transformation, +C singular values. +C +C ****************************************************************** +C +C .. Parameters .. + INTEGER IMAX, IMIN + PARAMETER ( IMAX = 1, IMIN = 2 ) + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) +C .. Scalar Arguments .. + INTEGER INFO, LDA, M, N, RANK + DOUBLE PRECISION RCOND, SVLMAX +C .. Array Arguments .. + INTEGER JPVT( * ) + DOUBLE PRECISION A( LDA, * ), DWORK( * ), SVAL( 3 ), TAU( * ) +C .. +C .. Local Scalars .. + INTEGER I, ISMAX, ISMIN, ITEMP, J, MN, PVT + DOUBLE PRECISION AII, C1, C2, S1, S2, SMAX, SMAXPR, SMIN, + $ SMINPR, TEMP, TEMP2, TOLZ +C .. +C .. External Functions .. + INTEGER IDAMAX + DOUBLE PRECISION DLAMCH, DNRM2 + EXTERNAL DLAMCH, DNRM2, IDAMAX +C .. External Subroutines .. + EXTERNAL DLAIC1, DLARF, DLARFG, DSCAL, DSWAP, XERBLA +C .. Intrinsic Functions .. + INTRINSIC ABS, MAX, MIN, SQRT +C .. +C .. Executable Statements .. +C +C Test the input scalar arguments. +C + INFO = 0 + IF( M.LT.0 ) THEN + INFO = -1 + ELSE IF( N.LT.0 ) THEN + INFO = -2 + ELSE IF( LDA.LT.MAX( 1, M ) ) THEN + INFO = -4 + ELSE IF( RCOND.LT.ZERO .OR. RCOND.GT.ONE ) THEN + INFO = -5 + ELSE IF( SVLMAX.LT.ZERO ) THEN + INFO = -6 + END IF +C + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'MB03OY', -INFO ) + RETURN + END IF +C +C Quick return if possible. +C + MN = MIN( M, N ) + IF( MN.EQ.0 ) THEN + RANK = 0 + SVAL( 1 ) = ZERO + SVAL( 2 ) = ZERO + SVAL( 3 ) = ZERO + RETURN + END IF +C + TOLZ = SQRT( DLAMCH( 'Epsilon' ) ) + ISMIN = 1 + ISMAX = ISMIN + N +C +C Initialize partial column norms and pivoting vector. The first n +C elements of DWORK store the exact column norms. The already used +C leading part is then overwritten by the condition estimator. +C + DO 10 I = 1, N + DWORK( I ) = DNRM2( M, A( 1, I ), 1 ) + DWORK( N+I ) = DWORK( I ) + JPVT( I ) = I + 10 CONTINUE +C +C Compute factorization and determine RANK using incremental +C condition estimation. +C + RANK = 0 +C + 20 CONTINUE + IF( RANK.LT.MN ) THEN + I = RANK + 1 +C +C Determine ith pivot column and swap if necessary. +C + PVT = ( I-1 ) + IDAMAX( N-I+1, DWORK( I ), 1 ) +C + IF( PVT.NE.I ) THEN + CALL DSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 ) + ITEMP = JPVT( PVT ) + JPVT( PVT ) = JPVT( I ) + JPVT( I ) = ITEMP + DWORK( PVT ) = DWORK( I ) + DWORK( N+PVT ) = DWORK( N+I ) + END IF +C +C Save A(I,I) and generate elementary reflector H(i). +C + IF( I.LT.M ) THEN + AII = A( I, I ) + CALL DLARFG( M-I+1, A( I, I ), A( I+1, I ), 1, TAU( I ) ) + ELSE + TAU( M ) = ZERO + END IF +C + IF( RANK.EQ.0 ) THEN +C +C Initialize; exit if matrix is zero (RANK = 0). +C + SMAX = ABS( A( 1, 1 ) ) + IF ( SMAX.EQ.ZERO ) THEN + SVAL( 1 ) = ZERO + SVAL( 2 ) = ZERO + SVAL( 3 ) = ZERO + RETURN + END IF + SMIN = SMAX + SMAXPR = SMAX + SMINPR = SMIN + C1 = ONE + C2 = ONE + ELSE +C +C One step of incremental condition estimation. +C + CALL DLAIC1( IMIN, RANK, DWORK( ISMIN ), SMIN, A( 1, I ), + $ A( I, I ), SMINPR, S1, C1 ) + CALL DLAIC1( IMAX, RANK, DWORK( ISMAX ), SMAX, A( 1, I ), + $ A( I, I ), SMAXPR, S2, C2 ) + END IF +C + IF( SVLMAX*RCOND.LE.SMAXPR ) THEN + IF( SVLMAX*RCOND.LE.SMINPR ) THEN + IF( SMAXPR*RCOND.LE.SMINPR ) THEN +C +C Continue factorization, as rank is at least RANK. +C + IF( I.LT.N ) THEN +C +C Apply H(i) to A(i:m,i+1:n) from the left. +C + AII = A( I, I ) + A( I, I ) = ONE + CALL DLARF( 'Left', M-I+1, N-I, A( I, I ), 1, + $ TAU( I ), A( I, I+1 ), LDA, + $ DWORK( 2*N+1 ) ) + A( I, I ) = AII + END IF +C +C Update partial column norms. +C + DO 30 J = I + 1, N + IF( DWORK( J ).NE.ZERO ) THEN + TEMP = ABS( A( I, J ) ) / DWORK( J ) + TEMP = MAX( ( ONE + TEMP )*( ONE - TEMP ), ZERO) + TEMP2 = TEMP*( DWORK( J ) / DWORK( N+J ) )**2 + IF( TEMP2.LE.TOLZ ) THEN + IF( M-I.GT.0 ) THEN + DWORK( J ) = DNRM2( M-I, A( I+1, J ), 1 ) + DWORK( N+J ) = DWORK( J ) + ELSE + DWORK( J ) = ZERO + DWORK( N+J ) = ZERO + END IF + ELSE + DWORK( J ) = DWORK( J )*SQRT( TEMP ) + END IF + END IF + 30 CONTINUE +C + DO 40 I = 1, RANK + DWORK( ISMIN+I-1 ) = S1*DWORK( ISMIN+I-1 ) + DWORK( ISMAX+I-1 ) = S2*DWORK( ISMAX+I-1 ) + 40 CONTINUE +C + DWORK( ISMIN+RANK ) = C1 + DWORK( ISMAX+RANK ) = C2 + SMIN = SMINPR + SMAX = SMAXPR + RANK = RANK + 1 + GO TO 20 + END IF + END IF + END IF + END IF +C +C Restore the changed part of the (RANK+1)-th column and set SVAL. +C + IF ( RANK.LT.N ) THEN + IF ( I.LT.M ) THEN + CALL DSCAL( M-I, -A( I, I )*TAU( I ), A( I+1, I ), 1 ) + A( I, I ) = AII + END IF + END IF + IF ( RANK.EQ.0 ) THEN + SMIN = ZERO + SMINPR = ZERO + END IF + SVAL( 1 ) = SMAX + SVAL( 2 ) = SMIN + SVAL( 3 ) = SMINPR +C + RETURN +C *** Last line of MB03OY *** + END Property changes on: trunk/octave-forge/main/control/devel/ctrbf/MB03OY.f ___________________________________________________________________ Added: svn:executable + * Added: trunk/octave-forge/main/control/devel/ctrbf/TB01UD.dat =================================================================== --- trunk/octave-forge/main/control/devel/ctrbf/TB01UD.dat (rev 0) +++ trunk/octave-forge/main/control/devel/ctrbf/TB01UD.dat 2011-10-21 17:05:35 UTC (rev 8814) @@ -0,0 +1,9 @@ + TB01UD EXAMPLE PROGRAM DATA + 3 2 2 0.0 I + -1.0 0.0 0.0 + -2.0 -2.0 -2.0 + -1.0 0.0 -3.0 + 1.0 0.0 0.0 + 0.0 2.0 1.0 + 0.0 2.0 1.0 + 1.0 0.0 0.0 Property changes on: trunk/octave-forge/main/control/devel/ctrbf/TB01UD.dat ___________________________________________________________________ Added: svn:executable + * Added: trunk/octave-forge/main/control/devel/ctrbf/TB01UD.f =================================================================== --- trunk/octave-forge/main/control/devel/ctrbf/TB01UD.f (rev 0) +++ trunk/octave-forge/main/control/devel/ctrbf/TB01UD.f 2011-10-21 17:05:35 UTC (rev 8814) @@ -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-2010 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 Property changes on: trunk/octave-forge/main/control/devel/ctrbf/TB01UD.f ___________________________________________________________________ Added: svn:executable + * Added: trunk/octave-forge/main/control/devel/ctrbf/TB01UD.res =================================================================== --- trunk/octave-forge/main/control/devel/ctrbf/TB01UD.res (rev 0) +++ trunk/octave-forge/main/control/devel/ctrbf/TB01UD.res 2011-10-21 17:05:35 UTC (rev 8814) @@ -0,0 +1,25 @@ + TB01UD EXAMPLE PROGRAM RESULTS + + The order of the controllable state-space representation = 2 + + The transformed state dynamics matrix of a controllable realization is + -3.0000 2.2361 + 0.0000 -1.0000 + + and the dimensions of its diagonal blocks are + 2 + + The transformed input/state matrix B of a controllable realization is + 0.0000 -2.2361 + 1.0000 0.0000 + + The transformed output/state matrix C of a controllable realization is + -2.2361 0.0000 + 0.0000 1.0000 + + The controllability index of the transformed system representation = 1 + + The similarity transformation matrix Z is + 0.0000 1.0000 0.0000 + -0.8944 0.0000 -0.4472 + -0.4472 0.0000 0.8944 Property changes on: trunk/octave-forge/main/control/devel/ctrbf/TB01UD.res ___________________________________________________________________ Added: svn:executable + * Added: trunk/octave-forge/main/control/devel/ctrbf/common.cc =================================================================== --- trunk/octave-forge/main/control/devel/ctrbf/common.cc (rev 0) +++ trunk/octave-forge/main/control/devel/ctrbf/common.cc 2011-10-21 17:05:35 UTC (rev 8814) @@ -0,0 +1,53 @@ +/* + +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 LTI Syncope. If not, see <http://www.gnu.org/licenses/>. + +Common code for oct-files. + +Author: Lukas Reichlin <luk...@gm...> +Created: April 2010 +Version: 0.1 + +*/ + + +int max (int a, int b) +{ + if (a > b) + return a; + else + return b; +} + +int max (int a, int b, int c) +{ + return max (max (a, b), c); +} + +int max (int a, int b, int c, int d) +{ + return max (max (a, b), max (c, d)); +} + +int min (int a, int b) +{ + if (a < b) + return a; + else + return b; +} Added: trunk/octave-forge/main/control/devel/ctrbf/makefile_ctrbf.m =================================================================== --- trunk/octave-forge/main/control/devel/ctrbf/makefile_ctrbf.m (rev 0) +++ trunk/octave-forge/main/control/devel/ctrbf/makefile_ctrbf.m 2011-10-21 17:05:35 UTC (rev 8814) @@ -0,0 +1,2 @@ +mkoctfile sltb01ud.cc \ + TB01UD.f MB01PD.f MB03OY.f MB01QD.f \ No newline at end of file Added: trunk/octave-forge/main/control/devel/ctrbf/sltb01ud.cc =================================================================== --- trunk/octave-forge/main/control/devel/ctrbf/sltb01ud.cc (rev 0) +++ trunk/octave-forge/main/control/devel/ctrbf/sltb01ud.cc 2011-10-21 17:05:35 UTC (rev 8814) @@ -0,0 +1,140 @@ +/* + +Copyright (C) 2011 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 LTI Syncope. If not, see <http://www.gnu.org/licenses/>. + +Orthogonal canonical form. +Uses SLICOT TB01UD by courtesy of NICONET e.V. +<http://www.slicot.org> + +Author: Lukas Reichlin <luk...@gm...> +Created: October 2011 +Version: 0.1 + +*/ + +#include <octave/oct.h> +#include <f77-fcn.h> +#include "common.cc" + +extern "C" +{ + int F77_FUNC (tb01ud, TB01UD) + (char& JOBZ, + int& N, int& M, int& P, + double* A, int& LDA, + double* B, int& LDB, + double* C, int& LDC, + int& NCONT, int& INDCON, + int* NBLK, + double* Z, int& LDZ, + double* TAU, + double& TOL, + int* IWORK, + double* DWORK, int& LDWORK, + int& INFO); +} + +DEFUN_DLD (sltb01ud, args, nargout, + "-*- ... [truncated message content] |
From: <par...@us...> - 2011-11-02 17:08:24
|
Revision: 8938 http://octave.svn.sourceforge.net/octave/?rev=8938&view=rev Author: paramaniac Date: 2011-11-02 17:08:18 +0000 (Wed, 02 Nov 2011) Log Message: ----------- control: update draft code Modified Paths: -------------- trunk/octave-forge/main/control/devel/ctrbf/test_ctrbf.m trunk/octave-forge/main/control/devel/ctrbf.m Modified: trunk/octave-forge/main/control/devel/ctrbf/test_ctrbf.m =================================================================== --- trunk/octave-forge/main/control/devel/ctrbf/test_ctrbf.m 2011-11-02 17:03:06 UTC (rev 8937) +++ trunk/octave-forge/main/control/devel/ctrbf/test_ctrbf.m 2011-11-02 17:08:18 UTC (rev 8938) @@ -12,6 +12,8 @@ [ac, bc, cc, z, ncont] = sltb01ud (a, b, c, 0.0) +z\a*z + a = [ 1 1 4 -2 ]; Modified: trunk/octave-forge/main/control/devel/ctrbf.m =================================================================== --- trunk/octave-forge/main/control/devel/ctrbf.m 2011-11-02 17:03:06 UTC (rev 8937) +++ trunk/octave-forge/main/control/devel/ctrbf.m 2011-11-02 17:08:18 UTC (rev 8938) @@ -1,18 +1,20 @@ -## Copyright (C) 2010 Benjamin Fernandez <ma...@be...> +## Copyright (C) 2010 Benjamin Fernandez +## Copyright (C) 2011 Lukas F. Reichlin ## -## This program is free software; you can redistribute it and/or modify +## 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 2 of the License, or +## the Free Software Foundation, either version 3 of the License, or ## (at your option) any later version. -## -## This program is distributed in the hope that it will be useful, +## +## 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 Octave; see the file COPYING. If not, see -## <http://www.gnu.org/licenses/>. +## along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @deftypefn{Function File} {[@var{Abar}, @var{Bbar}, @var{Cbar}, @var{T}, @var{K}] =} ctrbf (@var{A}, @var{B}, @var{C}) @@ -43,28 +45,53 @@ ## containing the number of controllable states. ## @end deftypefn -## Author: Benjamin Fernandez <benjas@benjas-laptop> +## Author: Benjamin Fernandez <ma...@be...> ## Created: 2010-04-30 +## Version: 0.1 -function [Abar, Bbar, Cbar, T, K] = ctrbf (A, B, C, TOL) +function [ac, bc, cc, z, ncont] = ctrbf (a, b, c, tol = []) if (nargin < 3 || nargin > 4) print_usage (); endif - if (nargin == 3) - TOL = length (A) * norm (A,1) * eps; + if (isempty (tol)) + tol = 0; # default tolerance + elseif (! is_real_scalar (tol)) + error ("ctrbf: tol must be a real scalar"); endif - Co = ctrb (A, B); - [nrc, ncc] = size (Co); - rco = rank (Co, TOL); - lr = nrc - rco; - [U, S, V] = svd (Co); - K = U(:, 1:rco); # Basis column space - T = U; # [B orth(B)] - Abar = T \ A * T; - Bbar = T \ B; - Cbar = C * T; + [ac, bc, cc, z, ncont] = sltb01ud (a, b, c, 0.0); -endfunction \ No newline at end of file +endfunction + +%!shared Ao, Bo, Co, Zo, Ae, Be, Ce, Ze, NCONT +%! A = [ -1.0 0.0 0.0 +%! -2.0 -2.0 -2.0 +%! -1.0 0.0 -3.0 ]; +%! +%! B = [ 1.0 0.0 0.0 +%! 0.0 2.0 1.0 ].'; +%! +%! C = [ 0.0 2.0 1.0 +%! 1.0 0.0 0.0 ]; +%! +%! [Ao, Bo, Co, Zo, NCONT] = ctrbf (A, B, C); +%! +%! Ae = [ -3.0000 2.2361 +%! 0.0000 -1.0000 ]; +%! +%! Be = [ 0.0000 -2.2361 +%! 1.0000 0.0000 ]; +%! +%! Ce = [ -2.2361 0.0000 +%! 0.0000 1.0000 ]; +%! +%! Ze = [ 0.0000 1.0000 0.0000 +%! -0.8944 0.0000 -0.4472 +%! -0.4472 0.0000 0.8944 ]; +%! +%!assert (Ao(1:NCONT, 1:NCONT), Ae, 1e-4); +%!assert (Bo(1:NCONT, :), Be, 1e-4); +%!assert (Co(:, 1:NCONT), Ce, 1e-4); +%!assert (Zo, Ze, 1e-4); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-11-05 16:22:20
|
Revision: 9002 http://octave.svn.sourceforge.net/octave/?rev=9002&view=rev Author: paramaniac Date: 2011-11-05 16:22:14 +0000 (Sat, 05 Nov 2011) Log Message: ----------- control: rename makefile Added Paths: ----------- trunk/octave-forge/main/control/devel/makefile_control.m Removed Paths: ------------- trunk/octave-forge/main/control/devel/makefile_all.m Deleted: trunk/octave-forge/main/control/devel/makefile_all.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_all.m 2011-11-05 16:19:38 UTC (rev 9001) +++ trunk/octave-forge/main/control/devel/makefile_all.m 2011-11-05 16:22:14 UTC (rev 9002) @@ -1,31 +0,0 @@ -## ============================================================================== -## Developer Makefile for OCT-files -## ============================================================================== -## USAGE: * fetch control from Octave-Forge by svn -## * add control/inst, control/src and control/devel to your Octave path -## * run makefile_* -## ============================================================================== - -homedir = pwd (); -develdir = fileparts (which ("makefile_all")); -srcdir = [develdir, "/../src"]; -cd (srcdir); - -makefile_chol -makefile_conversions -makefile_h2syn -makefile_hankel -makefile_helpers -makefile_hinfsyn -makefile_lqr -makefile_lyap -makefile_minreal -makefile_ncfsyn -makefile_norm -makefile_place -makefile_scale -makefile_staircase -makefile_tustin -makefile_zero - -cd (homedir); Copied: trunk/octave-forge/main/control/devel/makefile_control.m (from rev 9000, trunk/octave-forge/main/control/devel/makefile_all.m) =================================================================== --- trunk/octave-forge/main/control/devel/makefile_control.m (rev 0) +++ trunk/octave-forge/main/control/devel/makefile_control.m 2011-11-05 16:22:14 UTC (rev 9002) @@ -0,0 +1,31 @@ +## ============================================================================== +## Developer Makefile for OCT-files +## ============================================================================== +## USAGE: * fetch control from Octave-Forge by svn +## * add control/inst, control/src and control/devel to your Octave path +## * run makefile_control +## ============================================================================== + +homedir = pwd (); +develdir = fileparts (which ("makefile_control")); +srcdir = [develdir, "/../src"]; +cd (srcdir); + +makefile_chol +makefile_conversions +makefile_h2syn +makefile_hankel +makefile_helpers +makefile_hinfsyn +makefile_lqr +makefile_lyap +makefile_minreal +makefile_ncfsyn +makefile_norm +makefile_place +makefile_scale +makefile_staircase +makefile_tustin +makefile_zero + +cd (homedir); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-09-15 14:36:28
|
Revision: 11027 http://octave.svn.sourceforge.net/octave/?rev=11027&view=rev Author: paramaniac Date: 2012-09-15 14:36:20 +0000 (Sat, 15 Sep 2012) Log Message: ----------- control: add draft code for multiplot time responses Added Paths: ----------- trunk/octave-forge/main/control/devel/__time_response_2__.m trunk/octave-forge/main/control/devel/step2.m Added: trunk/octave-forge/main/control/devel/__time_response_2__.m =================================================================== --- trunk/octave-forge/main/control/devel/__time_response_2__.m (rev 0) +++ trunk/octave-forge/main/control/devel/__time_response_2__.m 2012-09-15 14:36:20 UTC (rev 11027) @@ -0,0 +1,363 @@ +## Copyright (C) 2009, 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 LTI Syncope. If not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## Common code for the time response functions step, impulse and initial. + +## Author: Lukas Reichlin <luk...@gm...> +## Created: October 2009 +## Version: 0.2 + +% function [y, t, x_arr] = __time_response_2__ (sys, resptype, plotflag, tfinal, dt, x0, sysname) +function [y, t, x_arr] = __time_response_2__ (resptype, args) + + sys_idx = cellfun (@isa, args, {"lti"}); # look for LTI models + sys_cell = cellfun (@ss, args(sys_idx)); # system must be proper + + if (! size_equal (sys_cell{:})) + error ("%s: models must have equal sizes", resptype); + endif + + tmp = cellfun (@is_real_matrix, args); + vec_idx = find (tmp); + n_vec = length (vec_idx) + + if (n_vec >= 1) + arg = args{vec_idx(1)}; + + endif + + ## extract tfinal/t, dt, x0 + + [tfinal, dt] = cellfun (@__sim_horizon__, sys_cell, {tfinal}, dt); + + tfinal = max (tfinal); + + % __sim_horizon__ (sys, tfinal, dt); + + + hier sim_horizon + + ct_idx = cellfun (@isct, sys_cell) + sys_dt_cell = sys_cell; + tmp = (@c2d, sys_cell(ct_idx), dt, {"zoh"}, "uniformoutput", false) + sys_dt_cell(ct_idx) = tmp; + + +endfunction +%{ + if (! isa (sys, "ss")) + sys = ss (sys); # sys must be proper + endif + + if (is_real_vector (tfinal) && length (tfinal) > 1) # time vector t passed + dt = tfinal(2) - tfinal(1); # assume that t is regularly spaced + tfinal = tfinal(end); + endif + + [A, B, C, D, tsam] = ssdata (sys); + + discrete = ! isct (sys); # static gains are treated as analog systems + tsam = abs (tsam); # use 1 second if tsam is unspecified (-1) + + if (discrete) + if (! isempty (dt)) + warning ("time_response: argument dt has no effect on sampling time of discrete system"); + endif + + dt = tsam; + endif + + [tfinal, dt] = __sim_horizon__ (A, discrete, tfinal, dt); + + if (! discrete) + sys = c2d (sys, dt, "zoh"); + endif + + [F, G] = ssdata (sys); # matrices C and D don't change + + n = rows (F); # number of states + m = columns (G); # number of inputs + p = rows (C); # number of outputs + + ## time vector + t = reshape (0 : dt : tfinal, [], 1); + l_t = length (t); +%} + +function [y, x_arr] = __initial_response__ (sys, sys_dt, t, x0) + + [A, B, C, D] = ssdata (sys); + [F, G] = ssdata (sys_dt); + + n = rows (F); # number of states + m = columns (G); # number of inputs + p = rows (C); # number of outputs + l_t = length (t); + + ## preallocate memory + y = zeros (l_t, p); + x_arr = zeros (l_t, n); + + ## initial conditions + x = reshape (x0, [], 1); # make sure that x is a column vector + + if (n != length (x0) || ! is_real_vector (x0)) + error ("initial: x0 must be a real vector with %d elements", n); + endif + + ## simulation + for k = 1 : l_t + y(k, :) = C * x; + x_arr(k, :) = x; + x = F * x; + endfor + +endfunction + + +function [y, x_arr] = __step_response__ (sys_dt, t) + + [F, G, C, D] = ssdata (sys_dt); + + n = rows (F); # number of states + m = columns (G); # number of inputs + p = rows (C); # number of outputs + l_t = length (t); + + ## preallocate memory + y = zeros (l_t, p, m); + x_arr = zeros (l_t, n, m); + + for j = 1 : m # for every input channel + ## initial conditions + x = zeros (n, 1); + u = zeros (m, 1); + u(j) = 1; + + ## simulation + for k = 1 : l_t + y(k, :, j) = C * x + D * u; + x_arr(k, :, j) = x; + x = F * x + G * u; + endfor + endfor + +endfunction + + +function [y, x_arr] = __impulse_response__ (sys, sys_dt, t) + + [A, B, C, D, dt] = ssdata (sys); + [F, G] = ssdata (sys_dt); + discrete = ! isct (sys); + + n = rows (F); # number of states + m = columns (G); # number of inputs + p = rows (C); # number of outputs + l_t = length (t); + + ## preallocate memory + y = zeros (l_t, p, m); + x_arr = zeros (l_t, n, m); + + for j = 1 : m # for every input channel + ## initial conditions + u = zeros (m, 1); + u(j) = 1; + + if (discrete) + x = zeros (n, 1); # zero by definition + y(1, :, j) = D * u / dt; + x_arr(1, :, j) = x; + x = G * u / dt; + else + x = B * u; # B, not G! + y(1, :, j) = C * x; + x_arr(1, :, j) = x; + x = F * x; + endif + + ## simulation + for k = 2 : l_t + y (k, :, j) = C * x; + x_arr(k, :, j) = x; + x = F * x; + endfor + endfor + + if (discrete) + y *= dt; + x_arr *= dt; + endif + +endfunction + +%{ + if (plotflag) # display plot + + ## TODO: Set correct titles, especially for multi-input systems + + stable = isstable (sys); + outname = get (sys, "outname"); + outname = __labels__ (outname, "y_"); + + if (strcmp (resptype, "initial")) + cols = 1; + else + cols = m; + endif + + if (discrete) # discrete system + for k = 1 : p + for j = 1 : cols + + subplot (p, cols, (k-1)*cols+j); + + if (stable) + stairs (t, [y(:, k, j), yfinal(k, j) * ones(l_t, 1)]); + else + stairs (t, y(:, k, j)); + endif + + grid ("on"); + + if (k == 1 && j == 1) + title (str); + endif + + if (j == 1) + ylabel (sprintf ("Amplitude %s", outname{k})); + endif + + endfor + endfor + + xlabel ("Time [s]"); + + else # continuous system + for k = 1 : p + for j = 1 : cols + + subplot (p, cols, (k-1)*cols+j); + + if (stable) + plot (t, [y(:, k, j), yfinal(k, j) * ones(l_t, 1)]); + else + plot (t, y(:, k, j)); + endif + + grid ("on"); + + if (k == 1 && j == 1) + title (str); + endif + + if (j == 1) + ylabel (sprintf ("Amplitude %s", outname{k})); + endif + + endfor + endfor + + xlabel ("Time [s]"); + + endif + endif + +endfunction +%} + + +% function [tfinal, dt] = __sim_horizon__ (A, discrete, tfinal, Ts) +function [tfinal, dt] = __sim_horizon__ (sys, tfinal, Ts) + + ## code based on __stepimp__.m of Kai P. Mueller and A. Scottedward Hodel + + TOL = 1.0e-10; # values below TOL are assumed to be zero + N_MIN = 50; # min number of points + N_MAX = 2000; # max number of points + N_DEF = 1000; # default number of points + T_DEF = 10; # default simulation time + + ev = pole (sys); + n = length (ev); + + if (discrete) + ## perform bilinear transformation on poles in z + for k = 1 : n + pol = ev(k); + if (abs (pol + 1) < TOL) + ev(k) = 0; + else + ev(k) = 2 / Ts * (pol - 1) / (pol + 1); + endif + endfor + endif + + ## remove poles near zero from eigenvalue array ev + nk = n; + for k = 1 : n + if (abs (real (ev(k))) < TOL) + ev(k) = 0; + nk -= 1; + endif + endfor + + if (nk == 0) + if (isempty (tfinal)) + tfinal = T_DEF; + endif + + if (! discrete) + dt = tfinal / N_DEF; + endif + else + ev = ev(find (ev)); + ev_max = max (abs (ev)); + + if (! discrete) + dt = 0.2 * pi / ev_max; + endif + + if (isempty (tfinal)) + ev_min = min (abs (real (ev))); + tfinal = 5.0 / ev_min; + + ## round up + yy = 10^(ceil (log10 (tfinal)) - 1); + tfinal = yy * ceil (tfinal / yy); + endif + + if (! discrete) + N = tfinal / dt; + + if (N < N_MIN) + dt = tfinal / N_MIN; + endif + + if (N > N_MAX) + dt = tfinal / N_MAX; + endif + endif + endif + + if (! isempty (Ts)) # catch case cont. system with dt specified + dt = Ts; + endif + +endfunction Added: trunk/octave-forge/main/control/devel/step2.m =================================================================== --- trunk/octave-forge/main/control/devel/step2.m (rev 0) +++ trunk/octave-forge/main/control/devel/step2.m 2012-09-15 14:36:20 UTC (rev 11027) @@ -0,0 +1,80 @@ +## Copyright (C) 2009 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 LTI Syncope. If not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn{Function File} {[@var{y}, @var{t}, @var{x}] =} step (@var{sys}) +## @deftypefnx{Function File} {[@var{y}, @var{t}, @var{x}] =} step (@var{sys}, @var{t}) +## @deftypefnx{Function File} {[@var{y}, @var{t}, @var{x}] =} step (@var{sys}, @var{tfinal}) +## @deftypefnx{Function File} {[@var{y}, @var{t}, @var{x}] =} step (@var{sys}, @var{tfinal}, @var{dt}) +## Step response of LTI system. +## If no output arguments are given, the response is printed on the screen. +## +## @strong{Inputs} +## @table @var +## @item sys +## LTI model. +## @item t +## Time vector. Should be evenly spaced. If not specified, it is calculated by +## the poles of the system to reflect adequately the response transients. +## @item tfinal +## Optional simulation horizon. If not specified, it is calculated by +## the poles of the system to reflect adequately the response transients. +## @item dt +## Optional sampling time. Be sure to choose it small enough to capture transient +## phenomena. If not specified, it is calculated by the poles of the system. +## @end table +## +## @strong{Outputs} +## @table @var +## @item y +## Output response array. Has as many rows as time samples (length of t) +## and as many columns as outputs. +## @item t +## Time row vector. +## @item x +## State trajectories array. Has @code{length (t)} rows and as many columns as states. +## @end table +## +## @seealso{impulse, initial, lsim} +## @end deftypefn + +## Author: Lukas Reichlin <luk...@gm...> +## Created: October 2009 +## Version: 0.1 + +% function [y_r, t_r, x_r] = step2 (sys, tfinal = [], dt = []) +function [y_r, t_r, x_r] = step2 (varargin) + + if (nargin == 0) + print_usage (); + endif + + %tmp = cellfun (@isa, varargin, {"lti"}); + %sys_idx = find (tmp); + + + %sys_names = arrayfun (@inputname, sys_idx); + + [y, t, x] = __time_response__ ("step", varargin, ! nargout); + + if (nargout) + y_r = y; + t_r = t; + x_r = x; + endif + +endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-12-03 06:54:54
|
Revision: 9254 http://octave.svn.sourceforge.net/octave/?rev=9254&view=rev Author: paramaniac Date: 2011-12-03 06:54:48 +0000 (Sat, 03 Dec 2011) Log Message: ----------- control: add test cases Modified Paths: -------------- trunk/octave-forge/main/control/devel/test_tf2dss.m Added Paths: ----------- trunk/octave-forge/main/control/devel/bug_minreal.m trunk/octave-forge/main/control/devel/bug_plotting.m Added: trunk/octave-forge/main/control/devel/bug_minreal.m =================================================================== --- trunk/octave-forge/main/control/devel/bug_minreal.m (rev 0) +++ trunk/octave-forge/main/control/devel/bug_minreal.m 2011-12-03 06:54:48 UTC (rev 9254) @@ -0,0 +1,20 @@ +P = ss (-1, 1, 1, 0) +Pi = inv (P) +minreal (Pi) + +%{ + ** On entry to TG01JD parameter number 20 had an illegal value +error: sltg01jd: exception encountered in Fortran subroutine tg01jd_ +error: called from: +error: /Users/lukas/control/inst/@ss/__minreal__.m at line 39, column 19 +error: /Users/lukas/control/inst/@lti/minreal.m at line 38, column 7 +%} + + +ss (inv (tf (Boeing707 ))) + +%{ +octave(211,0x7fff70ddbcc0) malloc: *** error for object 0x108d0e2d8: incorrect checksum for freed object - object was probably modified after being freed. +*** set a breakpoint in malloc_error_break to debug +panic: Abort trap -- stopping myself... +%} \ No newline at end of file Added: trunk/octave-forge/main/control/devel/bug_plotting.m =================================================================== --- trunk/octave-forge/main/control/devel/bug_plotting.m (rev 0) +++ trunk/octave-forge/main/control/devel/bug_plotting.m 2011-12-03 06:54:48 UTC (rev 9254) @@ -0,0 +1,7 @@ +s = tf ('s'); +G = 1/(s^2+1) + +[mag, pha] = bode (G) % this works OK + +figure +bode (G) % fltk and gnuplot both crash \ No newline at end of file Modified: trunk/octave-forge/main/control/devel/test_tf2dss.m =================================================================== --- trunk/octave-forge/main/control/devel/test_tf2dss.m 2011-12-02 14:50:14 UTC (rev 9253) +++ trunk/octave-forge/main/control/devel/test_tf2dss.m 2011-12-03 06:54:48 UTC (rev 9254) @@ -13,4 +13,41 @@ P = ss (sys) -ss (inv (tf (Boeing707 ))) \ No newline at end of file +ss (inv (tf (Boeing707 ))) + +%{ +octave:2> test_tf2dss + +Transfer function "sys" from input "u1" to output ... + + s^2 + 5 s + 6 + y1: ------------- + s^2 + 5 s + 7 + + s^2 + 8 s + 6 + y2: ------------- + s + 7 + +Transfer function "sys" from input "u2" to output ... + + y1: s + 2 + + s^2 + 3 s + 2 + y2: ------------- + s^2 + 5 s + 5 + +Transfer function "sys" from input "u3" to output ... + + s^4 + 2 s^3 - s^2 + 2 s + 2 + y1: --------------------------- + s^2 - s + 1 + + 3 s^2 + 4 s + y2: ----------- + s + 1 + +Continuous-time model. +octave(211,0x7fff70ddbcc0) malloc: *** error for object 0x108d0e2d8: incorrect checksum for freed object - object was probably modified after being freed. +*** set a breakpoint in malloc_error_break to debug +panic: Abort trap -- stopping myself... +%} \ No newline at end of file This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-12-15 12:44:37
|
Revision: 9388 http://octave.svn.sourceforge.net/octave/?rev=9388&view=rev Author: paramaniac Date: 2011-12-15 12:44:26 +0000 (Thu, 15 Dec 2011) Log Message: ----------- control: remove mac-specific compiler flags from developer makefiles (new macports octave-devel version) Modified Paths: -------------- trunk/octave-forge/main/control/devel/makefile_hankel.m trunk/octave-forge/main/control/devel/makefile_staircase.m Modified: trunk/octave-forge/main/control/devel/makefile_hankel.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_hankel.m 2011-12-15 10:42:14 UTC (rev 9387) +++ trunk/octave-forge/main/control/devel/makefile_hankel.m 2011-12-15 12:44:26 UTC (rev 9388) @@ -14,8 +14,7 @@ srcdir = [develdir, "/../src"]; cd (srcdir); -mkoctfile "-Wl,-framework" "-Wl,vecLib" \ - slab13ad.cc \ +mkoctfile slab13ad.cc \ AB13AD.f TB01ID.f TB01KD.f AB13AX.f MA02DD.f \ MB03UD.f TB01LD.f SB03OU.f MB03QX.f select.f \ SB03OT.f MB03QD.f MB04ND.f MB04OD.f MB03QY.f \ Modified: trunk/octave-forge/main/control/devel/makefile_staircase.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_staircase.m 2011-12-15 10:42:14 UTC (rev 9387) +++ trunk/octave-forge/main/control/devel/makefile_staircase.m 2011-12-15 12:44:26 UTC (rev 9388) @@ -15,18 +15,15 @@ cd (srcdir); ## staircase form using orthogonal transformations -mkoctfile "-Wl,-framework" "-Wl,vecLib" \ - slab01od.cc \ +mkoctfile slab01od.cc \ AB01OD.f AB01ND.f MB03OY.f MB01PD.f MB01QD.f ## controllability staircase form of descriptor state-space models -mkoctfile "-Wl,-framework" "-Wl,vecLib" \ - sltg01hd.cc \ +mkoctfile sltg01hd.cc \ TG01HD.f TG01HX.f ## observability staircase form of descriptor state-space models -mkoctfile "-Wl,-framework" "-Wl,vecLib" \ - sltg01id.cc \ +mkoctfile sltg01id.cc \ TG01ID.f TB01XD.f MA02CD.f AB07MD.f TG01HX.f \ MA02BD.f This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-12-15 16:00:36
|
Revision: 9391 http://octave.svn.sourceforge.net/octave/?rev=9391&view=rev Author: paramaniac Date: 2011-12-15 16:00:30 +0000 (Thu, 15 Dec 2011) Log Message: ----------- control: remove obsolete remarks Modified Paths: -------------- trunk/octave-forge/main/control/devel/makefile_hankel.m trunk/octave-forge/main/control/devel/makefile_staircase.m Modified: trunk/octave-forge/main/control/devel/makefile_hankel.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_hankel.m 2011-12-15 16:00:14 UTC (rev 9390) +++ trunk/octave-forge/main/control/devel/makefile_hankel.m 2011-12-15 16:00:30 UTC (rev 9391) @@ -4,9 +4,6 @@ ## USAGE: * fetch control from Octave-Forge by svn ## * add control/inst, control/src and control/devel to your Octave path ## * run makefile_* -## NOTES: * The option "-Wl,-framework" "-Wl,vecLib" is needed for MacPorts' -## octave-devel @3.3.52_1+gcc44 on MacOS X 10.6.4. However, this option -## breaks other platforms. See MacPorts Ticket #26640. ## ============================================================================== homedir = pwd (); Modified: trunk/octave-forge/main/control/devel/makefile_staircase.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_staircase.m 2011-12-15 16:00:14 UTC (rev 9390) +++ trunk/octave-forge/main/control/devel/makefile_staircase.m 2011-12-15 16:00:30 UTC (rev 9391) @@ -4,9 +4,6 @@ ## USAGE: * fetch control from Octave-Forge by svn ## * add control/inst, control/src and control/devel to your Octave path ## * run makefile_* -## NOTES: * The option "-Wl,-framework" "-Wl,vecLib" is needed for MacPorts' -## octave-devel @3.3.52_1+gcc44 on MacOS X 10.6.4. However, this option -## breaks other platforms. See MacPorts Ticket #26640. ## ============================================================================== homedir = pwd (); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-12-29 17:39:26
|
Revision: 9479 http://octave.svn.sourceforge.net/octave/?rev=9479&view=rev Author: paramaniac Date: 2011-12-29 17:39:19 +0000 (Thu, 29 Dec 2011) Log Message: ----------- control: update developer makefiles for new linking mode of mkoctfile 3.5.x (fingers crossed) Modified Paths: -------------- trunk/octave-forge/main/control/devel/makefile_chol.m trunk/octave-forge/main/control/devel/makefile_conversions.m trunk/octave-forge/main/control/devel/makefile_h2syn.m trunk/octave-forge/main/control/devel/makefile_hankel.m trunk/octave-forge/main/control/devel/makefile_hinfsyn.m trunk/octave-forge/main/control/devel/makefile_lqr.m trunk/octave-forge/main/control/devel/makefile_lyap.m trunk/octave-forge/main/control/devel/makefile_minreal.m trunk/octave-forge/main/control/devel/makefile_ncfsyn.m trunk/octave-forge/main/control/devel/makefile_norm.m trunk/octave-forge/main/control/devel/makefile_place.m trunk/octave-forge/main/control/devel/makefile_scale.m trunk/octave-forge/main/control/devel/makefile_staircase.m trunk/octave-forge/main/control/devel/makefile_tustin.m trunk/octave-forge/main/control/devel/makefile_zero.m Modified: trunk/octave-forge/main/control/devel/makefile_chol.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_chol.m 2011-12-29 13:39:07 UTC (rev 9478) +++ trunk/octave-forge/main/control/devel/makefile_chol.m 2011-12-29 17:39:19 UTC (rev 9479) @@ -14,12 +14,14 @@ mkoctfile slsb03od.cc \ SB03OD.f select.f SB03OU.f SB03OT.f MB04ND.f \ MB04OD.f SB03OR.f SB03OY.f SB04PX.f MB04NY.f \ - MB04OY.f SB03OV.f - + MB04OY.f SB03OV.f \ + "$(mkoctfile -p BLAS_LIBS)" + mkoctfile slsg03bd.cc \ SG03BD.f SG03BV.f SG03BU.f SG03BW.f SG03BX.f \ - SG03BY.f MB02UU.f MB02UV.f + SG03BY.f MB02UU.f MB02UV.f \ + "$(mkoctfile -p BLAS_LIBS)" system ("rm *.o"); cd (homedir); Modified: trunk/octave-forge/main/control/devel/makefile_conversions.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_conversions.m 2011-12-29 13:39:07 UTC (rev 9478) +++ trunk/octave-forge/main/control/devel/makefile_conversions.m 2011-12-29 17:39:19 UTC (rev 9479) @@ -15,16 +15,19 @@ mkoctfile sltb04bd.cc \ TB04BD.f MC01PY.f TB01ID.f TB01ZD.f MC01PD.f \ TB04BX.f MA02AD.f MB02RD.f MB01PD.f MB02SD.f \ - MB01QD.f + MB01QD.f \ + "$(mkoctfile -p BLAS_LIBS)" ## descriptor to regular state-space mkoctfile slsb10jd.cc \ - SB10JD.f + SB10JD.f \ + "$(mkoctfile -p BLAS_LIBS)" ## transfer function to state-space mkoctfile sltd04ad.cc \ TD04AD.f TD03AY.f TB01PD.f TB01XD.f AB07MD.f \ - TB01UD.f TB01ID.f MB01PD.f MB03OY.f MB01QD.f + TB01UD.f TB01ID.f MB01PD.f MB03OY.f MB01QD.f \ + "$(mkoctfile -p BLAS_LIBS)" system ("rm *.o"); cd (homedir); Modified: trunk/octave-forge/main/control/devel/makefile_h2syn.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_h2syn.m 2011-12-29 13:39:07 UTC (rev 9478) +++ trunk/octave-forge/main/control/devel/makefile_h2syn.m 2011-12-29 17:39:19 UTC (rev 9479) @@ -17,14 +17,16 @@ MB01SD.f SB02MS.f SB02MV.f SB02MW.f MA02AD.f \ SB02QD.f MB02PD.f SB03QX.f SB03QY.f MB01RX.f \ MB01RY.f SB03SX.f SB03SY.f select.f SB03MX.f \ - SB03MY.f MB01UD.f SB03MV.f SB03MW.f SB04PX.f + SB03MY.f MB01UD.f SB03MV.f SB03MW.f SB04PX.f \ + "$(mkoctfile -p BLAS_LIBS)" mkoctfile slsb10ed.cc \ SB10ED.f SB10SD.f SB10TD.f SB10PD.f MB01RX.f \ SB02SD.f SB02OD.f MB01RU.f SB02OU.f SB02OV.f \ SB02OW.f MB01RY.f SB02OY.f SB03SX.f SB03SY.f \ MA02ED.f select.f SB03MX.f SB02MR.f SB02MV.f \ - MB01UD.f SB03MV.f SB04PX.f + MB01UD.f SB03MV.f SB04PX.f \ + "$(mkoctfile -p BLAS_LIBS)" system ("rm *.o"); cd (homedir); Modified: trunk/octave-forge/main/control/devel/makefile_hankel.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_hankel.m 2011-12-29 13:39:07 UTC (rev 9478) +++ trunk/octave-forge/main/control/devel/makefile_hankel.m 2011-12-29 17:39:19 UTC (rev 9479) @@ -16,7 +16,8 @@ MB03UD.f TB01LD.f SB03OU.f MB03QX.f select.f \ SB03OT.f MB03QD.f MB04ND.f MB04OD.f MB03QY.f \ SB03OR.f SB03OY.f SB04PX.f MB04NY.f MB04OY.f \ - SB03OV.f + SB03OV.f \ + "$(mkoctfile -p BLAS_LIBS)" system ("rm *.o"); cd (homedir); Modified: trunk/octave-forge/main/control/devel/makefile_hinfsyn.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_hinfsyn.m 2011-12-29 13:39:07 UTC (rev 9478) +++ trunk/octave-forge/main/control/devel/makefile_hinfsyn.m 2011-12-29 17:39:19 UTC (rev 9479) @@ -17,14 +17,16 @@ SB02RU.f SB02MR.f MB01SD.f SB02MS.f SB02MV.f \ SB02MW.f SB02QD.f MB02PD.f SB03QX.f SB03QY.f \ MB01RY.f SB03SX.f SB03SY.f select.f SB03MX.f \ - SB03MY.f MB01UD.f SB03MV.f SB03MW.f SB04PX.f + SB03MY.f MB01UD.f SB03MV.f SB03MW.f SB04PX.f \ + "$(mkoctfile -p BLAS_LIBS)" mkoctfile slsb10dd.cc \ SB10DD.f MB01RU.f MB01RX.f SB02SD.f SB02OD.f \ MA02AD.f SB02OU.f SB02OV.f SB02OW.f MB01RY.f \ SB02OY.f SB03SX.f SB03SY.f MA02ED.f select.f \ SB03MX.f SB02MR.f SB02MV.f MB01UD.f SB03MV.f \ - SB04PX.f + SB04PX.f \ + "$(mkoctfile -p BLAS_LIBS)" system ("rm *.o"); cd (homedir); Modified: trunk/octave-forge/main/control/devel/makefile_lqr.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_lqr.m 2011-12-29 13:39:07 UTC (rev 9478) +++ trunk/octave-forge/main/control/devel/makefile_lqr.m 2011-12-29 17:39:19 UTC (rev 9479) @@ -13,11 +13,13 @@ mkoctfile slsb02od.cc \ SB02OD.f SB02OU.f SB02OV.f SB02OW.f SB02OY.f \ - SB02MR.f SB02MV.f + SB02MR.f SB02MV.f \ + "$(mkoctfile -p BLAS_LIBS)" mkoctfile slsg02ad.cc \ SG02AD.f SB02OU.f SB02OV.f SB02OW.f SB02OY.f \ - MB01SD.f MB02VD.f MB02PD.f MA02GD.f + MB01SD.f MB02VD.f MB02PD.f MA02GD.f \ + "$(mkoctfile -p BLAS_LIBS)" system ("rm *.o"); cd (homedir); Modified: trunk/octave-forge/main/control/devel/makefile_lyap.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_lyap.m 2011-12-29 13:39:07 UTC (rev 9478) +++ trunk/octave-forge/main/control/devel/makefile_lyap.m 2011-12-29 17:39:19 UTC (rev 9479) @@ -14,21 +14,25 @@ ## Lypunov mkoctfile slsb03md.cc \ SB03MD.f select.f SB03MX.f SB03MY.f MB01RD.f \ - SB03MV.f SB03MW.f SB04PX.f + SB03MV.f SB03MW.f SB04PX.f \ + "$(mkoctfile -p BLAS_LIBS)" ## Sylvester mkoctfile slsb04md.cc \ - SB04MD.f SB04MU.f SB04MY.f SB04MR.f SB04MW.f + SB04MD.f SB04MU.f SB04MY.f SB04MR.f SB04MW.f \ + "$(mkoctfile -p BLAS_LIBS)" mkoctfile slsb04qd.cc \ - SB04QD.f SB04QU.f SB04QY.f SB04MW.f SB04QR.f + SB04QD.f SB04QU.f SB04QY.f SB04MW.f SB04QR.f \ + "$(mkoctfile -p BLAS_LIBS)" ## Generalized Lyapunov mkoctfile slsg03ad.cc \ SG03AD.f MB01RW.f MB01RD.f SG03AX.f SG03AY.f \ - MB02UU.f MB02UV.f + MB02UU.f MB02UV.f \ + "$(mkoctfile -p BLAS_LIBS)" system ("rm *.o"); cd (homedir); Modified: trunk/octave-forge/main/control/devel/makefile_minreal.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_minreal.m 2011-12-29 13:39:07 UTC (rev 9478) +++ trunk/octave-forge/main/control/devel/makefile_minreal.m 2011-12-29 17:39:19 UTC (rev 9479) @@ -13,10 +13,12 @@ mkoctfile sltb01pd.cc \ TB01PD.f TB01XD.f TB01ID.f AB07MD.f TB01UD.f \ - MB03OY.f MB01PD.f MB01QD.f + MB03OY.f MB01PD.f MB01QD.f \ + "$(mkoctfile -p BLAS_LIBS)" mkoctfile sltg01jd.cc \ - TG01JD.f TG01AD.f TB01XD.f MA02CD.f TG01HX.f + TG01JD.f TG01AD.f TB01XD.f MA02CD.f TG01HX.f \ + "$(mkoctfile -p BLAS_LIBS)" system ("rm *.o"); cd (homedir); Modified: trunk/octave-forge/main/control/devel/makefile_ncfsyn.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_ncfsyn.m 2011-12-29 13:39:07 UTC (rev 9478) +++ trunk/octave-forge/main/control/devel/makefile_ncfsyn.m 2011-12-29 17:39:19 UTC (rev 9479) @@ -18,18 +18,22 @@ MB01RU.f SB02QD.f SB02MV.f SB02MW.f SB02MR.f \ MA02AD.f MB02PD.f MB01SD.f MB01UD.f SB03SY.f \ MB01RX.f SB03MX.f SB03SX.f MB01RY.f SB03QY.f \ - SB03QX.f SB03MY.f SB04PX.f SB03MV.f SB03MW.f + SB03QX.f SB03MY.f SB04PX.f SB03MV.f SB03MW.f \ + "$(mkoctfile -p BLAS_LIBS)" ## H-infinity loop shaping - discrete-time - strictly proper case mkoctfile slsb10kd.cc \ SB10KD.f SB02OD.f select.f SB02OY.f SB02OW.f \ - SB02OV.f SB02MV.f SB02OU.f SB02MR.f + SB02OV.f SB02MV.f SB02OU.f SB02MR.f \ + "$(mkoctfile -p BLAS_LIBS)" + ## H-infinity loop shaping - discrete-time - proper case mkoctfile slsb10zd.cc \ SB10ZD.f MA02AD.f SB02OD.f select.f MB01RX.f \ MB02VD.f SB02OY.f SB02OW.f SB02OV.f SB02OU.f \ - SB02MR.f MA02GD.f SB02MV.f + SB02MR.f MA02GD.f SB02MV.f \ + "$(mkoctfile -p BLAS_LIBS)" system ("rm *.o"); cd (homedir); Modified: trunk/octave-forge/main/control/devel/makefile_norm.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_norm.m 2011-12-29 13:39:07 UTC (rev 9478) +++ trunk/octave-forge/main/control/devel/makefile_norm.m 2011-12-29 17:39:19 UTC (rev 9479) @@ -16,7 +16,8 @@ AB13BD.f SB08DD.f SB03OU.f SB01FY.f TB01LD.f \ SB03OT.f MB04ND.f MB04OD.f MB03QX.f select.f \ SB03OR.f MB04OX.f MB03QD.f SB03OY.f MA02AD.f \ - MB03QY.f SB04PX.f MB04NY.f MB04OY.f SB03OV.f + MB03QY.f SB04PX.f MB04NY.f MB04OY.f SB03OV.f \ + "$(mkoctfile -p BLAS_LIBS)" ## L-inf norm mkoctfile slab13dd.cc \ @@ -25,7 +26,9 @@ MB03XP.f MB04DD.f MB04QB.f MB04TB.f MB03XU.f \ MB04TS.f UE01MD.f MB02RD.f MB02SD.f MB04QC.f \ MB04QF.f MB03YA.f MB03YD.f MB02RZ.f MB04QU.f \ - MB02SZ.f MB03YT.f + MB02SZ.f MB03YT.f \ + "$(mkoctfile -p BLAS_LIBS)" \ + "$(mkoctfile -p FLIBS)" system ("rm *.o"); cd (homedir); Modified: trunk/octave-forge/main/control/devel/makefile_place.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_place.m 2011-12-29 13:39:07 UTC (rev 9478) +++ trunk/octave-forge/main/control/devel/makefile_place.m 2011-12-29 17:39:19 UTC (rev 9479) @@ -13,7 +13,8 @@ mkoctfile slsb01bd.cc \ SB01BD.f MB03QD.f MB03QY.f SB01BX.f SB01BY.f \ - select.f + select.f \ + "$(mkoctfile -p BLAS_LIBS)" system ("rm *.o"); cd (homedir); Modified: trunk/octave-forge/main/control/devel/makefile_scale.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_scale.m 2011-12-29 13:39:07 UTC (rev 9478) +++ trunk/octave-forge/main/control/devel/makefile_scale.m 2011-12-29 17:39:19 UTC (rev 9479) @@ -13,11 +13,13 @@ ## scaling of state-space models mkoctfile sltb01id.cc \ - TB01ID.f + TB01ID.f \ + "$(mkoctfile -p BLAS_LIBS)" ## scaling of descriptor state-space models mkoctfile sltg01ad.cc \ - TG01AD.f + TG01AD.f \ + "$(mkoctfile -p BLAS_LIBS)" system ("rm *.o"); cd (homedir); Modified: trunk/octave-forge/main/control/devel/makefile_staircase.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_staircase.m 2011-12-29 13:39:07 UTC (rev 9478) +++ trunk/octave-forge/main/control/devel/makefile_staircase.m 2011-12-29 17:39:19 UTC (rev 9479) @@ -13,16 +13,19 @@ ## staircase form using orthogonal transformations mkoctfile slab01od.cc \ - AB01OD.f AB01ND.f MB03OY.f MB01PD.f MB01QD.f + AB01OD.f AB01ND.f MB03OY.f MB01PD.f MB01QD.f \ + "$(mkoctfile -p BLAS_LIBS)" ## controllability staircase form of descriptor state-space models mkoctfile sltg01hd.cc \ - TG01HD.f TG01HX.f + TG01HD.f TG01HX.f \ + "$(mkoctfile -p BLAS_LIBS)" ## observability staircase form of descriptor state-space models mkoctfile sltg01id.cc \ TG01ID.f TB01XD.f MA02CD.f AB07MD.f TG01HX.f \ - MA02BD.f + MA02BD.f \ + "$(mkoctfile -p BLAS_LIBS)" system ("rm *.o"); cd (homedir); \ No newline at end of file Modified: trunk/octave-forge/main/control/devel/makefile_tustin.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_tustin.m 2011-12-29 13:39:07 UTC (rev 9478) +++ trunk/octave-forge/main/control/devel/makefile_tustin.m 2011-12-29 17:39:19 UTC (rev 9479) @@ -13,7 +13,8 @@ ## bilinear transformation mkoctfile slab04md.cc \ - AB04MD.f + AB04MD.f \ + "$(mkoctfile -p BLAS_LIBS)" system ("rm *.o"); cd (homedir); Modified: trunk/octave-forge/main/control/devel/makefile_zero.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_zero.m 2011-12-29 13:39:07 UTC (rev 9478) +++ trunk/octave-forge/main/control/devel/makefile_zero.m 2011-12-29 17:39:19 UTC (rev 9479) @@ -13,16 +13,19 @@ ## transmission zeros of state-space models mkoctfile slab08nd.cc \ - AB08ND.f AB08NX.f TB01ID.f MB03OY.f MB03PY.f + AB08ND.f AB08NX.f TB01ID.f MB03OY.f MB03PY.f \ + "$(mkoctfile -p BLAS_LIBS)" ## transmission zeros of descriptor state-space models mkoctfile slag08bd.cc \ AG08BD.f AG08BY.f TG01AD.f TB01XD.f MA02CD.f \ - TG01FD.f MA02BD.f MB03OY.f + TG01FD.f MA02BD.f MB03OY.f \ + "$(mkoctfile -p BLAS_LIBS)" ## gain of descriptor state-space models mkoctfile sltg04bx.cc \ - TG04BX.f MB02RD.f MB02SD.f + TG04BX.f MB02RD.f MB02SD.f \ + "$(mkoctfile -p BLAS_LIBS)" system ("rm *.o"); cd (homedir); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-12-29 21:16:04
|
Revision: 9481 http://octave.svn.sourceforge.net/octave/?rev=9481&view=rev Author: paramaniac Date: 2011-12-29 21:15:57 +0000 (Thu, 29 Dec 2011) Log Message: ----------- control: update developer makefiles for new linking mode of mkoctfile 3.5.x (for real this time?) Modified Paths: -------------- trunk/octave-forge/main/control/devel/makefile_chol.m trunk/octave-forge/main/control/devel/makefile_conversions.m trunk/octave-forge/main/control/devel/makefile_h2syn.m trunk/octave-forge/main/control/devel/makefile_hankel.m trunk/octave-forge/main/control/devel/makefile_hinfsyn.m trunk/octave-forge/main/control/devel/makefile_lqr.m trunk/octave-forge/main/control/devel/makefile_lyap.m trunk/octave-forge/main/control/devel/makefile_minreal.m trunk/octave-forge/main/control/devel/makefile_ncfsyn.m trunk/octave-forge/main/control/devel/makefile_norm.m trunk/octave-forge/main/control/devel/makefile_place.m trunk/octave-forge/main/control/devel/makefile_scale.m trunk/octave-forge/main/control/devel/makefile_staircase.m trunk/octave-forge/main/control/devel/makefile_tustin.m trunk/octave-forge/main/control/devel/makefile_zero.m Modified: trunk/octave-forge/main/control/devel/makefile_chol.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_chol.m 2011-12-29 17:43:24 UTC (rev 9480) +++ trunk/octave-forge/main/control/devel/makefile_chol.m 2011-12-29 21:15:57 UTC (rev 9481) @@ -15,12 +15,13 @@ SB03OD.f select.f SB03OU.f SB03OT.f MB04ND.f \ MB04OD.f SB03OR.f SB03OY.f SB04PX.f MB04NY.f \ MB04OY.f SB03OV.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ "$(mkoctfile -p BLAS_LIBS)" - mkoctfile slsg03bd.cc \ SG03BD.f SG03BV.f SG03BU.f SG03BW.f SG03BX.f \ SG03BY.f MB02UU.f MB02UV.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ "$(mkoctfile -p BLAS_LIBS)" system ("rm *.o"); Modified: trunk/octave-forge/main/control/devel/makefile_conversions.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_conversions.m 2011-12-29 17:43:24 UTC (rev 9480) +++ trunk/octave-forge/main/control/devel/makefile_conversions.m 2011-12-29 21:15:57 UTC (rev 9481) @@ -16,17 +16,20 @@ TB04BD.f MC01PY.f TB01ID.f TB01ZD.f MC01PD.f \ TB04BX.f MA02AD.f MB02RD.f MB01PD.f MB02SD.f \ MB01QD.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ "$(mkoctfile -p BLAS_LIBS)" ## descriptor to regular state-space mkoctfile slsb10jd.cc \ SB10JD.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ "$(mkoctfile -p BLAS_LIBS)" ## transfer function to state-space mkoctfile sltd04ad.cc \ TD04AD.f TD03AY.f TB01PD.f TB01XD.f AB07MD.f \ TB01UD.f TB01ID.f MB01PD.f MB03OY.f MB01QD.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ "$(mkoctfile -p BLAS_LIBS)" system ("rm *.o"); Modified: trunk/octave-forge/main/control/devel/makefile_h2syn.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_h2syn.m 2011-12-29 17:43:24 UTC (rev 9480) +++ trunk/octave-forge/main/control/devel/makefile_h2syn.m 2011-12-29 21:15:57 UTC (rev 9481) @@ -18,6 +18,7 @@ SB02QD.f MB02PD.f SB03QX.f SB03QY.f MB01RX.f \ MB01RY.f SB03SX.f SB03SY.f select.f SB03MX.f \ SB03MY.f MB01UD.f SB03MV.f SB03MW.f SB04PX.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ "$(mkoctfile -p BLAS_LIBS)" mkoctfile slsb10ed.cc \ @@ -26,6 +27,7 @@ SB02OW.f MB01RY.f SB02OY.f SB03SX.f SB03SY.f \ MA02ED.f select.f SB03MX.f SB02MR.f SB02MV.f \ MB01UD.f SB03MV.f SB04PX.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ "$(mkoctfile -p BLAS_LIBS)" system ("rm *.o"); Modified: trunk/octave-forge/main/control/devel/makefile_hankel.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_hankel.m 2011-12-29 17:43:24 UTC (rev 9480) +++ trunk/octave-forge/main/control/devel/makefile_hankel.m 2011-12-29 21:15:57 UTC (rev 9481) @@ -17,6 +17,7 @@ SB03OT.f MB03QD.f MB04ND.f MB04OD.f MB03QY.f \ SB03OR.f SB03OY.f SB04PX.f MB04NY.f MB04OY.f \ SB03OV.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ "$(mkoctfile -p BLAS_LIBS)" system ("rm *.o"); Modified: trunk/octave-forge/main/control/devel/makefile_hinfsyn.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_hinfsyn.m 2011-12-29 17:43:24 UTC (rev 9480) +++ trunk/octave-forge/main/control/devel/makefile_hinfsyn.m 2011-12-29 21:15:57 UTC (rev 9481) @@ -18,6 +18,7 @@ SB02MW.f SB02QD.f MB02PD.f SB03QX.f SB03QY.f \ MB01RY.f SB03SX.f SB03SY.f select.f SB03MX.f \ SB03MY.f MB01UD.f SB03MV.f SB03MW.f SB04PX.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ "$(mkoctfile -p BLAS_LIBS)" mkoctfile slsb10dd.cc \ @@ -26,6 +27,7 @@ SB02OY.f SB03SX.f SB03SY.f MA02ED.f select.f \ SB03MX.f SB02MR.f SB02MV.f MB01UD.f SB03MV.f \ SB04PX.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ "$(mkoctfile -p BLAS_LIBS)" system ("rm *.o"); Modified: trunk/octave-forge/main/control/devel/makefile_lqr.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_lqr.m 2011-12-29 17:43:24 UTC (rev 9480) +++ trunk/octave-forge/main/control/devel/makefile_lqr.m 2011-12-29 21:15:57 UTC (rev 9481) @@ -14,11 +14,13 @@ mkoctfile slsb02od.cc \ SB02OD.f SB02OU.f SB02OV.f SB02OW.f SB02OY.f \ SB02MR.f SB02MV.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ "$(mkoctfile -p BLAS_LIBS)" mkoctfile slsg02ad.cc \ SG02AD.f SB02OU.f SB02OV.f SB02OW.f SB02OY.f \ MB01SD.f MB02VD.f MB02PD.f MA02GD.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ "$(mkoctfile -p BLAS_LIBS)" system ("rm *.o"); Modified: trunk/octave-forge/main/control/devel/makefile_lyap.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_lyap.m 2011-12-29 17:43:24 UTC (rev 9480) +++ trunk/octave-forge/main/control/devel/makefile_lyap.m 2011-12-29 21:15:57 UTC (rev 9481) @@ -15,16 +15,19 @@ mkoctfile slsb03md.cc \ SB03MD.f select.f SB03MX.f SB03MY.f MB01RD.f \ SB03MV.f SB03MW.f SB04PX.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ "$(mkoctfile -p BLAS_LIBS)" ## Sylvester mkoctfile slsb04md.cc \ SB04MD.f SB04MU.f SB04MY.f SB04MR.f SB04MW.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ "$(mkoctfile -p BLAS_LIBS)" mkoctfile slsb04qd.cc \ SB04QD.f SB04QU.f SB04QY.f SB04MW.f SB04QR.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ "$(mkoctfile -p BLAS_LIBS)" @@ -32,6 +35,7 @@ mkoctfile slsg03ad.cc \ SG03AD.f MB01RW.f MB01RD.f SG03AX.f SG03AY.f \ MB02UU.f MB02UV.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ "$(mkoctfile -p BLAS_LIBS)" system ("rm *.o"); Modified: trunk/octave-forge/main/control/devel/makefile_minreal.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_minreal.m 2011-12-29 17:43:24 UTC (rev 9480) +++ trunk/octave-forge/main/control/devel/makefile_minreal.m 2011-12-29 21:15:57 UTC (rev 9481) @@ -14,10 +14,12 @@ mkoctfile sltb01pd.cc \ TB01PD.f TB01XD.f TB01ID.f AB07MD.f TB01UD.f \ MB03OY.f MB01PD.f MB01QD.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ "$(mkoctfile -p BLAS_LIBS)" mkoctfile sltg01jd.cc \ TG01JD.f TG01AD.f TB01XD.f MA02CD.f TG01HX.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ "$(mkoctfile -p BLAS_LIBS)" system ("rm *.o"); Modified: trunk/octave-forge/main/control/devel/makefile_ncfsyn.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_ncfsyn.m 2011-12-29 17:43:24 UTC (rev 9480) +++ trunk/octave-forge/main/control/devel/makefile_ncfsyn.m 2011-12-29 21:15:57 UTC (rev 9481) @@ -19,20 +19,22 @@ MA02AD.f MB02PD.f MB01SD.f MB01UD.f SB03SY.f \ MB01RX.f SB03MX.f SB03SX.f MB01RY.f SB03QY.f \ SB03QX.f SB03MY.f SB04PX.f SB03MV.f SB03MW.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ "$(mkoctfile -p BLAS_LIBS)" ## H-infinity loop shaping - discrete-time - strictly proper case mkoctfile slsb10kd.cc \ SB10KD.f SB02OD.f select.f SB02OY.f SB02OW.f \ SB02OV.f SB02MV.f SB02OU.f SB02MR.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ "$(mkoctfile -p BLAS_LIBS)" - ## H-infinity loop shaping - discrete-time - proper case mkoctfile slsb10zd.cc \ SB10ZD.f MA02AD.f SB02OD.f select.f MB01RX.f \ MB02VD.f SB02OY.f SB02OW.f SB02OV.f SB02OU.f \ SB02MR.f MA02GD.f SB02MV.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ "$(mkoctfile -p BLAS_LIBS)" system ("rm *.o"); Modified: trunk/octave-forge/main/control/devel/makefile_norm.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_norm.m 2011-12-29 17:43:24 UTC (rev 9480) +++ trunk/octave-forge/main/control/devel/makefile_norm.m 2011-12-29 21:15:57 UTC (rev 9481) @@ -17,6 +17,7 @@ SB03OT.f MB04ND.f MB04OD.f MB03QX.f select.f \ SB03OR.f MB04OX.f MB03QD.f SB03OY.f MA02AD.f \ MB03QY.f SB04PX.f MB04NY.f MB04OY.f SB03OV.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ "$(mkoctfile -p BLAS_LIBS)" ## L-inf norm @@ -27,6 +28,7 @@ MB04TS.f UE01MD.f MB02RD.f MB02SD.f MB04QC.f \ MB04QF.f MB03YA.f MB03YD.f MB02RZ.f MB04QU.f \ MB02SZ.f MB03YT.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ "$(mkoctfile -p BLAS_LIBS)" \ "$(mkoctfile -p FLIBS)" Modified: trunk/octave-forge/main/control/devel/makefile_place.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_place.m 2011-12-29 17:43:24 UTC (rev 9480) +++ trunk/octave-forge/main/control/devel/makefile_place.m 2011-12-29 21:15:57 UTC (rev 9481) @@ -14,6 +14,7 @@ mkoctfile slsb01bd.cc \ SB01BD.f MB03QD.f MB03QY.f SB01BX.f SB01BY.f \ select.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ "$(mkoctfile -p BLAS_LIBS)" system ("rm *.o"); Modified: trunk/octave-forge/main/control/devel/makefile_scale.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_scale.m 2011-12-29 17:43:24 UTC (rev 9480) +++ trunk/octave-forge/main/control/devel/makefile_scale.m 2011-12-29 21:15:57 UTC (rev 9481) @@ -14,11 +14,13 @@ ## scaling of state-space models mkoctfile sltb01id.cc \ TB01ID.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ "$(mkoctfile -p BLAS_LIBS)" ## scaling of descriptor state-space models mkoctfile sltg01ad.cc \ TG01AD.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ "$(mkoctfile -p BLAS_LIBS)" system ("rm *.o"); Modified: trunk/octave-forge/main/control/devel/makefile_staircase.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_staircase.m 2011-12-29 17:43:24 UTC (rev 9480) +++ trunk/octave-forge/main/control/devel/makefile_staircase.m 2011-12-29 21:15:57 UTC (rev 9481) @@ -14,17 +14,20 @@ ## staircase form using orthogonal transformations mkoctfile slab01od.cc \ AB01OD.f AB01ND.f MB03OY.f MB01PD.f MB01QD.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ "$(mkoctfile -p BLAS_LIBS)" ## controllability staircase form of descriptor state-space models mkoctfile sltg01hd.cc \ TG01HD.f TG01HX.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ "$(mkoctfile -p BLAS_LIBS)" ## observability staircase form of descriptor state-space models mkoctfile sltg01id.cc \ TG01ID.f TB01XD.f MA02CD.f AB07MD.f TG01HX.f \ MA02BD.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ "$(mkoctfile -p BLAS_LIBS)" system ("rm *.o"); Modified: trunk/octave-forge/main/control/devel/makefile_tustin.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_tustin.m 2011-12-29 17:43:24 UTC (rev 9480) +++ trunk/octave-forge/main/control/devel/makefile_tustin.m 2011-12-29 21:15:57 UTC (rev 9481) @@ -14,6 +14,7 @@ ## bilinear transformation mkoctfile slab04md.cc \ AB04MD.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ "$(mkoctfile -p BLAS_LIBS)" system ("rm *.o"); Modified: trunk/octave-forge/main/control/devel/makefile_zero.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_zero.m 2011-12-29 17:43:24 UTC (rev 9480) +++ trunk/octave-forge/main/control/devel/makefile_zero.m 2011-12-29 21:15:57 UTC (rev 9481) @@ -14,17 +14,20 @@ ## transmission zeros of state-space models mkoctfile slab08nd.cc \ AB08ND.f AB08NX.f TB01ID.f MB03OY.f MB03PY.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ "$(mkoctfile -p BLAS_LIBS)" ## transmission zeros of descriptor state-space models mkoctfile slag08bd.cc \ AG08BD.f AG08BY.f TG01AD.f TB01XD.f MA02CD.f \ TG01FD.f MA02BD.f MB03OY.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ "$(mkoctfile -p BLAS_LIBS)" ## gain of descriptor state-space models mkoctfile sltg04bx.cc \ TG04BX.f MB02RD.f MB02SD.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ "$(mkoctfile -p BLAS_LIBS)" system ("rm *.o"); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2011-12-30 13:42:35
|
Revision: 9487 http://octave.svn.sourceforge.net/octave/?rev=9487&view=rev Author: paramaniac Date: 2011-12-30 13:42:29 +0000 (Fri, 30 Dec 2011) Log Message: ----------- control: finish ctrbf and obsvf Modified Paths: -------------- trunk/octave-forge/main/control/devel/ctrbf/makefile_ctrbf.m trunk/octave-forge/main/control/devel/ctrbf.m trunk/octave-forge/main/control/devel/makefile_staircase.m trunk/octave-forge/main/control/devel/obsvf.m Modified: trunk/octave-forge/main/control/devel/ctrbf/makefile_ctrbf.m =================================================================== --- trunk/octave-forge/main/control/devel/ctrbf/makefile_ctrbf.m 2011-12-29 23:37:26 UTC (rev 9486) +++ trunk/octave-forge/main/control/devel/ctrbf/makefile_ctrbf.m 2011-12-30 13:42:29 UTC (rev 9487) @@ -1,2 +1,4 @@ mkoctfile sltb01ud.cc \ - TB01UD.f MB01PD.f MB03OY.f MB01QD.f \ No newline at end of file + TB01UD.f MB01PD.f MB03OY.f MB01QD.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ + "$(mkoctfile -p BLAS_LIBS)" Modified: trunk/octave-forge/main/control/devel/ctrbf.m =================================================================== --- trunk/octave-forge/main/control/devel/ctrbf.m 2011-12-29 23:37:26 UTC (rev 9486) +++ trunk/octave-forge/main/control/devel/ctrbf.m 2011-12-30 13:42:29 UTC (rev 9487) @@ -17,7 +17,9 @@ ## along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn{Function File} {[@var{Abar}, @var{Bbar}, @var{Cbar}, @var{T}, @var{K}] =} ctrbf (@var{A}, @var{B}, @var{C}) +## @deftypefn{Function File} {[@var{sysbar}, @var{T}, @var{K}] =} ctrbf (@var{sys}) +## @deftypefnx{Function File} {[@var{sysbar}, @var{T}, @var{K}] =} ctrbf (@var{sys}, @var{tol}) +## @deftypefnx{Function File} {[@var{Abar}, @var{Bbar}, @var{Cbar}, @var{T}, @var{K}] =} ctrbf (@var{A}, @var{B}, @var{C}) ## @deftypefnx{Function File} {[@var{Abar}, @var{Bbar}, @var{Cbar}, @var{T}, @var{K}] =} ctrbf (@var{A}, @var{B}, @var{C}, @var{TOL}) ## If Co=ctrb(A,B) has rank r <= n = SIZE(A,1), then there is a ## similarity transformation Tc such that Tc = [t1 t2] where t1 @@ -49,11 +51,28 @@ ## Created: 2010-04-30 ## Version: 0.1 -function [ac, bc, cc, z, ncont] = ctrbf (a, b, c, tol = []) +function [ac, bc, cc, z, ncont] = ctrbf (a, b = [], c, tol = []) - if (nargin < 3 || nargin > 4) + if (nargin < 1 || nargin > 4) print_usage (); endif + + islti = isa (a, "lti"); + + if (islti) + if (nargin > 2) + print_usage (); + endif + sys = a; + tol = b; + [a, b, c] = ssdata (sys); + else + if (nargin < 3) + print_usage (); + endif + sys = ss (a, b, c); + [a, b, c] = ssdata (sys); + endif if (isempty (tol)) tol = 0; # default tolerance @@ -63,6 +82,12 @@ [ac, bc, cc, z, ncont] = sltb01ud (a, b, c, tol); + if (islti) + ac = set (sys, "a", ac, "b", bc, "c", cc, "scaled", false); + bc = z; + cc = ncont; + endif + endfunction %!shared Ao, Bo, Co, Zo, Ae, Be, Ce, Ze, NCONT Modified: trunk/octave-forge/main/control/devel/makefile_staircase.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_staircase.m 2011-12-29 23:37:26 UTC (rev 9486) +++ trunk/octave-forge/main/control/devel/makefile_staircase.m 2011-12-30 13:42:29 UTC (rev 9487) @@ -30,5 +30,11 @@ "$(mkoctfile -p LAPACK_LIBS)" \ "$(mkoctfile -p BLAS_LIBS)" +## controllable block Hessenberg realization +mkoctfile sltb01ud.cc \ + TB01UD.f MB01PD.f MB03OY.f MB01QD.f \ + "$(mkoctfile -p LAPACK_LIBS)" \ + "$(mkoctfile -p BLAS_LIBS)" + system ("rm *.o"); cd (homedir); \ No newline at end of file Modified: trunk/octave-forge/main/control/devel/obsvf.m =================================================================== --- trunk/octave-forge/main/control/devel/obsvf.m 2011-12-29 23:37:26 UTC (rev 9486) +++ trunk/octave-forge/main/control/devel/obsvf.m 2011-12-30 13:42:29 UTC (rev 9487) @@ -1,21 +1,25 @@ -## Copyright (C) 2010 Benjamin Fernandez +## Copyright (C) 2010 Benjamin Fernandez +## Copyright (C) 2011 Lukas F. Reichlin ## -## This program is free software; you can redistribute it and/or modify +## 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 2 of the License, or +## the Free Software Foundation, either version 3 of the License, or ## (at your option) any later version. -## -## This program is distributed in the hope that it will be useful, +## +## 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 Octave; see the file COPYING. If not, see -## <http://www.gnu.org/licenses/>. +## along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn{Function File} {[@var{Abar}, @var{Bbar}, @var{Cbar}, @var{T}, @var{K}] =} obsvf (@var{A}, @var{B}, @var{C}) +## @deftypefn{Function File} {[@var{sysbar}, @var{T}, @var{K}] =} obsvf (@var{sys}) +## @deftypefnx{Function File} {[@var{sysbar}, @var{T}, @var{K}] =} obsvf (@var{sys}, @var{tol}) +## @deftypefnx{Function File} {[@var{Abar}, @var{Bbar}, @var{Cbar}, @var{T}, @var{K}] =} obsvf (@var{A}, @var{B}, @var{C}) ## @deftypefnx{Function File} {[@var{Abar}, @var{Bbar}, @var{Cbar}, @var{T}, @var{K}] =} obsvf (@var{A}, @var{B}, @var{C}, @var{TOL}) ## If Ob=obsv(A,C) has rank r <= n = SIZE(A,1), then there is a ## similarity transformation Tc such that To = [t1;t2] where t1 is c @@ -45,26 +49,29 @@ ## Author: Benjamin Fernandez <benjas@benjas-laptop> ## Created: 2010-05-02 +## Version: 0.1 -function [Abar, Bbar, Cbar, T, K] = obsvf (A, B, C, TOL) +function [ac, bc, cc, z, ncont] = obsvf (a, b = [], c, tol = []) - if (nargin < 3 || nargin > 4) + if (nargin < 1 || nargin > 4) print_usage (); endif - - if (nargin == 3) - TOL = length (A) * norm (A, 1) * eps; + + if (isa (a, "lti")) + if (nargin > 2) + print_usage (); + endif + [ac, bc, cc] = ctrbf (a.', b); # [sysbar, z, ncont] = ctrbf (sys.', tol); + ac = ac.'; + z = ncont = []; + else + if (nargin < 3) + print_usage (); + endif + [ac, tmp, cc, z, ncont] = ctrbf (a.', c.', b.', tol); + ac = ac.'; + bc = cc.'; + cc = tmp.'; endif - Ob = obsv (A, C); - [nro, nco] = size (Ob); - rob = rank (Ob); - lr = nco - rob; - [U, S, V] = svd (Ob); - K = V(:, 1:rob); # Basis raw space - T = V; # [c; orth(c)]; - Abar = T \ A * T; - Bbar = T \ B; - Cbar = C * T; - -endfunction \ No newline at end of file +endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-02-22 18:50:27
|
Revision: 9657 http://octave.svn.sourceforge.net/octave/?rev=9657&view=rev Author: paramaniac Date: 2012-02-22 18:50:21 +0000 (Wed, 22 Feb 2012) Log Message: ----------- control: reorganize package makefile (4) Modified Paths: -------------- trunk/octave-forge/main/control/devel/makefile_chol.m trunk/octave-forge/main/control/devel/makefile_control.m Modified: trunk/octave-forge/main/control/devel/makefile_chol.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_chol.m 2012-02-22 18:34:38 UTC (rev 9656) +++ trunk/octave-forge/main/control/devel/makefile_chol.m 2012-02-22 18:50:21 UTC (rev 9657) @@ -11,18 +11,8 @@ srcdir = [develdir, "/../src"]; cd (srcdir); -mkoctfile slsb03od.cc \ - SB03OD.f select.f SB03OU.f SB03OT.f MB04ND.f \ - MB04OD.f SB03OR.f SB03OY.f SB04PX.f MB04NY.f \ - MB04OY.f SB03OV.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" +system ("rm -rf *.o slsb03od.oct slsg03bd.oct"); +system ("make slsb03od.oct slsg03bd.oct") -mkoctfile slsg03bd.cc \ - SG03BD.f SG03BV.f SG03BU.f SG03BW.f SG03BX.f \ - SG03BY.f MB02UU.f MB02UV.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -system ("rm *.o"); +system ("rm -rf *.o"); cd (homedir); Modified: trunk/octave-forge/main/control/devel/makefile_control.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_control.m 2012-02-22 18:34:38 UTC (rev 9656) +++ trunk/octave-forge/main/control/devel/makefile_control.m 2012-02-22 18:50:21 UTC (rev 9657) @@ -30,7 +30,7 @@ makefile_zero %} -system ("make clean"); +system ("make realclean"); system ("make -j4 all"); system ("rm *.o"); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-02-22 19:21:15
|
Revision: 9660 http://octave.svn.sourceforge.net/octave/?rev=9660&view=rev Author: paramaniac Date: 2012-02-22 19:21:08 +0000 (Wed, 22 Feb 2012) Log Message: ----------- control: touch up makefiles Modified Paths: -------------- trunk/octave-forge/main/control/devel/makefile_chol.m trunk/octave-forge/main/control/devel/makefile_conversions.m Modified: trunk/octave-forge/main/control/devel/makefile_chol.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_chol.m 2012-02-22 19:13:32 UTC (rev 9659) +++ trunk/octave-forge/main/control/devel/makefile_chol.m 2012-02-22 19:21:08 UTC (rev 9660) @@ -12,7 +12,14 @@ cd (srcdir); system ("rm -rf *.o slsb03od.oct slsg03bd.oct"); -system ("make slsb03od.oct slsg03bd.oct") +system ("make slsb03od.oct slsg03bd.oct"); system ("rm -rf *.o"); cd (homedir); + +%{ +system ("rm -rf *.o "); +system ("make "); + +system ("rm -rf *.o"); +%} \ No newline at end of file Modified: trunk/octave-forge/main/control/devel/makefile_conversions.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_conversions.m 2012-02-22 19:13:32 UTC (rev 9659) +++ trunk/octave-forge/main/control/devel/makefile_conversions.m 2012-02-22 19:21:08 UTC (rev 9660) @@ -12,25 +12,13 @@ cd (srcdir); ## state-space to transfer function -mkoctfile sltb04bd.cc \ - TB04BD.f MC01PY.f TB01ID.f TB01ZD.f MC01PD.f \ - TB04BX.f MA02AD.f MB02RD.f MB01PD.f MB02SD.f \ - MB01QD.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" ## descriptor to regular state-space -mkoctfile slsb10jd.cc \ - SB10JD.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" ## transfer function to state-space -mkoctfile sltd04ad.cc \ - TD04AD.f TD03AY.f TB01PD.f TB01XD.f AB07MD.f \ - TB01UD.f TB01ID.f MB01PD.f MB03OY.f MB01QD.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" -system ("rm *.o"); +system ("rm -rf *.o sltb04bd.oct slsb10jd.oct sltd04ad.oct"); +system ("make sltb04bd.oct slsb10jd.oct sltd04ad.oct"); + +system ("rm -rf *.o"); cd (homedir); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-02-26 12:42:32
|
Revision: 9698 http://octave.svn.sourceforge.net/octave/?rev=9698&view=rev Author: paramaniac Date: 2012-02-26 12:42:25 +0000 (Sun, 26 Feb 2012) Log Message: ----------- control: remove obsolete developer makefiles Modified Paths: -------------- trunk/octave-forge/main/control/devel/makefile_control.m Removed Paths: ------------- trunk/octave-forge/main/control/devel/makefile_chol.m trunk/octave-forge/main/control/devel/makefile_conred.m trunk/octave-forge/main/control/devel/makefile_conversions.m trunk/octave-forge/main/control/devel/makefile_h2syn.m trunk/octave-forge/main/control/devel/makefile_hankel.m trunk/octave-forge/main/control/devel/makefile_hinfsyn.m trunk/octave-forge/main/control/devel/makefile_lqr.m trunk/octave-forge/main/control/devel/makefile_lyap.m trunk/octave-forge/main/control/devel/makefile_minreal.m trunk/octave-forge/main/control/devel/makefile_modred.m trunk/octave-forge/main/control/devel/makefile_ncfsyn.m trunk/octave-forge/main/control/devel/makefile_norm.m trunk/octave-forge/main/control/devel/makefile_place.m trunk/octave-forge/main/control/devel/makefile_scale.m trunk/octave-forge/main/control/devel/makefile_staircase.m trunk/octave-forge/main/control/devel/makefile_tustin.m trunk/octave-forge/main/control/devel/makefile_zero.m Deleted: trunk/octave-forge/main/control/devel/makefile_chol.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_chol.m 2012-02-26 12:08:31 UTC (rev 9697) +++ trunk/octave-forge/main/control/devel/makefile_chol.m 2012-02-26 12:42:25 UTC (rev 9698) @@ -1,25 +0,0 @@ -## ============================================================================== -## Developer Makefile for OCT-files -## ============================================================================== -## USAGE: * fetch control from Octave-Forge by svn -## * add control/inst, control/src and control/devel to your Octave path -## * run makefile_* -## ============================================================================== - -homedir = pwd (); -develdir = fileparts (which ("makefile_chol")); -srcdir = [develdir, "/../src"]; -cd (srcdir); - -system ("rm -rf *.o slsb03od.oct slsg03bd.oct"); -system ("make slsb03od.oct slsg03bd.oct"); - -system ("rm -rf *.o"); -cd (homedir); - -%{ -system ("rm -rf *.o "); -system ("make "); - -system ("rm -rf *.o"); -%} \ No newline at end of file Deleted: trunk/octave-forge/main/control/devel/makefile_conred.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_conred.m 2012-02-26 12:08:31 UTC (rev 9697) +++ trunk/octave-forge/main/control/devel/makefile_conred.m 2012-02-26 12:42:25 UTC (rev 9698) @@ -1,44 +0,0 @@ -## ============================================================================== -## Developer Makefile for OCT-files -## ============================================================================== -## USAGE: * fetch control-devel from Octave-Forge by svn -## * add control-devel/inst, control-devel/src and control-devel/devel -## to your Octave path (by an .octaverc file) -## * run makefile_conred -## ============================================================================== - -homedir = pwd (); -develdir = fileparts (which ("makefile_conred")); -srcdir = [develdir, "/../src"]; -cd (srcdir); - -mkoctfile slsb16ad.cc \ - SB16AD.f TB01ID.f SB16AY.f TB01KD.f AB09IX.f \ - MB04OD.f MB01WD.f SB03OD.f MB03UD.f AB05PD.f \ - AB09DD.f AB07ND.f TB01LD.f AB05QD.f SB03OU.f \ - MA02AD.f MB03QX.f select.f MB01YD.f MB01ZD.f \ - SB03OT.f MB04OY.f MB03QD.f MB04ND.f MB03QY.f \ - SB03OR.f SB03OY.f SB04PX.f MB04NY.f SB03OV.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -mkoctfile slsb16bd.cc \ - SB16BD.f AB09AD.f AB09BD.f SB08GD.f SB08HD.f \ - TB01ID.f AB09AX.f MA02GD.f AB09BX.f TB01WD.f \ - MA02DD.f MB03UD.f select.f AB09DD.f SB03OU.f \ - MA02AD.f SB03OT.f MB04ND.f MB04OD.f SB03OR.f \ - SB03OY.f SB04PX.f MB04NY.f MB04OY.f SB03OV.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -mkoctfile slsb16cd.cc \ - SB16CD.f SB16CY.f AB09IX.f SB03OD.f MB02UD.f \ - AB09DD.f MA02AD.f MB03UD.f select.f SB03OU.f \ - MB01SD.f SB03OT.f MB04ND.f MB04OD.f SB03OR.f \ - SB03OY.f SB04PX.f MB04NY.f MB04OY.f SB03OV.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -system ("rm *.o"); -cd (homedir); - Modified: trunk/octave-forge/main/control/devel/makefile_control.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_control.m 2012-02-26 12:08:31 UTC (rev 9697) +++ trunk/octave-forge/main/control/devel/makefile_control.m 2012-02-26 12:42:25 UTC (rev 9698) @@ -11,26 +11,8 @@ srcdir = [develdir, "/../src"]; cd (srcdir); -%{ -makefile_chol -makefile_conversions -makefile_h2syn -makefile_hankel -makefile_helpers -makefile_hinfsyn -makefile_lqr -makefile_lyap -makefile_minreal -makefile_ncfsyn -makefile_norm -makefile_place -makefile_scale -makefile_staircase -makefile_tustin -makefile_zero -%} - -system ("make realclean"); +system ("make realclean"); # recompile slicotlibrary.a +## system ("make clean"); system ("make -j4 all"); system ("rm *.o"); Deleted: trunk/octave-forge/main/control/devel/makefile_conversions.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_conversions.m 2012-02-26 12:08:31 UTC (rev 9697) +++ trunk/octave-forge/main/control/devel/makefile_conversions.m 2012-02-26 12:42:25 UTC (rev 9698) @@ -1,24 +0,0 @@ -## ============================================================================== -## Developer Makefile for OCT-files -## ============================================================================== -## USAGE: * fetch control from Octave-Forge by svn -## * add control/inst, control/src and control/devel to your Octave path -## * run makefile_* -## ============================================================================== - -homedir = pwd (); -develdir = fileparts (which ("makefile_conversions")); -srcdir = [develdir, "/../src"]; -cd (srcdir); - -## state-space to transfer function - -## descriptor to regular state-space - -## transfer function to state-space - -system ("rm -rf *.o sltb04bd.oct slsb10jd.oct sltd04ad.oct"); -system ("make sltb04bd.oct slsb10jd.oct sltd04ad.oct"); - -system ("rm -rf *.o"); -cd (homedir); Deleted: trunk/octave-forge/main/control/devel/makefile_h2syn.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_h2syn.m 2012-02-26 12:08:31 UTC (rev 9697) +++ trunk/octave-forge/main/control/devel/makefile_h2syn.m 2012-02-26 12:42:25 UTC (rev 9698) @@ -1,34 +0,0 @@ -## ============================================================================== -## Developer Makefile for OCT-files -## ============================================================================== -## USAGE: * fetch control from Octave-Forge by svn -## * add control/inst, control/src and control/devel to your Octave path -## * run makefile_* -## ============================================================================== - -homedir = pwd (); -develdir = fileparts (which ("makefile_h2syn")); -srcdir = [develdir, "/../src"]; -cd (srcdir); - -mkoctfile slsb10hd.cc \ - SB10HD.f SB10UD.f SB10VD.f SB10WD.f SB02RD.f \ - MB01RU.f SB02SD.f MA02ED.f SB02RU.f SB02MR.f \ - MB01SD.f SB02MS.f SB02MV.f SB02MW.f MA02AD.f \ - SB02QD.f MB02PD.f SB03QX.f SB03QY.f MB01RX.f \ - MB01RY.f SB03SX.f SB03SY.f select.f SB03MX.f \ - SB03MY.f MB01UD.f SB03MV.f SB03MW.f SB04PX.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -mkoctfile slsb10ed.cc \ - SB10ED.f SB10SD.f SB10TD.f SB10PD.f MB01RX.f \ - SB02SD.f SB02OD.f MB01RU.f SB02OU.f SB02OV.f \ - SB02OW.f MB01RY.f SB02OY.f SB03SX.f SB03SY.f \ - MA02ED.f select.f SB03MX.f SB02MR.f SB02MV.f \ - MB01UD.f SB03MV.f SB04PX.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -system ("rm *.o"); -cd (homedir); Deleted: trunk/octave-forge/main/control/devel/makefile_hankel.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_hankel.m 2012-02-26 12:08:31 UTC (rev 9697) +++ trunk/octave-forge/main/control/devel/makefile_hankel.m 2012-02-26 12:42:25 UTC (rev 9698) @@ -1,24 +0,0 @@ -## ============================================================================== -## Developer Makefile for OCT-files -## ============================================================================== -## USAGE: * fetch control from Octave-Forge by svn -## * add control/inst, control/src and control/devel to your Octave path -## * run makefile_* -## ============================================================================== - -homedir = pwd (); -develdir = fileparts (which ("makefile_hankel")); -srcdir = [develdir, "/../src"]; -cd (srcdir); - -mkoctfile slab13ad.cc \ - AB13AD.f TB01ID.f TB01KD.f AB13AX.f MA02DD.f \ - MB03UD.f TB01LD.f SB03OU.f MB03QX.f select.f \ - SB03OT.f MB03QD.f MB04ND.f MB04OD.f MB03QY.f \ - SB03OR.f SB03OY.f SB04PX.f MB04NY.f MB04OY.f \ - SB03OV.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -system ("rm *.o"); -cd (homedir); Deleted: trunk/octave-forge/main/control/devel/makefile_hinfsyn.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_hinfsyn.m 2012-02-26 12:08:31 UTC (rev 9697) +++ trunk/octave-forge/main/control/devel/makefile_hinfsyn.m 2012-02-26 12:42:25 UTC (rev 9698) @@ -1,34 +0,0 @@ -## ============================================================================== -## Developer Makefile for OCT-files -## ============================================================================== -## USAGE: * fetch control from Octave-Forge by svn -## * add control/inst, control/src and control/devel to your Octave path -## * run makefile_* -## ============================================================================== - -homedir = pwd (); -develdir = fileparts (which ("makefile_hinfsyn")); -srcdir = [develdir, "/../src"]; -cd (srcdir); - -mkoctfile slsb10fd.cc \ - SB10FD.f SB10PD.f SB10QD.f SB10RD.f SB02RD.f \ - MB01RU.f MB01RX.f MA02AD.f SB02SD.f MA02ED.f \ - SB02RU.f SB02MR.f MB01SD.f SB02MS.f SB02MV.f \ - SB02MW.f SB02QD.f MB02PD.f SB03QX.f SB03QY.f \ - MB01RY.f SB03SX.f SB03SY.f select.f SB03MX.f \ - SB03MY.f MB01UD.f SB03MV.f SB03MW.f SB04PX.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -mkoctfile slsb10dd.cc \ - SB10DD.f MB01RU.f MB01RX.f SB02SD.f SB02OD.f \ - MA02AD.f SB02OU.f SB02OV.f SB02OW.f MB01RY.f \ - SB02OY.f SB03SX.f SB03SY.f MA02ED.f select.f \ - SB03MX.f SB02MR.f SB02MV.f MB01UD.f SB03MV.f \ - SB04PX.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -system ("rm *.o"); -cd (homedir); Deleted: trunk/octave-forge/main/control/devel/makefile_lqr.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_lqr.m 2012-02-26 12:08:31 UTC (rev 9697) +++ trunk/octave-forge/main/control/devel/makefile_lqr.m 2012-02-26 12:42:25 UTC (rev 9698) @@ -1,27 +0,0 @@ -## ============================================================================== -## Developer Makefile for OCT-files -## ============================================================================== -## USAGE: * fetch control from Octave-Forge by svn -## * add control/inst, control/src and control/devel to your Octave path -## * run makefile_* -## ============================================================================== - -homedir = pwd (); -develdir = fileparts (which ("makefile_lqr")); -srcdir = [develdir, "/../src"]; -cd (srcdir); - -mkoctfile slsb02od.cc \ - SB02OD.f SB02OU.f SB02OV.f SB02OW.f SB02OY.f \ - SB02MR.f SB02MV.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -mkoctfile slsg02ad.cc \ - SG02AD.f SB02OU.f SB02OV.f SB02OW.f SB02OY.f \ - MB01SD.f MB02VD.f MB02PD.f MA02GD.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -system ("rm *.o"); -cd (homedir); Deleted: trunk/octave-forge/main/control/devel/makefile_lyap.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_lyap.m 2012-02-26 12:08:31 UTC (rev 9697) +++ trunk/octave-forge/main/control/devel/makefile_lyap.m 2012-02-26 12:42:25 UTC (rev 9698) @@ -1,42 +0,0 @@ -## ============================================================================== -## Developer Makefile for OCT-files -## ============================================================================== -## USAGE: * fetch control from Octave-Forge by svn -## * add control/inst, control/src and control/devel to your Octave path -## * run makefile_* -## ============================================================================== - -homedir = pwd (); -develdir = fileparts (which ("makefile_lyap")); -srcdir = [develdir, "/../src"]; -cd (srcdir); - -## Lypunov -mkoctfile slsb03md.cc \ - SB03MD.f select.f SB03MX.f SB03MY.f MB01RD.f \ - SB03MV.f SB03MW.f SB04PX.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - - -## Sylvester -mkoctfile slsb04md.cc \ - SB04MD.f SB04MU.f SB04MY.f SB04MR.f SB04MW.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -mkoctfile slsb04qd.cc \ - SB04QD.f SB04QU.f SB04QY.f SB04MW.f SB04QR.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - - -## Generalized Lyapunov -mkoctfile slsg03ad.cc \ - SG03AD.f MB01RW.f MB01RD.f SG03AX.f SG03AY.f \ - MB02UU.f MB02UV.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -system ("rm *.o"); -cd (homedir); Deleted: trunk/octave-forge/main/control/devel/makefile_minreal.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_minreal.m 2012-02-26 12:08:31 UTC (rev 9697) +++ trunk/octave-forge/main/control/devel/makefile_minreal.m 2012-02-26 12:42:25 UTC (rev 9698) @@ -1,26 +0,0 @@ -## ============================================================================== -## Developer Makefile for OCT-files -## ============================================================================== -## USAGE: * fetch control from Octave-Forge by svn -## * add control/inst, control/src and control/devel to your Octave path -## * run makefile_* -## ============================================================================== - -homedir = pwd (); -develdir = fileparts (which ("makefile_minreal")); -srcdir = [develdir, "/../src"]; -cd (srcdir); - -mkoctfile sltb01pd.cc \ - TB01PD.f TB01XD.f TB01ID.f AB07MD.f TB01UD.f \ - MB03OY.f MB01PD.f MB01QD.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -mkoctfile sltg01jd.cc \ - TG01JD.f TG01AD.f TB01XD.f MA02CD.f TG01HX.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -system ("rm *.o"); -cd (homedir); Deleted: trunk/octave-forge/main/control/devel/makefile_modred.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_modred.m 2012-02-26 12:08:31 UTC (rev 9697) +++ trunk/octave-forge/main/control/devel/makefile_modred.m 2012-02-26 12:42:25 UTC (rev 9698) @@ -1,32 +0,0 @@ -## ============================================================================== -## Developer Makefile for OCT-files -## ============================================================================== -## USAGE: * fetch control-devel from Octave-Forge by svn -## * add control-devel/inst, control-devel/src and control-devel/devel -## to your Octave path (by an .octaverc file) -## * run makefile_modred -## ============================================================================== - -homedir = pwd (); -develdir = fileparts (which ("makefile_modred")); -srcdir = [develdir, "/../src"]; -cd (srcdir); - -mkoctfile slab09hd.cc \ - slicotlibrary.a \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -mkoctfile slab09id.cc \ - slicotlibrary.a \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -mkoctfile slab09jd.cc \ - slicotlibrary.a \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -system ("rm *.o"); -cd (homedir); - Deleted: trunk/octave-forge/main/control/devel/makefile_ncfsyn.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_ncfsyn.m 2012-02-26 12:08:31 UTC (rev 9697) +++ trunk/octave-forge/main/control/devel/makefile_ncfsyn.m 2012-02-26 12:42:25 UTC (rev 9698) @@ -1,41 +0,0 @@ -## ============================================================================== -## Developer Makefile for OCT-files -## ============================================================================== -## USAGE: * fetch control from Octave-Forge by svn -## * add control/inst, control/src and control/devel to your Octave path -## * run makefile_* -## ============================================================================== - -homedir = pwd (); -develdir = fileparts (which ("makefile_ncfsyn")); -srcdir = [develdir, "/../src"]; -cd (srcdir); - -## H-infinity loop shaping - continuous-time -mkoctfile slsb10id.cc \ - SB10ID.f SB02RD.f select.f SB10JD.f MB02VD.f \ - MA02GD.f SB02MS.f MA02ED.f SB02RU.f SB02SD.f \ - MB01RU.f SB02QD.f SB02MV.f SB02MW.f SB02MR.f \ - MA02AD.f MB02PD.f MB01SD.f MB01UD.f SB03SY.f \ - MB01RX.f SB03MX.f SB03SX.f MB01RY.f SB03QY.f \ - SB03QX.f SB03MY.f SB04PX.f SB03MV.f SB03MW.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -## H-infinity loop shaping - discrete-time - strictly proper case -mkoctfile slsb10kd.cc \ - SB10KD.f SB02OD.f select.f SB02OY.f SB02OW.f \ - SB02OV.f SB02MV.f SB02OU.f SB02MR.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -## H-infinity loop shaping - discrete-time - proper case -mkoctfile slsb10zd.cc \ - SB10ZD.f MA02AD.f SB02OD.f select.f MB01RX.f \ - MB02VD.f SB02OY.f SB02OW.f SB02OV.f SB02OU.f \ - SB02MR.f MA02GD.f SB02MV.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -system ("rm *.o"); -cd (homedir); Deleted: trunk/octave-forge/main/control/devel/makefile_norm.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_norm.m 2012-02-26 12:08:31 UTC (rev 9697) +++ trunk/octave-forge/main/control/devel/makefile_norm.m 2012-02-26 12:42:25 UTC (rev 9698) @@ -1,36 +0,0 @@ -## ============================================================================== -## Developer Makefile for OCT-files -## ============================================================================== -## USAGE: * fetch control from Octave-Forge by svn -## * add control/inst, control/src and control/devel to your Octave path -## * run makefile_* -## ============================================================================== - -homedir = pwd (); -develdir = fileparts (which ("makefile_norm")); -srcdir = [develdir, "/../src"]; -cd (srcdir); - -## H-2 norm -mkoctfile slab13bd.cc \ - AB13BD.f SB08DD.f SB03OU.f SB01FY.f TB01LD.f \ - SB03OT.f MB04ND.f MB04OD.f MB03QX.f select.f \ - SB03OR.f MB04OX.f MB03QD.f SB03OY.f MA02AD.f \ - MB03QY.f SB04PX.f MB04NY.f MB04OY.f SB03OV.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -## L-inf norm -mkoctfile slab13dd.cc \ - AB13DD.f MA02AD.f MB01SD.f MB03XD.f TB01ID.f \ - TG01AD.f TG01BD.f AB13DX.f MA01AD.f MA02ID.f \ - MB03XP.f MB04DD.f MB04QB.f MB04TB.f MB03XU.f \ - MB04TS.f UE01MD.f MB02RD.f MB02SD.f MB04QC.f \ - MB04QF.f MB03YA.f MB03YD.f MB02RZ.f MB04QU.f \ - MB02SZ.f MB03YT.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" \ - "$(mkoctfile -p FLIBS)" - -system ("rm *.o"); -cd (homedir); Deleted: trunk/octave-forge/main/control/devel/makefile_place.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_place.m 2012-02-26 12:08:31 UTC (rev 9697) +++ trunk/octave-forge/main/control/devel/makefile_place.m 2012-02-26 12:42:25 UTC (rev 9698) @@ -1,21 +0,0 @@ -## ============================================================================== -## Developer Makefile for OCT-files -## ============================================================================== -## USAGE: * fetch control from Octave-Forge by svn -## * add control/inst, control/src and control/devel to your Octave path -## * run makefile_* -## ============================================================================== - -homedir = pwd (); -develdir = fileparts (which ("makefile_place")); -srcdir = [develdir, "/../src"]; -cd (srcdir); - -mkoctfile slsb01bd.cc \ - SB01BD.f MB03QD.f MB03QY.f SB01BX.f SB01BY.f \ - select.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -system ("rm *.o"); -cd (homedir); Deleted: trunk/octave-forge/main/control/devel/makefile_scale.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_scale.m 2012-02-26 12:08:31 UTC (rev 9697) +++ trunk/octave-forge/main/control/devel/makefile_scale.m 2012-02-26 12:42:25 UTC (rev 9698) @@ -1,27 +0,0 @@ -## ============================================================================== -## Developer Makefile for OCT-files -## ============================================================================== -## USAGE: * fetch control from Octave-Forge by svn -## * add control/inst, control/src and control/devel to your Octave path -## * run makefile_* -## ============================================================================== - -homedir = pwd (); -develdir = fileparts (which ("makefile_scale")); -srcdir = [develdir, "/../src"]; -cd (srcdir); - -## scaling of state-space models -mkoctfile sltb01id.cc \ - TB01ID.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -## scaling of descriptor state-space models -mkoctfile sltg01ad.cc \ - TG01AD.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -system ("rm *.o"); -cd (homedir); Deleted: trunk/octave-forge/main/control/devel/makefile_staircase.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_staircase.m 2012-02-26 12:08:31 UTC (rev 9697) +++ trunk/octave-forge/main/control/devel/makefile_staircase.m 2012-02-26 12:42:25 UTC (rev 9698) @@ -1,40 +0,0 @@ -## ============================================================================== -## Developer Makefile for OCT-files -## ============================================================================== -## USAGE: * fetch control from Octave-Forge by svn -## * add control/inst, control/src and control/devel to your Octave path -## * run makefile_* -## ============================================================================== - -homedir = pwd (); -develdir = fileparts (which ("makefile_staircase")); -srcdir = [develdir, "/../src"]; -cd (srcdir); - -## staircase form using orthogonal transformations -mkoctfile slab01od.cc \ - AB01OD.f AB01ND.f MB03OY.f MB01PD.f MB01QD.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -## controllability staircase form of descriptor state-space models -mkoctfile sltg01hd.cc \ - TG01HD.f TG01HX.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -## observability staircase form of descriptor state-space models -mkoctfile sltg01id.cc \ - TG01ID.f TB01XD.f MA02CD.f AB07MD.f TG01HX.f \ - MA02BD.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -## controllable block Hessenberg realization -mkoctfile sltb01ud.cc \ - TB01UD.f MB01PD.f MB03OY.f MB01QD.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -system ("rm *.o"); -cd (homedir); \ No newline at end of file Deleted: trunk/octave-forge/main/control/devel/makefile_tustin.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_tustin.m 2012-02-26 12:08:31 UTC (rev 9697) +++ trunk/octave-forge/main/control/devel/makefile_tustin.m 2012-02-26 12:42:25 UTC (rev 9698) @@ -1,21 +0,0 @@ -## ============================================================================== -## Developer Makefile for OCT-files -## ============================================================================== -## USAGE: * fetch control from Octave-Forge by svn -## * add control/inst, control/src and control/devel to your Octave path -## * run makefile_* -## ============================================================================== - -homedir = pwd (); -develdir = fileparts (which ("makefile_tustin")); -srcdir = [develdir, "/../src"]; -cd (srcdir); - -## bilinear transformation -mkoctfile slab04md.cc \ - AB04MD.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -system ("rm *.o"); -cd (homedir); Deleted: trunk/octave-forge/main/control/devel/makefile_zero.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_zero.m 2012-02-26 12:08:31 UTC (rev 9697) +++ trunk/octave-forge/main/control/devel/makefile_zero.m 2012-02-26 12:42:25 UTC (rev 9698) @@ -1,34 +0,0 @@ -## ============================================================================== -## Developer Makefile for OCT-files -## ============================================================================== -## USAGE: * fetch control from Octave-Forge by svn -## * add control/inst, control/src and control/devel to your Octave path -## * run makefile_* -## ============================================================================== - -homedir = pwd (); -develdir = fileparts (which ("makefile_zero")); -srcdir = [develdir, "/../src"]; -cd (srcdir); - -## transmission zeros of state-space models -mkoctfile slab08nd.cc \ - AB08ND.f AB08NX.f TB01ID.f MB03OY.f MB03PY.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -## transmission zeros of descriptor state-space models -mkoctfile slag08bd.cc \ - AG08BD.f AG08BY.f TG01AD.f TB01XD.f MA02CD.f \ - TG01FD.f MA02BD.f MB03OY.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -## gain of descriptor state-space models -mkoctfile sltg04bx.cc \ - TG04BX.f MB02RD.f MB02SD.f \ - "$(mkoctfile -p LAPACK_LIBS)" \ - "$(mkoctfile -p BLAS_LIBS)" - -system ("rm *.o"); -cd (homedir); This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-04-24 06:58:37
|
Revision: 10318 http://octave.svn.sourceforge.net/octave/?rev=10318&view=rev Author: paramaniac Date: 2012-04-24 06:58:30 +0000 (Tue, 24 Apr 2012) Log Message: ----------- control: draft code for James B. Rawlings Modified Paths: -------------- trunk/octave-forge/main/control/devel/lqe.m Added Paths: ----------- trunk/octave-forge/main/control/devel/dlqe.m Added: trunk/octave-forge/main/control/devel/dlqe.m =================================================================== --- trunk/octave-forge/main/control/devel/dlqe.m (rev 0) +++ trunk/octave-forge/main/control/devel/dlqe.m 2012-04-24 06:58:30 UTC (rev 10318) @@ -0,0 +1,22 @@ +function [l, p, z, e] = dlqe (a, g, c, q, r, s = []) + + if (nargin < 5 || nargin > 6) + print_usage (); + endif + + if (isempty (s)) + [~, p, e] = dlqr (a.', c.', g*q*g.', r); + else + [~, p, e] = dlqr (a.', c.', g*q*g.', r, g*s); + endif + + ## k computed by dlqr would be + ## k = (r + c*p*c.') \ (c*p*a.' + s.') + ## such that l = a \ k.' + ## what about the s term? + + l = p*c.' / (c*p*c.' + r); + ## z = p - p*c.' / (c*p*c.' + r) * c*p; + z = p - l*c*p; + +endfunction \ No newline at end of file Modified: trunk/octave-forge/main/control/devel/lqe.m =================================================================== --- trunk/octave-forge/main/control/devel/lqe.m 2012-04-23 20:22:00 UTC (rev 10317) +++ trunk/octave-forge/main/control/devel/lqe.m 2012-04-24 06:58:30 UTC (rev 10318) @@ -1,15 +1,19 @@ -function [g, x, l] = lqe (a, c, q, r = [], s = [], e = []) +function [l, p, e] = lqe (a, g, c, q = [], r = [], s = []) if (nargin < 3 || nargin > 6) print_usage (); endif if (isa (a, "lti")) - [g, x, l] = lqr (a.', c, q, r); # lqr (sys.', q, r, s) + [l, p, e] = lqr (a.', g, c, q); # lqe (sys.', q, r, s), g=I, works like lqr (sys.', q, r, s).' + elseif (isempty (g)) + [l, p, e] = lqr (a.', c.', q, r, s); # lqe (a, [], c, q, r, s), g=I, works like lqr (a.', c.', q, r, s).' + elseif (isempty (s)) + [l, p, e] = lqr (a.', c.', g*q*g.', r); else - [g, x, l] = lqr (a.', c.', q, r, s, e.'); + [l, p, e] = lqr (a.', c.', g*q*g.', r, g*s); endif - - g = g.' + l = l.'; + endfunction \ No newline at end of file This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-05-19 22:16:03
|
Revision: 10465 http://octave.svn.sourceforge.net/octave/?rev=10465&view=rev Author: paramaniac Date: 2012-05-19 22:15:57 +0000 (Sat, 19 May 2012) Log Message: ----------- control: work on docstrings Modified Paths: -------------- trunk/octave-forge/main/control/devel/dlqe.m trunk/octave-forge/main/control/devel/lqe.m Modified: trunk/octave-forge/main/control/devel/dlqe.m =================================================================== --- trunk/octave-forge/main/control/devel/dlqe.m 2012-05-19 08:39:04 UTC (rev 10464) +++ trunk/octave-forge/main/control/devel/dlqe.m 2012-05-19 22:15:57 UTC (rev 10465) @@ -1,3 +1,77 @@ +## Copyright (C) 2012 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 LTI Syncope. If not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{g}, @var{x}, @var{l}] =} lqr (@var{sys}, @var{q}, @var{r}) +## @deftypefnx {Function File} {[@var{g}, @var{x}, @var{l}] =} lqr (@var{sys}, @var{q}, @var{r}, @var{s}) +## @deftypefnx {Function File} {[@var{g}, @var{x}, @var{l}] =} lqr (@var{a}, @var{b}, @var{q}, @var{r}) +## @deftypefnx {Function File} {[@var{g}, @var{x}, @var{l}] =} lqr (@var{a}, @var{b}, @var{q}, @var{r}, @var{s}) +## @deftypefnx {Function File} {[@var{g}, @var{x}, @var{l}] =} lqr (@var{a}, @var{b}, @var{q}, @var{r}, @var{[]}, @var{e}) +## @deftypefnx {Function File} {[@var{g}, @var{x}, @var{l}] =} lqr (@var{a}, @var{b}, @var{q}, @var{r}, @var{s}, @var{e}) +## Linear-quadratic regulator. +## +## @strong{Inputs} +## @table @var +## @item sys +## Continuous or discrete-time LTI model. +## @item a +## State transition matrix of continuous-time system. +## @item b +## Input matrix of continuous-time system. +## @item q +## State weighting matrix. +## @item r +## Input weighting matrix. +## @item s +## Optional cross term matrix. If @var{s} is not specified, a zero matrix is assumed. +## @item e +## Optional descriptor matrix. If @var{e} is not specified, an identity matrix is assumed. +## @end table +## +## @strong{Outputs} +## @table @var +## @item g +## State feedback matrix. +## @item x +## Unique stabilizing solution of the continuous-time Riccati equation. +## @item l +## Closed-loop poles. +## @end table +## +## @strong{Equations} +## @example +## @group +## . +## x = A x + B u, x(0) = x0 +## +## inf +## J(x0) = INT (x' Q x + u' R u + 2 x' S u) dt +## 0 +## +## L = eig (A - B*G) +## @end group +## @end example +## @seealso{care, dare, dlqr} +## @end deftypefn + +## Author: Lukas Reichlin <luk...@gm...> +## Created: April 2012 +## Version: 0.1 + function [l, p, z, e] = dlqe (a, g, c, q, r, s = []) if (nargin < 5 || nargin > 6) @@ -19,4 +93,4 @@ ## z = p - p*c.' / (c*p*c.' + r) * c*p; z = p - l*c*p; -endfunction \ No newline at end of file +endfunction Modified: trunk/octave-forge/main/control/devel/lqe.m =================================================================== --- trunk/octave-forge/main/control/devel/lqe.m 2012-05-19 08:39:04 UTC (rev 10464) +++ trunk/octave-forge/main/control/devel/lqe.m 2012-05-19 22:15:57 UTC (rev 10465) @@ -1,3 +1,77 @@ +## Copyright (C) 2012 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 LTI Syncope. If not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{g}, @var{x}, @var{l}] =} lqe (@var{sys}, @var{q}, @var{r}) +## @deftypefnx {Function File} {[@var{g}, @var{x}, @var{l}] =} lqe (@var{sys}, @var{q}, @var{r}, @var{s}) +## @deftypefnx {Function File} {[@var{g}, @var{x}, @var{l}] =} lqe (@var{a}, @var{g}, @var{c}, @var{q}, @var{r}) +## @deftypefnx {Function File} {[@var{g}, @var{x}, @var{l}] =} lqe (@var{a}, @var{g}, @var{c}, @var{q}, @var{r}, @var{s}) +## @deftypefnx {Function File} {[@var{g}, @var{x}, @var{l}] =} lqe (@var{a}, @var{[]}, @var{c}, @var{q}, @var{r}) +## @deftypefnx {Function File} {[@var{g}, @var{x}, @var{l}] =} lqe (@var{a}, @var{[]}, @var{c}, @var{q}, @var{r}, @var{s}) +## Linear-quadratic estimator. +## +## @strong{Inputs} +## @table @var +## @item sys +## Continuous or discrete-time LTI model. +## @item a +## State transition matrix of continuous-time system. +## @item b +## Input matrix of continuous-time system. +## @item q +## State weighting matrix. +## @item r +## Input weighting matrix. +## @item s +## Optional cross term matrix. If @var{s} is not specified, a zero matrix is assumed. +## @item e +## Optional descriptor matrix. If @var{e} is not specified, an identity matrix is assumed. +## @end table +## +## @strong{Outputs} +## @table @var +## @item l +## Observer gain matrix. +## @item p +## Unique stabilizing solution of the continuous-time Riccati equation. +## @item e +## Closed-loop poles. +## @end table +## +## @strong{Equations} +## @example +## @group +## . +## x = A x + B u, x(0) = x0 +## +## inf +## J(x0) = INT (x' Q x + u' R u + 2 x' S u) dt +## 0 +## +## L = eig (A - B*G) +## @end group +## @end example +## @seealso{care, dare, dlqr} +## @end deftypefn + +## Author: Lukas Reichlin <luk...@gm...> +## Created: April 2012 +## Version: 0.1 + function [l, p, e] = lqe (a, g, c, q = [], r = [], s = []) if (nargin < 3 || nargin > 6) @@ -16,4 +90,4 @@ l = l.'; -endfunction \ No newline at end of file +endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-06-05 06:48:13
|
Revision: 10563 http://octave.svn.sourceforge.net/octave/?rev=10563&view=rev Author: paramaniac Date: 2012-06-05 06:48:07 +0000 (Tue, 05 Jun 2012) Log Message: ----------- control: touch up lqe doc Modified Paths: -------------- trunk/octave-forge/main/control/devel/dlqe.m trunk/octave-forge/main/control/devel/lqe.m Modified: trunk/octave-forge/main/control/devel/dlqe.m =================================================================== --- trunk/octave-forge/main/control/devel/dlqe.m 2012-06-04 20:12:13 UTC (rev 10562) +++ trunk/octave-forge/main/control/devel/dlqe.m 2012-06-05 06:48:07 UTC (rev 10563) @@ -16,13 +16,11 @@ ## along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn {Function File} {[@var{g}, @var{x}, @var{l}] =} lqr (@var{sys}, @var{q}, @var{r}) -## @deftypefnx {Function File} {[@var{g}, @var{x}, @var{l}] =} lqr (@var{sys}, @var{q}, @var{r}, @var{s}) -## @deftypefnx {Function File} {[@var{g}, @var{x}, @var{l}] =} lqr (@var{a}, @var{b}, @var{q}, @var{r}) -## @deftypefnx {Function File} {[@var{g}, @var{x}, @var{l}] =} lqr (@var{a}, @var{b}, @var{q}, @var{r}, @var{s}) -## @deftypefnx {Function File} {[@var{g}, @var{x}, @var{l}] =} lqr (@var{a}, @var{b}, @var{q}, @var{r}, @var{[]}, @var{e}) -## @deftypefnx {Function File} {[@var{g}, @var{x}, @var{l}] =} lqr (@var{a}, @var{b}, @var{q}, @var{r}, @var{s}, @var{e}) -## Linear-quadratic regulator. +## @deftypefn {Function File} {[@var{l}, @var{p}, @var{z}, @var{e}] =} lqe (@var{a}, @var{g}, @var{c}, @var{q}, @var{r}) +## @deftypefnx {Function File} {[@var{l}, @var{p}, @var{z}, @var{e}] =} lqe (@var{a}, @var{g}, @var{c}, @var{q}, @var{r}, @var{s}) +## @deftypefnx {Function File} {[@var{l}, @var{p}, @var{z}, @var{e}] =} lqe (@var{a}, @var{[]}, @var{c}, @var{q}, @var{r}) +## @deftypefnx {Function File} {[@var{l}, @var{p}, @var{z}, @var{e}] =} lqe (@var{a}, @var{[]}, @var{c}, @var{q}, @var{r}, @var{s}) +## Linear-quadratic estimator for discrete-time systems. ## ## @strong{Inputs} ## @table @var @@ -78,7 +76,9 @@ print_usage (); endif - if (isempty (s)) + if (isempty (g)) + [~, p, e] = dlqr (a.', c.', q, r, s); # dlqe (a, [], c, q, r, s), g=I + elseif (isempty (s)) [~, p, e] = dlqr (a.', c.', g*q*g.', r); else [~, p, e] = dlqr (a.', c.', g*q*g.', r, g*s); Modified: trunk/octave-forge/main/control/devel/lqe.m =================================================================== --- trunk/octave-forge/main/control/devel/lqe.m 2012-06-04 20:12:13 UTC (rev 10562) +++ trunk/octave-forge/main/control/devel/lqe.m 2012-06-05 06:48:07 UTC (rev 10563) @@ -16,12 +16,12 @@ ## along with LTI Syncope. If not, see <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn {Function File} {[@var{g}, @var{x}, @var{l}] =} lqe (@var{sys}, @var{q}, @var{r}) -## @deftypefnx {Function File} {[@var{g}, @var{x}, @var{l}] =} lqe (@var{sys}, @var{q}, @var{r}, @var{s}) -## @deftypefnx {Function File} {[@var{g}, @var{x}, @var{l}] =} lqe (@var{a}, @var{g}, @var{c}, @var{q}, @var{r}) -## @deftypefnx {Function File} {[@var{g}, @var{x}, @var{l}] =} lqe (@var{a}, @var{g}, @var{c}, @var{q}, @var{r}, @var{s}) -## @deftypefnx {Function File} {[@var{g}, @var{x}, @var{l}] =} lqe (@var{a}, @var{[]}, @var{c}, @var{q}, @var{r}) -## @deftypefnx {Function File} {[@var{g}, @var{x}, @var{l}] =} lqe (@var{a}, @var{[]}, @var{c}, @var{q}, @var{r}, @var{s}) +## @deftypefn {Function File} {[@var{l}, @var{p}, @var{e}] =} lqe (@var{sys}, @var{q}, @var{r}) +## @deftypefnx {Function File} {[@var{l}, @var{p}, @var{e}] =} lqe (@var{sys}, @var{q}, @var{r}, @var{s}) +## @deftypefnx {Function File} {[@var{l}, @var{p}, @var{e}] =} lqe (@var{a}, @var{g}, @var{c}, @var{q}, @var{r}) +## @deftypefnx {Function File} {[@var{l}, @var{p}, @var{e}] =} lqe (@var{a}, @var{g}, @var{c}, @var{q}, @var{r}, @var{s}) +## @deftypefnx {Function File} {[@var{l}, @var{p}, @var{e}] =} lqe (@var{a}, @var{[]}, @var{c}, @var{q}, @var{r}) +## @deftypefnx {Function File} {[@var{l}, @var{p}, @var{e}] =} lqe (@var{a}, @var{[]}, @var{c}, @var{q}, @var{r}, @var{s}) ## Linear-quadratic estimator. ## ## @strong{Inputs} @@ -79,7 +79,7 @@ endif if (isa (a, "lti")) - [l, p, e] = lqr (a.', g, c, q); # lqe (sys.', q, r, s), g=I, works like lqr (sys.', q, r, s).' + [l, p, e] = lqr (a.', g, c, q); # lqe (sys, q, r, s), g=I, works like lqr (sys.', q, r, s).' elseif (isempty (g)) [l, p, e] = lqr (a.', c.', q, r, s); # lqe (a, [], c, q, r, s), g=I, works like lqr (a.', c.', q, r, s).' elseif (isempty (s)) This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-06-08 17:47:50
|
Revision: 10600 http://octave.svn.sourceforge.net/octave/?rev=10600&view=rev Author: paramaniac Date: 2012-06-08 17:47:43 +0000 (Fri, 08 Jun 2012) Log Message: ----------- control: remove cruft Removed Paths: ------------- trunk/octave-forge/main/control/devel/makefile_helpers.m trunk/octave-forge/main/control/devel/minreal.m trunk/octave-forge/main/control/devel/test_ctrbf.m trunk/octave-forge/main/control/devel/tf2ss_draft.m Deleted: trunk/octave-forge/main/control/devel/makefile_helpers.m =================================================================== --- trunk/octave-forge/main/control/devel/makefile_helpers.m 2012-06-08 16:45:39 UTC (rev 10599) +++ trunk/octave-forge/main/control/devel/makefile_helpers.m 2012-06-08 17:47:43 UTC (rev 10600) @@ -1,23 +0,0 @@ -## ============================================================================== -## Developer Makefile for OCT-files -## ============================================================================== -## USAGE: * fetch control from Octave-Forge by svn -## * add control/inst, control/src and control/devel to your Octave path -## * run makefile_* -## ============================================================================== - -homedir = pwd (); -develdir = fileparts (which ("makefile_helpers")); -srcdir = [develdir, "/../src"]; -cd (srcdir); - -mkoctfile is_real_scalar.cc - -mkoctfile is_real_vector.cc - -mkoctfile is_real_matrix.cc - -mkoctfile is_real_square_matrix.cc - -system ("rm *.o"); -cd (homedir); \ No newline at end of file Deleted: trunk/octave-forge/main/control/devel/minreal.m =================================================================== --- trunk/octave-forge/main/control/devel/minreal.m 2012-06-08 16:45:39 UTC (rev 10599) +++ trunk/octave-forge/main/control/devel/minreal.m 2012-06-08 17:47:43 UTC (rev 10600) @@ -1,45 +0,0 @@ -## Copyright (C) 2010 Benjamin Fernandez <ma...@be...> -## -## This program 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 2 of the License, or -## (at your option) any later version. -## -## This program 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 Octave; see the file COPYING. If not, see -## <http://www.gnu.org/licenses/>. - -## -*- texinfo -*- -## @deftypefn{Function File} {[@var{Amin}, @var{Bmin}, @var{Cmin}] =} ctrbf (@var{A}, @var{B}, @var{C}) -## @deftypefnx{Function File} {[@var{Amin}, @var{Bmin}, @var{Cmin}] =} ctrbf (@var{A}, @var{B}, @var{C}, @var{TOL}) -## Minimal realization of the system (A,B,C). -## If the system is controlable and observable, the syste es minimal. -## If the system is not controlable or/and observable, a new system created -## with the controlabe and observable subspace is equivalent -## which means that Cco(sI-Aco)^(-1)Bco = C(sI-A)^(-1)B. -## @end deftypefn - -## Author: Benjamin Fernandez <ma...@be...> -## Created: 2010-07-10 - -function [Amin,Bmin,Cmin] = minreal(A,B,C,TOL) - if(nargin<3 || nargin>4) - print_usage(); - endif - n = length(A); - if(nargin == 3) - TOL = n*norm(A,1)*eps; - endif - [Abar,Bbar,Cbar,T,K,Ac,Bc,Cc,Cnc,lc] = ctrbf(A,B,C); - [Abar,Bbar,Cbar,T,K,Amin,Bmin,Cmin,Cno,lo] = obsvf(Ac,Bc,Cc); - l = lc+lo; - disp('States reduced:'); - disp(l); - -endfunction - Deleted: trunk/octave-forge/main/control/devel/test_ctrbf.m =================================================================== --- trunk/octave-forge/main/control/devel/test_ctrbf.m 2012-06-08 16:45:39 UTC (rev 10599) +++ trunk/octave-forge/main/control/devel/test_ctrbf.m 2012-06-08 17:47:43 UTC (rev 10600) @@ -1,84 +0,0 @@ -a = [ -1.0 0.0 0.0 - -2.0 -2.0 -2.0 - -1.0 0.0 -3.0 ]; - -b = [ 1.0 0.0 0.0 - 0.0 2.0 1.0 ].'; - - -c = [ 0.0 2.0 1.0 - 1.0 0.0 0.0 ]; - - -[ac, bc, cc, z, ncont] = sltb01ud (a, b, c, 0.0) - -z\a*z - -a = [ 1 1 - 4 -2 ]; - -b = [ 1 -1 - 1 -1 ]; - -c = [ 1 0 - 0 1 ]; - -[ac, bc, cc, z, ncont] = sltb01ud (a, b, c, 0.0) - -%{ - The transformed state dynamics matrix of a controllable realization is - -3.0000 2.2361 - 0.0000 -1.0000 - - and the dimensions of its diagonal blocks are - 2 - - The transformed input/state matrix B of a controllable realization is - 0.0000 -2.2361 - 1.0000 0.0000 - - The transformed output/state matrix C of a controllable realization is - -2.2361 0.0000 - 0.0000 1.0000 - - The controllability index of the transformed system representation = 1 - - The similarity transformation matrix Z is - 0.0000 1.0000 0.0000 - -0.8944 0.0000 -0.4472 - -0.4472 0.0000 0.8944 -%} -%{ -A = - 1 1 - 4 -2 - -B = - 1 -1 - 1 -1 - -C = - 1 0 - 0 1 -and locate the uncontrollable mode. - -[Abar,Bbar,Cbar,T,k]=ctrbf(A,B,C) - -Abar = - -3.0000 0 - -3.0000 2.0000 - -Bbar = - 0.0000 0.0000 - 1.4142 -1.4142 - -Cbar = - -0.7071 0.7071 - 0.7071 0.7071 - -T = - -0.7071 0.7071 - 0.7071 0.7071 -k = - 1 0 -%} \ No newline at end of file Deleted: trunk/octave-forge/main/control/devel/tf2ss_draft.m =================================================================== --- trunk/octave-forge/main/control/devel/tf2ss_draft.m 2012-06-08 16:45:39 UTC (rev 10599) +++ trunk/octave-forge/main/control/devel/tf2ss_draft.m 2012-06-08 17:47:43 UTC (rev 10600) @@ -1,76 +0,0 @@ -function retsys = tf2ss_draft () - -num = {[1, 5, 7], [1]; [1, 7], [1, 5, 5]}; -den = {[1, 5, 6], [1, 2]; [1, 8, 6], [1, 3, 2]}; -sys = tf (num, den); - -%{ -sys = tf (1, [1, 0]) -sys = tf (1, [1, 1]) - -sys = tf (1, conv ([1, 1, 1], [1, 4, 6, 4, 1])) - -sys = tf () -sys = tf ("s") -%} - -sys = tf (WestlandLynx); -tol = sqrt (eps) - - [p, m] = size (sys); - [num, den] = tfdata (sys); - - numc = cell (p, m); - denc = cell (p, 1); - - ## multiply all denominators in a row and - ## update each numerator accordingly - for i = 1 : p - denc(i) = __conv__ (den{i,:}); - for j = 1 : m - idx = setdiff (1:m, j); - numc(i,j) = __conv__ (num{i,j}, den{i,idx}); - endfor - endfor - - len_numc = cellfun (@length, numc); - len_denc = cellfun (@length, denc); - - ## tfpoly ensures that there are no leading zeros - tmp = len_numc > repmat (len_denc, 1, m); - if (any (tmp(:))) - error ("tf: tf2ss: system must be proper"); - endif - - max_len_denc = max (len_denc(:)); - ucoeff = zeros (p, m, max_len_denc); - dcoeff = zeros (p, max_len_denc); - index = len_denc-1; - - for i = 1 : p - len = len_denc(i); - dcoeff(i, 1:len) = denc{i}; - for j = 1 : m - ucoeff(i, j, len-len_numc(i,j)+1 : len) = numc{i,j}; - endfor - endfor -index, prod (index), eps*prod (index) -min (sqrt (eps), eps*prod (index)) - [a, b, c, d] = sltd04ad (ucoeff, dcoeff, index, tol); - - retsys = ss (a, b, c, d); - -endfunction - - -function vec = __conv__ (vec, varargin) - - if (nargin == 1) - return; - else - for k = 1 : nargin-1 - vec = conv (vec, varargin{k}); - endfor - endif - -endfunction This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |
From: <par...@us...> - 2012-06-12 05:12:37
|
Revision: 10619 http://octave.svn.sourceforge.net/octave/?rev=10619&view=rev Author: paramaniac Date: 2012-06-12 05:12:31 +0000 (Tue, 12 Jun 2012) Log Message: ----------- control: add the sage worksheet I used to compute mimo transfer function interconnections Added Paths: ----------- trunk/octave-forge/main/control/devel/mimo_transfer_function_sage_worksheet.zip Property Changed: ---------------- trunk/octave-forge/main/control/devel/ Property changes on: trunk/octave-forge/main/control/devel ___________________________________________________________________ Deleted: svn:ignore - *.zip Added: trunk/octave-forge/main/control/devel/mimo_transfer_function_sage_worksheet.zip =================================================================== (Binary files differ) Property changes on: trunk/octave-forge/main/control/devel/mimo_transfer_function_sage_worksheet.zip ___________________________________________________________________ Added: svn:mime-type + application/octet-stream This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |