001:       SUBROUTINE DTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
002:      $                   LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK,
003:      $                   INFO )
004: *
005: *  -- LAPACK routine (version 3.2) --
006: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
007: *     November 2006
008: *
009: *     Modified to call DLACN2 in place of DLACON, 5 Feb 03, SJH.
010: *
011: *     .. Scalar Arguments ..
012:       CHARACTER          HOWMNY, JOB
013:       INTEGER            INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
014: *     ..
015: *     .. Array Arguments ..
016:       LOGICAL            SELECT( * )
017:       INTEGER            IWORK( * )
018:       DOUBLE PRECISION   S( * ), SEP( * ), T( LDT, * ), VL( LDVL, * ),
019:      $                   VR( LDVR, * ), WORK( LDWORK, * )
020: *     ..
021: *
022: *  Purpose
023: *  =======
024: *
025: *  DTRSNA estimates reciprocal condition numbers for specified
026: *  eigenvalues and/or right eigenvectors of a real upper
027: *  quasi-triangular matrix T (or of any matrix Q*T*Q**T with Q
028: *  orthogonal).
029: *
030: *  T must be in Schur canonical form (as returned by DHSEQR), that is,
031: *  block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
032: *  2-by-2 diagonal block has its diagonal elements equal and its
033: *  off-diagonal elements of opposite sign.
034: *
035: *  Arguments
036: *  =========
037: *
038: *  JOB     (input) CHARACTER*1
039: *          Specifies whether condition numbers are required for
040: *          eigenvalues (S) or eigenvectors (SEP):
041: *          = 'E': for eigenvalues only (S);
042: *          = 'V': for eigenvectors only (SEP);
043: *          = 'B': for both eigenvalues and eigenvectors (S and SEP).
044: *
045: *  HOWMNY  (input) CHARACTER*1
046: *          = 'A': compute condition numbers for all eigenpairs;
047: *          = 'S': compute condition numbers for selected eigenpairs
048: *                 specified by the array SELECT.
049: *
050: *  SELECT  (input) LOGICAL array, dimension (N)
051: *          If HOWMNY = 'S', SELECT specifies the eigenpairs for which
052: *          condition numbers are required. To select condition numbers
053: *          for the eigenpair corresponding to a real eigenvalue w(j),
054: *          SELECT(j) must be set to .TRUE.. To select condition numbers
055: *          corresponding to a complex conjugate pair of eigenvalues w(j)
056: *          and w(j+1), either SELECT(j) or SELECT(j+1) or both, must be
057: *          set to .TRUE..
058: *          If HOWMNY = 'A', SELECT is not referenced.
059: *
060: *  N       (input) INTEGER
061: *          The order of the matrix T. N >= 0.
062: *
063: *  T       (input) DOUBLE PRECISION array, dimension (LDT,N)
064: *          The upper quasi-triangular matrix T, in Schur canonical form.
065: *
066: *  LDT     (input) INTEGER
067: *          The leading dimension of the array T. LDT >= max(1,N).
068: *
069: *  VL      (input) DOUBLE PRECISION array, dimension (LDVL,M)
070: *          If JOB = 'E' or 'B', VL must contain left eigenvectors of T
071: *          (or of any Q*T*Q**T with Q orthogonal), corresponding to the
072: *          eigenpairs specified by HOWMNY and SELECT. The eigenvectors
073: *          must be stored in consecutive columns of VL, as returned by
074: *          DHSEIN or DTREVC.
075: *          If JOB = 'V', VL is not referenced.
076: *
077: *  LDVL    (input) INTEGER
078: *          The leading dimension of the array VL.
079: *          LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N.
080: *
081: *  VR      (input) DOUBLE PRECISION array, dimension (LDVR,M)
082: *          If JOB = 'E' or 'B', VR must contain right eigenvectors of T
083: *          (or of any Q*T*Q**T with Q orthogonal), corresponding to the
084: *          eigenpairs specified by HOWMNY and SELECT. The eigenvectors
085: *          must be stored in consecutive columns of VR, as returned by
086: *          DHSEIN or DTREVC.
087: *          If JOB = 'V', VR is not referenced.
088: *
089: *  LDVR    (input) INTEGER
090: *          The leading dimension of the array VR.
091: *          LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N.
092: *
093: *  S       (output) DOUBLE PRECISION array, dimension (MM)
094: *          If JOB = 'E' or 'B', the reciprocal condition numbers of the
095: *          selected eigenvalues, stored in consecutive elements of the
096: *          array. For a complex conjugate pair of eigenvalues two
097: *          consecutive elements of S are set to the same value. Thus
098: *          S(j), SEP(j), and the j-th columns of VL and VR all
099: *          correspond to the same eigenpair (but not in general the
100: *          j-th eigenpair, unless all eigenpairs are selected).
101: *          If JOB = 'V', S is not referenced.
102: *
103: *  SEP     (output) DOUBLE PRECISION array, dimension (MM)
104: *          If JOB = 'V' or 'B', the estimated reciprocal condition
105: *          numbers of the selected eigenvectors, stored in consecutive
106: *          elements of the array. For a complex eigenvector two
107: *          consecutive elements of SEP are set to the same value. If
108: *          the eigenvalues cannot be reordered to compute SEP(j), SEP(j)
109: *          is set to 0; this can only occur when the true value would be
110: *          very small anyway.
111: *          If JOB = 'E', SEP is not referenced.
112: *
113: *  MM      (input) INTEGER
114: *          The number of elements in the arrays S (if JOB = 'E' or 'B')
115: *           and/or SEP (if JOB = 'V' or 'B'). MM >= M.
116: *
117: *  M       (output) INTEGER
118: *          The number of elements of the arrays S and/or SEP actually
119: *          used to store the estimated condition numbers.
120: *          If HOWMNY = 'A', M is set to N.
121: *
122: *  WORK    (workspace) DOUBLE PRECISION array, dimension (LDWORK,N+6)
123: *          If JOB = 'E', WORK is not referenced.
124: *
125: *  LDWORK  (input) INTEGER
126: *          The leading dimension of the array WORK.
127: *          LDWORK >= 1; and if JOB = 'V' or 'B', LDWORK >= N.
128: *
129: *  IWORK   (workspace) INTEGER array, dimension (2*(N-1))
130: *          If JOB = 'E', IWORK is not referenced.
131: *
132: *  INFO    (output) INTEGER
133: *          = 0: successful exit
134: *          < 0: if INFO = -i, the i-th argument had an illegal value
135: *
136: *  Further Details
137: *  ===============
138: *
139: *  The reciprocal of the condition number of an eigenvalue lambda is
140: *  defined as
141: *
142: *          S(lambda) = |v'*u| / (norm(u)*norm(v))
143: *
144: *  where u and v are the right and left eigenvectors of T corresponding
145: *  to lambda; v' denotes the conjugate-transpose of v, and norm(u)
146: *  denotes the Euclidean norm. These reciprocal condition numbers always
147: *  lie between zero (very badly conditioned) and one (very well
148: *  conditioned). If n = 1, S(lambda) is defined to be 1.
149: *
150: *  An approximate error bound for a computed eigenvalue W(i) is given by
151: *
152: *                      EPS * norm(T) / S(i)
153: *
154: *  where EPS is the machine precision.
155: *
156: *  The reciprocal of the condition number of the right eigenvector u
157: *  corresponding to lambda is defined as follows. Suppose
158: *
159: *              T = ( lambda  c  )
160: *                  (   0    T22 )
161: *
162: *  Then the reciprocal condition number is
163: *
164: *          SEP( lambda, T22 ) = sigma-min( T22 - lambda*I )
165: *
166: *  where sigma-min denotes the smallest singular value. We approximate
167: *  the smallest singular value by the reciprocal of an estimate of the
168: *  one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is
169: *  defined to be abs(T(1,1)).
170: *
171: *  An approximate error bound for a computed right eigenvector VR(i)
172: *  is given by
173: *
174: *                      EPS * norm(T) / SEP(i)
175: *
176: *  =====================================================================
177: *
178: *     .. Parameters ..
179:       DOUBLE PRECISION   ZERO, ONE, TWO
180:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
181: *     ..
182: *     .. Local Scalars ..
183:       LOGICAL            PAIR, SOMCON, WANTBH, WANTS, WANTSP
184:       INTEGER            I, IERR, IFST, ILST, J, K, KASE, KS, N2, NN
185:       DOUBLE PRECISION   BIGNUM, COND, CS, DELTA, DUMM, EPS, EST, LNRM,
186:      $                   MU, PROD, PROD1, PROD2, RNRM, SCALE, SMLNUM, SN
187: *     ..
188: *     .. Local Arrays ..
189:       INTEGER            ISAVE( 3 )
190:       DOUBLE PRECISION   DUMMY( 1 )
191: *     ..
192: *     .. External Functions ..
193:       LOGICAL            LSAME
194:       DOUBLE PRECISION   DDOT, DLAMCH, DLAPY2, DNRM2
195:       EXTERNAL           LSAME, DDOT, DLAMCH, DLAPY2, DNRM2
196: *     ..
197: *     .. External Subroutines ..
198:       EXTERNAL           DLACN2, DLACPY, DLAQTR, DTREXC, XERBLA
199: *     ..
200: *     .. Intrinsic Functions ..
201:       INTRINSIC          ABS, MAX, SQRT
202: *     ..
203: *     .. Executable Statements ..
204: *
205: *     Decode and test the input parameters
206: *
207:       WANTBH = LSAME( JOB, 'B' )
208:       WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
209:       WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH
210: *
211:       SOMCON = LSAME( HOWMNY, 'S' )
212: *
213:       INFO = 0
214:       IF( .NOT.WANTS .AND. .NOT.WANTSP ) THEN
215:          INFO = -1
216:       ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
217:          INFO = -2
218:       ELSE IF( N.LT.0 ) THEN
219:          INFO = -4
220:       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
221:          INFO = -6
222:       ELSE IF( LDVL.LT.1 .OR. ( WANTS .AND. LDVL.LT.N ) ) THEN
223:          INFO = -8
224:       ELSE IF( LDVR.LT.1 .OR. ( WANTS .AND. LDVR.LT.N ) ) THEN
225:          INFO = -10
226:       ELSE
227: *
228: *        Set M to the number of eigenpairs for which condition numbers
229: *        are required, and test MM.
230: *
231:          IF( SOMCON ) THEN
232:             M = 0
233:             PAIR = .FALSE.
234:             DO 10 K = 1, N
235:                IF( PAIR ) THEN
236:                   PAIR = .FALSE.
237:                ELSE
238:                   IF( K.LT.N ) THEN
239:                      IF( T( K+1, K ).EQ.ZERO ) THEN
240:                         IF( SELECT( K ) )
241:      $                     M = M + 1
242:                      ELSE
243:                         PAIR = .TRUE.
244:                         IF( SELECT( K ) .OR. SELECT( K+1 ) )
245:      $                     M = M + 2
246:                      END IF
247:                   ELSE
248:                      IF( SELECT( N ) )
249:      $                  M = M + 1
250:                   END IF
251:                END IF
252:    10       CONTINUE
253:          ELSE
254:             M = N
255:          END IF
256: *
257:          IF( MM.LT.M ) THEN
258:             INFO = -13
259:          ELSE IF( LDWORK.LT.1 .OR. ( WANTSP .AND. LDWORK.LT.N ) ) THEN
260:             INFO = -16
261:          END IF
262:       END IF
263:       IF( INFO.NE.0 ) THEN
264:          CALL XERBLA( 'DTRSNA', -INFO )
265:          RETURN
266:       END IF
267: *
268: *     Quick return if possible
269: *
270:       IF( N.EQ.0 )
271:      $   RETURN
272: *
273:       IF( N.EQ.1 ) THEN
274:          IF( SOMCON ) THEN
275:             IF( .NOT.SELECT( 1 ) )
276:      $         RETURN
277:          END IF
278:          IF( WANTS )
279:      $      S( 1 ) = ONE
280:          IF( WANTSP )
281:      $      SEP( 1 ) = ABS( T( 1, 1 ) )
282:          RETURN
283:       END IF
284: *
285: *     Get machine constants
286: *
287:       EPS = DLAMCH( 'P' )
288:       SMLNUM = DLAMCH( 'S' ) / EPS
289:       BIGNUM = ONE / SMLNUM
290:       CALL DLABAD( SMLNUM, BIGNUM )
291: *
292:       KS = 0
293:       PAIR = .FALSE.
294:       DO 60 K = 1, N
295: *
296: *        Determine whether T(k,k) begins a 1-by-1 or 2-by-2 block.
297: *
298:          IF( PAIR ) THEN
299:             PAIR = .FALSE.
300:             GO TO 60
301:          ELSE
302:             IF( K.LT.N )
303:      $         PAIR = T( K+1, K ).NE.ZERO
304:          END IF
305: *
306: *        Determine whether condition numbers are required for the k-th
307: *        eigenpair.
308: *
309:          IF( SOMCON ) THEN
310:             IF( PAIR ) THEN
311:                IF( .NOT.SELECT( K ) .AND. .NOT.SELECT( K+1 ) )
312:      $            GO TO 60
313:             ELSE
314:                IF( .NOT.SELECT( K ) )
315:      $            GO TO 60
316:             END IF
317:          END IF
318: *
319:          KS = KS + 1
320: *
321:          IF( WANTS ) THEN
322: *
323: *           Compute the reciprocal condition number of the k-th
324: *           eigenvalue.
325: *
326:             IF( .NOT.PAIR ) THEN
327: *
328: *              Real eigenvalue.
329: *
330:                PROD = DDOT( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
331:                RNRM = DNRM2( N, VR( 1, KS ), 1 )
332:                LNRM = DNRM2( N, VL( 1, KS ), 1 )
333:                S( KS ) = ABS( PROD ) / ( RNRM*LNRM )
334:             ELSE
335: *
336: *              Complex eigenvalue.
337: *
338:                PROD1 = DDOT( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
339:                PROD1 = PROD1 + DDOT( N, VR( 1, KS+1 ), 1, VL( 1, KS+1 ),
340:      $                 1 )
341:                PROD2 = DDOT( N, VL( 1, KS ), 1, VR( 1, KS+1 ), 1 )
342:                PROD2 = PROD2 - DDOT( N, VL( 1, KS+1 ), 1, VR( 1, KS ),
343:      $                 1 )
344:                RNRM = DLAPY2( DNRM2( N, VR( 1, KS ), 1 ),
345:      $                DNRM2( N, VR( 1, KS+1 ), 1 ) )
346:                LNRM = DLAPY2( DNRM2( N, VL( 1, KS ), 1 ),
347:      $                DNRM2( N, VL( 1, KS+1 ), 1 ) )
348:                COND = DLAPY2( PROD1, PROD2 ) / ( RNRM*LNRM )
349:                S( KS ) = COND
350:                S( KS+1 ) = COND
351:             END IF
352:          END IF
353: *
354:          IF( WANTSP ) THEN
355: *
356: *           Estimate the reciprocal condition number of the k-th
357: *           eigenvector.
358: *
359: *           Copy the matrix T to the array WORK and swap the diagonal
360: *           block beginning at T(k,k) to the (1,1) position.
361: *
362:             CALL DLACPY( 'Full', N, N, T, LDT, WORK, LDWORK )
363:             IFST = K
364:             ILST = 1
365:             CALL DTREXC( 'No Q', N, WORK, LDWORK, DUMMY, 1, IFST, ILST,
366:      $                   WORK( 1, N+1 ), IERR )
367: *
368:             IF( IERR.EQ.1 .OR. IERR.EQ.2 ) THEN
369: *
370: *              Could not swap because blocks not well separated
371: *
372:                SCALE = ONE
373:                EST = BIGNUM
374:             ELSE
375: *
376: *              Reordering successful
377: *
378:                IF( WORK( 2, 1 ).EQ.ZERO ) THEN
379: *
380: *                 Form C = T22 - lambda*I in WORK(2:N,2:N).
381: *
382:                   DO 20 I = 2, N
383:                      WORK( I, I ) = WORK( I, I ) - WORK( 1, 1 )
384:    20             CONTINUE
385:                   N2 = 1
386:                   NN = N - 1
387:                ELSE
388: *
389: *                 Triangularize the 2 by 2 block by unitary
390: *                 transformation U = [  cs   i*ss ]
391: *                                    [ i*ss   cs  ].
392: *                 such that the (1,1) position of WORK is complex
393: *                 eigenvalue lambda with positive imaginary part. (2,2)
394: *                 position of WORK is the complex eigenvalue lambda
395: *                 with negative imaginary  part.
396: *
397:                   MU = SQRT( ABS( WORK( 1, 2 ) ) )*
398:      $                 SQRT( ABS( WORK( 2, 1 ) ) )
399:                   DELTA = DLAPY2( MU, WORK( 2, 1 ) )
400:                   CS = MU / DELTA
401:                   SN = -WORK( 2, 1 ) / DELTA
402: *
403: *                 Form
404: *
405: *                 C' = WORK(2:N,2:N) + i*[rwork(1) ..... rwork(n-1) ]
406: *                                        [   mu                     ]
407: *                                        [         ..               ]
408: *                                        [             ..           ]
409: *                                        [                  mu      ]
410: *                 where C' is conjugate transpose of complex matrix C,
411: *                 and RWORK is stored starting in the N+1-st column of
412: *                 WORK.
413: *
414:                   DO 30 J = 3, N
415:                      WORK( 2, J ) = CS*WORK( 2, J )
416:                      WORK( J, J ) = WORK( J, J ) - WORK( 1, 1 )
417:    30             CONTINUE
418:                   WORK( 2, 2 ) = ZERO
419: *
420:                   WORK( 1, N+1 ) = TWO*MU
421:                   DO 40 I = 2, N - 1
422:                      WORK( I, N+1 ) = SN*WORK( 1, I+1 )
423:    40             CONTINUE
424:                   N2 = 2
425:                   NN = 2*( N-1 )
426:                END IF
427: *
428: *              Estimate norm(inv(C'))
429: *
430:                EST = ZERO
431:                KASE = 0
432:    50          CONTINUE
433:                CALL DLACN2( NN, WORK( 1, N+2 ), WORK( 1, N+4 ), IWORK,
434:      $                      EST, KASE, ISAVE )
435:                IF( KASE.NE.0 ) THEN
436:                   IF( KASE.EQ.1 ) THEN
437:                      IF( N2.EQ.1 ) THEN
438: *
439: *                       Real eigenvalue: solve C'*x = scale*c.
440: *
441:                         CALL DLAQTR( .TRUE., .TRUE., N-1, WORK( 2, 2 ),
442:      $                               LDWORK, DUMMY, DUMM, SCALE,
443:      $                               WORK( 1, N+4 ), WORK( 1, N+6 ),
444:      $                               IERR )
445:                      ELSE
446: *
447: *                       Complex eigenvalue: solve
448: *                       C'*(p+iq) = scale*(c+id) in real arithmetic.
449: *
450:                         CALL DLAQTR( .TRUE., .FALSE., N-1, WORK( 2, 2 ),
451:      $                               LDWORK, WORK( 1, N+1 ), MU, SCALE,
452:      $                               WORK( 1, N+4 ), WORK( 1, N+6 ),
453:      $                               IERR )
454:                      END IF
455:                   ELSE
456:                      IF( N2.EQ.1 ) THEN
457: *
458: *                       Real eigenvalue: solve C*x = scale*c.
459: *
460:                         CALL DLAQTR( .FALSE., .TRUE., N-1, WORK( 2, 2 ),
461:      $                               LDWORK, DUMMY, DUMM, SCALE,
462:      $                               WORK( 1, N+4 ), WORK( 1, N+6 ),
463:      $                               IERR )
464:                      ELSE
465: *
466: *                       Complex eigenvalue: solve
467: *                       C*(p+iq) = scale*(c+id) in real arithmetic.
468: *
469:                         CALL DLAQTR( .FALSE., .FALSE., N-1,
470:      $                               WORK( 2, 2 ), LDWORK,
471:      $                               WORK( 1, N+1 ), MU, SCALE,
472:      $                               WORK( 1, N+4 ), WORK( 1, N+6 ),
473:      $                               IERR )
474: *
475:                      END IF
476:                   END IF
477: *
478:                   GO TO 50
479:                END IF
480:             END IF
481: *
482:             SEP( KS ) = SCALE / MAX( EST, SMLNUM )
483:             IF( PAIR )
484:      $         SEP( KS+1 ) = SEP( KS )
485:          END IF
486: *
487:          IF( PAIR )
488:      $      KS = KS + 1
489: *
490:    60 CONTINUE
491:       RETURN
492: *
493: *     End of DTRSNA
494: *
495:       END
496: