1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
|
SUBROUTINE IB01OD( CTRL, NOBR, L, SV, N, TOL, IWARN, INFO )
C
C RELEASE 4.0, WGS COPYRIGHT 2000.
C
C PURPOSE
C
C To estimate the system order, based on the singular values of the
C relevant part of the triangular factor of the concatenated block
C Hankel matrices.
C
C ARGUMENTS
C
C Mode Parameters
C
C CTRL CHARACTER*1
C Specifies whether or not the user's confirmation of the
C system order estimate is desired, as follows:
C = 'C': user's confirmation;
C = 'N': no confirmation.
C If CTRL = 'C', a reverse communication routine, IB01OY,
C is called, and, after inspecting the singular values and
C system order estimate, n, the user may accept n or set
C a new value.
C IB01OY is not called by the routine if CTRL = 'N'.
C
C Input/Output Parameters
C
C NOBR (input) INTEGER
C The number of block rows, s, in the processed input and
C output block Hankel matrices. NOBR > 0.
C
C L (input) INTEGER
C The number of system outputs. L > 0.
C
C SV (input) DOUBLE PRECISION array, dimension ( L*NOBR )
C The singular values of the relevant part of the triangular
C factor from the QR factorization of the concatenated block
C Hankel matrices.
C
C N (output) INTEGER
C The estimated order of the system.
C
C Tolerances
C
C TOL DOUBLE PRECISION
C Absolute tolerance used for determining an estimate of
C the system order. If TOL >= 0, the estimate is
C indicated by the index of the last singular value greater
C than or equal to TOL. (Singular values less than TOL
C are considered as zero.) When TOL = 0, an internally
C computed default value, TOL = NOBR*EPS*SV(1), is used,
C where SV(1) is the maximal singular value, and EPS is
C the relative machine precision (see LAPACK Library routine
C DLAMCH). When TOL < 0, the estimate is indicated by the
C index of the singular value that has the largest
C logarithmic gap to its successor.
C
C Warning Indicator
C
C IWARN INTEGER
C = 0: no warning;
C = 3: all singular values were exactly zero, hence N = 0.
C (Both input and output were identically zero.)
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 singular values are compared to the given, or default TOL, and
C the estimated order n is returned, possibly after user's
C confirmation.
C
C CONTRIBUTOR
C
C V. Sima, Research Institute for Informatics, Bucharest, Aug. 1999.
C
C REVISIONS
C
C August 2000.
C
C KEYWORDS
C
C Identification methods, multivariable systems, singular value
C decomposition.
C
C ******************************************************************
C
C .. Parameters ..
DOUBLE PRECISION ZERO
PARAMETER ( ZERO = 0.0D0 )
C .. Scalar Arguments ..
DOUBLE PRECISION TOL
INTEGER INFO, IWARN, L, N, NOBR
CHARACTER CTRL
C .. Array Arguments ..
DOUBLE PRECISION SV(*)
C .. Local Scalars ..
DOUBLE PRECISION GAP, RNRM, TOLL
INTEGER I, IERR, LNOBR
LOGICAL CONTRL
C .. External Functions ..
LOGICAL LSAME
DOUBLE PRECISION DLAMCH
EXTERNAL DLAMCH, LSAME
C .. External Subroutines ..
EXTERNAL IB01OY, XERBLA
C .. Intrinsic Functions ..
INTRINSIC DBLE, LOG10
C ..
C .. Executable Statements ..
C
C Check the scalar input parameters.
C
CONTRL = LSAME( CTRL, 'C' )
LNOBR = L*NOBR
IWARN = 0
INFO = 0
IF( .NOT.( CONTRL .OR. LSAME( CTRL, 'N' ) ) ) THEN
INFO = -1
ELSE IF( NOBR.LE.0 ) THEN
INFO = -2
ELSE IF( L.LE.0 ) THEN
INFO = -3
END IF
C
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'IB01OD', -INFO )
RETURN
END IF
C
C Set TOL if necessay.
C
TOLL = TOL
IF ( TOLL.EQ.ZERO)
$ TOLL = DLAMCH( 'Precision' )*SV(1)*DBLE( NOBR )
C
C Obtain the system order.
C
N = 0
IF ( SV(1).NE.ZERO ) THEN
N = NOBR
IF ( TOLL.GE.ZERO) THEN
C
C Estimate n based on the tolerance TOLL.
C
DO 10 I = 1, NOBR - 1
IF ( SV(I+1).LT.TOLL ) THEN
N = I
GO TO 30
END IF
10 CONTINUE
ELSE
C
C Estimate n based on the largest logarithmic gap between
C two consecutive singular values.
C
GAP = ZERO
DO 20 I = 1, NOBR - 1
RNRM = SV(I+1)
IF ( RNRM.NE.ZERO ) THEN
RNRM = LOG10( SV(I) ) - LOG10( RNRM )
IF ( RNRM.GT.GAP ) THEN
GAP = RNRM
N = I
END IF
ELSE
IF ( GAP.EQ.ZERO )
$ N = I
GO TO 30
END IF
20 CONTINUE
END IF
END IF
C
30 CONTINUE
IF ( N.EQ.0 ) THEN
C
C Return with N = 0 if all singular values are zero.
C
IWARN = 3
RETURN
END IF
C
IF ( CONTRL ) THEN
C
C Ask confirmation of the system order.
C
CALL IB01OY( LNOBR, NOBR-1, N, SV, IERR )
END IF
RETURN
C
C *** Last line of IB01OD ***
END
|
