LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
cget37.f
Go to the documentation of this file.
1*> \brief \b CGET37
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8* Definition:
9* ===========
10*
11* SUBROUTINE CGET37( RMAX, LMAX, NINFO, KNT, NIN )
12*
13* .. Scalar Arguments ..
14* INTEGER KNT, NIN
15* ..
16* .. Array Arguments ..
17* INTEGER LMAX( 3 ), NINFO( 3 )
18* REAL RMAX( 3 )
19* ..
20*
21*
22*> \par Purpose:
23* =============
24*>
25*> \verbatim
26*>
27*> CGET37 tests CTRSNA, a routine for estimating condition numbers of
28*> eigenvalues and/or right eigenvectors of a matrix.
29*>
30*> The test matrices are read from a file with logical unit number NIN.
31*> \endverbatim
32*
33* Arguments:
34* ==========
35*
36*> \param[out] RMAX
37*> \verbatim
38*> RMAX is REAL array, dimension (3)
39*> Value of the largest test ratio.
40*> RMAX(1) = largest ratio comparing different calls to CTRSNA
41*> RMAX(2) = largest error in reciprocal condition
42*> numbers taking their conditioning into account
43*> RMAX(3) = largest error in reciprocal condition
44*> numbers not taking their conditioning into
45*> account (may be larger than RMAX(2))
46*> \endverbatim
47*>
48*> \param[out] LMAX
49*> \verbatim
50*> LMAX is INTEGER array, dimension (3)
51*> LMAX(i) is example number where largest test ratio
52*> RMAX(i) is achieved. Also:
53*> If CGEHRD returns INFO nonzero on example i, LMAX(1)=i
54*> If CHSEQR returns INFO nonzero on example i, LMAX(2)=i
55*> If CTRSNA returns INFO nonzero on example i, LMAX(3)=i
56*> \endverbatim
57*>
58*> \param[out] NINFO
59*> \verbatim
60*> NINFO is INTEGER array, dimension (3)
61*> NINFO(1) = No. of times CGEHRD returned INFO nonzero
62*> NINFO(2) = No. of times CHSEQR returned INFO nonzero
63*> NINFO(3) = No. of times CTRSNA returned INFO nonzero
64*> \endverbatim
65*>
66*> \param[out] KNT
67*> \verbatim
68*> KNT is INTEGER
69*> Total number of examples tested.
70*> \endverbatim
71*>
72*> \param[in] NIN
73*> \verbatim
74*> NIN is INTEGER
75*> Input logical unit number
76*> \endverbatim
77*
78* Authors:
79* ========
80*
81*> \author Univ. of Tennessee
82*> \author Univ. of California Berkeley
83*> \author Univ. of Colorado Denver
84*> \author NAG Ltd.
85*
86*> \ingroup complex_eig
87*
88* =====================================================================
89 SUBROUTINE cget37( RMAX, LMAX, NINFO, KNT, NIN )
90*
91* -- LAPACK test routine --
92* -- LAPACK is a software package provided by Univ. of Tennessee, --
93* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
94*
95* .. Scalar Arguments ..
96 INTEGER KNT, NIN
97* ..
98* .. Array Arguments ..
99 INTEGER LMAX( 3 ), NINFO( 3 )
100 REAL RMAX( 3 )
101* ..
102*
103* =====================================================================
104*
105* .. Parameters ..
106 REAL ZERO, ONE, TWO
107 parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
108 REAL EPSIN
109 parameter( epsin = 5.9605e-8 )
110 INTEGER LDT, LWORK
111 parameter( ldt = 20, lwork = 2*ldt*( 10+ldt ) )
112* ..
113* .. Local Scalars ..
114 INTEGER I, ICMP, INFO, ISCL, ISRT, J, KMIN, M, N
115 REAL BIGNUM, EPS, SMLNUM, TNRM, TOL, TOLIN, V,
116 $ VCMIN, VMAX, VMIN, VMUL
117* ..
118* .. Local Arrays ..
119 LOGICAL SELECT( LDT )
120 INTEGER LCMP( 3 )
121 REAL 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 CDUM( 1 ), LE( LDT, LDT ), RE( LDT, LDT ),
126 $ T( LDT, LDT ), TMP( LDT, LDT ), W( LDT ),
127 $ WORK( LWORK ), WTMP( LDT )
128* ..
129* .. External Functions ..
130 REAL CLANGE, SLAMCH
131 EXTERNAL clange, slamch
132* ..
133* .. External Subroutines ..
134 EXTERNAL ccopy, cgehrd, chseqr, clacpy, csscal, ctrevc,
135 $ ctrsna, scopy, sscal
136* ..
137* .. Intrinsic Functions ..
138 INTRINSIC aimag, max, real, sqrt
139* ..
140* .. Executable Statements ..
141*
142 eps = slamch( 'P' )
143 smlnum = slamch( 'S' ) / eps
144 bignum = one / smlnum
145*
146* EPSIN = 2**(-24) = precision to which input data computed
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* Read input data until N=0. Assume input eigenvalues are sorted
164* lexicographically (increasing by real part if ISRT = 0,
165* increasing by imaginary part if ISRT = 1)
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 = clange( 'M', n, n, tmp, ldt, rwork )
178 DO 260 iscl = 1, 3
179*
180* Scale input matrix
181*
182 knt = knt + 1
183 CALL clacpy( 'F', n, n, tmp, ldt, t, ldt )
184 vmul = val( iscl )
185 DO 40 i = 1, n
186 CALL csscal( n, vmul, t( 1, i ), 1 )
187 40 CONTINUE
188 IF( tnrm.EQ.zero )
189 $ vmul = one
190*
191* Compute eigenvalues and eigenvectors
192*
193 CALL cgehrd( 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* Compute Schur form
207*
208 CALL chseqr( '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* Compute eigenvectors
217*
218 DO 70 i = 1, n
219 SELECT( i ) = .true.
220 70 CONTINUE
221 CALL ctrevc( 'B', 'A', SELECT, n, t, ldt, le, ldt, re, ldt, n,
222 $ m, work, rwork, info )
223*
224* Compute condition numbers
225*
226 CALL ctrsna( '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* Sort eigenvalues and condition numbers lexicographically
235* to compare with inputs
236*
237 CALL ccopy( n, w, 1, wtmp, 1 )
238 IF( isrt.EQ.0 ) THEN
239*
240* Sort by increasing real part
241*
242 DO 80 i = 1, n
243 wsrt( i ) = real( w( i ) )
244 80 CONTINUE
245 ELSE
246*
247* Sort by increasing imaginary part
248*
249 DO 90 i = 1, n
250 wsrt( i ) = aimag( w( i ) )
251 90 CONTINUE
252 END IF
253 CALL scopy( n, s, 1, stmp, 1 )
254 CALL scopy( n, sep, 1, septmp, 1 )
255 CALL sscal( 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 = real( 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* Compare condition numbers for eigenvalues
279* taking their condition numbers into account
280*
281 v = max( two*real( 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* Compare condition numbers for eigenvectors
316* taking their condition numbers into account
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* Compare condition numbers for eigenvalues
350* without taking their condition numbers into account
351*
352 DO 140 i = 1, n
353 IF( sin( i ).LE.real( 2*n )*eps .AND. stmp( i ).LE.
354 $ real( 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* Compare condition numbers for eigenvectors
375* without taking their condition numbers into account
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* Compute eigenvalue condition numbers only and compare
399*
400 vmax = zero
401 dum( 1 ) = -one
402 CALL scopy( n, dum, 0, stmp, 1 )
403 CALL scopy( n, dum, 0, septmp, 1 )
404 CALL ctrsna( '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* Compute eigenvector condition numbers only and compare
419*
420 CALL scopy( n, dum, 0, stmp, 1 )
421 CALL scopy( n, dum, 0, septmp, 1 )
422 CALL ctrsna( '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* Compute all condition numbers using SELECT and compare
437*
438 DO 180 i = 1, n
439 SELECT( i ) = .true.
440 180 CONTINUE
441 CALL scopy( n, dum, 0, stmp, 1 )
442 CALL scopy( n, dum, 0, septmp, 1 )
443 CALL ctrsna( '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* Compute eigenvalue condition numbers using SELECT and compare
458*
459 CALL scopy( n, dum, 0, stmp, 1 )
460 CALL scopy( n, dum, 0, septmp, 1 )
461 CALL ctrsna( '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* Compute eigenvector condition numbers using SELECT and compare
476*
477 CALL scopy( n, dum, 0, stmp, 1 )
478 CALL scopy( n, dum, 0, septmp, 1 )
479 CALL ctrsna( '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* Select second and next to last eigenvalues
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 ccopy( n, re( 1, 2 ), 1, re( 1, 1 ), 1 )
509 CALL ccopy( 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 ccopy( n, re( 1, n-1 ), 1, re( 1, 2 ), 1 )
516 CALL ccopy( n, le( 1, n-1 ), 1, le( 1, 2 ), 1 )
517 END IF
518*
519* Compute all selected condition numbers
520*
521 CALL scopy( icmp, dum, 0, stmp, 1 )
522 CALL scopy( icmp, dum, 0, septmp, 1 )
523 CALL ctrsna( '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* Compute selected eigenvalue condition numbers
539*
540 CALL scopy( icmp, dum, 0, stmp, 1 )
541 CALL scopy( icmp, dum, 0, septmp, 1 )
542 CALL ctrsna( '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* Compute selected eigenvector condition numbers
558*
559 CALL scopy( icmp, dum, 0, stmp, 1 )
560 CALL scopy( icmp, dum, 0, septmp, 1 )
561 CALL ctrsna( '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* End of CGET37
584*
585 END
cget37
subroutine cget37(rmax, lmax, ninfo, knt, nin)
CGET37
Definition cget37.f:90
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
ccopy
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
cgehrd
subroutine cgehrd(n, ilo, ihi, a, lda, tau, work, lwork, info)
CGEHRD
Definition cgehrd.f:167
chseqr
subroutine chseqr(job, compz, n, ilo, ihi, h, ldh, w, z, ldz, work, lwork, info)
CHSEQR
Definition chseqr.f:299
clacpy
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:103
csscal
subroutine csscal(n, sa, cx, incx)
CSSCAL
Definition csscal.f:78
sscal
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
ctrevc
subroutine ctrevc(side, howmny, select, n, t, ldt, vl, ldvl, vr, ldvr, mm, m, work, rwork, info)
CTREVC
Definition ctrevc.f:218
ctrsna
subroutine ctrsna(job, howmny, select, n, t, ldt, vl, ldvl, vr, ldvr, s, sep, mm, m, work, ldwork, rwork, info)
CTRSNA
Definition ctrsna.f:249