Annotation of rpl/lapack/lapack/dtrsna.f, revision 1.9
1.9 ! bertrand 1: *> \brief \b DTRSNA
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DTRSNA + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrsna.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrsna.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrsna.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
! 22: * LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK,
! 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: * INTEGER IWORK( * )
! 32: * DOUBLE PRECISION S( * ), SEP( * ), T( LDT, * ), VL( LDVL, * ),
! 33: * $ VR( LDVR, * ), WORK( LDWORK, * )
! 34: * ..
! 35: *
! 36: *
! 37: *> \par Purpose:
! 38: * =============
! 39: *>
! 40: *> \verbatim
! 41: *>
! 42: *> DTRSNA estimates reciprocal condition numbers for specified
! 43: *> eigenvalues and/or right eigenvectors of a real upper
! 44: *> quasi-triangular matrix T (or of any matrix Q*T*Q**T with Q
! 45: *> orthogonal).
! 46: *>
! 47: *> T must be in Schur canonical form (as returned by DHSEQR), that is,
! 48: *> block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
! 49: *> 2-by-2 diagonal block has its diagonal elements equal and its
! 50: *> off-diagonal elements of opposite sign.
! 51: *> \endverbatim
! 52: *
! 53: * Arguments:
! 54: * ==========
! 55: *
! 56: *> \param[in] JOB
! 57: *> \verbatim
! 58: *> JOB is CHARACTER*1
! 59: *> Specifies whether condition numbers are required for
! 60: *> eigenvalues (S) or eigenvectors (SEP):
! 61: *> = 'E': for eigenvalues only (S);
! 62: *> = 'V': for eigenvectors only (SEP);
! 63: *> = 'B': for both eigenvalues and eigenvectors (S and SEP).
! 64: *> \endverbatim
! 65: *>
! 66: *> \param[in] HOWMNY
! 67: *> \verbatim
! 68: *> HOWMNY is CHARACTER*1
! 69: *> = 'A': compute condition numbers for all eigenpairs;
! 70: *> = 'S': compute condition numbers for selected eigenpairs
! 71: *> specified by the array SELECT.
! 72: *> \endverbatim
! 73: *>
! 74: *> \param[in] SELECT
! 75: *> \verbatim
! 76: *> SELECT is LOGICAL array, dimension (N)
! 77: *> If HOWMNY = 'S', SELECT specifies the eigenpairs for which
! 78: *> condition numbers are required. To select condition numbers
! 79: *> for the eigenpair corresponding to a real eigenvalue w(j),
! 80: *> SELECT(j) must be set to .TRUE.. To select condition numbers
! 81: *> corresponding to a complex conjugate pair of eigenvalues w(j)
! 82: *> and w(j+1), either SELECT(j) or SELECT(j+1) or both, must be
! 83: *> set to .TRUE..
! 84: *> If HOWMNY = 'A', SELECT is not referenced.
! 85: *> \endverbatim
! 86: *>
! 87: *> \param[in] N
! 88: *> \verbatim
! 89: *> N is INTEGER
! 90: *> The order of the matrix T. N >= 0.
! 91: *> \endverbatim
! 92: *>
! 93: *> \param[in] T
! 94: *> \verbatim
! 95: *> T is DOUBLE PRECISION array, dimension (LDT,N)
! 96: *> The upper quasi-triangular matrix T, in Schur canonical form.
! 97: *> \endverbatim
! 98: *>
! 99: *> \param[in] LDT
! 100: *> \verbatim
! 101: *> LDT is INTEGER
! 102: *> The leading dimension of the array T. LDT >= max(1,N).
! 103: *> \endverbatim
! 104: *>
! 105: *> \param[in] VL
! 106: *> \verbatim
! 107: *> VL is DOUBLE PRECISION array, dimension (LDVL,M)
! 108: *> If JOB = 'E' or 'B', VL must contain left eigenvectors of T
! 109: *> (or of any Q*T*Q**T with Q orthogonal), corresponding to the
! 110: *> eigenpairs specified by HOWMNY and SELECT. The eigenvectors
! 111: *> must be stored in consecutive columns of VL, as returned by
! 112: *> DHSEIN or DTREVC.
! 113: *> If JOB = 'V', VL is not referenced.
! 114: *> \endverbatim
! 115: *>
! 116: *> \param[in] LDVL
! 117: *> \verbatim
! 118: *> LDVL is INTEGER
! 119: *> The leading dimension of the array VL.
! 120: *> LDVL >= 1; and if JOB = 'E' or 'B', LDVL >= N.
! 121: *> \endverbatim
! 122: *>
! 123: *> \param[in] VR
! 124: *> \verbatim
! 125: *> VR is DOUBLE PRECISION array, dimension (LDVR,M)
! 126: *> If JOB = 'E' or 'B', VR must contain right eigenvectors of T
! 127: *> (or of any Q*T*Q**T with Q orthogonal), corresponding to the
! 128: *> eigenpairs specified by HOWMNY and SELECT. The eigenvectors
! 129: *> must be stored in consecutive columns of VR, as returned by
! 130: *> DHSEIN or DTREVC.
! 131: *> If JOB = 'V', VR is not referenced.
! 132: *> \endverbatim
! 133: *>
! 134: *> \param[in] LDVR
! 135: *> \verbatim
! 136: *> LDVR is INTEGER
! 137: *> The leading dimension of the array VR.
! 138: *> LDVR >= 1; and if JOB = 'E' or 'B', LDVR >= N.
! 139: *> \endverbatim
! 140: *>
! 141: *> \param[out] S
! 142: *> \verbatim
! 143: *> S is DOUBLE PRECISION array, dimension (MM)
! 144: *> If JOB = 'E' or 'B', the reciprocal condition numbers of the
! 145: *> selected eigenvalues, stored in consecutive elements of the
! 146: *> array. For a complex conjugate pair of eigenvalues two
! 147: *> consecutive elements of S are set to the same value. Thus
! 148: *> S(j), SEP(j), and the j-th columns of VL and VR all
! 149: *> correspond to the same eigenpair (but not in general the
! 150: *> j-th eigenpair, unless all eigenpairs are selected).
! 151: *> If JOB = 'V', S is not referenced.
! 152: *> \endverbatim
! 153: *>
! 154: *> \param[out] SEP
! 155: *> \verbatim
! 156: *> SEP is DOUBLE PRECISION array, dimension (MM)
! 157: *> If JOB = 'V' or 'B', the estimated reciprocal condition
! 158: *> numbers of the selected eigenvectors, stored in consecutive
! 159: *> elements of the array. For a complex eigenvector two
! 160: *> consecutive elements of SEP are set to the same value. If
! 161: *> the eigenvalues cannot be reordered to compute SEP(j), SEP(j)
! 162: *> is set to 0; this can only occur when the true value would be
! 163: *> very small anyway.
! 164: *> If JOB = 'E', SEP is not referenced.
! 165: *> \endverbatim
! 166: *>
! 167: *> \param[in] MM
! 168: *> \verbatim
! 169: *> MM is INTEGER
! 170: *> The number of elements in the arrays S (if JOB = 'E' or 'B')
! 171: *> and/or SEP (if JOB = 'V' or 'B'). MM >= M.
! 172: *> \endverbatim
! 173: *>
! 174: *> \param[out] M
! 175: *> \verbatim
! 176: *> M is INTEGER
! 177: *> The number of elements of the arrays S and/or SEP actually
! 178: *> used to store the estimated condition numbers.
! 179: *> If HOWMNY = 'A', M is set to N.
! 180: *> \endverbatim
! 181: *>
! 182: *> \param[out] WORK
! 183: *> \verbatim
! 184: *> WORK is DOUBLE PRECISION array, dimension (LDWORK,N+6)
! 185: *> If JOB = 'E', WORK is not referenced.
! 186: *> \endverbatim
! 187: *>
! 188: *> \param[in] LDWORK
! 189: *> \verbatim
! 190: *> LDWORK is INTEGER
! 191: *> The leading dimension of the array WORK.
! 192: *> LDWORK >= 1; and if JOB = 'V' or 'B', LDWORK >= N.
! 193: *> \endverbatim
! 194: *>
! 195: *> \param[out] IWORK
! 196: *> \verbatim
! 197: *> IWORK is INTEGER array, dimension (2*(N-1))
! 198: *> If JOB = 'E', IWORK is not referenced.
! 199: *> \endverbatim
! 200: *>
! 201: *> \param[out] INFO
! 202: *> \verbatim
! 203: *> INFO is INTEGER
! 204: *> = 0: successful exit
! 205: *> < 0: if INFO = -i, the i-th argument had an illegal value
! 206: *> \endverbatim
! 207: *
! 208: * Authors:
! 209: * ========
! 210: *
! 211: *> \author Univ. of Tennessee
! 212: *> \author Univ. of California Berkeley
! 213: *> \author Univ. of Colorado Denver
! 214: *> \author NAG Ltd.
! 215: *
! 216: *> \date November 2011
! 217: *
! 218: *> \ingroup doubleOTHERcomputational
! 219: *
! 220: *> \par Further Details:
! 221: * =====================
! 222: *>
! 223: *> \verbatim
! 224: *>
! 225: *> The reciprocal of the condition number of an eigenvalue lambda is
! 226: *> defined as
! 227: *>
! 228: *> S(lambda) = |v**T*u| / (norm(u)*norm(v))
! 229: *>
! 230: *> where u and v are the right and left eigenvectors of T corresponding
! 231: *> to lambda; v**T denotes the transpose of v, and norm(u)
! 232: *> denotes the Euclidean norm. These reciprocal condition numbers always
! 233: *> lie between zero (very badly conditioned) and one (very well
! 234: *> conditioned). If n = 1, S(lambda) is defined to be 1.
! 235: *>
! 236: *> An approximate error bound for a computed eigenvalue W(i) is given by
! 237: *>
! 238: *> EPS * norm(T) / S(i)
! 239: *>
! 240: *> where EPS is the machine precision.
! 241: *>
! 242: *> The reciprocal of the condition number of the right eigenvector u
! 243: *> corresponding to lambda is defined as follows. Suppose
! 244: *>
! 245: *> T = ( lambda c )
! 246: *> ( 0 T22 )
! 247: *>
! 248: *> Then the reciprocal condition number is
! 249: *>
! 250: *> SEP( lambda, T22 ) = sigma-min( T22 - lambda*I )
! 251: *>
! 252: *> where sigma-min denotes the smallest singular value. We approximate
! 253: *> the smallest singular value by the reciprocal of an estimate of the
! 254: *> one-norm of the inverse of T22 - lambda*I. If n = 1, SEP(1) is
! 255: *> defined to be abs(T(1,1)).
! 256: *>
! 257: *> An approximate error bound for a computed right eigenvector VR(i)
! 258: *> is given by
! 259: *>
! 260: *> EPS * norm(T) / SEP(i)
! 261: *> \endverbatim
! 262: *>
! 263: * =====================================================================
1.1 bertrand 264: SUBROUTINE DTRSNA( JOB, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
265: $ LDVR, S, SEP, MM, M, WORK, LDWORK, IWORK,
266: $ INFO )
267: *
1.9 ! bertrand 268: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 269: * -- LAPACK is a software package provided by Univ. of Tennessee, --
270: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 271: * November 2011
1.1 bertrand 272: *
273: * .. Scalar Arguments ..
274: CHARACTER HOWMNY, JOB
275: INTEGER INFO, LDT, LDVL, LDVR, LDWORK, M, MM, N
276: * ..
277: * .. Array Arguments ..
278: LOGICAL SELECT( * )
279: INTEGER IWORK( * )
280: DOUBLE PRECISION S( * ), SEP( * ), T( LDT, * ), VL( LDVL, * ),
281: $ VR( LDVR, * ), WORK( LDWORK, * )
282: * ..
283: *
284: * =====================================================================
285: *
286: * .. Parameters ..
287: DOUBLE PRECISION ZERO, ONE, TWO
288: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
289: * ..
290: * .. Local Scalars ..
291: LOGICAL PAIR, SOMCON, WANTBH, WANTS, WANTSP
292: INTEGER I, IERR, IFST, ILST, J, K, KASE, KS, N2, NN
293: DOUBLE PRECISION BIGNUM, COND, CS, DELTA, DUMM, EPS, EST, LNRM,
294: $ MU, PROD, PROD1, PROD2, RNRM, SCALE, SMLNUM, SN
295: * ..
296: * .. Local Arrays ..
297: INTEGER ISAVE( 3 )
298: DOUBLE PRECISION DUMMY( 1 )
299: * ..
300: * .. External Functions ..
301: LOGICAL LSAME
302: DOUBLE PRECISION DDOT, DLAMCH, DLAPY2, DNRM2
303: EXTERNAL LSAME, DDOT, DLAMCH, DLAPY2, DNRM2
304: * ..
305: * .. External Subroutines ..
306: EXTERNAL DLACN2, DLACPY, DLAQTR, DTREXC, XERBLA
307: * ..
308: * .. Intrinsic Functions ..
309: INTRINSIC ABS, MAX, SQRT
310: * ..
311: * .. Executable Statements ..
312: *
313: * Decode and test the input parameters
314: *
315: WANTBH = LSAME( JOB, 'B' )
316: WANTS = LSAME( JOB, 'E' ) .OR. WANTBH
317: WANTSP = LSAME( JOB, 'V' ) .OR. WANTBH
318: *
319: SOMCON = LSAME( HOWMNY, 'S' )
320: *
321: INFO = 0
322: IF( .NOT.WANTS .AND. .NOT.WANTSP ) THEN
323: INFO = -1
324: ELSE IF( .NOT.LSAME( HOWMNY, 'A' ) .AND. .NOT.SOMCON ) THEN
325: INFO = -2
326: ELSE IF( N.LT.0 ) THEN
327: INFO = -4
328: ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
329: INFO = -6
330: ELSE IF( LDVL.LT.1 .OR. ( WANTS .AND. LDVL.LT.N ) ) THEN
331: INFO = -8
332: ELSE IF( LDVR.LT.1 .OR. ( WANTS .AND. LDVR.LT.N ) ) THEN
333: INFO = -10
334: ELSE
335: *
336: * Set M to the number of eigenpairs for which condition numbers
337: * are required, and test MM.
338: *
339: IF( SOMCON ) THEN
340: M = 0
341: PAIR = .FALSE.
342: DO 10 K = 1, N
343: IF( PAIR ) THEN
344: PAIR = .FALSE.
345: ELSE
346: IF( K.LT.N ) THEN
347: IF( T( K+1, K ).EQ.ZERO ) THEN
348: IF( SELECT( K ) )
349: $ M = M + 1
350: ELSE
351: PAIR = .TRUE.
352: IF( SELECT( K ) .OR. SELECT( K+1 ) )
353: $ M = M + 2
354: END IF
355: ELSE
356: IF( SELECT( N ) )
357: $ M = M + 1
358: END IF
359: END IF
360: 10 CONTINUE
361: ELSE
362: M = N
363: END IF
364: *
365: IF( MM.LT.M ) THEN
366: INFO = -13
367: ELSE IF( LDWORK.LT.1 .OR. ( WANTSP .AND. LDWORK.LT.N ) ) THEN
368: INFO = -16
369: END IF
370: END IF
371: IF( INFO.NE.0 ) THEN
372: CALL XERBLA( 'DTRSNA', -INFO )
373: RETURN
374: END IF
375: *
376: * Quick return if possible
377: *
378: IF( N.EQ.0 )
379: $ RETURN
380: *
381: IF( N.EQ.1 ) THEN
382: IF( SOMCON ) THEN
383: IF( .NOT.SELECT( 1 ) )
384: $ RETURN
385: END IF
386: IF( WANTS )
387: $ S( 1 ) = ONE
388: IF( WANTSP )
389: $ SEP( 1 ) = ABS( T( 1, 1 ) )
390: RETURN
391: END IF
392: *
393: * Get machine constants
394: *
395: EPS = DLAMCH( 'P' )
396: SMLNUM = DLAMCH( 'S' ) / EPS
397: BIGNUM = ONE / SMLNUM
398: CALL DLABAD( SMLNUM, BIGNUM )
399: *
400: KS = 0
401: PAIR = .FALSE.
402: DO 60 K = 1, N
403: *
404: * Determine whether T(k,k) begins a 1-by-1 or 2-by-2 block.
405: *
406: IF( PAIR ) THEN
407: PAIR = .FALSE.
408: GO TO 60
409: ELSE
410: IF( K.LT.N )
411: $ PAIR = T( K+1, K ).NE.ZERO
412: END IF
413: *
414: * Determine whether condition numbers are required for the k-th
415: * eigenpair.
416: *
417: IF( SOMCON ) THEN
418: IF( PAIR ) THEN
419: IF( .NOT.SELECT( K ) .AND. .NOT.SELECT( K+1 ) )
420: $ GO TO 60
421: ELSE
422: IF( .NOT.SELECT( K ) )
423: $ GO TO 60
424: END IF
425: END IF
426: *
427: KS = KS + 1
428: *
429: IF( WANTS ) THEN
430: *
431: * Compute the reciprocal condition number of the k-th
432: * eigenvalue.
433: *
434: IF( .NOT.PAIR ) THEN
435: *
436: * Real eigenvalue.
437: *
438: PROD = DDOT( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
439: RNRM = DNRM2( N, VR( 1, KS ), 1 )
440: LNRM = DNRM2( N, VL( 1, KS ), 1 )
441: S( KS ) = ABS( PROD ) / ( RNRM*LNRM )
442: ELSE
443: *
444: * Complex eigenvalue.
445: *
446: PROD1 = DDOT( N, VR( 1, KS ), 1, VL( 1, KS ), 1 )
447: PROD1 = PROD1 + DDOT( N, VR( 1, KS+1 ), 1, VL( 1, KS+1 ),
448: $ 1 )
449: PROD2 = DDOT( N, VL( 1, KS ), 1, VR( 1, KS+1 ), 1 )
450: PROD2 = PROD2 - DDOT( N, VL( 1, KS+1 ), 1, VR( 1, KS ),
451: $ 1 )
452: RNRM = DLAPY2( DNRM2( N, VR( 1, KS ), 1 ),
453: $ DNRM2( N, VR( 1, KS+1 ), 1 ) )
454: LNRM = DLAPY2( DNRM2( N, VL( 1, KS ), 1 ),
455: $ DNRM2( N, VL( 1, KS+1 ), 1 ) )
456: COND = DLAPY2( PROD1, PROD2 ) / ( RNRM*LNRM )
457: S( KS ) = COND
458: S( KS+1 ) = COND
459: END IF
460: END IF
461: *
462: IF( WANTSP ) THEN
463: *
464: * Estimate the reciprocal condition number of the k-th
465: * eigenvector.
466: *
467: * Copy the matrix T to the array WORK and swap the diagonal
468: * block beginning at T(k,k) to the (1,1) position.
469: *
470: CALL DLACPY( 'Full', N, N, T, LDT, WORK, LDWORK )
471: IFST = K
472: ILST = 1
473: CALL DTREXC( 'No Q', N, WORK, LDWORK, DUMMY, 1, IFST, ILST,
474: $ WORK( 1, N+1 ), IERR )
475: *
476: IF( IERR.EQ.1 .OR. IERR.EQ.2 ) THEN
477: *
478: * Could not swap because blocks not well separated
479: *
480: SCALE = ONE
481: EST = BIGNUM
482: ELSE
483: *
484: * Reordering successful
485: *
486: IF( WORK( 2, 1 ).EQ.ZERO ) THEN
487: *
488: * Form C = T22 - lambda*I in WORK(2:N,2:N).
489: *
490: DO 20 I = 2, N
491: WORK( I, I ) = WORK( I, I ) - WORK( 1, 1 )
492: 20 CONTINUE
493: N2 = 1
494: NN = N - 1
495: ELSE
496: *
497: * Triangularize the 2 by 2 block by unitary
498: * transformation U = [ cs i*ss ]
499: * [ i*ss cs ].
500: * such that the (1,1) position of WORK is complex
501: * eigenvalue lambda with positive imaginary part. (2,2)
502: * position of WORK is the complex eigenvalue lambda
503: * with negative imaginary part.
504: *
505: MU = SQRT( ABS( WORK( 1, 2 ) ) )*
506: $ SQRT( ABS( WORK( 2, 1 ) ) )
507: DELTA = DLAPY2( MU, WORK( 2, 1 ) )
508: CS = MU / DELTA
509: SN = -WORK( 2, 1 ) / DELTA
510: *
511: * Form
512: *
1.8 bertrand 513: * C**T = WORK(2:N,2:N) + i*[rwork(1) ..... rwork(n-1) ]
514: * [ mu ]
515: * [ .. ]
516: * [ .. ]
517: * [ mu ]
518: * where C**T is transpose of matrix C,
1.1 bertrand 519: * and RWORK is stored starting in the N+1-st column of
520: * WORK.
521: *
522: DO 30 J = 3, N
523: WORK( 2, J ) = CS*WORK( 2, J )
524: WORK( J, J ) = WORK( J, J ) - WORK( 1, 1 )
525: 30 CONTINUE
526: WORK( 2, 2 ) = ZERO
527: *
528: WORK( 1, N+1 ) = TWO*MU
529: DO 40 I = 2, N - 1
530: WORK( I, N+1 ) = SN*WORK( 1, I+1 )
531: 40 CONTINUE
532: N2 = 2
533: NN = 2*( N-1 )
534: END IF
535: *
1.8 bertrand 536: * Estimate norm(inv(C**T))
1.1 bertrand 537: *
538: EST = ZERO
539: KASE = 0
540: 50 CONTINUE
541: CALL DLACN2( NN, WORK( 1, N+2 ), WORK( 1, N+4 ), IWORK,
542: $ EST, KASE, ISAVE )
543: IF( KASE.NE.0 ) THEN
544: IF( KASE.EQ.1 ) THEN
545: IF( N2.EQ.1 ) THEN
546: *
1.8 bertrand 547: * Real eigenvalue: solve C**T*x = scale*c.
1.1 bertrand 548: *
549: CALL DLAQTR( .TRUE., .TRUE., N-1, WORK( 2, 2 ),
550: $ LDWORK, DUMMY, DUMM, SCALE,
551: $ WORK( 1, N+4 ), WORK( 1, N+6 ),
552: $ IERR )
553: ELSE
554: *
555: * Complex eigenvalue: solve
1.8 bertrand 556: * C**T*(p+iq) = scale*(c+id) in real arithmetic.
1.1 bertrand 557: *
558: CALL DLAQTR( .TRUE., .FALSE., N-1, WORK( 2, 2 ),
559: $ LDWORK, WORK( 1, N+1 ), MU, SCALE,
560: $ WORK( 1, N+4 ), WORK( 1, N+6 ),
561: $ IERR )
562: END IF
563: ELSE
564: IF( N2.EQ.1 ) THEN
565: *
566: * Real eigenvalue: solve C*x = scale*c.
567: *
568: CALL DLAQTR( .FALSE., .TRUE., N-1, WORK( 2, 2 ),
569: $ LDWORK, DUMMY, DUMM, SCALE,
570: $ WORK( 1, N+4 ), WORK( 1, N+6 ),
571: $ IERR )
572: ELSE
573: *
574: * Complex eigenvalue: solve
575: * C*(p+iq) = scale*(c+id) in real arithmetic.
576: *
577: CALL DLAQTR( .FALSE., .FALSE., N-1,
578: $ WORK( 2, 2 ), LDWORK,
579: $ WORK( 1, N+1 ), MU, SCALE,
580: $ WORK( 1, N+4 ), WORK( 1, N+6 ),
581: $ IERR )
582: *
583: END IF
584: END IF
585: *
586: GO TO 50
587: END IF
588: END IF
589: *
590: SEP( KS ) = SCALE / MAX( EST, SMLNUM )
591: IF( PAIR )
592: $ SEP( KS+1 ) = SEP( KS )
593: END IF
594: *
595: IF( PAIR )
596: $ KS = KS + 1
597: *
598: 60 CONTINUE
599: RETURN
600: *
601: * End of DTRSNA
602: *
603: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtrsna.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack