90
91
92
93
94
95
96 INTEGER KNT, NIN
97
98
99 INTEGER LMAX( 3 ), NINFO( 3 )
100 DOUBLE PRECISION RMAX( 3 )
101
102
103
104
105
106 DOUBLE PRECISION ZERO, ONE, TWO
107 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
108 DOUBLE PRECISION EPSIN
109 parameter( epsin = 5.9605d-8 )
110 INTEGER LDT, LWORK
111 parameter( ldt = 20, lwork = 2*ldt*( 10+ldt ) )
112
113
114 INTEGER I, ICMP, INFO, ISCL, ISRT, J, KMIN, M, N
115 DOUBLE PRECISION BIGNUM, EPS, SMLNUM, TNRM, TOL, TOLIN, V,
116 $ VCMIN, VMAX, VMIN, VMUL
117
118
119 LOGICAL SELECT( LDT )
120 INTEGER LCMP( 3 )
121 DOUBLE PRECISION DUM( 1 ), RWORK( 2*LDT ), S( LDT ), SEP( LDT ),
122 $ SEPIN( LDT ), SEPTMP( LDT ), SIN( LDT ),
123 $ STMP( LDT ), VAL( 3 ), WIIN( LDT ),
124 $ WRIN( LDT ), WSRT( LDT )
125 COMPLEX*16 CDUM( 1 ), LE( LDT, LDT ), RE( LDT, LDT ),
126 $ T( LDT, LDT ), TMP( LDT, LDT ), W( LDT ),
127 $ WORK( LWORK ), WTMP( LDT )
128
129
130 DOUBLE PRECISION DLAMCH, ZLANGE
132
133
136
137
138 INTRINSIC dble, dimag, max, sqrt
139
140
141
143 smlnum =
dlamch(
'S' ) / eps
144 bignum = one / smlnum
145
146
147
148 eps = max( eps, epsin )
149 rmax( 1 ) = zero
150 rmax( 2 ) = zero
151 rmax( 3 ) = zero
152 lmax( 1 ) = 0
153 lmax( 2 ) = 0
154 lmax( 3 ) = 0
155 knt = 0
156 ninfo( 1 ) = 0
157 ninfo( 2 ) = 0
158 ninfo( 3 ) = 0
159 val( 1 ) = sqrt( smlnum )
160 val( 2 ) = one
161 val( 3 ) = sqrt( bignum )
162
163
164
165
166
167 10 CONTINUE
168 READ( nin, fmt = * )n, isrt
169 IF( n.EQ.0 )
170 $ RETURN
171 DO 20 i = 1, n
172 READ( nin, fmt = * )( tmp( i, j ), j = 1, n )
173 20 CONTINUE
174 DO 30 i = 1, n
175 READ( nin, fmt = * )wrin( i ), wiin( i ), sin( i ), sepin( i )
176 30 CONTINUE
177 tnrm =
zlange(
'M', n, n, tmp, ldt, rwork )
178 DO 260 iscl = 1, 3
179
180
181
182 knt = knt + 1
183 CALL zlacpy(
'F', n, n, tmp, ldt, t, ldt )
184 vmul = val( iscl )
185 DO 40 i = 1, n
186 CALL zdscal( n, vmul, t( 1, i ), 1 )
187 40 CONTINUE
188 IF( tnrm.EQ.zero )
189 $ vmul = one
190
191
192
193 CALL zgehrd( n, 1, n, t, ldt, work( 1 ), work( n+1 ), lwork-n,
194 $ info )
195 IF( info.NE.0 ) THEN
196 lmax( 1 ) = knt
197 ninfo( 1 ) = ninfo( 1 ) + 1
198 GO TO 260
199 END IF
200 DO 60 j = 1, n - 2
201 DO 50 i = j + 2, n
202 t( i, j ) = zero
203 50 CONTINUE
204 60 CONTINUE
205
206
207
208 CALL zhseqr(
'S',
'N', n, 1, n, t, ldt, w, cdum, 1, work,
209 $ lwork, info )
210 IF( info.NE.0 ) THEN
211 lmax( 2 ) = knt
212 ninfo( 2 ) = ninfo( 2 ) + 1
213 GO TO 260
214 END IF
215
216
217
218 DO 70 i = 1, n
219 SELECT( i ) = .true.
220 70 CONTINUE
221 CALL ztrevc(
'B',
'A',
SELECT, n, t, ldt, le, ldt, re, ldt, n,
222 $ m, work, rwork, info )
223
224
225
226 CALL ztrsna(
'B',
'A',
SELECT, n, t, ldt, le, ldt, re, ldt, s,
227 $ sep, n, m, work, n, rwork, info )
228 IF( info.NE.0 ) THEN
229 lmax( 3 ) = knt
230 ninfo( 3 ) = ninfo( 3 ) + 1
231 GO TO 260
232 END IF
233
234
235
236
237 CALL zcopy( n, w, 1, wtmp, 1 )
238 IF( isrt.EQ.0 ) THEN
239
240
241
242 DO 80 i = 1, n
243 wsrt( i ) = dble( w( i ) )
244 80 CONTINUE
245 ELSE
246
247
248
249 DO 90 i = 1, n
250 wsrt( i ) = dimag( w( i ) )
251 90 CONTINUE
252 END IF
253 CALL dcopy( n, s, 1, stmp, 1 )
254 CALL dcopy( n, sep, 1, septmp, 1 )
255 CALL dscal( n, one / vmul, septmp, 1 )
256 DO 110 i = 1, n - 1
257 kmin = i
258 vmin = wsrt( i )
259 DO 100 j = i + 1, n
260 IF( wsrt( j ).LT.vmin ) THEN
261 kmin = j
262 vmin = wsrt( j )
263 END IF
264 100 CONTINUE
265 wsrt( kmin ) = wsrt( i )
266 wsrt( i ) = vmin
267 vcmin = dble( wtmp( i ) )
268 wtmp( i ) = w( kmin )
269 wtmp( kmin ) = vcmin
270 vmin = stmp( kmin )
271 stmp( kmin ) = stmp( i )
272 stmp( i ) = vmin
273 vmin = septmp( kmin )
274 septmp( kmin ) = septmp( i )
275 septmp( i ) = vmin
276 110 CONTINUE
277
278
279
280
281 v = max( two*dble( n )*eps*tnrm, smlnum )
282 IF( tnrm.EQ.zero )
283 $ v = one
284 DO 120 i = 1, n
285 IF( v.GT.septmp( i ) ) THEN
286 tol = one
287 ELSE
288 tol = v / septmp( i )
289 END IF
290 IF( v.GT.sepin( i ) ) THEN
291 tolin = one
292 ELSE
293 tolin = v / sepin( i )
294 END IF
295 tol = max( tol, smlnum / eps )
296 tolin = max( tolin, smlnum / eps )
297 IF( eps*( sin( i )-tolin ).GT.stmp( i )+tol ) THEN
298 vmax = one / eps
299 ELSE IF( sin( i )-tolin.GT.stmp( i )+tol ) THEN
300 vmax = ( sin( i )-tolin ) / ( stmp( i )+tol )
301 ELSE IF( sin( i )+tolin.LT.eps*( stmp( i )-tol ) ) THEN
302 vmax = one / eps
303 ELSE IF( sin( i )+tolin.LT.stmp( i )-tol ) THEN
304 vmax = ( stmp( i )-tol ) / ( sin( i )+tolin )
305 ELSE
306 vmax = one
307 END IF
308 IF( vmax.GT.rmax( 2 ) ) THEN
309 rmax( 2 ) = vmax
310 IF( ninfo( 2 ).EQ.0 )
311 $ lmax( 2 ) = knt
312 END IF
313 120 CONTINUE
314
315
316
317
318 DO 130 i = 1, n
319 IF( v.GT.septmp( i )*stmp( i ) ) THEN
320 tol = septmp( i )
321 ELSE
322 tol = v / stmp( i )
323 END IF
324 IF( v.GT.sepin( i )*sin( i ) ) THEN
325 tolin = sepin( i )
326 ELSE
327 tolin = v / sin( i )
328 END IF
329 tol = max( tol, smlnum / eps )
330 tolin = max( tolin, smlnum / eps )
331 IF( eps*( sepin( i )-tolin ).GT.septmp( i )+tol ) THEN
332 vmax = one / eps
333 ELSE IF( sepin( i )-tolin.GT.septmp( i )+tol ) THEN
334 vmax = ( sepin( i )-tolin ) / ( septmp( i )+tol )
335 ELSE IF( sepin( i )+tolin.LT.eps*( septmp( i )-tol ) ) THEN
336 vmax = one / eps
337 ELSE IF( sepin( i )+tolin.LT.septmp( i )-tol ) THEN
338 vmax = ( septmp( i )-tol ) / ( sepin( i )+tolin )
339 ELSE
340 vmax = one
341 END IF
342 IF( vmax.GT.rmax( 2 ) ) THEN
343 rmax( 2 ) = vmax
344 IF( ninfo( 2 ).EQ.0 )
345 $ lmax( 2 ) = knt
346 END IF
347 130 CONTINUE
348
349
350
351
352 DO 140 i = 1, n
353 IF( sin( i ).LE.dble( 2*n )*eps .AND. stmp( i ).LE.
354 $ dble( 2*n )*eps ) THEN
355 vmax = one
356 ELSE IF( eps*sin( i ).GT.stmp( i ) ) THEN
357 vmax = one / eps
358 ELSE IF( sin( i ).GT.stmp( i ) ) THEN
359 vmax = sin( i ) / stmp( i )
360 ELSE IF( sin( i ).LT.eps*stmp( i ) ) THEN
361 vmax = one / eps
362 ELSE IF( sin( i ).LT.stmp( i ) ) THEN
363 vmax = stmp( i ) / sin( i )
364 ELSE
365 vmax = one
366 END IF
367 IF( vmax.GT.rmax( 3 ) ) THEN
368 rmax( 3 ) = vmax
369 IF( ninfo( 3 ).EQ.0 )
370 $ lmax( 3 ) = knt
371 END IF
372 140 CONTINUE
373
374
375
376
377 DO 150 i = 1, n
378 IF( sepin( i ).LE.v .AND. septmp( i ).LE.v ) THEN
379 vmax = one
380 ELSE IF( eps*sepin( i ).GT.septmp( i ) ) THEN
381 vmax = one / eps
382 ELSE IF( sepin( i ).GT.septmp( i ) ) THEN
383 vmax = sepin( i ) / septmp( i )
384 ELSE IF( sepin( i ).LT.eps*septmp( i ) ) THEN
385 vmax = one / eps
386 ELSE IF( sepin( i ).LT.septmp( i ) ) THEN
387 vmax = septmp( i ) / sepin( i )
388 ELSE
389 vmax = one
390 END IF
391 IF( vmax.GT.rmax( 3 ) ) THEN
392 rmax( 3 ) = vmax
393 IF( ninfo( 3 ).EQ.0 )
394 $ lmax( 3 ) = knt
395 END IF
396 150 CONTINUE
397
398
399
400 vmax = zero
401 dum( 1 ) = -one
402 CALL dcopy( n, dum, 0, stmp, 1 )
403 CALL dcopy( n, dum, 0, septmp, 1 )
404 CALL ztrsna(
'E',
'A',
SELECT, n, t, ldt, le, ldt, re, ldt,
405 $ stmp, septmp, n, m, work, n, rwork, info )
406 IF( info.NE.0 ) THEN
407 lmax( 3 ) = knt
408 ninfo( 3 ) = ninfo( 3 ) + 1
409 GO TO 260
410 END IF
411 DO 160 i = 1, n
412 IF( stmp( i ).NE.s( i ) )
413 $ vmax = one / eps
414 IF( septmp( i ).NE.dum( 1 ) )
415 $ vmax = one / eps
416 160 CONTINUE
417
418
419
420 CALL dcopy( n, dum, 0, stmp, 1 )
421 CALL dcopy( n, dum, 0, septmp, 1 )
422 CALL ztrsna(
'V',
'A',
SELECT, n, t, ldt, le, ldt, re, ldt,
423 $ stmp, septmp, n, m, work, n, rwork, info )
424 IF( info.NE.0 ) THEN
425 lmax( 3 ) = knt
426 ninfo( 3 ) = ninfo( 3 ) + 1
427 GO TO 260
428 END IF
429 DO 170 i = 1, n
430 IF( stmp( i ).NE.dum( 1 ) )
431 $ vmax = one / eps
432 IF( septmp( i ).NE.sep( i ) )
433 $ vmax = one / eps
434 170 CONTINUE
435
436
437
438 DO 180 i = 1, n
439 SELECT( i ) = .true.
440 180 CONTINUE
441 CALL dcopy( n, dum, 0, stmp, 1 )
442 CALL dcopy( n, dum, 0, septmp, 1 )
443 CALL ztrsna(
'B',
'S',
SELECT, n, t, ldt, le, ldt, re, ldt,
444 $ stmp, septmp, n, m, work, n, rwork, info )
445 IF( info.NE.0 ) THEN
446 lmax( 3 ) = knt
447 ninfo( 3 ) = ninfo( 3 ) + 1
448 GO TO 260
449 END IF
450 DO 190 i = 1, n
451 IF( septmp( i ).NE.sep( i ) )
452 $ vmax = one / eps
453 IF( stmp( i ).NE.s( i ) )
454 $ vmax = one / eps
455 190 CONTINUE
456
457
458
459 CALL dcopy( n, dum, 0, stmp, 1 )
460 CALL dcopy( n, dum, 0, septmp, 1 )
461 CALL ztrsna(
'E',
'S',
SELECT, n, t, ldt, le, ldt, re, ldt,
462 $ stmp, septmp, n, m, work, n, rwork, info )
463 IF( info.NE.0 ) THEN
464 lmax( 3 ) = knt
465 ninfo( 3 ) = ninfo( 3 ) + 1
466 GO TO 260
467 END IF
468 DO 200 i = 1, n
469 IF( stmp( i ).NE.s( i ) )
470 $ vmax = one / eps
471 IF( septmp( i ).NE.dum( 1 ) )
472 $ vmax = one / eps
473 200 CONTINUE
474
475
476
477 CALL dcopy( n, dum, 0, stmp, 1 )
478 CALL dcopy( n, dum, 0, septmp, 1 )
479 CALL ztrsna(
'V',
'S',
SELECT, n, t, ldt, le, ldt, re, ldt,
480 $ stmp, septmp, n, m, work, n, rwork, info )
481 IF( info.NE.0 ) THEN
482 lmax( 3 ) = knt
483 ninfo( 3 ) = ninfo( 3 ) + 1
484 GO TO 260
485 END IF
486 DO 210 i = 1, n
487 IF( stmp( i ).NE.dum( 1 ) )
488 $ vmax = one / eps
489 IF( septmp( i ).NE.sep( i ) )
490 $ vmax = one / eps
491 210 CONTINUE
492 IF( vmax.GT.rmax( 1 ) ) THEN
493 rmax( 1 ) = vmax
494 IF( ninfo( 1 ).EQ.0 )
495 $ lmax( 1 ) = knt
496 END IF
497
498
499
500 DO 220 i = 1, n
501 SELECT( i ) = .false.
502 220 CONTINUE
503 icmp = 0
504 IF( n.GT.1 ) THEN
505 icmp = 1
506 lcmp( 1 ) = 2
507 SELECT( 2 ) = .true.
508 CALL zcopy( n, re( 1, 2 ), 1, re( 1, 1 ), 1 )
509 CALL zcopy( n, le( 1, 2 ), 1, le( 1, 1 ), 1 )
510 END IF
511 IF( n.GT.3 ) THEN
512 icmp = 2
513 lcmp( 2 ) = n - 1
514 SELECT( n-1 ) = .true.
515 CALL zcopy( n, re( 1, n-1 ), 1, re( 1, 2 ), 1 )
516 CALL zcopy( n, le( 1, n-1 ), 1, le( 1, 2 ), 1 )
517 END IF
518
519
520
521 CALL dcopy( icmp, dum, 0, stmp, 1 )
522 CALL dcopy( icmp, dum, 0, septmp, 1 )
523 CALL ztrsna(
'B',
'S',
SELECT, n, t, ldt, le, ldt, re, ldt,
524 $ stmp, septmp, n, m, work, n, rwork, info )
525 IF( info.NE.0 ) THEN
526 lmax( 3 ) = knt
527 ninfo( 3 ) = ninfo( 3 ) + 1
528 GO TO 260
529 END IF
530 DO 230 i = 1, icmp
531 j = lcmp( i )
532 IF( septmp( i ).NE.sep( j ) )
533 $ vmax = one / eps
534 IF( stmp( i ).NE.s( j ) )
535 $ vmax = one / eps
536 230 CONTINUE
537
538
539
540 CALL dcopy( icmp, dum, 0, stmp, 1 )
541 CALL dcopy( icmp, dum, 0, septmp, 1 )
542 CALL ztrsna(
'E',
'S',
SELECT, n, t, ldt, le, ldt, re, ldt,
543 $ stmp, septmp, n, m, work, n, rwork, info )
544 IF( info.NE.0 ) THEN
545 lmax( 3 ) = knt
546 ninfo( 3 ) = ninfo( 3 ) + 1
547 GO TO 260
548 END IF
549 DO 240 i = 1, icmp
550 j = lcmp( i )
551 IF( stmp( i ).NE.s( j ) )
552 $ vmax = one / eps
553 IF( septmp( i ).NE.dum( 1 ) )
554 $ vmax = one / eps
555 240 CONTINUE
556
557
558
559 CALL dcopy( icmp, dum, 0, stmp, 1 )
560 CALL dcopy( icmp, dum, 0, septmp, 1 )
561 CALL ztrsna(
'V',
'S',
SELECT, n, t, ldt, le, ldt, re, ldt,
562 $ stmp, septmp, n, m, work, n, rwork, info )
563 IF( info.NE.0 ) THEN
564 lmax( 3 ) = knt
565 ninfo( 3 ) = ninfo( 3 ) + 1
566 GO TO 260
567 END IF
568 DO 250 i = 1, icmp
569 j = lcmp( i )
570 IF( stmp( i ).NE.dum( 1 ) )
571 $ vmax = one / eps
572 IF( septmp( i ).NE.sep( j ) )
573 $ vmax = one / eps
574 250 CONTINUE
575 IF( vmax.GT.rmax( 1 ) ) THEN
576 rmax( 1 ) = vmax
577 IF( ninfo( 1 ).EQ.0 )
578 $ lmax( 1 ) = knt
579 END IF
580 260 CONTINUE
581 GO TO 10
582
583
584
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
subroutine zgehrd(n, ilo, ihi, a, lda, tau, work, lwork, info)
ZGEHRD
subroutine zhseqr(job, compz, n, ilo, ihi, h, ldh, w, z, ldz, work, lwork, info)
ZHSEQR
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
double precision function zlange(norm, m, n, a, lda, work)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine dscal(n, da, dx, incx)
DSCAL
subroutine zdscal(n, da, zx, incx)
ZDSCAL
subroutine ztrevc(side, howmny, select, n, t, ldt, vl, ldvl, vr, ldvr, mm, m, work, rwork, info)
ZTREVC
subroutine ztrsna(job, howmny, select, n, t, ldt, vl, ldvl, vr, ldvr, s, sep, mm, m, work, ldwork, rwork, info)
ZTRSNA