001:       SUBROUTINE ZTGSNA( JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL,
002:      $                   LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK,
003:      $                   IWORK, 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: *     .. Scalar Arguments ..
011:       CHARACTER          HOWMNY, JOB
012:       INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, M, MM, N
013: *     ..
014: *     .. Array Arguments ..
015:       LOGICAL            SELECT( * )
016:       INTEGER            IWORK( * )
017:       DOUBLE PRECISION   DIF( * ), S( * )
018:       COMPLEX*16         A( LDA, * ), B( LDB, * ), VL( LDVL, * ),
019:      $                   VR( LDVR, * ), WORK( * )
020: *     ..
021: *
022: *  Purpose
023: *  =======
024: *
025: *  ZTGSNA estimates reciprocal condition numbers for specified
026: *  eigenvalues and/or eigenvectors of a matrix pair (A, B).
027: *
028: *  (A, B) must be in generalized Schur canonical form, that is, A and
029: *  B are both upper triangular.
030: *
031: *  Arguments
032: *  =========
033: *
034: *  JOB     (input) CHARACTER*1
035: *          Specifies whether condition numbers are required for
036: *          eigenvalues (S) or eigenvectors (DIF):
037: *          = 'E': for eigenvalues only (S);
038: *          = 'V': for eigenvectors only (DIF);
039: *          = 'B': for both eigenvalues and eigenvectors (S and DIF).
040: *
041: *  HOWMNY  (input) CHARACTER*1
042: *          = 'A': compute condition numbers for all eigenpairs;
043: *          = 'S': compute condition numbers for selected eigenpairs
044: *                 specified by the array SELECT.
045: *
046: *  SELECT  (input) LOGICAL array, dimension (N)
047: *          If HOWMNY = 'S', SELECT specifies the eigenpairs for which
048: *          condition numbers are required. To select condition numbers
049: *          for the corresponding j-th eigenvalue and/or eigenvector,
050: *          SELECT(j) must be set to .TRUE..
051: *          If HOWMNY = 'A', SELECT is not referenced.
052: *
053: *  N       (input) INTEGER
054: *          The order of the square matrix pair (A, B). N >= 0.
055: *
056: *  A       (input) COMPLEX*16 array, dimension (LDA,N)
057: *          The upper triangular matrix A in the pair (A,B).
058: *
059: *  LDA     (input) INTEGER
060: *          The leading dimension of the array A. LDA >= max(1,N).
061: *
062: *  B       (input) COMPLEX*16 array, dimension (LDB,N)
063: *          The upper triangular matrix B in the pair (A, B).
064: *
065: *  LDB     (input) INTEGER
066: *          The leading dimension of the array B. LDB >= max(1,N).
067: *
068: *  VL      (input) COMPLEX*16 array, dimension (LDVL,M)
069: *          IF JOB = 'E' or 'B', VL must contain left eigenvectors of
070: *          (A, B), corresponding to the eigenpairs specified by HOWMNY
071: *          and SELECT.  The eigenvectors must be stored in consecutive
072: *          columns of VL, as returned by ZTGEVC.
073: *          If JOB = 'V', VL is not referenced.
074: *
075: *  LDVL    (input) INTEGER
076: *          The leading dimension of the array VL. LDVL >= 1; and
077: *          If JOB = 'E' or 'B', LDVL >= N.
078: *
079: *  VR      (input) COMPLEX*16 array, dimension (LDVR,M)
080: *          IF JOB = 'E' or 'B', VR must contain right eigenvectors of
081: *          (A, B), corresponding to the eigenpairs specified by HOWMNY
082: *          and SELECT.  The eigenvectors must be stored in consecutive
083: *          columns of VR, as returned by ZTGEVC.
084: *          If JOB = 'V', VR is not referenced.
085: *
086: *  LDVR    (input) INTEGER
087: *          The leading dimension of the array VR. LDVR >= 1;
088: *          If JOB = 'E' or 'B', LDVR >= N.
089: *
090: *  S       (output) DOUBLE PRECISION array, dimension (MM)
091: *          If JOB = 'E' or 'B', the reciprocal condition numbers of the
092: *          selected eigenvalues, stored in consecutive elements of the
093: *          array.
094: *          If JOB = 'V', S is not referenced.
095: *
096: *  DIF     (output) DOUBLE PRECISION array, dimension (MM)
097: *          If JOB = 'V' or 'B', the estimated reciprocal condition
098: *          numbers of the selected eigenvectors, stored in consecutive
099: *          elements of the array.
100: *          If the eigenvalues cannot be reordered to compute DIF(j),
101: *          DIF(j) is set to 0; this can only occur when the true value
102: *          would be very small anyway.
103: *          For each eigenvalue/vector specified by SELECT, DIF stores
104: *          a Frobenius norm-based estimate of Difl.
105: *          If JOB = 'E', DIF is not referenced.
106: *
107: *  MM      (input) INTEGER
108: *          The number of elements in the arrays S and DIF. MM >= M.
109: *
110: *  M       (output) INTEGER
111: *          The number of elements of the arrays S and DIF used to store
112: *          the specified condition numbers; for each selected eigenvalue
113: *          one element is used. If HOWMNY = 'A', M is set to N.
114: *
115: *  WORK    (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
116: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
117: *
118: *  LWORK  (input) INTEGER
119: *          The dimension of the array WORK. LWORK >= max(1,N).
120: *          If JOB = 'V' or 'B', LWORK >= max(1,2*N*N).
121: *
122: *  IWORK   (workspace) INTEGER array, dimension (N+2)
123: *          If JOB = 'E', IWORK is not referenced.
124: *
125: *  INFO    (output) INTEGER
126: *          = 0: Successful exit
127: *          < 0: If INFO = -i, the i-th argument had an illegal value
128: *
129: *  Further Details
130: *  ===============
131: *
132: *  The reciprocal of the condition number of the i-th generalized
133: *  eigenvalue w = (a, b) is defined as
134: *
135: *          S(I) = (|v'Au|**2 + |v'Bu|**2)**(1/2) / (norm(u)*norm(v))
136: *
137: *  where u and v are the right and left eigenvectors of (A, B)
138: *  corresponding to w; |z| denotes the absolute value of the complex
139: *  number, and norm(u) denotes the 2-norm of the vector u. The pair
140: *  (a, b) corresponds to an eigenvalue w = a/b (= v'Au/v'Bu) of the
141: *  matrix pair (A, B). If both a and b equal zero, then (A,B) is
142: *  singular and S(I) = -1 is returned.
143: *
144: *  An approximate error bound on the chordal distance between the i-th
145: *  computed generalized eigenvalue w and the corresponding exact
146: *  eigenvalue lambda is
147: *
148: *          chord(w, lambda) <=   EPS * norm(A, B) / S(I),
149: *
150: *  where EPS is the machine precision.
151: *
152: *  The reciprocal of the condition number of the right eigenvector u
153: *  and left eigenvector v corresponding to the generalized eigenvalue w
154: *  is defined as follows. Suppose
155: *
156: *                   (A, B) = ( a   *  ) ( b  *  )  1
157: *                            ( 0  A22 ),( 0 B22 )  n-1
158: *                              1  n-1     1 n-1
159: *
160: *  Then the reciprocal condition number DIF(I) is
161: *
162: *          Difl[(a, b), (A22, B22)]  = sigma-min( Zl )
163: *
164: *  where sigma-min(Zl) denotes the smallest singular value of
165: *
166: *         Zl = [ kron(a, In-1) -kron(1, A22) ]
167: *              [ kron(b, In-1) -kron(1, B22) ].
168: *
169: *  Here In-1 is the identity matrix of size n-1 and X' is the conjugate
170: *  transpose of X. kron(X, Y) is the Kronecker product between the
171: *  matrices X and Y.
172: *
173: *  We approximate the smallest singular value of Zl with an upper
174: *  bound. This is done by ZLATDF.
175: *
176: *  An approximate error bound for a computed eigenvector VL(i) or
177: *  VR(i) is given by
178: *
179: *                      EPS * norm(A, B) / DIF(i).
180: *
181: *  See ref. [2-3] for more details and further references.
182: *
183: *  Based on contributions by
184: *     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
185: *     Umea University, S-901 87 Umea, Sweden.
186: *
187: *  References
188: *  ==========
189: *
190: *  [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
191: *      Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
192: *      M.S. Moonen et al (eds), Linear Algebra for Large Scale and
193: *      Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
194: *
195: *  [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
196: *      Eigenvalues of a Regular Matrix Pair (A, B) and Condition
197: *      Estimation: Theory, Algorithms and Software, Report
198: *      UMINF - 94.04, Department of Computing Science, Umea University,
199: *      S-901 87 Umea, Sweden, 1994. Also as LAPACK Working Note 87.
200: *      To appear in Numerical Algorithms, 1996.
201: *
202: *  [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
203: *      for Solving the Generalized Sylvester Equation and Estimating the
204: *      Separation between Regular Matrix Pairs, Report UMINF - 93.23,
205: *      Department of Computing Science, Umea University, S-901 87 Umea,
206: *      Sweden, December 1993, Revised April 1994, Also as LAPACK Working
207: *      Note 75.
208: *      To appear in ACM Trans. on Math. Software, Vol 22, No 1, 1996.
209: *
210: *  =====================================================================
211: *
212: *     .. Parameters ..
213:       DOUBLE PRECISION   ZERO, ONE
214:       INTEGER            IDIFJB
215:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, IDIFJB = 3 )
216: *     ..
217: *     .. Local Scalars ..
218:       LOGICAL            LQUERY, SOMCON, WANTBH, WANTDF, WANTS
219:       INTEGER            I, IERR, IFST, ILST, K, KS, LWMIN, N1, N2
220:       DOUBLE PRECISION   BIGNUM, COND, EPS, LNRM, RNRM, SCALE, SMLNUM
221:       COMPLEX*16         YHAX, YHBX
222: *     ..
223: *     .. Local Arrays ..
224:       COMPLEX*16         DUMMY( 1 ), DUMMY1( 1 )
225: *     ..
226: *     .. External Functions ..
227:       LOGICAL            LSAME
228:       DOUBLE PRECISION   DLAMCH, DLAPY2, DZNRM2
229:       COMPLEX*16         ZDOTC
230:       EXTERNAL           LSAME, DLAMCH, DLAPY2, DZNRM2, ZDOTC
231: *     ..
232: *     .. External Subroutines ..
233:       EXTERNAL           DLABAD, XERBLA, ZGEMV, ZLACPY, ZTGEXC, ZTGSYL
234: *     ..
235: *     .. Intrinsic Functions ..
236:       INTRINSIC          ABS, DCMPLX, MAX
237: *     ..
238: *     .. Executable Statements ..
239: *
240: *     Decode and test the input parameters
241: *
242:       WANTBH = LSAME( JOB, 'B' )
243:       WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
244:       WANTDF = LSAME( JOB, 'V' ) .OR. WANTBH
245: *
246:       SOMCON = LSAME( HOWMNY, 'S' )
247: *
248:       INFO = 0
249:       LQUERY = ( LWORK.EQ.-1 )
250: *
251:       IF( .NOT.WANTS .AND. .NOT.WANTDF ) THEN
252:          INFO = -1
253:       ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
254:          INFO = -2
255:       ELSE IF( N.LT.0 ) THEN
256:          INFO = -4
257:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
258:          INFO = -6
259:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
260:          INFO = -8
261:       ELSE IF( WANTS .AND. LDVL.LT.N ) THEN
262:          INFO = -10
263:       ELSE IF( WANTS .AND. LDVR.LT.N ) THEN
264:          INFO = -12
265:       ELSE
266: *
267: *        Set M to the number of eigenpairs for which condition numbers
268: *        are required, and test MM.
269: *
270:          IF( SOMCON ) THEN
271:             M = 0
272:             DO 10 K = 1, N
273:                IF( SELECT( K ) )
274:      $            M = M + 1
275:    10       CONTINUE
276:          ELSE
277:             M = N
278:          END IF
279: *
280:          IF( N.EQ.0 ) THEN
281:             LWMIN = 1
282:          ELSE IF( LSAME( JOB, 'V' ) .OR. LSAME( JOB, 'B' ) ) THEN
283:             LWMIN = 2*N*N
284:          ELSE
285:             LWMIN = N
286:          END IF
287:          WORK( 1 ) = LWMIN
288: *
289:          IF( MM.LT.M ) THEN
290:             INFO = -15
291:          ELSE IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
292:             INFO = -18
293:          END IF
294:       END IF
295: *
296:       IF( INFO.NE.0 ) THEN
297:          CALL XERBLA( 'ZTGSNA', -INFO )
298:          RETURN
299:       ELSE IF( LQUERY ) THEN
300:          RETURN
301:       END IF
302: *
303: *     Quick return if possible
304: *
305:       IF( N.EQ.0 )
306:      $   RETURN
307: *
308: *     Get machine constants
309: *
310:       EPS = DLAMCH( 'P' )
311:       SMLNUM = DLAMCH( 'S' ) / EPS
312:       BIGNUM = ONE / SMLNUM
313:       CALL DLABAD( SMLNUM, BIGNUM )
314:       KS = 0
315:       DO 20 K = 1, N
316: *
317: *        Determine whether condition numbers are required for the k-th
318: *        eigenpair.
319: *
320:          IF( SOMCON ) THEN
321:             IF( .NOT.SELECT( K ) )
322:      $         GO TO 20
323:          END IF
324: *
325:          KS = KS + 1
326: *
327:          IF( WANTS ) THEN
328: *
329: *           Compute the reciprocal condition number of the k-th
330: *           eigenvalue.
331: *
332:             RNRM = DZNRM2( N, VR( 1, KS ), 1 )
333:             LNRM = DZNRM2( N, VL( 1, KS ), 1 )
334:             CALL ZGEMV( 'N', N, N, DCMPLX( ONE, ZERO ), A, LDA,
335:      $                  VR( 1, KS ), 1, DCMPLX( ZERO, ZERO ), WORK, 1 )
336:             YHAX = ZDOTC( N, WORK, 1, VL( 1, KS ), 1 )
337:             CALL ZGEMV( 'N', N, N, DCMPLX( ONE, ZERO ), B, LDB,
338:      $                  VR( 1, KS ), 1, DCMPLX( ZERO, ZERO ), WORK, 1 )
339:             YHBX = ZDOTC( N, WORK, 1, VL( 1, KS ), 1 )
340:             COND = DLAPY2( ABS( YHAX ), ABS( YHBX ) )
341:             IF( COND.EQ.ZERO ) THEN
342:                S( KS ) = -ONE
343:             ELSE
344:                S( KS ) = COND / ( RNRM*LNRM )
345:             END IF
346:          END IF
347: *
348:          IF( WANTDF ) THEN
349:             IF( N.EQ.1 ) THEN
350:                DIF( KS ) = DLAPY2( ABS( A( 1, 1 ) ), ABS( B( 1, 1 ) ) )
351:             ELSE
352: *
353: *              Estimate the reciprocal condition number of the k-th
354: *              eigenvectors.
355: *
356: *              Copy the matrix (A, B) to the array WORK and move the
357: *              (k,k)th pair to the (1,1) position.
358: *
359:                CALL ZLACPY( 'Full', N, N, A, LDA, WORK, N )
360:                CALL ZLACPY( 'Full', N, N, B, LDB, WORK( N*N+1 ), N )
361:                IFST = K
362:                ILST = 1
363: *
364:                CALL ZTGEXC( .FALSE., .FALSE., N, WORK, N, WORK( N*N+1 ),
365:      $                      N, DUMMY, 1, DUMMY1, 1, IFST, ILST, IERR )
366: *
367:                IF( IERR.GT.0 ) THEN
368: *
369: *                 Ill-conditioned problem - swap rejected.
370: *
371:                   DIF( KS ) = ZERO
372:                ELSE
373: *
374: *                 Reordering successful, solve generalized Sylvester
375: *                 equation for R and L,
376: *                            A22 * R - L * A11 = A12
377: *                            B22 * R - L * B11 = B12,
378: *                 and compute estimate of Difl[(A11,B11), (A22, B22)].
379: *
380:                   N1 = 1
381:                   N2 = N - N1
382:                   I = N*N + 1
383:                   CALL ZTGSYL( 'N', IDIFJB, N2, N1, WORK( N*N1+N1+1 ),
384:      $                         N, WORK, N, WORK( N1+1 ), N,
385:      $                         WORK( N*N1+N1+I ), N, WORK( I ), N,
386:      $                         WORK( N1+I ), N, SCALE, DIF( KS ), DUMMY,
387:      $                         1, IWORK, IERR )
388:                END IF
389:             END IF
390:          END IF
391: *
392:    20 CONTINUE
393:       WORK( 1 ) = LWMIN
394:       RETURN
395: *
396: *     End of ZTGSNA
397: *
398:       END
399: