88
89
90
91
92
93
94 INTEGER KNT, LMAX, NIN
95 DOUBLE PRECISION RMAX
96
97
98 INTEGER NINFO( 3 )
99
100
101
102
103
104 DOUBLE PRECISION ZERO, ONE
105 parameter( zero = 0.0d0, one = 1.0d0 )
106 INTEGER LDT, LWORK
107 parameter( ldt = 10, lwork = 2*ldt*ldt )
108
109
110 INTEGER I, IFST, IFST1, IFST2, IFSTSV, ILST, ILST1,
111 $ ILST2, ILSTSV, INFO1, INFO2, J, LOC, N
112 DOUBLE PRECISION EPS, RES
113
114
115 DOUBLE PRECISION Q( LDT, LDT ), RESULT( 2 ), T1( LDT, LDT ),
116 $ T2( LDT, LDT ), TMP( LDT, LDT ), WORK( LWORK )
117
118
119 DOUBLE PRECISION DLAMCH
121
122
124
125
126 INTRINSIC abs, sign
127
128
129
131 rmax = zero
132 lmax = 0
133 knt = 0
134 ninfo( 1 ) = 0
135 ninfo( 2 ) = 0
136 ninfo( 3 ) = 0
137
138
139
140 10 CONTINUE
141 READ( nin, fmt = * )n, ifst, ilst
142 IF( n.EQ.0 )
143 $ RETURN
144 knt = knt + 1
145 DO 20 i = 1, n
146 READ( nin, fmt = * )( tmp( i, j ), j = 1, n )
147 20 CONTINUE
148 CALL dlacpy(
'F', n, n, tmp, ldt, t1, ldt )
149 CALL dlacpy(
'F', n, n, tmp, ldt, t2, ldt )
150 ifstsv = ifst
151 ilstsv = ilst
152 ifst1 = ifst
153 ilst1 = ilst
154 ifst2 = ifst
155 ilst2 = ilst
156 res = zero
157
158
159
160 CALL dlaset(
'Full', n, n, zero, one, q, ldt )
161 CALL dtrexc(
'N', n, t1, ldt, q, ldt, ifst1, ilst1, work, info1 )
162 DO 40 i = 1, n
163 DO 30 j = 1, n
164 IF( i.EQ.j .AND. q( i, j ).NE.one )
165 $ res = res + one / eps
166 IF( i.NE.j .AND. q( i, j ).NE.zero )
167 $ res = res + one / eps
168 30 CONTINUE
169 40 CONTINUE
170
171
172
173 CALL dlaset(
'Full', n, n, zero, one, q, ldt )
174 CALL dtrexc(
'V', n, t2, ldt, q, ldt, ifst2, ilst2, work, info2 )
175
176
177
178 DO 60 i = 1, n
179 DO 50 j = 1, n
180 IF( t1( i, j ).NE.t2( i, j ) )
181 $ res = res + one / eps
182 50 CONTINUE
183 60 CONTINUE
184 IF( ifst1.NE.ifst2 )
185 $ res = res + one / eps
186 IF( ilst1.NE.ilst2 )
187 $ res = res + one / eps
188 IF( info1.NE.info2 )
189 $ res = res + one / eps
190
191
192
193 IF( info2.NE.0 ) THEN
194 ninfo( info2 ) = ninfo( info2 ) + 1
195 ELSE
196 IF( abs( ifst2-ifstsv ).GT.1 )
197 $ res = res + one / eps
198 IF( abs( ilst2-ilstsv ).GT.1 )
199 $ res = res + one / eps
200 END IF
201
202
203
204 CALL dhst01( n, 1, n, tmp, ldt, t2, ldt, q, ldt, work, lwork,
205 $ result )
206 res = res + result( 1 ) + result( 2 )
207
208
209
210 loc = 1
211 70 CONTINUE
212 IF( t2( loc+1, loc ).NE.zero ) THEN
213
214
215
216 IF( t2( loc, loc+1 ).EQ.zero .OR. t2( loc, loc ).NE.
217 $ t2( loc+1, loc+1 ) .OR. sign( one, t2( loc, loc+1 ) ).EQ.
218 $ sign( one, t2( loc+1, loc ) ) )res = res + one / eps
219 DO 80 i = loc + 2, n
220 IF( t2( i, loc ).NE.zero )
221 $ res = res + one / res
222 IF( t2( i, loc+1 ).NE.zero )
223 $ res = res + one / res
224 80 CONTINUE
225 loc = loc + 2
226 ELSE
227
228
229
230 DO 90 i = loc + 1, n
231 IF( t2( i, loc ).NE.zero )
232 $ res = res + one / res
233 90 CONTINUE
234 loc = loc + 1
235 END IF
236 IF( loc.LT.n )
237 $ GO TO 70
238 IF( res.GT.rmax ) THEN
239 rmax = res
240 lmax = knt
241 END IF
242 GO TO 10
243
244
245
subroutine dhst01(n, ilo, ihi, a, lda, h, ldh, q, ldq, work, lwork, result)
DHST01
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
double precision function dlamch(cmach)
DLAMCH
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
subroutine dtrexc(compq, n, t, ldt, q, ldq, ifst, ilst, work, info)
DTREXC