Menu
▾
▴
[Octave-cvsupdate] SF.net SVN: octave:[8255] trunk/octave-forge/main/control/devel
|
From: <par...@us...> - 2011-05-09 21:34:52
|
Revision: 8255
http://octave.svn.sourceforge.net/octave/?rev=8255&view=rev
Author: paramaniac
Date: 2011-05-09 21:34:44 +0000 (Mon, 09 May 2011)
Log Message:
-----------
control: add draft code for prescale command
Added Paths:
-----------
trunk/octave-forge/main/control/devel/scale/
trunk/octave-forge/main/control/devel/scale/TB01ID.dat
trunk/octave-forge/main/control/devel/scale/TB01ID.f
trunk/octave-forge/main/control/devel/scale/TB01ID.html
trunk/octave-forge/main/control/devel/scale/TB01ID.res
trunk/octave-forge/main/control/devel/scale/common.cc
trunk/octave-forge/main/control/devel/scale/makefile_scale.m
trunk/octave-forge/main/control/devel/scale/scale.m
trunk/octave-forge/main/control/devel/scale/sltb01id.cc
Property Changed:
----------------
trunk/octave-forge/main/control/devel/
trunk/octave-forge/main/control/devel/ss2tf/
Property changes on: trunk/octave-forge/main/control/devel
___________________________________________________________________
Added: svn:ignore
+ makefiles_2011_05_05_2158_no_veclib_args.zip
Property changes on: trunk/octave-forge/main/control/devel/scale
___________________________________________________________________
Added: svn:ignore
+ TB01ID.oct
sltb01id.oct
Added: trunk/octave-forge/main/control/devel/scale/TB01ID.dat
===================================================================
--- trunk/octave-forge/main/control/devel/scale/TB01ID.dat (rev 0)
+++ trunk/octave-forge/main/control/devel/scale/TB01ID.dat 2011-05-09 21:34:44 UTC (rev 8255)
@@ -0,0 +1,17 @@
+ TB01ID EXAMPLE PROGRAM DATA
+ 5 2 5 A 0.0
+ 0.0 1.0000e+000 0.0 0.0 0.0
+ -1.5800e+006 -1.2570e+003 0.0 0.0 0.0
+ 3.5410e+014 0.0 -1.4340e+003 0.0 -5.3300e+011
+ 0.0 0.0 0.0 0.0 1.0000e+000
+ 0.0 0.0 0.0 -1.8630e+004 -1.4820e+000
+ 0.0 0.0
+ 1.1030e+002 0.0
+ 0.0 0.0
+ 0.0 0.0
+ 0.0 8.3330e-003
+ 1.0000e+000 0.0 0.0 0.0 0.0
+ 0.0 0.0 1.0000e+000 0.0 0.0
+ 0.0 0.0 0.0 1.0000e+000 0.0
+ 6.6640e-001 0.0 -6.2000e-013 0.0 0.0
+ 0.0 0.0 -1.0000e-003 1.8960e+006 1.5080e+002
Added: trunk/octave-forge/main/control/devel/scale/TB01ID.f
===================================================================
--- trunk/octave-forge/main/control/devel/scale/TB01ID.f (rev 0)
+++ trunk/octave-forge/main/control/devel/scale/TB01ID.f 2011-05-09 21:34:44 UTC (rev 8255)
@@ -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-2009 NICONET e.V.
+C
+C This program is free software: you can redistribute it and/or
+C modify it under the terms of the GNU General Public License as
+C published by the Free Software Foundation, either version 2 of
+C the License, or (at your option) any later version.
+C
+C This program is distributed in the hope that it will be useful,
+C but WITHOUT ANY WARRANTY; without even the implied warranty of
+C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+C GNU General Public License for more details.
+C
+C You should have received a copy of the GNU General Public License
+C along with this program. If not, see
+C <http://www.gnu.org/licenses/>.
+C
+C PURPOSE
+C
+C To 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
Added: trunk/octave-forge/main/control/devel/scale/TB01ID.html
===================================================================
--- trunk/octave-forge/main/control/devel/scale/TB01ID.html (rev 0)
+++ trunk/octave-forge/main/control/devel/scale/TB01ID.html 2011-05-09 21:34:44 UTC (rev 8255)
@@ -0,0 +1,312 @@
+<HTML>
+<HEAD><TITLE>TB01ID - SLICOT Library Routine Documentation</TITLE>
+</HEAD>
+<BODY>
+
+<H2><A Name="TB01ID">TB01ID</A></H2>
+<H3>
+Balancing a system matrix corresponding to a triplet (A,B,C)
+</H3>
+<A HREF ="#Specification"><B>[Specification]</B></A>
+<A HREF ="#Arguments"><B>[Arguments]</B></A>
+<A HREF ="#Method"><B>[Method]</B></A>
+<A HREF ="#References"><B>[References]</B></A>
+<A HREF ="#Comments"><B>[Comments]</B></A>
+<A HREF ="#Example"><B>[Example]</B></A>
+
+<P>
+<B><FONT SIZE="+1">Purpose</FONT></B>
+<PRE>
+ To reduce the 1-norm of a system matrix
+
+ S = ( A B )
+ ( C 0 )
+
+ corresponding to the triple (A,B,C), by balancing. This involves
+ a diagonal similarity transformation inv(D)*A*D applied
+ iteratively to A to make the rows and columns of
+ -1
+ diag(D,I) * S * diag(D,I)
+
+ as close in norm as possible.
+
+ The balancing can be performed optionally on the following
+ particular system matrices
+
+ S = A, S = ( A B ) or S = ( A )
+ ( C )
+
+</PRE>
+<A name="Specification"><B><FONT SIZE="+1">Specification</FONT></B></A>
+<PRE>
+ SUBROUTINE TB01ID( JOB, N, M, P, MAXRED, A, LDA, B, LDB, C, LDC,
+ $ SCALE, INFO )
+C .. Scalar Arguments ..
+ CHARACTER JOB
+ INTEGER INFO, LDA, LDB, LDC, M, N, P
+ DOUBLE PRECISION MAXRED
+C .. Array Arguments ..
+ DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ),
+ $ SCALE( * )
+
+</PRE>
+<A name="Arguments"><B><FONT SIZE="+1">Arguments</FONT></B></A>
+<P>
+
+<B>Mode Parameters</B>
+<PRE>
+ JOB CHARACTER*1
+ Indicates which matrices are involved in balancing, as
+ follows:
+ = 'A': All matrices are involved in balancing;
+ = 'B': B and A matrices are involved in balancing;
+ = 'C': C and A matrices are involved in balancing;
+ = 'N': B and C matrices are not involved in balancing.
+
+</PRE>
+<B>Input/Output Parameters</B>
+<PRE>
+ N (input) INTEGER
+ The order of the matrix A, the number of rows of matrix B
+ and the number of columns of matrix C.
+ N represents the dimension of the state vector. N >= 0.
+
+ M (input) INTEGER.
+ The number of columns of matrix B.
+ M represents the dimension of input vector. M >= 0.
+
+ P (input) INTEGER.
+ The number of rows of matrix C.
+ P represents the dimension of output vector. P >= 0.
+
+ MAXRED (input/output) DOUBLE PRECISION
+ On entry, the maximum allowed reduction in the 1-norm of
+ S (in an iteration) if zero rows or columns are
+ encountered.
+ If MAXRED > 0.0, MAXRED must be larger than one (to enable
+ the norm reduction).
+ If MAXRED <= 0.0, then the value 10.0 for MAXRED is
+ used.
+ On exit, if the 1-norm of the given matrix S is non-zero,
+ the ratio between the 1-norm of the given matrix and the
+ 1-norm of the balanced matrix.
+
+ A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
+ On entry, the leading N-by-N part of this array must
+ contain the system state matrix A.
+ On exit, the leading N-by-N part of this array contains
+ the balanced matrix inv(D)*A*D.
+
+ LDA INTEGER
+ The leading dimension of the array A. LDA >= max(1,N).
+
+ B (input/output) DOUBLE PRECISION array, dimension (LDB,M)
+ On entry, if M > 0, the leading N-by-M part of this array
+ must contain the system input matrix B.
+ On exit, if M > 0, the leading N-by-M part of this array
+ contains the balanced matrix inv(D)*B.
+ The array B is not referenced if M = 0.
+
+ LDB INTEGER
+ The leading dimension of the array B.
+ LDB >= MAX(1,N) if M > 0.
+ LDB >= 1 if M = 0.
+
+ C (input/output) DOUBLE PRECISION array, dimension (LDC,N)
+ On entry, if P > 0, the leading P-by-N part of this array
+ must contain the system output matrix C.
+ On exit, if P > 0, the leading P-by-N part of this array
+ contains the balanced matrix C*D.
+ The array C is not referenced if P = 0.
+
+ LDC INTEGER
+ The leading dimension of the array C. LDC >= MAX(1,P).
+
+ SCALE (output) DOUBLE PRECISION array, dimension (N)
+ The scaling factors applied to S. If D(j) is the scaling
+ factor applied to row and column j, then SCALE(j) = D(j),
+ for j = 1,...,N.
+
+</PRE>
+<B>Error Indicator</B>
+<PRE>
+ INFO INTEGER
+ = 0: successful exit.
+ < 0: if INFO = -i, the i-th argument had an illegal
+ value.
+
+</PRE>
+<A name="Method"><B><FONT SIZE="+1">Method</FONT></B></A>
+<PRE>
+ Balancing consists of applying a diagonal similarity
+ transformation
+ -1
+ diag(D,I) * S * diag(D,I)
+
+ to make the 1-norms of each row of the first N rows of S and its
+ corresponding column nearly equal.
+
+ Information about the diagonal matrix D is returned in the vector
+ SCALE.
+
+</PRE>
+<A name="References"><B><FONT SIZE="+1">References</FONT></B></A>
+<PRE>
+ [1] Anderson, E., Bai, Z., Bischof, C., Demmel, J., Dongarra, J.,
+ Du Croz, J., Greenbaum, A., Hammarling, S., McKenney, A.,
+ Ostrouchov, S., and Sorensen, D.
+ LAPACK Users' Guide: Second Edition.
+ SIAM, Philadelphia, 1995.
+
+</PRE>
+<A name="Numerical Aspects"><B><FONT SIZE="+1">Numerical Aspects</FONT></B></A>
+<PRE>
+ None.
+
+</PRE>
+
+<A name="Comments"><B><FONT SIZE="+1">Further Comments</FONT></B></A>
+<PRE>
+ None
+</PRE>
+
+<A name="Example"><B><FONT SIZE="+1">Example</FONT></B></A>
+<P>
+<B>Program Text</B>
+<PRE>
+* TB01ID EXAMPLE PROGRAM TEXT.
+* Copyright (c) 2002-2009 NICONET e.V.
+*
+* .. Parameters ..
+ INTEGER NIN, NOUT
+ PARAMETER ( NIN = 5, NOUT = 6 )
+ INTEGER NMAX, MMAX, PMAX
+ PARAMETER ( NMAX = 20, MMAX = 20, PMAX = 20 )
+ INTEGER LDA, LDB, LDC
+ PARAMETER ( LDA = NMAX, LDB = NMAX, LDC = PMAX )
+* .. Local Scalars ..
+ CHARACTER*1 JOB
+ INTEGER I, INFO, J, M, N, P
+ DOUBLE PRECISION MAXRED
+* .. Local Arrays ..
+ DOUBLE PRECISION A(LDA,NMAX), B(LDB,MMAX), C(LDC,NMAX),
+ $ SCALE(NMAX)
+* .. External Subroutines ..
+ EXTERNAL TB01ID, UD01MD
+* .. Executable Statements ..
+*
+ WRITE ( NOUT, FMT = 99999 )
+* Skip the heading in the data file and read the data.
+ READ ( NIN, FMT = '()' )
+ READ ( NIN, FMT = * ) N, M, P, JOB, MAXRED
+ IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
+ WRITE ( NOUT, FMT = 99993 ) N
+ ELSE
+ READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,N )
+ IF ( M.LT.0 .OR. M.GT.MMAX ) THEN
+ WRITE ( NOUT, FMT = 99992 ) M
+ ELSE
+ READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,M ), I = 1,N )
+ IF ( P.LT.0 .OR. P.GT.MMAX ) THEN
+ WRITE ( NOUT, FMT = 99991 ) P
+ ELSE
+ READ ( NIN, FMT = * ) ( ( C(I,J), J = 1,N ), I = 1,P )
+* Balance system matrix S.
+ CALL TB01ID( JOB, N, M, P, MAXRED, A, LDA, B, LDB, C,
+ $ LDC, SCALE, INFO )
+*
+ IF ( INFO.NE.0 ) THEN
+ WRITE ( NOUT, FMT = 99998 ) INFO
+ ELSE
+ CALL UD01MD( N, N, 5, NOUT, A, LDA,
+ $ 'The balanced matrix A', INFO )
+ IF ( M.GT.0 )
+ $ CALL UD01MD( N, M, 5, NOUT, B, LDB,
+ $ 'The balanced matrix B', INFO )
+ IF ( P.GT.0 )
+ $ CALL UD01MD( P, N, 5, NOUT, C, LDC,
+ $ 'The balanced matrix C', INFO )
+ CALL UD01MD( 1, N, 5, NOUT, SCALE, 1,
+ $ 'The scaling vector SCALE', INFO )
+ WRITE ( NOUT, FMT = 99994 ) MAXRED
+ END IF
+ END IF
+ END IF
+ END IF
+ STOP
+*
+99999 FORMAT (' TB01ID EXAMPLE PROGRAM RESULTS',/1X)
+99998 FORMAT (' INFO on exit from TB01ID = ',I2)
+99994 FORMAT (/' MAXRED is ',E13.4)
+99993 FORMAT (/' N is out of range.',/' N = ',I5)
+99992 FORMAT (/' M is out of range.',/' M = ',I5)
+99991 FORMAT (/' P is out of range.',/' P = ',I5)
+ END
+</PRE>
+<B>Program Data</B>
+<PRE>
+ TB01ID EXAMPLE PROGRAM DATA
+ 5 2 5 A 0.0
+ 0.0 1.0000e+000 0.0 0.0 0.0
+ -1.5800e+006 -1.2570e+003 0.0 0.0 0.0
+ 3.5410e+014 0.0 -1.4340e+003 0.0 -5.3300e+011
+ 0.0 0.0 0.0 0.0 1.0000e+000
+ 0.0 0.0 0.0 -1.8630e+004 -1.4820e+000
+ 0.0 0.0
+ 1.1030e+002 0.0
+ 0.0 0.0
+ 0.0 0.0
+ 0.0 8.3330e-003
+ 1.0000e+000 0.0 0.0 0.0 0.0
+ 0.0 0.0 1.0000e+000 0.0 0.0
+ 0.0 0.0 0.0 1.0000e+000 0.0
+ 6.6640e-001 0.0 -6.2000e-013 0.0 0.0
+ 0.0 0.0 -1.0000e-003 1.8960e+006 1.5080e+002
+</PRE>
+<B>Program Results</B>
+<PRE>
+ TB01ID EXAMPLE PROGRAM RESULTS
+
+ The balanced matrix A ( 5X 5)
+
+ 1 2 3 4 5
+ 1 0.0000000D+00 0.1000000D+05 0.0000000D+00 0.0000000D+00 0.0000000D+00
+ 2 -0.1580000D+03 -0.1257000D+04 0.0000000D+00 0.0000000D+00 0.0000000D+00
+ 3 0.3541000D+05 0.0000000D+00 -0.1434000D+04 0.0000000D+00 -0.5330000D+03
+ 4 0.0000000D+00 0.0000000D+00 0.0000000D+00 0.0000000D+00 0.1000000D+03
+ 5 0.0000000D+00 0.0000000D+00 0.0000000D+00 -0.1863000D+03 -0.1482000D+01
+
+ The balanced matrix B ( 5X 2)
+
+ 1 2
+ 1 0.0000000D+00 0.0000000D+00
+ 2 0.1103000D+04 0.0000000D+00
+ 3 0.0000000D+00 0.0000000D+00
+ 4 0.0000000D+00 0.0000000D+00
+ 5 0.0000000D+00 0.8333000D+02
+
+ The balanced matrix C ( 5X 5)
+
+ 1 2 3 4 5
+ 1 0.1000000D-04 0.0000000D+00 0.0000000D+00 0.0000000D+00 0.0000000D+00
+ 2 0.0000000D+00 0.0000000D+00 0.1000000D+06 0.0000000D+00 0.0000000D+00
+ 3 0.0000000D+00 0.0000000D+00 0.0000000D+00 0.1000000D-05 0.0000000D+00
+ 4 0.6664000D-05 0.0000000D+00 -0.6200000D-07 0.0000000D+00 0.0000000D+00
+ 5 0.0000000D+00 0.0000000D+00 -0.1000000D+03 0.1896000D+01 0.1508000D-01
+
+ The scaling vector SCALE ( 1X 5)
+
+ 1 2 3 4 5
+ 1 0.1000000D-04 0.1000000D+00 0.1000000D+06 0.1000000D-05 0.1000000D-03
+
+
+ MAXRED is 0.3488E+10
+</PRE>
+
+<HR>
+<p>
+Click <a href=../../SLICOT/arc/T/TB/TB01/TB01ID.tar.gz><B>here</B></a> to get a compressed (gzip) tar file containing the source code
+of the routine, the example program, data, documentation, and related files.
+</p>
+<A HREF=../libindex.html><B>Return to index</B></A></BODY>
+</HTML>
Added: trunk/octave-forge/main/control/devel/scale/TB01ID.res
===================================================================
--- trunk/octave-forge/main/control/devel/scale/TB01ID.res (rev 0)
+++ trunk/octave-forge/main/control/devel/scale/TB01ID.res 2011-05-09 21:34:44 UTC (rev 8255)
@@ -0,0 +1,36 @@
+ TB01ID EXAMPLE PROGRAM RESULTS
+
+ The balanced matrix A ( 5X 5)
+
+ 1 2 3 4 5
+ 1 0.0000000D+00 0.1000000D+05 0.0000000D+00 0.0000000D+00 0.0000000D+00
+ 2 -0.1580000D+03 -0.1257000D+04 0.0000000D+00 0.0000000D+00 0.0000000D+00
+ 3 0.3541000D+05 0.0000000D+00 -0.1434000D+04 0.0000000D+00 -0.5330000D+03
+ 4 0.0000000D+00 0.0000000D+00 0.0000000D+00 0.0000000D+00 0.1000000D+03
+ 5 0.0000000D+00 0.0000000D+00 0.0000000D+00 -0.1863000D+03 -0.1482000D+01
+
+ The balanced matrix B ( 5X 2)
+
+ 1 2
+ 1 0.0000000D+00 0.0000000D+00
+ 2 0.1103000D+04 0.0000000D+00
+ 3 0.0000000D+00 0.0000000D+00
+ 4 0.0000000D+00 0.0000000D+00
+ 5 0.0000000D+00 0.8333000D+02
+
+ The balanced matrix C ( 5X 5)
+
+ 1 2 3 4 5
+ 1 0.1000000D-04 0.0000000D+00 0.0000000D+00 0.0000000D+00 0.0000000D+00
+ 2 0.0000000D+00 0.0000000D+00 0.1000000D+06 0.0000000D+00 0.0000000D+00
+ 3 0.0000000D+00 0.0000000D+00 0.0000000D+00 0.1000000D-05 0.0000000D+00
+ 4 0.6664000D-05 0.0000000D+00 -0.6200000D-07 0.0000000D+00 0.0000000D+00
+ 5 0.0000000D+00 0.0000000D+00 -0.1000000D+03 0.1896000D+01 0.1508000D-01
+
+ The scaling vector SCALE ( 1X 5)
+
+ 1 2 3 4 5
+ 1 0.1000000D-04 0.1000000D+00 0.1000000D+06 0.1000000D-05 0.1000000D-03
+
+
+ MAXRED is 0.3488E+10
Added: trunk/octave-forge/main/control/devel/scale/common.cc
===================================================================
--- trunk/octave-forge/main/control/devel/scale/common.cc (rev 0)
+++ trunk/octave-forge/main/control/devel/scale/common.cc 2011-05-09 21:34:44 UTC (rev 8255)
@@ -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/scale/makefile_scale.m
===================================================================
--- trunk/octave-forge/main/control/devel/scale/makefile_scale.m (rev 0)
+++ trunk/octave-forge/main/control/devel/scale/makefile_scale.m 2011-05-09 21:34:44 UTC (rev 8255)
@@ -0,0 +1,5 @@
+## scaling of state-space models
+mkoctfile "-Wl,-framework" "-Wl,vecLib" \
+ sltb01id.cc \
+ TB01ID.f
+
Added: trunk/octave-forge/main/control/devel/scale/scale.m
===================================================================
--- trunk/octave-forge/main/control/devel/scale/scale.m (rev 0)
+++ trunk/octave-forge/main/control/devel/scale/scale.m 2011-05-09 21:34:44 UTC (rev 8255)
@@ -0,0 +1,68 @@
+
+a = [
+ 0.0 1.0000e+000 0.0 0.0 0.0
+ -1.5800e+006 -1.2570e+003 0.0 0.0 0.0
+ 3.5410e+014 0.0 -1.4340e+003 0.0 -5.3300e+011
+ 0.0 0.0 0.0 0.0 1.0000e+000
+ 0.0 0.0 0.0 -1.8630e+004 -1.4820e+000
+];
+
+b = [
+ 0.0 0.0
+ 1.1030e+002 0.0
+ 0.0 0.0
+ 0.0 0.0
+ 0.0 8.3330e-003
+];
+
+c = [
+ 1.0000e+000 0.0 0.0 0.0 0.0
+ 0.0 0.0 1.0000e+000 0.0 0.0
+ 0.0 0.0 0.0 1.0000e+000 0.0
+ 6.6640e-001 0.0 -6.2000e-013 0.0 0.0
+ 0.0 0.0 -1.0000e-003 1.8960e+006 1.5080e+002
+];
+
+[as, bs, cs, maxred, scale] = sltb01id (a, b, c, 0);
+
+as, bs, cs, maxred, scale
+
+
+%{
+ TB01ID EXAMPLE PROGRAM RESULTS
+
+ The balanced matrix A ( 5X 5)
+
+ 1 2 3 4 5
+ 1 0.0000000D+00 0.1000000D+05 0.0000000D+00 0.0000000D+00 0.0000000D+00
+ 2 -0.1580000D+03 -0.1257000D+04 0.0000000D+00 0.0000000D+00 0.0000000D+00
+ 3 0.3541000D+05 0.0000000D+00 -0.1434000D+04 0.0000000D+00 -0.5330000D+03
+ 4 0.0000000D+00 0.0000000D+00 0.0000000D+00 0.0000000D+00 0.1000000D+03
+ 5 0.0000000D+00 0.0000000D+00 0.0000000D+00 -0.1863000D+03 -0.1482000D+01
+
+ The balanced matrix B ( 5X 2)
+
+ 1 2
+ 1 0.0000000D+00 0.0000000D+00
+ 2 0.1103000D+04 0.0000000D+00
+ 3 0.0000000D+00 0.0000000D+00
+ 4 0.0000000D+00 0.0000000D+00
+ 5 0.0000000D+00 0.8333000D+02
+
+ The balanced matrix C ( 5X 5)
+
+ 1 2 3 4 5
+ 1 0.1000000D-04 0.0000000D+00 0.0000000D+00 0.0000000D+00 0.0000000D+00
+ 2 0.0000000D+00 0.0000000D+00 0.1000000D+06 0.0000000D+00 0.0000000D+00
+ 3 0.0000000D+00 0.0000000D+00 0.0000000D+00 0.1000000D-05 0.0000000D+00
+ 4 0.6664000D-05 0.0000000D+00 -0.6200000D-07 0.0000000D+00 0.0000000D+00
+ 5 0.0000000D+00 0.0000000D+00 -0.1000000D+03 0.1896000D+01 0.1508000D-01
+
+ The scaling vector SCALE ( 1X 5)
+
+ 1 2 3 4 5
+ 1 0.1000000D-04 0.1000000D+00 0.1000000D+06 0.1000000D-05 0.1000000D-03
+
+
+ MAXRED is 0.3488E+10
+%}
\ No newline at end of file
Added: trunk/octave-forge/main/control/devel/scale/sltb01id.cc
===================================================================
--- trunk/octave-forge/main/control/devel/scale/sltb01id.cc (rev 0)
+++ trunk/octave-forge/main/control/devel/scale/sltb01id.cc 2011-05-09 21:34:44 UTC (rev 8255)
@@ -0,0 +1,114 @@
+/*
+
+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/>.
+
+Balance state-space model.
+Uses SLICOT TB01ID by courtesy of NICONET e.V.
+<http://www.slicot.org>
+
+Author: Lukas Reichlin <luk...@gm...>
+Created: May 2011
+Version: 0.1
+
+*/
+
+#include <octave/oct.h>
+#include <f77-fcn.h>
+#include "common.cc"
+
+extern "C"
+{
+ int F77_FUNC (tb01id, TB01ID)
+ (char& JOB,
+ int& N, int& M, int& P,
+ double& MAXRED,
+ double* A, int& LDA,
+ double* B, int& LDB,
+ double* C, int& LDC,
+ double* SCALE,
+ int& INFO);
+}
+
+DEFUN_DLD (sltb01id, args, nargout,
+ "-*- texinfo -*-\n\
+Slicot TB01ID Release 5.0\n\
+No argument checking.\n\
+For internal use only.")
+{
+ int nargin = args.length ();
+ octave_value_list retval;
+
+ if (nargin != 4)
+ {
+ print_usage ();
+ }
+ else
+ {
+ // arguments in
+ char job = 'A';
+
+ Matrix a = args(0).matrix_value ();
+ Matrix b = args(1).matrix_value ();
+ Matrix c = args(2).matrix_value ();
+ double maxred = args(3).double_value ();
+
+ int n = a.rows (); // n: number of states
+ int m = b.columns (); // m: number of inputs
+ int p = c.rows (); // p: number of outputs
+
+ int lda = max (1, n);
+ int ldb = max (1, n);
+ int ldc = max (1, p);
+
+
+ // arguments out
+ ColumnVector scale (n);
+
+
+ // error indicators
+ int info = 0;
+
+
+ // SLICOT routine TB01ID
+ F77_XFCN (tb01id, TB01ID,
+ (job,
+ n, m, p,
+ maxred,
+ a.fortran_vec (), lda,
+ b.fortran_vec (), ldb,
+ c.fortran_vec (), ldc,
+ scale.fortran_vec (),
+ info));
+
+ if (f77_exception_encountered)
+ error ("ss: prescale: sltb01id: exception in SLICOT subroutine TB01ID");
+
+ if (info != 0)
+ error ("ss: prescale: sltb01id: TB01ID returned info = %d", info);
+
+
+ // return values
+ retval(0) = a;
+ retval(1) = b;
+ retval(2) = c;
+ retval(3) = octave_value (maxred);
+ retval(4) = scale;
+ }
+
+ return retval;
+}
Property changes on: trunk/octave-forge/main/control/devel/ss2tf
___________________________________________________________________
Added: svn:ignore
+ sltb04bd.oct
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
