github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/internal/testdata/dlasqtest/dlasq6.f (about) 1 *> \brief \b DLASQ6 computes one dqd transform in ping-pong form. Used by sbdsqr and sstegr. 2 * 3 * =========== DOCUMENTATION =========== 4 * 5 * Online html documentation available at 6 * http://www.netlib.org/lapack/explore-html/ 7 * 8 *> \htmlonly 9 *> Download DLASQ6 + dependencies 10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq6.f"> 11 *> [TGZ]</a> 12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq6.f"> 13 *> [ZIP]</a> 14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq6.f"> 15 *> [TXT]</a> 16 *> \endhtmlonly 17 * 18 * Definition: 19 * =========== 20 * 21 * SUBROUTINE DLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, 22 * DNM1, DNM2 ) 23 * 24 * .. Scalar Arguments .. 25 * INTEGER I0, N0, PP 26 * DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DNM1, DNM2 27 * .. 28 * .. Array Arguments .. 29 * DOUBLE PRECISION Z( * ) 30 * .. 31 * 32 * 33 *> \par Purpose: 34 * ============= 35 *> 36 *> \verbatim 37 *> 38 *> DLASQ6 computes one dqd (shift equal to zero) transform in 39 *> ping-pong form, with protection against underflow and overflow. 40 *> \endverbatim 41 * 42 * Arguments: 43 * ========== 44 * 45 *> \param[in] I0 46 *> \verbatim 47 *> I0 is INTEGER 48 *> First index. 49 *> \endverbatim 50 *> 51 *> \param[in] N0 52 *> \verbatim 53 *> N0 is INTEGER 54 *> Last index. 55 *> \endverbatim 56 *> 57 *> \param[in] Z 58 *> \verbatim 59 *> Z is DOUBLE PRECISION array, dimension ( 4*N ) 60 *> Z holds the qd array. EMIN is stored in Z(4*N0) to avoid 61 *> an extra argument. 62 *> \endverbatim 63 *> 64 *> \param[in] PP 65 *> \verbatim 66 *> PP is INTEGER 67 *> PP=0 for ping, PP=1 for pong. 68 *> \endverbatim 69 *> 70 *> \param[out] DMIN 71 *> \verbatim 72 *> DMIN is DOUBLE PRECISION 73 *> Minimum value of d. 74 *> \endverbatim 75 *> 76 *> \param[out] DMIN1 77 *> \verbatim 78 *> DMIN1 is DOUBLE PRECISION 79 *> Minimum value of d, excluding D( N0 ). 80 *> \endverbatim 81 *> 82 *> \param[out] DMIN2 83 *> \verbatim 84 *> DMIN2 is DOUBLE PRECISION 85 *> Minimum value of d, excluding D( N0 ) and D( N0-1 ). 86 *> \endverbatim 87 *> 88 *> \param[out] DN 89 *> \verbatim 90 *> DN is DOUBLE PRECISION 91 *> d(N0), the last value of d. 92 *> \endverbatim 93 *> 94 *> \param[out] DNM1 95 *> \verbatim 96 *> DNM1 is DOUBLE PRECISION 97 *> d(N0-1). 98 *> \endverbatim 99 *> 100 *> \param[out] DNM2 101 *> \verbatim 102 *> DNM2 is DOUBLE PRECISION 103 *> d(N0-2). 104 *> \endverbatim 105 * 106 * Authors: 107 * ======== 108 * 109 *> \author Univ. of Tennessee 110 *> \author Univ. of California Berkeley 111 *> \author Univ. of Colorado Denver 112 *> \author NAG Ltd. 113 * 114 *> \date September 2012 115 * 116 *> \ingroup auxOTHERcomputational 117 * 118 * ===================================================================== 119 SUBROUTINE DLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, 120 $ DNM1, DNM2 ) 121 * 122 * -- LAPACK computational routine (version 3.4.2) -- 123 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 124 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 125 * September 2012 126 * 127 * .. Scalar Arguments .. 128 INTEGER I0, N0, PP 129 DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DNM1, DNM2 130 * .. 131 * .. Array Arguments .. 132 DOUBLE PRECISION Z( * ) 133 * .. 134 * 135 * ===================================================================== 136 * 137 * .. Parameter .. 138 DOUBLE PRECISION ZERO 139 PARAMETER ( ZERO = 0.0D0 ) 140 * .. 141 * .. Local Scalars .. 142 INTEGER J4, J4P2 143 DOUBLE PRECISION D, EMIN, SAFMIN, TEMP 144 * .. 145 * .. External Function .. 146 DOUBLE PRECISION DLAMCH 147 EXTERNAL DLAMCH 148 * .. 149 * .. Intrinsic Functions .. 150 INTRINSIC MIN 151 * .. 152 * .. Executable Statements .. 153 * 154 IF( ( N0-I0-1 ).LE.0 ) 155 $ RETURN 156 * 157 158 print *, "In dlasq6" 159 STOP 160 161 SAFMIN = DLAMCH( 'Safe minimum' ) 162 J4 = 4*I0 + PP - 3 163 EMIN = Z( J4+4 ) 164 D = Z( J4 ) 165 DMIN = D 166 * 167 IF( PP.EQ.0 ) THEN 168 DO 10 J4 = 4*I0, 4*( N0-3 ), 4 169 Z( J4-2 ) = D + Z( J4-1 ) 170 IF( Z( J4-2 ).EQ.ZERO ) THEN 171 Z( J4 ) = ZERO 172 D = Z( J4+1 ) 173 DMIN = D 174 EMIN = ZERO 175 ELSE IF( SAFMIN*Z( J4+1 ).LT.Z( J4-2 ) .AND. 176 $ SAFMIN*Z( J4-2 ).LT.Z( J4+1 ) ) THEN 177 TEMP = Z( J4+1 ) / Z( J4-2 ) 178 Z( J4 ) = Z( J4-1 )*TEMP 179 D = D*TEMP 180 ELSE 181 Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) ) 182 D = Z( J4+1 )*( D / Z( J4-2 ) ) 183 END IF 184 DMIN = MIN( DMIN, D ) 185 EMIN = MIN( EMIN, Z( J4 ) ) 186 10 CONTINUE 187 ELSE 188 DO 20 J4 = 4*I0, 4*( N0-3 ), 4 189 Z( J4-3 ) = D + Z( J4 ) 190 IF( Z( J4-3 ).EQ.ZERO ) THEN 191 Z( J4-1 ) = ZERO 192 D = Z( J4+2 ) 193 DMIN = D 194 EMIN = ZERO 195 ELSE IF( SAFMIN*Z( J4+2 ).LT.Z( J4-3 ) .AND. 196 $ SAFMIN*Z( J4-3 ).LT.Z( J4+2 ) ) THEN 197 TEMP = Z( J4+2 ) / Z( J4-3 ) 198 Z( J4-1 ) = Z( J4 )*TEMP 199 D = D*TEMP 200 ELSE 201 Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) ) 202 D = Z( J4+2 )*( D / Z( J4-3 ) ) 203 END IF 204 DMIN = MIN( DMIN, D ) 205 EMIN = MIN( EMIN, Z( J4-1 ) ) 206 20 CONTINUE 207 END IF 208 * 209 * Unroll last two steps. 210 * 211 DNM2 = D 212 DMIN2 = DMIN 213 J4 = 4*( N0-2 ) - PP 214 J4P2 = J4 + 2*PP - 1 215 Z( J4-2 ) = DNM2 + Z( J4P2 ) 216 IF( Z( J4-2 ).EQ.ZERO ) THEN 217 Z( J4 ) = ZERO 218 DNM1 = Z( J4P2+2 ) 219 DMIN = DNM1 220 EMIN = ZERO 221 ELSE IF( SAFMIN*Z( J4P2+2 ).LT.Z( J4-2 ) .AND. 222 $ SAFMIN*Z( J4-2 ).LT.Z( J4P2+2 ) ) THEN 223 TEMP = Z( J4P2+2 ) / Z( J4-2 ) 224 Z( J4 ) = Z( J4P2 )*TEMP 225 DNM1 = DNM2*TEMP 226 ELSE 227 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 228 DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) 229 END IF 230 DMIN = MIN( DMIN, DNM1 ) 231 * 232 DMIN1 = DMIN 233 J4 = J4 + 4 234 J4P2 = J4 + 2*PP - 1 235 Z( J4-2 ) = DNM1 + Z( J4P2 ) 236 IF( Z( J4-2 ).EQ.ZERO ) THEN 237 Z( J4 ) = ZERO 238 DN = Z( J4P2+2 ) 239 DMIN = DN 240 EMIN = ZERO 241 ELSE IF( SAFMIN*Z( J4P2+2 ).LT.Z( J4-2 ) .AND. 242 $ SAFMIN*Z( J4-2 ).LT.Z( J4P2+2 ) ) THEN 243 TEMP = Z( J4P2+2 ) / Z( J4-2 ) 244 Z( J4 ) = Z( J4P2 )*TEMP 245 DN = DNM1*TEMP 246 ELSE 247 Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) 248 DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) 249 END IF 250 DMIN = MIN( DMIN, DN ) 251 * 252 Z( J4+2 ) = DN 253 Z( 4*N0-PP ) = EMIN 254 RETURN 255 * 256 * End of DLASQ6 257 * 258 END