Annotation of rpl/lapack/lapack/dlasd8.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DLASD8( ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDDIFR,
! 2: $ DSIGMA, WORK, INFO )
! 3: *
! 4: * -- LAPACK auxiliary 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: * October 2006
! 8: *
! 9: * .. Scalar Arguments ..
! 10: INTEGER ICOMPQ, INFO, K, LDDIFR
! 11: * ..
! 12: * .. Array Arguments ..
! 13: DOUBLE PRECISION D( * ), DIFL( * ), DIFR( LDDIFR, * ),
! 14: $ DSIGMA( * ), VF( * ), VL( * ), WORK( * ),
! 15: $ Z( * )
! 16: * ..
! 17: *
! 18: * Purpose
! 19: * =======
! 20: *
! 21: * DLASD8 finds the square roots of the roots of the secular equation,
! 22: * as defined by the values in DSIGMA and Z. It makes the appropriate
! 23: * calls to DLASD4, and stores, for each element in D, the distance
! 24: * to its two nearest poles (elements in DSIGMA). It also updates
! 25: * the arrays VF and VL, the first and last components of all the
! 26: * right singular vectors of the original bidiagonal matrix.
! 27: *
! 28: * DLASD8 is called from DLASD6.
! 29: *
! 30: * Arguments
! 31: * =========
! 32: *
! 33: * ICOMPQ (input) INTEGER
! 34: * Specifies whether singular vectors are to be computed in
! 35: * factored form in the calling routine:
! 36: * = 0: Compute singular values only.
! 37: * = 1: Compute singular vectors in factored form as well.
! 38: *
! 39: * K (input) INTEGER
! 40: * The number of terms in the rational function to be solved
! 41: * by DLASD4. K >= 1.
! 42: *
! 43: * D (output) DOUBLE PRECISION array, dimension ( K )
! 44: * On output, D contains the updated singular values.
! 45: *
! 46: * Z (input/output) DOUBLE PRECISION array, dimension ( K )
! 47: * On entry, the first K elements of this array contain the
! 48: * components of the deflation-adjusted updating row vector.
! 49: * On exit, Z is updated.
! 50: *
! 51: * VF (input/output) DOUBLE PRECISION array, dimension ( K )
! 52: * On entry, VF contains information passed through DBEDE8.
! 53: * On exit, VF contains the first K components of the first
! 54: * components of all right singular vectors of the bidiagonal
! 55: * matrix.
! 56: *
! 57: * VL (input/output) DOUBLE PRECISION array, dimension ( K )
! 58: * On entry, VL contains information passed through DBEDE8.
! 59: * On exit, VL contains the first K components of the last
! 60: * components of all right singular vectors of the bidiagonal
! 61: * matrix.
! 62: *
! 63: * DIFL (output) DOUBLE PRECISION array, dimension ( K )
! 64: * On exit, DIFL(I) = D(I) - DSIGMA(I).
! 65: *
! 66: * DIFR (output) DOUBLE PRECISION array,
! 67: * dimension ( LDDIFR, 2 ) if ICOMPQ = 1 and
! 68: * dimension ( K ) if ICOMPQ = 0.
! 69: * On exit, DIFR(I,1) = D(I) - DSIGMA(I+1), DIFR(K,1) is not
! 70: * defined and will not be referenced.
! 71: *
! 72: * If ICOMPQ = 1, DIFR(1:K,2) is an array containing the
! 73: * normalizing factors for the right singular vector matrix.
! 74: *
! 75: * LDDIFR (input) INTEGER
! 76: * The leading dimension of DIFR, must be at least K.
! 77: *
! 78: * DSIGMA (input/output) DOUBLE PRECISION array, dimension ( K )
! 79: * On entry, the first K elements of this array contain the old
! 80: * roots of the deflated updating problem. These are the poles
! 81: * of the secular equation.
! 82: * On exit, the elements of DSIGMA may be very slightly altered
! 83: * in value.
! 84: *
! 85: * WORK (workspace) DOUBLE PRECISION array, dimension at least 3 * K
! 86: *
! 87: * INFO (output) INTEGER
! 88: * = 0: successful exit.
! 89: * < 0: if INFO = -i, the i-th argument had an illegal value.
! 90: * > 0: if INFO = 1, an singular value did not converge
! 91: *
! 92: * Further Details
! 93: * ===============
! 94: *
! 95: * Based on contributions by
! 96: * Ming Gu and Huan Ren, Computer Science Division, University of
! 97: * California at Berkeley, USA
! 98: *
! 99: * =====================================================================
! 100: *
! 101: * .. Parameters ..
! 102: DOUBLE PRECISION ONE
! 103: PARAMETER ( ONE = 1.0D+0 )
! 104: * ..
! 105: * .. Local Scalars ..
! 106: INTEGER I, IWK1, IWK2, IWK2I, IWK3, IWK3I, J
! 107: DOUBLE PRECISION DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, RHO, TEMP
! 108: * ..
! 109: * .. External Subroutines ..
! 110: EXTERNAL DCOPY, DLASCL, DLASD4, DLASET, XERBLA
! 111: * ..
! 112: * .. External Functions ..
! 113: DOUBLE PRECISION DDOT, DLAMC3, DNRM2
! 114: EXTERNAL DDOT, DLAMC3, DNRM2
! 115: * ..
! 116: * .. Intrinsic Functions ..
! 117: INTRINSIC ABS, SIGN, SQRT
! 118: * ..
! 119: * .. Executable Statements ..
! 120: *
! 121: * Test the input parameters.
! 122: *
! 123: INFO = 0
! 124: *
! 125: IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
! 126: INFO = -1
! 127: ELSE IF( K.LT.1 ) THEN
! 128: INFO = -2
! 129: ELSE IF( LDDIFR.LT.K ) THEN
! 130: INFO = -9
! 131: END IF
! 132: IF( INFO.NE.0 ) THEN
! 133: CALL XERBLA( 'DLASD8', -INFO )
! 134: RETURN
! 135: END IF
! 136: *
! 137: * Quick return if possible
! 138: *
! 139: IF( K.EQ.1 ) THEN
! 140: D( 1 ) = ABS( Z( 1 ) )
! 141: DIFL( 1 ) = D( 1 )
! 142: IF( ICOMPQ.EQ.1 ) THEN
! 143: DIFL( 2 ) = ONE
! 144: DIFR( 1, 2 ) = ONE
! 145: END IF
! 146: RETURN
! 147: END IF
! 148: *
! 149: * Modify values DSIGMA(i) to make sure all DSIGMA(i)-DSIGMA(j) can
! 150: * be computed with high relative accuracy (barring over/underflow).
! 151: * This is a problem on machines without a guard digit in
! 152: * add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2).
! 153: * The following code replaces DSIGMA(I) by 2*DSIGMA(I)-DSIGMA(I),
! 154: * which on any of these machines zeros out the bottommost
! 155: * bit of DSIGMA(I) if it is 1; this makes the subsequent
! 156: * subtractions DSIGMA(I)-DSIGMA(J) unproblematic when cancellation
! 157: * occurs. On binary machines with a guard digit (almost all
! 158: * machines) it does not change DSIGMA(I) at all. On hexadecimal
! 159: * and decimal machines with a guard digit, it slightly
! 160: * changes the bottommost bits of DSIGMA(I). It does not account
! 161: * for hexadecimal or decimal machines without guard digits
! 162: * (we know of none). We use a subroutine call to compute
! 163: * 2*DLAMBDA(I) to prevent optimizing compilers from eliminating
! 164: * this code.
! 165: *
! 166: DO 10 I = 1, K
! 167: DSIGMA( I ) = DLAMC3( DSIGMA( I ), DSIGMA( I ) ) - DSIGMA( I )
! 168: 10 CONTINUE
! 169: *
! 170: * Book keeping.
! 171: *
! 172: IWK1 = 1
! 173: IWK2 = IWK1 + K
! 174: IWK3 = IWK2 + K
! 175: IWK2I = IWK2 - 1
! 176: IWK3I = IWK3 - 1
! 177: *
! 178: * Normalize Z.
! 179: *
! 180: RHO = DNRM2( K, Z, 1 )
! 181: CALL DLASCL( 'G', 0, 0, RHO, ONE, K, 1, Z, K, INFO )
! 182: RHO = RHO*RHO
! 183: *
! 184: * Initialize WORK(IWK3).
! 185: *
! 186: CALL DLASET( 'A', K, 1, ONE, ONE, WORK( IWK3 ), K )
! 187: *
! 188: * Compute the updated singular values, the arrays DIFL, DIFR,
! 189: * and the updated Z.
! 190: *
! 191: DO 40 J = 1, K
! 192: CALL DLASD4( K, J, DSIGMA, Z, WORK( IWK1 ), RHO, D( J ),
! 193: $ WORK( IWK2 ), INFO )
! 194: *
! 195: * If the root finder fails, the computation is terminated.
! 196: *
! 197: IF( INFO.NE.0 ) THEN
! 198: RETURN
! 199: END IF
! 200: WORK( IWK3I+J ) = WORK( IWK3I+J )*WORK( J )*WORK( IWK2I+J )
! 201: DIFL( J ) = -WORK( J )
! 202: DIFR( J, 1 ) = -WORK( J+1 )
! 203: DO 20 I = 1, J - 1
! 204: WORK( IWK3I+I ) = WORK( IWK3I+I )*WORK( I )*
! 205: $ WORK( IWK2I+I ) / ( DSIGMA( I )-
! 206: $ DSIGMA( J ) ) / ( DSIGMA( I )+
! 207: $ DSIGMA( J ) )
! 208: 20 CONTINUE
! 209: DO 30 I = J + 1, K
! 210: WORK( IWK3I+I ) = WORK( IWK3I+I )*WORK( I )*
! 211: $ WORK( IWK2I+I ) / ( DSIGMA( I )-
! 212: $ DSIGMA( J ) ) / ( DSIGMA( I )+
! 213: $ DSIGMA( J ) )
! 214: 30 CONTINUE
! 215: 40 CONTINUE
! 216: *
! 217: * Compute updated Z.
! 218: *
! 219: DO 50 I = 1, K
! 220: Z( I ) = SIGN( SQRT( ABS( WORK( IWK3I+I ) ) ), Z( I ) )
! 221: 50 CONTINUE
! 222: *
! 223: * Update VF and VL.
! 224: *
! 225: DO 80 J = 1, K
! 226: DIFLJ = DIFL( J )
! 227: DJ = D( J )
! 228: DSIGJ = -DSIGMA( J )
! 229: IF( J.LT.K ) THEN
! 230: DIFRJ = -DIFR( J, 1 )
! 231: DSIGJP = -DSIGMA( J+1 )
! 232: END IF
! 233: WORK( J ) = -Z( J ) / DIFLJ / ( DSIGMA( J )+DJ )
! 234: DO 60 I = 1, J - 1
! 235: WORK( I ) = Z( I ) / ( DLAMC3( DSIGMA( I ), DSIGJ )-DIFLJ )
! 236: $ / ( DSIGMA( I )+DJ )
! 237: 60 CONTINUE
! 238: DO 70 I = J + 1, K
! 239: WORK( I ) = Z( I ) / ( DLAMC3( DSIGMA( I ), DSIGJP )+DIFRJ )
! 240: $ / ( DSIGMA( I )+DJ )
! 241: 70 CONTINUE
! 242: TEMP = DNRM2( K, WORK, 1 )
! 243: WORK( IWK2I+J ) = DDOT( K, WORK, 1, VF, 1 ) / TEMP
! 244: WORK( IWK3I+J ) = DDOT( K, WORK, 1, VL, 1 ) / TEMP
! 245: IF( ICOMPQ.EQ.1 ) THEN
! 246: DIFR( J, 2 ) = TEMP
! 247: END IF
! 248: 80 CONTINUE
! 249: *
! 250: CALL DCOPY( K, WORK( IWK2 ), 1, VF, 1 )
! 251: CALL DCOPY( K, WORK( IWK3 ), 1, VL, 1 )
! 252: *
! 253: RETURN
! 254: *
! 255: * End of DLASD8
! 256: *
! 257: END
! 258:
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlasd8.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack