OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
ctrsna.f
Go to the documentation of this file.
1*> \brief \b CTRSNA
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CTRSNA + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctrsna.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctrsna.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctrsna.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
22* LDVR, S, SEP, MM, M, WORK, LDWORK, RWORK,
23* INFO )
24*
25* .. Scalar Arguments ..
26* CHARACTER HOWMNY, JOB
27* INTEGER INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
28* ..
29* .. Array Arguments ..
30* LOGICAL SELECT( * )
31* REAL RWORK( * ), S( * ), SEP( * )
32* COMPLEX T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
33* $ WORK( LDWORK, * )
34* ..
35*
36*
37*> \par Purpose:
38* =============
39*>
40*> \verbatim
41*>
42*> CTRSNA estimates reciprocal condition numbers for specified
43*> eigenvalues and/or right eigenvectors of a complex upper triangular
44*> matrix T (or of any matrix Q*T*Q**H with Q unitary).
45*> \endverbatim
46*
47* Arguments:
48* ==========
49*
50*> \param[in] JOB
51*> \verbatim
52*> JOB is CHARACTER*1
53*> Specifies whether condition numbers are required for
54*> eigenvalues (S) or eigenvectors (SEP):
55*> = 'E': for eigenvalues only (S);
56*> = 'V': for eigenvectors only (SEP);
57*> = 'B': for both eigenvalues and eigenvectors (S and SEP).
58*> \endverbatim
59*>
60*> \param[in] HOWMNY
61*> \verbatim
62*> HOWMNY is CHARACTER*1
63*> = 'A': compute condition numbers for all eigenpairs;
64*> = 'S': compute condition numbers for selected eigenpairs
65*> specified by the array SELECT.
66*> \endverbatim
67*>
68*> \param[in] SELECT
69*> \verbatim
70*> SELECT is LOGICAL array, dimension (N)
71*> If HOWMNY = 'S', SELECT specifies the eigenpairs for which
72*> condition numbers are required. To select condition numbers
73*> for the j-th eigenpair, SELECT(j) must be set to .TRUE..
74*> If HOWMNY = 'A', SELECT is not referenced.
75*> \endverbatim
76*>
77*> \param[in] N
78*> \verbatim
79*> N is INTEGER
80*> The order of the matrix T. N >= 0.
81*> \endverbatim
82*>
83*> \param[in] T
84*> \verbatim
85*> T is COMPLEX array, dimension (LDT,N)
86*> The upper triangular matrix T.
87*> \endverbatim
88*>
89*> \param[in] LDT
90*> \verbatim
91*> LDT is INTEGER
92*> The leading dimension of the array T. LDT >= max(1,N).
93*> \endverbatim
94*>
95*> \param[in] VL
96*> \verbatim
97*> VL is COMPLEX array, dimension (LDVL,M)
98*> If JOB = 'E' or 'B', VL must contain left eigenvectors of T
99*> (or of any Q*T*Q**H with Q unitary), corresponding to the
100*> eigenpairs specified by HOWMNY and SELECT. The eigenvectors
101*> must be stored in consecutive columns of VL, as returned by
102*> CHSEIN or CTREVC.
103*> If JOB = 'V', VL is not referenced.
104*> \endverbatim
105*>
106*> \param[in] LDVL
107*> \verbatim
108*> LDVL is INTEGER
109*> The leading dimension of the array VL.
110*> LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N.
111*> \endverbatim
112*>
113*> \param[in] VR
114*> \verbatim
115*> VR is COMPLEX array, dimension (LDVR,M)
116*> If JOB = 'E' or 'B', VR must contain right eigenvectors of T
117*> (or of any Q*T*Q**H with Q unitary), corresponding to the
118*> eigenpairs specified by HOWMNY and SELECT. The eigenvectors
119*> must be stored in consecutive columns of VR, as returned by
120*> CHSEIN or CTREVC.
121*> If JOB = 'V', VR is not referenced.
122*> \endverbatim
123*>
124*> \param[in] LDVR
125*> \verbatim
126*> LDVR is INTEGER
127*> The leading dimension of the array VR.
128*> LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N.
129*> \endverbatim
130*>
131*> \param[out] S
132*> \verbatim
133*> S is REAL array, dimension (MM)
134*> If JOB = 'E' or 'B', the reciprocal condition numbers of the
135*> selected eigenvalues, stored in consecutive elements of the
136*> array. Thus S(j), SEP(j), and the j-th columns of VL and VR
137*> all correspond to the same eigenpair (but not in general the
138*> j-th eigenpair, unless all eigenpairs are selected).
139*> If JOB = 'V', S is not referenced.
140*> \endverbatim
141*>
142*> \param[out] SEP
143*> \verbatim
144*> SEP is REAL array, dimension (MM)
145*> If JOB = 'V' or 'B', the estimated reciprocal condition
146*> numbers of the selected eigenvectors, stored in consecutive
147*> elements of the array.
148*> If JOB = 'E', SEP is not referenced.
149*> \endverbatim
150*>
151*> \param[in] MM
152*> \verbatim
153*> MM is INTEGER
154*> The number of elements in the arrays S (if JOB = 'E' or 'B')
155*> and/or SEP (if JOB = 'V' or 'B'). MM >= M.
156*> \endverbatim
157*>
158*> \param[out] M
159*> \verbatim
160*> M is INTEGER
161*> The number of elements of the arrays S and/or SEP actually
162*> used to store the estimated condition numbers.
163*> If HOWMNY = 'A', M is set to N.
164*> \endverbatim
165*>
166*> \param[out] WORK
167*> \verbatim
168*> WORK is COMPLEX array, dimension (LDWORK,N+6)
169*> If JOB = 'E', WORK is not referenced.
170*> \endverbatim
171*>
172*> \param[in] LDWORK
173*> \verbatim
174*> LDWORK is INTEGER
175*> The leading dimension of the array WORK.
176*> LDWORK >= 1; and if JOB = 'V' or 'B', LDWORK >= N.
177*> \endverbatim
178*>
179*> \param[out] RWORK
180*> \verbatim
181*> RWORK is REAL array, dimension (N)
182*> If JOB = 'E', RWORK is not referenced.
183*> \endverbatim
184*>
185*> \param[out] INFO
186*> \verbatim
187*> INFO is INTEGER
188*> = 0: successful exit
189*> < 0: if INFO = -i, the i-th argument had an illegal value
190*> \endverbatim
191*
192* Authors:
193* ========
194*
195*> \author Univ. of Tennessee
196*> \author Univ. of California Berkeley
197*> \author Univ. of Colorado Denver
198*> \author NAG Ltd.
199*
200*> \ingroup complexOTHERcomputational
201*
202*> \par Further Details:
203* =====================
204*>
205*> \verbatim
206*>
207*> The reciprocal of the condition number of an eigenvalue lambda is
208*> defined as
209*>
210*> S(lambda) = |v**H*u| / (norm(u)*norm(v))
211*>
212*> where u and v are the right and left eigenvectors of T corresponding
213*> to lambda; v**H denotes the conjugate transpose of v, and norm(u)
214*> denotes the Euclidean norm. These reciprocal condition numbers always
215*> lie between zero (very badly conditioned) and one (very well
216*> conditioned). If n = 1, S(lambda) is defined to be 1.
217*>
218*> An approximate error bound for a computed eigenvalue W(i) is given by
219*>
220*> EPS * norm(T) / S(i)
221*>
222*> where EPS is the machine precision.
223*>
224*> The reciprocal of the condition number of the right eigenvector u
225*> corresponding to lambda is defined as follows. Suppose
226*>
227*> T = ( lambda c )
228*> ( 0 T22 )
229*>
230*> Then the reciprocal condition number is
231*>
232*> SEP( lambda, T22 ) = sigma-min( T22 - lambda*I )
233*>
234*> where sigma-min denotes the smallest singular value. We approximate
235*> the smallest singular value by the reciprocal of an estimate of the
236*> one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is
237*> defined to be abs(T(1,1)).
238*>
239*> An approximate error bound for a computed right eigenvector VR(i)
240*> is given by
241*>
242*> EPS * norm(T) / SEP(i)
243*> \endverbatim
244*>
245* =====================================================================
246 SUBROUTINE ctrsna( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
247 $ LDVR, S, SEP, MM, M, WORK, LDWORK, RWORK,
248 $ INFO )
249*
250* -- LAPACK computational routine --
251* -- LAPACK is a software package provided by Univ. of Tennessee, --
252* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
253*
254* .. Scalar Arguments ..
255 CHARACTER HOWMNY, JOB
256 INTEGER INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
257* ..
258* .. Array Arguments ..
259 LOGICAL SELECT( * )
260 REAL RWORK( * ), S( * ), SEP( * )
261 COMPLEX T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
262 $ work( ldwork, * )
263* ..
264*
265* =====================================================================
266*
267* .. Parameters ..
268 REAL ZERO, ONE
269 PARAMETER ( ZERO = 0.0e+0, one = 1.0+0 )
270* ..
271* .. Local Scalars ..
272 LOGICAL SOMCON, WANTBH, WANTS, WANTSP
273 CHARACTER NORMIN
274 INTEGER I, IERR, IX, J, K, KASE, KS
275 REAL BIGNUM, EPS, EST, LNRM, RNRM, SCALE, SMLNUM,
276 $ xnorm
277 COMPLEX CDUM, PROD
278* ..
279* .. Local Arrays ..
280 INTEGER ISAVE( 3 )
281 COMPLEX DUMMY( 1 )
282* ..
283* .. External Functions ..
284 LOGICAL LSAME
285 INTEGER ICAMAX
286 REAL SCNRM2, SLAMCH
287 COMPLEX CDOTC
288 EXTERNAL lsame, icamax, scnrm2, slamch, cdotc
289* ..
290* .. External Subroutines ..
291 EXTERNAL clacn2, clacpy, clatrs, csrscl, ctrexc, slabad,
292 $ xerbla
293* ..
294* .. Intrinsic Functions ..
295 INTRINSIC abs, aimag, max, real
296* ..
297* .. Statement Functions ..
298 REAL CABS1
299* ..
300* .. Statement Function definitions ..
301 cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( cdum ) )
302* ..
303* .. Executable Statements ..
304*
305* Decode and test the input parameters
306*
307 wantbh = lsame( job, 'B' )
308 wants = lsame( job, 'E' ) .OR. wantbh
309 wantsp = lsame( job, 'V' ) .OR. wantbh
310*
311 somcon = lsame( howmny, 'S' )
312*
313* Set M to the number of eigenpairs for which condition numbers are
314* to be computed.
315*
316 IF( somcon ) THEN
317 m = 0
318 DO 10 j = 1, n
319 IF( SELECT( j ) )
320 $ m = m + 1
321 10 CONTINUE
322 ELSE
323 m = n
324 END IF
325*
326 info = 0
327 IF( .NOT.wants .AND. .NOT.wantsp ) THEN
328 info = -1
329 ELSE IF( .NOT.lsame( howmny, 'A' ) .AND. .NOT.somcon ) THEN
330 info = -2
331 ELSE IF( n.LT.0 ) THEN
332 info = -4
333 ELSE IF( ldt.LT.max( 1, n ) ) THEN
334 info = -6
335 ELSE IF( ldvl.LT.1 .OR. ( wants .AND. ldvl.LT.n ) ) THEN
336 info = -8
337 ELSE IF( ldvr.LT.1 .OR. ( wants .AND. ldvr.LT.n ) ) THEN
338 info = -10
339 ELSE IF( mm.LT.m ) THEN
340 info = -13
341 ELSE IF( ldwork.LT.1 .OR. ( wantsp .AND. ldwork.LT.n ) ) THEN
342 info = -16
343 END IF
344 IF( info.NE.0 ) THEN
345 CALL xerbla( 'CTRSNA', -info )
346 RETURN
347 END IF
348*
349* Quick return if possible
350*
351 IF( n.EQ.0 )
352 $ RETURN
353*
354 IF( n.EQ.1 ) THEN
355 IF( somcon ) THEN
356 IF( .NOT.SELECT( 1 ) )
357 $ RETURN
358 END IF
359 IF( wants )
360 $ s( 1 ) = one
361 IF( wantsp )
362 $ sep( 1 ) = abs( t( 1, 1 ) )
363 RETURN
364 END IF
365*
366* Get machine constants
367*
368 eps = slamch( 'P' )
369 smlnum = slamch( 'S' ) / eps
370 bignum = one / smlnum
371 CALL slabad( smlnum, bignum )
372*
373 ks = 1
374 DO 50 k = 1, n
375*
376 IF( somcon ) THEN
377 IF( .NOT.SELECT( k ) )
378 $ GO TO 50
379 END IF
380*
381 IF( wants ) THEN
382*
383* Compute the reciprocal condition number of the k-th
384* eigenvalue.
385*
386 prod = cdotc( n, vr( 1, ks ), 1, vl( 1, ks ), 1 )
387 rnrm = scnrm2( n, vr( 1, ks ), 1 )
388 lnrm = scnrm2( n, vl( 1, ks ), 1 )
389 s( ks ) = abs( prod ) / ( rnrm*lnrm )
390*
391 END IF
392*
393 IF( wantsp ) THEN
394*
395* Estimate the reciprocal condition number of the k-th
396* eigenvector.
397*
398* Copy the matrix T to the array WORK and swap the k-th
399* diagonal element to the (1,1) position.
400*
401 CALL clacpy( 'Full', n, n, t, ldt, work, ldwork )
402 CALL ctrexc( 'No Q', n, work, ldwork, dummy, 1, k, 1, ierr )
403*
404* Form C = T22 - lambda*I in WORK(2:N,2:N).
405*
406 DO 20 i = 2, n
407 work( i, i ) = work( i, i ) - work( 1, 1 )
408 20 CONTINUE
409*
410* Estimate a lower bound for the 1-norm of inv(C**H). The 1st
411* and (N+1)th columns of WORK are used to store work vectors.
412*
413 sep( ks ) = zero
414 est = zero
415 kase = 0
416 normin = 'N'
417 30 CONTINUE
418 CALL clacn2( n-1, work( 1, n+1 ), work, est, kase, isave )
419*
420 IF( kase.NE.0 ) THEN
421 IF( kase.EQ.1 ) THEN
422*
423* Solve C**H*x = scale*b
424*
425 CALL clatrs( 'Upper', 'Conjugate transpose',
426 $ 'Nonunit', normin, n-1, work( 2, 2 ),
427 $ ldwork, work, scale, rwork, ierr )
428 ELSE
429*
430* Solve C*x = scale*b
431*
432 CALL clatrs( 'Upper', 'No transpose', 'Nonunit',
433 $ normin, n-1, work( 2, 2 ), ldwork, work,
434 $ scale, rwork, ierr )
435 END IF
436 normin = 'Y'
437 IF( scale.NE.one ) THEN
438*
439* Multiply by 1/SCALE if doing so will not cause
440* overflow.
441*
442 ix = icamax( n-1, work, 1 )
443 xnorm = cabs1( work( ix, 1 ) )
444 IF( scale.LT.xnorm*smlnum .OR. scale.EQ.zero )
445 $ GO TO 40
446 CALL csrscl( n, scale, work, 1 )
447 END IF
448 GO TO 30
449 END IF
450*
451 sep( ks ) = one / max( est, smlnum )
452 END IF
453*
454 40 CONTINUE
455 ks = ks + 1
456 50 CONTINUE
457 RETURN
458*
459* End of CTRSNA
460*
461 END
slabad
subroutine slabad(small, large)
SLABAD
Definition slabad.f:74
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
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
csrscl
subroutine csrscl(n, sa, sx, incx)
CSRSCL multiplies a vector by the reciprocal of a real scalar.
Definition csrscl.f:84
clatrs
subroutine clatrs(uplo, trans, diag, normin, n, a, lda, x, scale, cnorm, info)
CLATRS solves a triangular system of equations with the scale factor set to prevent overflow.
Definition clatrs.f:239
clacn2
subroutine clacn2(n, v, x, est, kase, isave)
CLACN2 estimates the 1-norm of a square matrix, using reverse communication for evaluating matrix-vec...
Definition clacn2.f:133
ctrexc
subroutine ctrexc(compq, n, t, ldt, q, ldq, ifst, ilst, info)
CTREXC
Definition ctrexc.f:126
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
max
#define max(a, b)
Definition macros.h:21