Annotation of rpl/lapack/lapack/dlaed9.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DLAED9( K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W,
! 2: $ S, LDS, INFO )
! 3: *
! 4: * -- LAPACK routine (version 3.2) --
! 5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 7: * November 2006
! 8: *
! 9: * .. Scalar Arguments ..
! 10: INTEGER INFO, K, KSTART, KSTOP, LDQ, LDS, N
! 11: DOUBLE PRECISION RHO
! 12: * ..
! 13: * .. Array Arguments ..
! 14: DOUBLE PRECISION D( * ), DLAMDA( * ), Q( LDQ, * ), S( LDS, * ),
! 15: $ W( * )
! 16: * ..
! 17: *
! 18: * Purpose
! 19: * =======
! 20: *
! 21: * DLAED9 finds the roots of the secular equation, as defined by the
! 22: * values in D, Z, and RHO, between KSTART and KSTOP. It makes the
! 23: * appropriate calls to DLAED4 and then stores the new matrix of
! 24: * eigenvectors for use in calculating the next level of Z vectors.
! 25: *
! 26: * Arguments
! 27: * =========
! 28: *
! 29: * K (input) INTEGER
! 30: * The number of terms in the rational function to be solved by
! 31: * DLAED4. K >= 0.
! 32: *
! 33: * KSTART (input) INTEGER
! 34: * KSTOP (input) INTEGER
! 35: * The updated eigenvalues Lambda(I), KSTART <= I <= KSTOP
! 36: * are to be computed. 1 <= KSTART <= KSTOP <= K.
! 37: *
! 38: * N (input) INTEGER
! 39: * The number of rows and columns in the Q matrix.
! 40: * N >= K (delation may result in N > K).
! 41: *
! 42: * D (output) DOUBLE PRECISION array, dimension (N)
! 43: * D(I) contains the updated eigenvalues
! 44: * for KSTART <= I <= KSTOP.
! 45: *
! 46: * Q (workspace) DOUBLE PRECISION array, dimension (LDQ,N)
! 47: *
! 48: * LDQ (input) INTEGER
! 49: * The leading dimension of the array Q. LDQ >= max( 1, N ).
! 50: *
! 51: * RHO (input) DOUBLE PRECISION
! 52: * The value of the parameter in the rank one update equation.
! 53: * RHO >= 0 required.
! 54: *
! 55: * DLAMDA (input) DOUBLE PRECISION array, dimension (K)
! 56: * The first K elements of this array contain the old roots
! 57: * of the deflated updating problem. These are the poles
! 58: * of the secular equation.
! 59: *
! 60: * W (input) DOUBLE PRECISION array, dimension (K)
! 61: * The first K elements of this array contain the components
! 62: * of the deflation-adjusted updating vector.
! 63: *
! 64: * S (output) DOUBLE PRECISION array, dimension (LDS, K)
! 65: * Will contain the eigenvectors of the repaired matrix which
! 66: * will be stored for subsequent Z vector calculation and
! 67: * multiplied by the previously accumulated eigenvectors
! 68: * to update the system.
! 69: *
! 70: * LDS (input) INTEGER
! 71: * The leading dimension of S. LDS >= max( 1, K ).
! 72: *
! 73: * INFO (output) INTEGER
! 74: * = 0: successful exit.
! 75: * < 0: if INFO = -i, the i-th argument had an illegal value.
! 76: * > 0: if INFO = 1, an eigenvalue did not converge
! 77: *
! 78: * Further Details
! 79: * ===============
! 80: *
! 81: * Based on contributions by
! 82: * Jeff Rutter, Computer Science Division, University of California
! 83: * at Berkeley, USA
! 84: *
! 85: * =====================================================================
! 86: *
! 87: * .. Local Scalars ..
! 88: INTEGER I, J
! 89: DOUBLE PRECISION TEMP
! 90: * ..
! 91: * .. External Functions ..
! 92: DOUBLE PRECISION DLAMC3, DNRM2
! 93: EXTERNAL DLAMC3, DNRM2
! 94: * ..
! 95: * .. External Subroutines ..
! 96: EXTERNAL DCOPY, DLAED4, XERBLA
! 97: * ..
! 98: * .. Intrinsic Functions ..
! 99: INTRINSIC MAX, SIGN, SQRT
! 100: * ..
! 101: * .. Executable Statements ..
! 102: *
! 103: * Test the input parameters.
! 104: *
! 105: INFO = 0
! 106: *
! 107: IF( K.LT.0 ) THEN
! 108: INFO = -1
! 109: ELSE IF( KSTART.LT.1 .OR. KSTART.GT.MAX( 1, K ) ) THEN
! 110: INFO = -2
! 111: ELSE IF( MAX( 1, KSTOP ).LT.KSTART .OR. KSTOP.GT.MAX( 1, K ) )
! 112: $ THEN
! 113: INFO = -3
! 114: ELSE IF( N.LT.K ) THEN
! 115: INFO = -4
! 116: ELSE IF( LDQ.LT.MAX( 1, K ) ) THEN
! 117: INFO = -7
! 118: ELSE IF( LDS.LT.MAX( 1, K ) ) THEN
! 119: INFO = -12
! 120: END IF
! 121: IF( INFO.NE.0 ) THEN
! 122: CALL XERBLA( 'DLAED9', -INFO )
! 123: RETURN
! 124: END IF
! 125: *
! 126: * Quick return if possible
! 127: *
! 128: IF( K.EQ.0 )
! 129: $ RETURN
! 130: *
! 131: * Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can
! 132: * be computed with high relative accuracy (barring over/underflow).
! 133: * This is a problem on machines without a guard digit in
! 134: * add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
! 135: * The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I),
! 136: * which on any of these machines zeros out the bottommost
! 137: * bit of DLAMDA(I) if it is 1; this makes the subsequent
! 138: * subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation
! 139: * occurs. On binary machines with a guard digit (almost all
! 140: * machines) it does not change DLAMDA(I) at all. On hexadecimal
! 141: * and decimal machines with a guard digit, it slightly
! 142: * changes the bottommost bits of DLAMDA(I). It does not account
! 143: * for hexadecimal or decimal machines without guard digits
! 144: * (we know of none). We use a subroutine call to compute
! 145: * 2*DLAMBDA(I) to prevent optimizing compilers from eliminating
! 146: * this code.
! 147: *
! 148: DO 10 I = 1, N
! 149: DLAMDA( I ) = DLAMC3( DLAMDA( I ), DLAMDA( I ) ) - DLAMDA( I )
! 150: 10 CONTINUE
! 151: *
! 152: DO 20 J = KSTART, KSTOP
! 153: CALL DLAED4( K, J, DLAMDA, W, Q( 1, J ), RHO, D( J ), INFO )
! 154: *
! 155: * If the zero finder fails, the computation is terminated.
! 156: *
! 157: IF( INFO.NE.0 )
! 158: $ GO TO 120
! 159: 20 CONTINUE
! 160: *
! 161: IF( K.EQ.1 .OR. K.EQ.2 ) THEN
! 162: DO 40 I = 1, K
! 163: DO 30 J = 1, K
! 164: S( J, I ) = Q( J, I )
! 165: 30 CONTINUE
! 166: 40 CONTINUE
! 167: GO TO 120
! 168: END IF
! 169: *
! 170: * Compute updated W.
! 171: *
! 172: CALL DCOPY( K, W, 1, S, 1 )
! 173: *
! 174: * Initialize W(I) = Q(I,I)
! 175: *
! 176: CALL DCOPY( K, Q, LDQ+1, W, 1 )
! 177: DO 70 J = 1, K
! 178: DO 50 I = 1, J - 1
! 179: W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
! 180: 50 CONTINUE
! 181: DO 60 I = J + 1, K
! 182: W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) )
! 183: 60 CONTINUE
! 184: 70 CONTINUE
! 185: DO 80 I = 1, K
! 186: W( I ) = SIGN( SQRT( -W( I ) ), S( I, 1 ) )
! 187: 80 CONTINUE
! 188: *
! 189: * Compute eigenvectors of the modified rank-1 modification.
! 190: *
! 191: DO 110 J = 1, K
! 192: DO 90 I = 1, K
! 193: Q( I, J ) = W( I ) / Q( I, J )
! 194: 90 CONTINUE
! 195: TEMP = DNRM2( K, Q( 1, J ), 1 )
! 196: DO 100 I = 1, K
! 197: S( I, J ) = Q( I, J ) / TEMP
! 198: 100 CONTINUE
! 199: 110 CONTINUE
! 200: *
! 201: 120 CONTINUE
! 202: RETURN
! 203: *
! 204: * End of DLAED9
! 205: *
! 206: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaed9.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack