1 SUBROUTINE DLASQ5( I0, N0, Z, PP, TAU, DMIN, DMIN1, DMIN2, DN, 2 $ DNM1, DNM2, IEEE ) 3* 4* -- LAPACK auxiliary routine (version 3.0) -- 5* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 6* Courant Institute, Argonne National Lab, and Rice University 7* May 17, 2000 8* 9* .. Scalar Arguments .. 10 LOGICAL IEEE 11 INTEGER I0, N0, PP 12 DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU 13* .. 14* .. Array Arguments .. 15 DOUBLE PRECISION Z( * ) 16* .. 17* 18* Purpose 19* ======= 20* 21* DLASQ5 computes one dqds transform in ping-pong form, one 22* version for IEEE machines another for non IEEE machines. 23* 24* Arguments 25* ========= 26* 27* I0 (input) INTEGER 28* First index. 29* 30* N0 (input) INTEGER 31* Last index. 32* 33* Z (input) DOUBLE PRECISION array, dimension ( 4*N ) 34* Z holds the qd array. EMIN is stored in Z(4*N0) to avoid 35* an extra argument. 36* 37* PP (input) INTEGER 38* PP=0 for ping, PP=1 for pong. 39* 40* TAU (input) DOUBLE PRECISION 41* This is the shift. 42* 43* DMIN (output) DOUBLE PRECISION 44* Minimum value of d. 45* 46* DMIN1 (output) DOUBLE PRECISION 47* Minimum value of d, excluding D( N0 ). 48* 49* DMIN2 (output) DOUBLE PRECISION 50* Minimum value of d, excluding D( N0 ) and D( N0-1 ). 51* 52* DN (output) DOUBLE PRECISION 53* d(N0), the last value of d. 54* 55* DNM1 (output) DOUBLE PRECISION 56* d(N0-1). 57* 58* DNM2 (output) DOUBLE PRECISION 59* d(N0-2). 60* 61* IEEE (input) LOGICAL 62* Flag for IEEE or non IEEE arithmetic. 63* 64* ===================================================================== 65* 66* .. Parameter .. 67 DOUBLE PRECISION ZERO 68 PARAMETER ( ZERO = 0.0D0 ) 69* .. 70* .. Local Scalars .. 71 INTEGER J4, J4P2 72 DOUBLE PRECISION D, EMIN, TEMP 73* .. 74* .. Intrinsic Functions .. 75 INTRINSIC MIN 76* .. 77* .. Executable Statements .. 78* 79 IF( ( N0-I0-1 ).LE.0 ) 80 $ RETURN 81* 82 J4 = 4*I0 + PP - 3 83 EMIN = Z( J4+4 ) 84 D = Z( J4 ) - TAU 85 DMIN = D 86 DMIN1 = -Z( J4 ) 87* 88 IF( IEEE ) THEN 89* 90* Code for IEEE arithmetic. 91* 92 IF( PP.EQ.0 ) THEN 93 DO 10 J4 = 4*I0, 4*( N0-3 ), 4 94 Z( J4-2 ) = D + Z( J4-1 ) 95 TEMP = Z( J4+1 ) / Z( J4-2 ) 96 D = D*TEMP - TAU 97 DMIN = MIN( DMIN, D ) 98 Z( J4 ) = Z( J4-1 )*TEMP 99 EMIN = MIN( Z( J4 ), EMIN ) 100 10 CONTINUE 101 ELSE 102 DO 20 J4 = 4*I0, 4*( N0-3 ), 4 103 Z( J4-3 ) = D + Z( J4 ) 104 TEMP = Z( J4+2 ) / Z( J4-3 ) 105 D = D*TEMP - TAU 106 DMIN = MIN( DMIN, D ) 107 Z( J4-1 ) = Z( J4 )*TEMP 108 EMIN = MIN( Z( J4-1 ), EMIN ) 109 20 CONTINUE 110 END IF 111* 112* Unroll last two steps. 113* 114 DNM2 = D 115 DMIN2 = DMIN 116 J4 = 4*( N0-2 ) - PP 117 J4P2 = J4 + 2*PP - 1 118 Z( J4-2 ) = DNM2 + Z( J4P2 ) 119 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 120 DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU 121 DMIN = MIN( DMIN, DNM1 ) 122* 123 DMIN1 = DMIN 124 J4 = J4 + 4 125 J4P2 = J4 + 2*PP - 1 126 Z( J4-2 ) = DNM1 + Z( J4P2 ) 127 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 128 DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU 129 DMIN = MIN( DMIN, DN ) 130* 131 ELSE 132* 133* Code for non IEEE arithmetic. 134* 135 IF( PP.EQ.0 ) THEN 136 DO 30 J4 = 4*I0, 4*( N0-3 ), 4 137 Z( J4-2 ) = D + Z( J4-1 ) 138 IF( D.LT.ZERO ) THEN 139 RETURN 140 ELSE 141 Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) ) 142 D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU 143 END IF 144 DMIN = MIN( DMIN, D ) 145 EMIN = MIN( EMIN, Z( J4 ) ) 146 30 CONTINUE 147 ELSE 148 DO 40 J4 = 4*I0, 4*( N0-3 ), 4 149 Z( J4-3 ) = D + Z( J4 ) 150 IF( D.LT.ZERO ) THEN 151 RETURN 152 ELSE 153 Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) ) 154 D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU 155 END IF 156 DMIN = MIN( DMIN, D ) 157 EMIN = MIN( EMIN, Z( J4-1 ) ) 158 40 CONTINUE 159 END IF 160* 161* Unroll last two steps. 162* 163 DNM2 = D 164 DMIN2 = DMIN 165 J4 = 4*( N0-2 ) - PP 166 J4P2 = J4 + 2*PP - 1 167 Z( J4-2 ) = DNM2 + Z( J4P2 ) 168 IF( DNM2.LT.ZERO ) THEN 169 RETURN 170 ELSE 171 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 172 DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU 173 END IF 174 DMIN = MIN( DMIN, DNM1 ) 175* 176 DMIN1 = DMIN 177 J4 = J4 + 4 178 J4P2 = J4 + 2*PP - 1 179 Z( J4-2 ) = DNM1 + Z( J4P2 ) 180 IF( DNM1.LT.ZERO ) THEN 181 RETURN 182 ELSE 183 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 184 DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU 185 END IF 186 DMIN = MIN( DMIN, DN ) 187* 188 END IF 189* 190 Z( J4+2 ) = DN 191 Z( 4*N0-PP ) = EMIN 192 RETURN 193* 194* End of DLASQ5 195* 196 END 197