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

◆ cget38()

subroutine cget38 ( real, dimension( 3 )  rmax,
integer, dimension( 3 )  lmax,
integer, dimension( 3 )  ninfo,
integer  knt,
integer  nin 
)

CGET38

Purpose:
 CGET38 tests CTRSEN, a routine for estimating condition numbers of a
 cluster of eigenvalues and/or its associated right invariant subspace

 The test matrices are read from a file with logical unit number NIN.
Parameters
[out]RMAX
          RMAX is REAL array, dimension (3)
          Values of the largest test ratios.
          RMAX(1) = largest residuals from CHST01 or comparing
                    different calls to CTRSEN
          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 CGEHRD returns INFO nonzero on example i, LMAX(1)=i
          If CHSEQR returns INFO nonzero on example i, LMAX(2)=i
          If CTRSEN returns INFO nonzero on example i, LMAX(3)=i
[out]NINFO
          NINFO is INTEGER array, dimension (3)
          NINFO(1) = No. of times CGEHRD returned INFO nonzero
          NINFO(2) = No. of times CHSEQR returned INFO nonzero
          NINFO(3) = No. of times CTRSEN 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 90 of file cget38.f.

91*
92* -- LAPACK test routine --
93* -- LAPACK is a software package provided by Univ. of Tennessee, --
94* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
95*
96* .. Scalar Arguments ..
97 INTEGER KNT, NIN
98* ..
99* .. Array Arguments ..
100 INTEGER LMAX( 3 ), NINFO( 3 )
101 REAL RMAX( 3 )
102* ..
103*
104* =====================================================================
105*
106* .. Parameters ..
107 INTEGER LDT, LWORK
108 parameter( ldt = 20, lwork = 2*ldt*( 10+ldt ) )
109 REAL ZERO, ONE, TWO
110 parameter( zero = 0.0e+0, one = 1.0e+0, two = 2.0e+0 )
111 REAL EPSIN
112 parameter( epsin = 5.9605e-8 )
113 COMPLEX CZERO
114 parameter( czero = ( 0.0e+0, 0.0e+0 ) )
115* ..
116* .. Local Scalars ..
117 INTEGER I, INFO, ISCL, ISRT, ITMP, J, KMIN, M, N, NDIM
118 REAL BIGNUM, EPS, S, SEP, SEPIN, SEPTMP, SIN,
119 $ SMLNUM, STMP, TNRM, TOL, TOLIN, V, VMAX, VMIN,
120 $ VMUL
121* ..
122* .. Local Arrays ..
123 LOGICAL SELECT( LDT )
124 INTEGER IPNT( LDT ), ISELEC( LDT )
125 REAL RESULT( 2 ), RWORK( LDT ), VAL( 3 ),
126 $ WSRT( LDT )
127 COMPLEX Q( LDT, LDT ), QSAV( LDT, LDT ),
128 $ QTMP( LDT, LDT ), T( LDT, LDT ),
129 $ TMP( LDT, LDT ), TSAV( LDT, LDT ),
130 $ TSAV1( LDT, LDT ), TTMP( LDT, LDT ), W( LDT ),
131 $ WORK( LWORK ), WTMP( LDT )
132* ..
133* .. External Functions ..
134 REAL CLANGE, SLAMCH
135 EXTERNAL clange, slamch
136* ..
137* .. External Subroutines ..
138 EXTERNAL cgehrd, chseqr, chst01, clacpy, csscal, ctrsen,
139 $ cunghr
140* ..
141* .. Intrinsic Functions ..
142 INTRINSIC aimag, max, real, sqrt
143* ..
144* .. Executable Statements ..
145*
146 eps = slamch( 'P' )
147 smlnum = slamch( 'S' ) / eps
148 bignum = one / smlnum
149*
150* EPSIN = 2**(-24) = precision to which input data computed
151*
152 eps = max( eps, epsin )
153 rmax( 1 ) = zero
154 rmax( 2 ) = zero
155 rmax( 3 ) = zero
156 lmax( 1 ) = 0
157 lmax( 2 ) = 0
158 lmax( 3 ) = 0
159 knt = 0
160 ninfo( 1 ) = 0
161 ninfo( 2 ) = 0
162 ninfo( 3 ) = 0
163 val( 1 ) = sqrt( smlnum )
164 val( 2 ) = one
165 val( 3 ) = sqrt( sqrt( bignum ) )
166*
167* Read input data until N=0. Assume input eigenvalues are sorted
168* lexicographically (increasing by real part, then decreasing by
169* imaginary part)
170*
171 10 CONTINUE
172 READ( nin, fmt = * )n, ndim, isrt
173 IF( n.EQ.0 )
174 $ RETURN
175 READ( nin, fmt = * )( iselec( i ), i = 1, ndim )
176 DO 20 i = 1, n
177 READ( nin, fmt = * )( tmp( i, j ), j = 1, n )
178 20 CONTINUE
179 READ( nin, fmt = * )sin, sepin
180*
181 tnrm = clange( 'M', n, n, tmp, ldt, rwork )
182 DO 200 iscl = 1, 3
183*
184* Scale input matrix
185*
186 knt = knt + 1
187 CALL clacpy( 'F', n, n, tmp, ldt, t, ldt )
188 vmul = val( iscl )
189 DO 30 i = 1, n
190 CALL csscal( n, vmul, t( 1, i ), 1 )
191 30 CONTINUE
192 IF( tnrm.EQ.zero )
193 $ vmul = one
194 CALL clacpy( 'F', n, n, t, ldt, tsav, ldt )
195*
196* Compute Schur form
197*
198 CALL cgehrd( 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 200
204 END IF
205*
206* Generate unitary matrix
207*
208 CALL clacpy( 'L', n, n, t, ldt, q, ldt )
209 CALL cunghr( n, 1, n, q, ldt, work( 1 ), work( n+1 ), lwork-n,
210 $ info )
211*
212* Compute Schur form
213*
214 DO 50 j = 1, n - 2
215 DO 40 i = j + 2, n
216 t( i, j ) = czero
217 40 CONTINUE
218 50 CONTINUE
219 CALL chseqr( 'S', 'V', n, 1, n, t, ldt, w, q, ldt, work, lwork,
220 $ info )
221 IF( info.NE.0 ) THEN
222 lmax( 2 ) = knt
223 ninfo( 2 ) = ninfo( 2 ) + 1
224 GO TO 200
225 END IF
226*
227* Sort, select eigenvalues
228*
229 DO 60 i = 1, n
230 ipnt( i ) = i
231 SELECT( i ) = .false.
232 60 CONTINUE
233 IF( isrt.EQ.0 ) THEN
234 DO 70 i = 1, n
235 wsrt( i ) = real( w( i ) )
236 70 CONTINUE
237 ELSE
238 DO 80 i = 1, n
239 wsrt( i ) = aimag( w( i ) )
240 80 CONTINUE
241 END IF
242 DO 100 i = 1, n - 1
243 kmin = i
244 vmin = wsrt( i )
245 DO 90 j = i + 1, n
246 IF( wsrt( j ).LT.vmin ) THEN
247 kmin = j
248 vmin = wsrt( j )
249 END IF
250 90 CONTINUE
251 wsrt( kmin ) = wsrt( i )
252 wsrt( i ) = vmin
253 itmp = ipnt( i )
254 ipnt( i ) = ipnt( kmin )
255 ipnt( kmin ) = itmp
256 100 CONTINUE
257 DO 110 i = 1, ndim
258 SELECT( ipnt( iselec( i ) ) ) = .true.
259 110 CONTINUE
260*
261* Compute condition numbers
262*
263 CALL clacpy( 'F', n, n, q, ldt, qsav, ldt )
264 CALL clacpy( 'F', n, n, t, ldt, tsav1, ldt )
265 CALL ctrsen( 'B', 'V', SELECT, n, t, ldt, q, ldt, wtmp, m, s,
266 $ sep, work, lwork, info )
267 IF( info.NE.0 ) THEN
268 lmax( 3 ) = knt
269 ninfo( 3 ) = ninfo( 3 ) + 1
270 GO TO 200
271 END IF
272 septmp = sep / vmul
273 stmp = s
274*
275* Compute residuals
276*
277 CALL chst01( n, 1, n, tsav, ldt, t, ldt, q, ldt, work, lwork,
278 $ rwork, result )
279 vmax = max( result( 1 ), result( 2 ) )
280 IF( vmax.GT.rmax( 1 ) ) THEN
281 rmax( 1 ) = vmax
282 IF( ninfo( 1 ).EQ.0 )
283 $ lmax( 1 ) = knt
284 END IF
285*
286* Compare condition number for eigenvalue cluster
287* taking its condition number into account
288*
289 v = max( two*real( n )*eps*tnrm, smlnum )
290 IF( tnrm.EQ.zero )
291 $ v = one
292 IF( v.GT.septmp ) THEN
293 tol = one
294 ELSE
295 tol = v / septmp
296 END IF
297 IF( v.GT.sepin ) THEN
298 tolin = one
299 ELSE
300 tolin = v / sepin
301 END IF
302 tol = max( tol, smlnum / eps )
303 tolin = max( tolin, smlnum / eps )
304 IF( eps*( sin-tolin ).GT.stmp+tol ) THEN
305 vmax = one / eps
306 ELSE IF( sin-tolin.GT.stmp+tol ) THEN
307 vmax = ( sin-tolin ) / ( stmp+tol )
308 ELSE IF( sin+tolin.LT.eps*( stmp-tol ) ) THEN
309 vmax = one / eps
310 ELSE IF( sin+tolin.LT.stmp-tol ) THEN
311 vmax = ( stmp-tol ) / ( sin+tolin )
312 ELSE
313 vmax = one
314 END IF
315 IF( vmax.GT.rmax( 2 ) ) THEN
316 rmax( 2 ) = vmax
317 IF( ninfo( 2 ).EQ.0 )
318 $ lmax( 2 ) = knt
319 END IF
320*
321* Compare condition numbers for invariant subspace
322* taking its condition number into account
323*
324 IF( v.GT.septmp*stmp ) THEN
325 tol = septmp
326 ELSE
327 tol = v / stmp
328 END IF
329 IF( v.GT.sepin*sin ) THEN
330 tolin = sepin
331 ELSE
332 tolin = v / sin
333 END IF
334 tol = max( tol, smlnum / eps )
335 tolin = max( tolin, smlnum / eps )
336 IF( eps*( sepin-tolin ).GT.septmp+tol ) THEN
337 vmax = one / eps
338 ELSE IF( sepin-tolin.GT.septmp+tol ) THEN
339 vmax = ( sepin-tolin ) / ( septmp+tol )
340 ELSE IF( sepin+tolin.LT.eps*( septmp-tol ) ) THEN
341 vmax = one / eps
342 ELSE IF( sepin+tolin.LT.septmp-tol ) THEN
343 vmax = ( septmp-tol ) / ( sepin+tolin )
344 ELSE
345 vmax = one
346 END IF
347 IF( vmax.GT.rmax( 2 ) ) THEN
348 rmax( 2 ) = vmax
349 IF( ninfo( 2 ).EQ.0 )
350 $ lmax( 2 ) = knt
351 END IF
352*
353* Compare condition number for eigenvalue cluster
354* without taking its condition number into account
355*
356 IF( sin.LE.real( 2*n )*eps .AND. stmp.LE.real( 2*n )*eps ) THEN
357 vmax = one
358 ELSE IF( eps*sin.GT.stmp ) THEN
359 vmax = one / eps
360 ELSE IF( sin.GT.stmp ) THEN
361 vmax = sin / stmp
362 ELSE IF( sin.LT.eps*stmp ) THEN
363 vmax = one / eps
364 ELSE IF( sin.LT.stmp ) THEN
365 vmax = stmp / sin
366 ELSE
367 vmax = one
368 END IF
369 IF( vmax.GT.rmax( 3 ) ) THEN
370 rmax( 3 ) = vmax
371 IF( ninfo( 3 ).EQ.0 )
372 $ lmax( 3 ) = knt
373 END IF
374*
375* Compare condition numbers for invariant subspace
376* without taking its condition number into account
377*
378 IF( sepin.LE.v .AND. septmp.LE.v ) THEN
379 vmax = one
380 ELSE IF( eps*sepin.GT.septmp ) THEN
381 vmax = one / eps
382 ELSE IF( sepin.GT.septmp ) THEN
383 vmax = sepin / septmp
384 ELSE IF( sepin.LT.eps*septmp ) THEN
385 vmax = one / eps
386 ELSE IF( sepin.LT.septmp ) THEN
387 vmax = septmp / sepin
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*
397* Compute eigenvalue condition number only and compare
398* Update Q
399*
400 vmax = zero
401 CALL clacpy( 'F', n, n, tsav1, ldt, ttmp, ldt )
402 CALL clacpy( 'F', n, n, qsav, ldt, qtmp, ldt )
403 septmp = -one
404 stmp = -one
405 CALL ctrsen( 'E', 'V', SELECT, n, ttmp, ldt, qtmp, ldt, wtmp,
406 $ m, stmp, septmp, work, lwork, info )
407 IF( info.NE.0 ) THEN
408 lmax( 3 ) = knt
409 ninfo( 3 ) = ninfo( 3 ) + 1
410 GO TO 200
411 END IF
412 IF( s.NE.stmp )
413 $ vmax = one / eps
414 IF( -one.NE.septmp )
415 $ vmax = one / eps
416 DO 130 i = 1, n
417 DO 120 j = 1, n
418 IF( ttmp( i, j ).NE.t( i, j ) )
419 $ vmax = one / eps
420 IF( qtmp( i, j ).NE.q( i, j ) )
421 $ vmax = one / eps
422 120 CONTINUE
423 130 CONTINUE
424*
425* Compute invariant subspace condition number only and compare
426* Update Q
427*
428 CALL clacpy( 'F', n, n, tsav1, ldt, ttmp, ldt )
429 CALL clacpy( 'F', n, n, qsav, ldt, qtmp, ldt )
430 septmp = -one
431 stmp = -one
432 CALL ctrsen( 'V', 'V', SELECT, n, ttmp, ldt, qtmp, ldt, wtmp,
433 $ m, stmp, septmp, work, lwork, info )
434 IF( info.NE.0 ) THEN
435 lmax( 3 ) = knt
436 ninfo( 3 ) = ninfo( 3 ) + 1
437 GO TO 200
438 END IF
439 IF( -one.NE.stmp )
440 $ vmax = one / eps
441 IF( sep.NE.septmp )
442 $ vmax = one / eps
443 DO 150 i = 1, n
444 DO 140 j = 1, n
445 IF( ttmp( i, j ).NE.t( i, j ) )
446 $ vmax = one / eps
447 IF( qtmp( i, j ).NE.q( i, j ) )
448 $ vmax = one / eps
449 140 CONTINUE
450 150 CONTINUE
451*
452* Compute eigenvalue condition number only and compare
453* Do not update Q
454*
455 CALL clacpy( 'F', n, n, tsav1, ldt, ttmp, ldt )
456 CALL clacpy( 'F', n, n, qsav, ldt, qtmp, ldt )
457 septmp = -one
458 stmp = -one
459 CALL ctrsen( 'E', 'N', SELECT, n, ttmp, ldt, qtmp, ldt, wtmp,
460 $ m, stmp, septmp, work, lwork, info )
461 IF( info.NE.0 ) THEN
462 lmax( 3 ) = knt
463 ninfo( 3 ) = ninfo( 3 ) + 1
464 GO TO 200
465 END IF
466 IF( s.NE.stmp )
467 $ vmax = one / eps
468 IF( -one.NE.septmp )
469 $ vmax = one / eps
470 DO 170 i = 1, n
471 DO 160 j = 1, n
472 IF( ttmp( i, j ).NE.t( i, j ) )
473 $ vmax = one / eps
474 IF( qtmp( i, j ).NE.qsav( i, j ) )
475 $ vmax = one / eps
476 160 CONTINUE
477 170 CONTINUE
478*
479* Compute invariant subspace condition number only and compare
480* Do not update Q
481*
482 CALL clacpy( 'F', n, n, tsav1, ldt, ttmp, ldt )
483 CALL clacpy( 'F', n, n, qsav, ldt, qtmp, ldt )
484 septmp = -one
485 stmp = -one
486 CALL ctrsen( 'V', 'N', SELECT, n, ttmp, ldt, qtmp, ldt, wtmp,
487 $ m, stmp, septmp, work, lwork, info )
488 IF( info.NE.0 ) THEN
489 lmax( 3 ) = knt
490 ninfo( 3 ) = ninfo( 3 ) + 1
491 GO TO 200
492 END IF
493 IF( -one.NE.stmp )
494 $ vmax = one / eps
495 IF( sep.NE.septmp )
496 $ vmax = one / eps
497 DO 190 i = 1, n
498 DO 180 j = 1, n
499 IF( ttmp( i, j ).NE.t( i, j ) )
500 $ vmax = one / eps
501 IF( qtmp( i, j ).NE.qsav( i, j ) )
502 $ vmax = one / eps
503 180 CONTINUE
504 190 CONTINUE
505 IF( vmax.GT.rmax( 1 ) ) THEN
506 rmax( 1 ) = vmax
507 IF( ninfo( 1 ).EQ.0 )
508 $ lmax( 1 ) = knt
509 END IF
510 200 CONTINUE
511 GO TO 10
512*
513* End of CGET38
514*
chst01
subroutine chst01(n, ilo, ihi, a, lda, h, ldh, q, ldq, work, lwork, rwork, result)
CHST01
Definition chst01.f:140
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
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
clange
real function clange(norm, m, n, a, lda, work)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition clange.f:115
csscal
subroutine csscal(n, sa, cx, incx)
CSSCAL
Definition csscal.f:78
ctrsen
subroutine ctrsen(job, compq, select, n, t, ldt, q, ldq, w, m, s, sep, work, lwork, info)
CTRSEN
Definition ctrsen.f:264
cunghr
subroutine cunghr(n, ilo, ihi, a, lda, tau, work, lwork, info)
CUNGHR
Definition cunghr.f:126
Here is the call graph for this function:
Here is the caller graph for this function: