SRC/dlasv2.f

   1 *> \brief \b DLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DLASV2 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasv2.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasv2.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasv2.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DLASV2( F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       DOUBLE PRECISION   CSL, CSR, F, G, H, SNL, SNR, SSMAX, SSMIN
  25 *       ..
  26 *
  27 *
  28 *> \par Purpose:
  29 *  =============
  30 *>
  31 *> \verbatim
  32 *>
  33 *> DLASV2 computes the singular value decomposition of a 2-by-2
  34 *> triangular matrix
  35 *>    [  F   G  ]
  36 *>    [  0   H  ].
  37 *> On return, abs(SSMAX) is the larger singular value, abs(SSMIN) is the
  38 *> smaller singular value, and (CSL,SNL) and (CSR,SNR) are the left and
  39 *> right singular vectors for abs(SSMAX), giving the decomposition
  40 *>
  41 *>    [ CSL  SNL ] [  F   G  ] [ CSR -SNR ]  =  [ SSMAX   0   ]
  42 *>    [-SNL  CSL ] [  0   H  ] [ SNR  CSR ]     [  0    SSMIN ].
  43 *> \endverbatim
  44 *
  45 *  Arguments:
  46 *  ==========
  47 *
  48 *> \param[in] F
  49 *> \verbatim
  50 *>          F is DOUBLE PRECISION
  51 *>          The (1,1) element of the 2-by-2 matrix.
  52 *> \endverbatim
  53 *>
  54 *> \param[in] G
  55 *> \verbatim
  56 *>          G is DOUBLE PRECISION
  57 *>          The (1,2) element of the 2-by-2 matrix.
  58 *> \endverbatim
  59 *>
  60 *> \param[in] H
  61 *> \verbatim
  62 *>          H is DOUBLE PRECISION
  63 *>          The (2,2) element of the 2-by-2 matrix.
  64 *> \endverbatim
  65 *>
  66 *> \param[out] SSMIN
  67 *> \verbatim
  68 *>          SSMIN is DOUBLE PRECISION
  69 *>          abs(SSMIN) is the smaller singular value.
  70 *> \endverbatim
  71 *>
  72 *> \param[out] SSMAX
  73 *> \verbatim
  74 *>          SSMAX is DOUBLE PRECISION
  75 *>          abs(SSMAX) is the larger singular value.
  76 *> \endverbatim
  77 *>
  78 *> \param[out] SNL
  79 *> \verbatim
  80 *>          SNL is DOUBLE PRECISION
  81 *> \endverbatim
  82 *>
  83 *> \param[out] CSL
  84 *> \verbatim
  85 *>          CSL is DOUBLE PRECISION
  86 *>          The vector (CSL, SNL) is a unit left singular vector for the
  87 *>          singular value abs(SSMAX).
  88 *> \endverbatim
  89 *>
  90 *> \param[out] SNR
  91 *> \verbatim
  92 *>          SNR is DOUBLE PRECISION
  93 *> \endverbatim
  94 *>
  95 *> \param[out] CSR
  96 *> \verbatim
  97 *>          CSR is DOUBLE PRECISION
  98 *>          The vector (CSR, SNR) is a unit right singular vector for the
  99 *>          singular value abs(SSMAX).
 100 *> \endverbatim
 101 *
 102 *  Authors:
 103 *  ========
 104 *
 105 *> \author Univ. of Tennessee
 106 *> \author Univ. of California Berkeley
 107 *> \author Univ. of Colorado Denver
 108 *> \author NAG Ltd.
 109 *
 110 *> \date September 2012
 111 *
 112 *> \ingroup auxOTHERauxiliary
 113 *
 114 *> \par Further Details:
 115 *  =====================
 116 *>
 117 *> \verbatim
 118 *>
 119 *>  Any input parameter may be aliased with any output parameter.
 120 *>
 121 *>  Barring over/underflow and assuming a guard digit in subtraction, all
 122 *>  output quantities are correct to within a few units in the last
 123 *>  place (ulps).
 124 *>
 125 *>  In IEEE arithmetic, the code works correctly if one matrix element is
 126 *>  infinite.
 127 *>
 128 *>  Overflow will not occur unless the largest singular value itself
 129 *>  overflows or is within a few ulps of overflow. (On machines with
 130 *>  partial overflow, like the Cray, overflow may occur if the largest
 131 *>  singular value is within a factor of 2 of overflow.)
 132 *>
 133 *>  Underflow is harmless if underflow is gradual. Otherwise, results
 134 *>  may correspond to a matrix modified by perturbations of size near
 135 *>  the underflow threshold.
 136 *> \endverbatim
 137 *>
 138 *  =====================================================================
 139       SUBROUTINE DLASV2( F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL )
 140 *
 141 *  -- LAPACK auxiliary routine (version 3.4.2) --
 142 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 143 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 144 *     September 2012
 145 *
 146 *     .. Scalar Arguments ..
 147       DOUBLE PRECISION   CSL, CSR, F, G, H, SNL, SNR, SSMAX, SSMIN
 148 *     ..
 149 *
 150 * =====================================================================
 151 *
 152 *     .. Parameters ..
 153       DOUBLE PRECISION   ZERO
 154       PARAMETER          ( ZERO = 0.0D0 )
 155       DOUBLE PRECISION   HALF
 156       PARAMETER          ( HALF = 0.5D0 )
 157       DOUBLE PRECISION   ONE
 158       PARAMETER          ( ONE = 1.0D0 )
 159       DOUBLE PRECISION   TWO
 160       PARAMETER          ( TWO = 2.0D0 )
 161       DOUBLE PRECISION   FOUR
 162       PARAMETER          ( FOUR = 4.0D0 )
 163 *     ..
 164 *     .. Local Scalars ..
 165       LOGICAL            GASMAL, SWAP
 166       INTEGER            PMAX
 167       DOUBLE PRECISION   A, CLT, CRT, D, FA, FT, GA, GT, HA, HT, L, M,
 168      $                   MM, R, S, SLT, SRT, T, TEMP, TSIGN, TT
 169 *     ..
 170 *     .. Intrinsic Functions ..
 171       INTRINSIC          ABS, SIGN, SQRT
 172 *     ..
 173 *     .. External Functions ..
 174       DOUBLE PRECISION   DLAMCH
 175       EXTERNAL           DLAMCH
 176 *     ..
 177 *     .. Executable Statements ..
 178 *
 179       FT = F
 180       FA = ABS( FT )
 181       HT = H
 182       HA = ABS( H )
 183 *
 184 *     PMAX points to the maximum absolute element of matrix
 185 *       PMAX = 1 if F largest in absolute values
 186 *       PMAX = 2 if G largest in absolute values
 187 *       PMAX = 3 if H largest in absolute values
 188 *
 189       PMAX = 1
 190       SWAP = ( HA.GT.FA )
 191       IF( SWAP ) THEN
 192          PMAX = 3
 193          TEMP = FT
 194          FT = HT
 195          HT = TEMP
 196          TEMP = FA
 197          FA = HA
 198          HA = TEMP
 199 *
 200 *        Now FA .ge. HA
 201 *
 202       END IF
 203       GT = G
 204       GA = ABS( GT )
 205       IF( GA.EQ.ZERO ) THEN
 206 *
 207 *        Diagonal matrix
 208 *
 209          SSMIN = HA
 210          SSMAX = FA
 211          CLT = ONE
 212          CRT = ONE
 213          SLT = ZERO
 214          SRT = ZERO
 215       ELSE
 216          GASMAL = .TRUE.
 217          IF( GA.GT.FA ) THEN
 218             PMAX = 2
 219             IF( ( FA / GA ).LT.DLAMCH( 'EPS' ) ) THEN
 220 *
 221 *              Case of very large GA
 222 *
 223                GASMAL = .FALSE.
 224                SSMAX = GA
 225                IF( HA.GT.ONE ) THEN
 226                   SSMIN = FA / ( GA / HA )
 227                ELSE
 228                   SSMIN = ( FA / GA )*HA
 229                END IF
 230                CLT = ONE
 231                SLT = HT / GT
 232                SRT = ONE
 233                CRT = FT / GT
 234             END IF
 235          END IF
 236          IF( GASMAL ) THEN
 237 *
 238 *           Normal case
 239 *
 240             D = FA - HA
 241             IF( D.EQ.FA ) THEN
 242 *
 243 *              Copes with infinite F or H
 244 *
 245                L = ONE
 246             ELSE
 247                L = D / FA
 248             END IF
 249 *
 250 *           Note that 0 .le. L .le. 1
 251 *
 252             M = GT / FT
 253 *
 254 *           Note that abs(M) .le. 1/macheps
 255 *
 256             T = TWO - L
 257 *
 258 *           Note that T .ge. 1
 259 *
 260             MM = M*M
 261             TT = T*T
 262             S = SQRT( TT+MM )
 263 *
 264 *           Note that 1 .le. S .le. 1 + 1/macheps
 265 *
 266             IF( L.EQ.ZERO ) THEN
 267                R = ABS( M )
 268             ELSE
 269                R = SQRT( L*L+MM )
 270             END IF
 271 *
 272 *           Note that 0 .le. R .le. 1 + 1/macheps
 273 *
 274             A = HALF*( S+R )
 275 *
 276 *           Note that 1 .le. A .le. 1 + abs(M)
 277 *
 278             SSMIN = HA / A
 279             SSMAX = FA*A
 280             IF( MM.EQ.ZERO ) THEN
 281 *
 282 *              Note that M is very tiny
 283 *
 284                IF( L.EQ.ZERO ) THEN
 285                   T = SIGN( TWO, FT )*SIGN( ONE, GT )
 286                ELSE
 287                   T = GT / SIGN( D, FT ) + M / T
 288                END IF
 289             ELSE
 290                T = ( M / ( S+T )+M / ( R+L ) )*( ONE+A )
 291             END IF
 292             L = SQRT( T*T+FOUR )
 293             CRT = TWO / L
 294             SRT = T / L
 295             CLT = ( CRT+SRT*M ) / A
 296             SLT = ( HT / FT )*SRT / A
 297          END IF
 298       END IF
 299       IF( SWAP ) THEN
 300          CSL = SRT
 301          SNL = CRT
 302          CSR = SLT
 303          SNR = CLT
 304       ELSE
 305          CSL = CLT
 306          SNL = SLT
 307          CSR = CRT
 308          SNR = SRT
 309       END IF
 310 *
 311 *     Correct signs of SSMAX and SSMIN
 312 *
 313       IF( PMAX.EQ.1 )
 314      $   TSIGN = SIGN( ONE, CSR )*SIGN( ONE, CSL )*SIGN( ONE, F )
 315       IF( PMAX.EQ.2 )
 316      $   TSIGN = SIGN( ONE, SNR )*SIGN( ONE, CSL )*SIGN( ONE, G )
 317       IF( PMAX.EQ.3 )
 318      $   TSIGN = SIGN( ONE, SNR )*SIGN( ONE, SNL )*SIGN( ONE, H )
 319       SSMAX = SIGN( SSMAX, TSIGN )
 320       SSMIN = SIGN( SSMIN, TSIGN*SIGN( ONE, F )*SIGN( ONE, H ) )
 321       RETURN
 322 *
 323 *     End of DLASV2
 324 *
 325       END