88
89
90
91
92
93
94 INTEGER KNT, LMAX, NIN
95 REAL RMAX
96
97
98 INTEGER NINFO( 3 )
99
100
101
102
103
104 REAL ZERO, ONE
105 parameter( zero = 0.0e0, one = 1.0e0 )
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 REAL EPS, RES
113
114
115 REAL Q( LDT, LDT ), RESULT( 2 ), T1( LDT, LDT ),
116 $ T2( LDT, LDT ), TMP( LDT, LDT ), WORK( LWORK )
117
118
119 REAL SLAMCH
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 slacpy(
'F', n, n, tmp, ldt, t1, ldt )
149 CALL slacpy(
'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 slaset(
'Full', n, n, zero, one, q, ldt )
161 CALL strexc(
'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 slaset(
'Full', n, n, zero, one, q, ldt )
174 CALL strexc(
'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 shst01( 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 slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
real function slamch(cmach)
SLAMCH
subroutine slaset(uplo, m, n, alpha, beta, a, lda)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
subroutine strexc(compq, n, t, ldt, q, ldq, ifst, ilst, work, info)
STREXC
subroutine shst01(n, ilo, ihi, a, lda, h, ldh, q, ldq, work, lwork, result)
SHST01