LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ sget37()

subroutine sget37 ( real, dimension( 3 )  RMAX,
integer, dimension( 3 )  LMAX,
integer, dimension( 3 )  NINFO,
integer  KNT,
integer  NIN 
)

SGET37

Purpose:
 SGET37 tests STRSNA, a routine for estimating condition numbers of
 eigenvalues and/or right eigenvectors of a matrix.

 The test matrices are read from a file with logical unit number NIN.
Parameters
[out]RMAX
          RMAX is REAL array, dimension (3)
          Value of the largest test ratio.
          RMAX(1) = largest ratio comparing different calls to STRSNA
          RMAX(2) = largest error in reciprocal condition
                    numbers taking their conditioning into account
          RMAX(3) = largest error in reciprocal condition
                    numbers not taking their conditioning into
                    account (may be larger than RMAX(2))
[out]LMAX
          LMAX is INTEGER array, dimension (3)
          LMAX(i) is example number where largest test ratio
          RMAX(i) is achieved. Also:
          If SGEHRD returns INFO nonzero on example i, LMAX(1)=i
          If SHSEQR returns INFO nonzero on example i, LMAX(2)=i
          If STRSNA returns INFO nonzero on example i, LMAX(3)=i
[out]NINFO
          NINFO is INTEGER array, dimension (3)
          NINFO(1) = No. of times SGEHRD returned INFO nonzero
          NINFO(2) = No. of times SHSEQR returned INFO nonzero
          NINFO(3) = No. of times STRSNA returned INFO nonzero
[out]KNT
          KNT is INTEGER
          Total number of examples tested.
[in]NIN
          NIN is INTEGER
          Input logical unit number
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 89 of file sget37.f.

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, IFND, INFO, ISCL, J, KMIN, M, N
115 REAL BIGNUM, EPS, SMLNUM, TNRM, TOL, TOLIN, V,
116 $ VIMIN, VMAX, VMUL, VRMIN
117* ..
118* .. Local Arrays ..
119 LOGICAL SELECT( LDT )
120 INTEGER IWORK( 2*LDT ), LCMP( 3 )
121 REAL DUM( 1 ), LE( LDT, LDT ), RE( LDT, LDT ),
122 $ S( LDT ), SEP( LDT ), SEPIN( LDT ),
123 $ SEPTMP( LDT ), SIN( LDT ), STMP( LDT ),
124 $ T( LDT, LDT ), TMP( LDT, LDT ), VAL( 3 ),
125 $ WI( LDT ), WIIN( LDT ), WITMP( LDT ),
126 $ WORK( LWORK ), WR( LDT ), WRIN( LDT ),
127 $ WRTMP( LDT )
128* ..
129* .. External Functions ..
130 REAL SLAMCH, SLANGE
131 EXTERNAL slamch, slange
132* ..
133* .. External Subroutines ..
134 EXTERNAL scopy, sgehrd, shseqr, slabad, slacpy, sscal,
135 $ strevc, strsna
136* ..
137* .. Intrinsic Functions ..
138 INTRINSIC max, real, sqrt
139* ..
140* .. Executable Statements ..
141*
142 eps = slamch( 'P' )
143 smlnum = slamch( 'S' ) / eps
144 bignum = one / smlnum
145 CALL slabad( smlnum, bignum )
146*
147* EPSIN = 2**(-24) = precision to which input data computed
148*
149 eps = max( eps, epsin )
150 rmax( 1 ) = zero
151 rmax( 2 ) = zero
152 rmax( 3 ) = zero
153 lmax( 1 ) = 0
154 lmax( 2 ) = 0
155 lmax( 3 ) = 0
156 knt = 0
157 ninfo( 1 ) = 0
158 ninfo( 2 ) = 0
159 ninfo( 3 ) = 0
160*
161 val( 1 ) = sqrt( smlnum )
162 val( 2 ) = one
163 val( 3 ) = sqrt( bignum )
164*
165* Read input data until N=0. Assume input eigenvalues are sorted
166* lexicographically (increasing by real part, then decreasing by
167* imaginary part)
168*
169 10 CONTINUE
170 READ( nin, fmt = * )n
171 IF( n.EQ.0 )
172 $ RETURN
173 DO 20 i = 1, n
174 READ( nin, fmt = * )( tmp( i, j ), j = 1, n )
175 20 CONTINUE
176 DO 30 i = 1, n
177 READ( nin, fmt = * )wrin( i ), wiin( i ), sin( i ), sepin( i )
178 30 CONTINUE
179 tnrm = slange( 'M', n, n, tmp, ldt, work )
180*
181* Begin test
182*
183 DO 240 iscl = 1, 3
184*
185* Scale input matrix
186*
187 knt = knt + 1
188 CALL slacpy( 'F', n, n, tmp, ldt, t, ldt )
189 vmul = val( iscl )
190 DO 40 i = 1, n
191 CALL sscal( n, vmul, t( 1, i ), 1 )
192 40 CONTINUE
193 IF( tnrm.EQ.zero )
194 $ vmul = one
195*
196* Compute eigenvalues and eigenvectors
197*
198 CALL sgehrd( n, 1, n, t, ldt, work( 1 ), work( n+1 ), lwork-n,
199 $ info )
200 IF( info.NE.0 ) THEN
201 lmax( 1 ) = knt
202 ninfo( 1 ) = ninfo( 1 ) + 1
203 GO TO 240
204 END IF
205 DO 60 j = 1, n - 2
206 DO 50 i = j + 2, n
207 t( i, j ) = zero
208 50 CONTINUE
209 60 CONTINUE
210*
211* Compute Schur form
212*
213 CALL shseqr( 'S', 'N', n, 1, n, t, ldt, wr, wi, dum, 1, work,
214 $ lwork, info )
215 IF( info.NE.0 ) THEN
216 lmax( 2 ) = knt
217 ninfo( 2 ) = ninfo( 2 ) + 1
218 GO TO 240
219 END IF
220*
221* Compute eigenvectors
222*
223 CALL strevc( 'Both', 'All', SELECT, n, t, ldt, le, ldt, re,
224 $ ldt, n, m, work, info )
225*
226* Compute condition numbers
227*
228 CALL strsna( 'Both', 'All', SELECT, n, t, ldt, le, ldt, re,
229 $ ldt, s, sep, n, m, work, n, iwork, info )
230 IF( info.NE.0 ) THEN
231 lmax( 3 ) = knt
232 ninfo( 3 ) = ninfo( 3 ) + 1
233 GO TO 240
234 END IF
235*
236* Sort eigenvalues and condition numbers lexicographically
237* to compare with inputs
238*
239 CALL scopy( n, wr, 1, wrtmp, 1 )
240 CALL scopy( n, wi, 1, witmp, 1 )
241 CALL scopy( n, s, 1, stmp, 1 )
242 CALL scopy( n, sep, 1, septmp, 1 )
243 CALL sscal( n, one / vmul, septmp, 1 )
244 DO 80 i = 1, n - 1
245 kmin = i
246 vrmin = wrtmp( i )
247 vimin = witmp( i )
248 DO 70 j = i + 1, n
249 IF( wrtmp( j ).LT.vrmin ) THEN
250 kmin = j
251 vrmin = wrtmp( j )
252 vimin = witmp( j )
253 END IF
254 70 CONTINUE
255 wrtmp( kmin ) = wrtmp( i )
256 witmp( kmin ) = witmp( i )
257 wrtmp( i ) = vrmin
258 witmp( i ) = vimin
259 vrmin = stmp( kmin )
260 stmp( kmin ) = stmp( i )
261 stmp( i ) = vrmin
262 vrmin = septmp( kmin )
263 septmp( kmin ) = septmp( i )
264 septmp( i ) = vrmin
265 80 CONTINUE
266*
267* Compare condition numbers for eigenvalues
268* taking their condition numbers into account
269*
270 v = max( two*real( n )*eps*tnrm, smlnum )
271 IF( tnrm.EQ.zero )
272 $ v = one
273 DO 90 i = 1, n
274 IF( v.GT.septmp( i ) ) THEN
275 tol = one
276 ELSE
277 tol = v / septmp( i )
278 END IF
279 IF( v.GT.sepin( i ) ) THEN
280 tolin = one
281 ELSE
282 tolin = v / sepin( i )
283 END IF
284 tol = max( tol, smlnum / eps )
285 tolin = max( tolin, smlnum / eps )
286 IF( eps*( sin( i )-tolin ).GT.stmp( i )+tol ) THEN
287 vmax = one / eps
288 ELSE IF( sin( i )-tolin.GT.stmp( i )+tol ) THEN
289 vmax = ( sin( i )-tolin ) / ( stmp( i )+tol )
290 ELSE IF( sin( i )+tolin.LT.eps*( stmp( i )-tol ) ) THEN
291 vmax = one / eps
292 ELSE IF( sin( i )+tolin.LT.stmp( i )-tol ) THEN
293 vmax = ( stmp( i )-tol ) / ( sin( i )+tolin )
294 ELSE
295 vmax = one
296 END IF
297 IF( vmax.GT.rmax( 2 ) ) THEN
298 rmax( 2 ) = vmax
299 IF( ninfo( 2 ).EQ.0 )
300 $ lmax( 2 ) = knt
301 END IF
302 90 CONTINUE
303*
304* Compare condition numbers for eigenvectors
305* taking their condition numbers into account
306*
307 DO 100 i = 1, n
308 IF( v.GT.septmp( i )*stmp( i ) ) THEN
309 tol = septmp( i )
310 ELSE
311 tol = v / stmp( i )
312 END IF
313 IF( v.GT.sepin( i )*sin( i ) ) THEN
314 tolin = sepin( i )
315 ELSE
316 tolin = v / sin( i )
317 END IF
318 tol = max( tol, smlnum / eps )
319 tolin = max( tolin, smlnum / eps )
320 IF( eps*( sepin( i )-tolin ).GT.septmp( i )+tol ) THEN
321 vmax = one / eps
322 ELSE IF( sepin( i )-tolin.GT.septmp( i )+tol ) THEN
323 vmax = ( sepin( i )-tolin ) / ( septmp( i )+tol )
324 ELSE IF( sepin( i )+tolin.LT.eps*( septmp( i )-tol ) ) THEN
325 vmax = one / eps
326 ELSE IF( sepin( i )+tolin.LT.septmp( i )-tol ) THEN
327 vmax = ( septmp( i )-tol ) / ( sepin( i )+tolin )
328 ELSE
329 vmax = one
330 END IF
331 IF( vmax.GT.rmax( 2 ) ) THEN
332 rmax( 2 ) = vmax
333 IF( ninfo( 2 ).EQ.0 )
334 $ lmax( 2 ) = knt
335 END IF
336 100 CONTINUE
337*
338* Compare condition numbers for eigenvalues
339* without taking their condition numbers into account
340*
341 DO 110 i = 1, n
342 IF( sin( i ).LE.real( 2*n )*eps .AND. stmp( i ).LE.
343 $ real( 2*n )*eps ) THEN
344 vmax = one
345 ELSE IF( eps*sin( i ).GT.stmp( i ) ) THEN
346 vmax = one / eps
347 ELSE IF( sin( i ).GT.stmp( i ) ) THEN
348 vmax = sin( i ) / stmp( i )
349 ELSE IF( sin( i ).LT.eps*stmp( i ) ) THEN
350 vmax = one / eps
351 ELSE IF( sin( i ).LT.stmp( i ) ) THEN
352 vmax = stmp( i ) / sin( i )
353 ELSE
354 vmax = one
355 END IF
356 IF( vmax.GT.rmax( 3 ) ) THEN
357 rmax( 3 ) = vmax
358 IF( ninfo( 3 ).EQ.0 )
359 $ lmax( 3 ) = knt
360 END IF
361 110 CONTINUE
362*
363* Compare condition numbers for eigenvectors
364* without taking their condition numbers into account
365*
366 DO 120 i = 1, n
367 IF( sepin( i ).LE.v .AND. septmp( i ).LE.v ) THEN
368 vmax = one
369 ELSE IF( eps*sepin( i ).GT.septmp( i ) ) THEN
370 vmax = one / eps
371 ELSE IF( sepin( i ).GT.septmp( i ) ) THEN
372 vmax = sepin( i ) / septmp( i )
373 ELSE IF( sepin( i ).LT.eps*septmp( i ) ) THEN
374 vmax = one / eps
375 ELSE IF( sepin( i ).LT.septmp( i ) ) THEN
376 vmax = septmp( i ) / sepin( i )
377 ELSE
378 vmax = one
379 END IF
380 IF( vmax.GT.rmax( 3 ) ) THEN
381 rmax( 3 ) = vmax
382 IF( ninfo( 3 ).EQ.0 )
383 $ lmax( 3 ) = knt
384 END IF
385 120 CONTINUE
386*
387* Compute eigenvalue condition numbers only and compare
388*
389 vmax = zero
390 dum( 1 ) = -one
391 CALL scopy( n, dum, 0, stmp, 1 )
392 CALL scopy( n, dum, 0, septmp, 1 )
393 CALL strsna( 'Eigcond', 'All', SELECT, n, t, ldt, le, ldt, re,
394 $ ldt, stmp, septmp, n, m, work, n, iwork, info )
395 IF( info.NE.0 ) THEN
396 lmax( 3 ) = knt
397 ninfo( 3 ) = ninfo( 3 ) + 1
398 GO TO 240
399 END IF
400 DO 130 i = 1, n
401 IF( stmp( i ).NE.s( i ) )
402 $ vmax = one / eps
403 IF( septmp( i ).NE.dum( 1 ) )
404 $ vmax = one / eps
405 130 CONTINUE
406*
407* Compute eigenvector condition numbers only and compare
408*
409 CALL scopy( n, dum, 0, stmp, 1 )
410 CALL scopy( n, dum, 0, septmp, 1 )
411 CALL strsna( 'Veccond', 'All', SELECT, n, t, ldt, le, ldt, re,
412 $ ldt, stmp, septmp, n, m, work, n, iwork, info )
413 IF( info.NE.0 ) THEN
414 lmax( 3 ) = knt
415 ninfo( 3 ) = ninfo( 3 ) + 1
416 GO TO 240
417 END IF
418 DO 140 i = 1, n
419 IF( stmp( i ).NE.dum( 1 ) )
420 $ vmax = one / eps
421 IF( septmp( i ).NE.sep( i ) )
422 $ vmax = one / eps
423 140 CONTINUE
424*
425* Compute all condition numbers using SELECT and compare
426*
427 DO 150 i = 1, n
428 SELECT( i ) = .true.
429 150 CONTINUE
430 CALL scopy( n, dum, 0, stmp, 1 )
431 CALL scopy( n, dum, 0, septmp, 1 )
432 CALL strsna( 'Bothcond', 'Some', SELECT, n, t, ldt, le, ldt,
433 $ re, ldt, stmp, septmp, n, m, work, n, iwork,
434 $ info )
435 IF( info.NE.0 ) THEN
436 lmax( 3 ) = knt
437 ninfo( 3 ) = ninfo( 3 ) + 1
438 GO TO 240
439 END IF
440 DO 160 i = 1, n
441 IF( septmp( i ).NE.sep( i ) )
442 $ vmax = one / eps
443 IF( stmp( i ).NE.s( i ) )
444 $ vmax = one / eps
445 160 CONTINUE
446*
447* Compute eigenvalue condition numbers using SELECT and compare
448*
449 CALL scopy( n, dum, 0, stmp, 1 )
450 CALL scopy( n, dum, 0, septmp, 1 )
451 CALL strsna( 'Eigcond', 'Some', SELECT, n, t, ldt, le, ldt, re,
452 $ ldt, stmp, septmp, n, m, work, n, iwork, info )
453 IF( info.NE.0 ) THEN
454 lmax( 3 ) = knt
455 ninfo( 3 ) = ninfo( 3 ) + 1
456 GO TO 240
457 END IF
458 DO 170 i = 1, n
459 IF( stmp( i ).NE.s( i ) )
460 $ vmax = one / eps
461 IF( septmp( i ).NE.dum( 1 ) )
462 $ vmax = one / eps
463 170 CONTINUE
464*
465* Compute eigenvector condition numbers using SELECT and compare
466*
467 CALL scopy( n, dum, 0, stmp, 1 )
468 CALL scopy( n, dum, 0, septmp, 1 )
469 CALL strsna( 'Veccond', 'Some', SELECT, n, t, ldt, le, ldt, re,
470 $ ldt, stmp, septmp, n, m, work, n, iwork, info )
471 IF( info.NE.0 ) THEN
472 lmax( 3 ) = knt
473 ninfo( 3 ) = ninfo( 3 ) + 1
474 GO TO 240
475 END IF
476 DO 180 i = 1, n
477 IF( stmp( i ).NE.dum( 1 ) )
478 $ vmax = one / eps
479 IF( septmp( i ).NE.sep( i ) )
480 $ vmax = one / eps
481 180 CONTINUE
482 IF( vmax.GT.rmax( 1 ) ) THEN
483 rmax( 1 ) = vmax
484 IF( ninfo( 1 ).EQ.0 )
485 $ lmax( 1 ) = knt
486 END IF
487*
488* Select first real and first complex eigenvalue
489*
490 IF( wi( 1 ).EQ.zero ) THEN
491 lcmp( 1 ) = 1
492 ifnd = 0
493 DO 190 i = 2, n
494 IF( ifnd.EQ.1 .OR. wi( i ).EQ.zero ) THEN
495 SELECT( i ) = .false.
496 ELSE
497 ifnd = 1
498 lcmp( 2 ) = i
499 lcmp( 3 ) = i + 1
500 CALL scopy( n, re( 1, i ), 1, re( 1, 2 ), 1 )
501 CALL scopy( n, re( 1, i+1 ), 1, re( 1, 3 ), 1 )
502 CALL scopy( n, le( 1, i ), 1, le( 1, 2 ), 1 )
503 CALL scopy( n, le( 1, i+1 ), 1, le( 1, 3 ), 1 )
504 END IF
505 190 CONTINUE
506 IF( ifnd.EQ.0 ) THEN
507 icmp = 1
508 ELSE
509 icmp = 3
510 END IF
511 ELSE
512 lcmp( 1 ) = 1
513 lcmp( 2 ) = 2
514 ifnd = 0
515 DO 200 i = 3, n
516 IF( ifnd.EQ.1 .OR. wi( i ).NE.zero ) THEN
517 SELECT( i ) = .false.
518 ELSE
519 lcmp( 3 ) = i
520 ifnd = 1
521 CALL scopy( n, re( 1, i ), 1, re( 1, 3 ), 1 )
522 CALL scopy( n, le( 1, i ), 1, le( 1, 3 ), 1 )
523 END IF
524 200 CONTINUE
525 IF( ifnd.EQ.0 ) THEN
526 icmp = 2
527 ELSE
528 icmp = 3
529 END IF
530 END IF
531*
532* Compute all selected condition numbers
533*
534 CALL scopy( icmp, dum, 0, stmp, 1 )
535 CALL scopy( icmp, dum, 0, septmp, 1 )
536 CALL strsna( 'Bothcond', 'Some', SELECT, n, t, ldt, le, ldt,
537 $ re, ldt, stmp, septmp, n, m, work, n, iwork,
538 $ info )
539 IF( info.NE.0 ) THEN
540 lmax( 3 ) = knt
541 ninfo( 3 ) = ninfo( 3 ) + 1
542 GO TO 240
543 END IF
544 DO 210 i = 1, icmp
545 j = lcmp( i )
546 IF( septmp( i ).NE.sep( j ) )
547 $ vmax = one / eps
548 IF( stmp( i ).NE.s( j ) )
549 $ vmax = one / eps
550 210 CONTINUE
551*
552* Compute selected eigenvalue condition numbers
553*
554 CALL scopy( icmp, dum, 0, stmp, 1 )
555 CALL scopy( icmp, dum, 0, septmp, 1 )
556 CALL strsna( 'Eigcond', 'Some', SELECT, n, t, ldt, le, ldt, re,
557 $ ldt, stmp, septmp, n, m, work, n, iwork, info )
558 IF( info.NE.0 ) THEN
559 lmax( 3 ) = knt
560 ninfo( 3 ) = ninfo( 3 ) + 1
561 GO TO 240
562 END IF
563 DO 220 i = 1, icmp
564 j = lcmp( i )
565 IF( stmp( i ).NE.s( j ) )
566 $ vmax = one / eps
567 IF( septmp( i ).NE.dum( 1 ) )
568 $ vmax = one / eps
569 220 CONTINUE
570*
571* Compute selected eigenvector condition numbers
572*
573 CALL scopy( icmp, dum, 0, stmp, 1 )
574 CALL scopy( icmp, dum, 0, septmp, 1 )
575 CALL strsna( 'Veccond', 'Some', SELECT, n, t, ldt, le, ldt, re,
576 $ ldt, stmp, septmp, n, m, work, n, iwork, info )
577 IF( info.NE.0 ) THEN
578 lmax( 3 ) = knt
579 ninfo( 3 ) = ninfo( 3 ) + 1
580 GO TO 240
581 END IF
582 DO 230 i = 1, icmp
583 j = lcmp( i )
584 IF( stmp( i ).NE.dum( 1 ) )
585 $ vmax = one / eps
586 IF( septmp( i ).NE.sep( j ) )
587 $ vmax = one / eps
588 230 CONTINUE
589 IF( vmax.GT.rmax( 1 ) ) THEN
590 rmax( 1 ) = vmax
591 IF( ninfo( 1 ).EQ.0 )
592 $ lmax( 1 ) = knt
593 END IF
594 240 CONTINUE
595 GO TO 10
596*
597* End of SGET37
598*
slabad
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
slacpy
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
slange
real function slange(NORM, M, N, A, LDA, WORK)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: slange.f:114
sgehrd
subroutine sgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
SGEHRD
Definition: sgehrd.f:167
strevc
subroutine strevc(SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, MM, M, WORK, INFO)
STREVC
Definition: strevc.f:222
strsna
subroutine strsna(JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK, INFO)
STRSNA
Definition: strsna.f:265
shseqr
subroutine shseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z, LDZ, WORK, LWORK, INFO)
SHSEQR
Definition: shseqr.f:316
scopy
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
sscal
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: