85
86
87
88
89
90
91 INTEGER KNT, LMAX, NIN, NINFO
92 DOUBLE PRECISION RMAX
93
94
95
96
97
98 DOUBLE PRECISION ZERO, ONE
99 parameter( zero = 0.0d+0, one = 1.0d+0 )
100 COMPLEX*16 CZERO, CONE
101 parameter( czero = ( 0.0d+0, 0.0d+0 ),
102 $ cone = ( 1.0d+0, 0.0d+0 ) )
103 INTEGER LDT, LWORK
104 parameter( ldt = 10, lwork = 2*ldt*ldt )
105
106
107 INTEGER I, IFST, ILST, INFO1, INFO2, J, N
108 DOUBLE PRECISION EPS, RES
109 COMPLEX*16 CTEMP
110
111
112 DOUBLE PRECISION RESULT( 2 ), RWORK( LDT )
113 COMPLEX*16 DIAG( LDT ), Q( LDT, LDT ), T1( LDT, LDT ),
114 $ T2( LDT, LDT ), TMP( LDT, LDT ), WORK( LWORK )
115
116
117 DOUBLE PRECISION DLAMCH
119
120
122
123
124
126 rmax = zero
127 lmax = 0
128 knt = 0
129 ninfo = 0
130
131
132
133 10 CONTINUE
134 READ( nin, fmt = * )n, ifst, ilst
135 IF( n.EQ.0 )
136 $ RETURN
137 knt = knt + 1
138 DO 20 i = 1, n
139 READ( nin, fmt = * )( tmp( i, j ), j = 1, n )
140 20 CONTINUE
141 CALL zlacpy(
'F', n, n, tmp, ldt, t1, ldt )
142 CALL zlacpy(
'F', n, n, tmp, ldt, t2, ldt )
143 res = zero
144
145
146
147 CALL zlaset(
'Full', n, n, czero, cone, q, ldt )
148 CALL ztrexc(
'N', n, t1, ldt, q, ldt, ifst, ilst, info1 )
149 DO 40 i = 1, n
150 DO 30 j = 1, n
151 IF( i.EQ.j .AND. q( i, j ).NE.cone )
152 $ res = res + one / eps
153 IF( i.NE.j .AND. q( i, j ).NE.czero )
154 $ res = res + one / eps
155 30 CONTINUE
156 40 CONTINUE
157
158
159
160 CALL zlaset(
'Full', n, n, czero, cone, q, ldt )
161 CALL ztrexc(
'V', n, t2, ldt, q, ldt, ifst, ilst, info2 )
162
163
164
165 DO 60 i = 1, n
166 DO 50 j = 1, n
167 IF( t1( i, j ).NE.t2( i, j ) )
168 $ res = res + one / eps
169 50 CONTINUE
170 60 CONTINUE
171 IF( info1.NE.0 .OR. info2.NE.0 )
172 $ ninfo = ninfo + 1
173 IF( info1.NE.info2 )
174 $ res = res + one / eps
175
176
177
178 CALL zcopy( n, tmp, ldt+1, diag, 1 )
179 IF( ifst.LT.ilst ) THEN
180 DO 70 i = ifst + 1, ilst
181 ctemp = diag( i )
182 diag( i ) = diag( i-1 )
183 diag( i-1 ) = ctemp
184 70 CONTINUE
185 ELSE IF( ifst.GT.ilst ) THEN
186 DO 80 i = ifst - 1, ilst, -1
187 ctemp = diag( i+1 )
188 diag( i+1 ) = diag( i )
189 diag( i ) = ctemp
190 80 CONTINUE
191 END IF
192 DO 90 i = 1, n
193 IF( t2( i, i ).NE.diag( i ) )
194 $ res = res + one / eps
195 90 CONTINUE
196
197
198
199 CALL zhst01( n, 1, n, tmp, ldt, t2, ldt, q, ldt, work, lwork,
200 $ rwork, result )
201 res = res + result( 1 ) + result( 2 )
202
203
204
205 DO 110 j = 1, n - 1
206 DO 100 i = j + 1, n
207 IF( t2( i, j ).NE.czero )
208 $ res = res + one / eps
209 100 CONTINUE
210 110 CONTINUE
211 IF( res.GT.rmax ) THEN
212 rmax = res
213 lmax = knt
214 END IF
215 GO TO 10
216
217
218
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
double precision function dlamch(cmach)
DLAMCH
subroutine zlaset(uplo, m, n, alpha, beta, a, lda)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
subroutine ztrexc(compq, n, t, ldt, q, ldq, ifst, ilst, info)
ZTREXC
subroutine zhst01(n, ilo, ihi, a, lda, h, ldh, q, ldq, work, lwork, rwork, result)
ZHST01