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