LAPACK: SRC/slasv2.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SLASV2( F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL )
00002 *
00003 *  -- LAPACK auxiliary routine (version 3.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       REAL               CSL, CSR, F, G, H, SNL, SNR, SSMAX, SSMIN
00010 *     ..
00011 *
00012 *  Purpose
00013 *  =======
00014 *
00015 *  SLASV2 computes the singular value decomposition of a 2-by-2
00016 *  triangular matrix
00017 *     [  F   G  ]
00018 *     [  0   H  ].
00019 *  On return, abs(SSMAX) is the larger singular value, abs(SSMIN) is the
00020 *  smaller singular value, and (CSL,SNL) and (CSR,SNR) are the left and
00021 *  right singular vectors for abs(SSMAX), giving the decomposition
00022 *
00023 *     [ CSL  SNL ] [  F   G  ] [ CSR -SNR ]  =  [ SSMAX   0   ]
00024 *     [-SNL  CSL ] [  0   H  ] [ SNR  CSR ]     [  0    SSMIN ].
00025 *
00026 *  Arguments
00027 *  =========
00028 *
00029 *  F       (input) REAL
00030 *          The (1,1) element of the 2-by-2 matrix.
00031 *
00032 *  G       (input) REAL
00033 *          The (1,2) element of the 2-by-2 matrix.
00034 *
00035 *  H       (input) REAL
00036 *          The (2,2) element of the 2-by-2 matrix.
00037 *
00038 *  SSMIN   (output) REAL
00039 *          abs(SSMIN) is the smaller singular value.
00040 *
00041 *  SSMAX   (output) REAL
00042 *          abs(SSMAX) is the larger singular value.
00043 *
00044 *  SNL     (output) REAL
00045 *  CSL     (output) REAL
00046 *          The vector (CSL, SNL) is a unit left singular vector for the
00047 *          singular value abs(SSMAX).
00048 *
00049 *  SNR     (output) REAL
00050 *  CSR     (output) REAL
00051 *          The vector (CSR, SNR) is a unit right singular vector for the
00052 *          singular value abs(SSMAX).
00053 *
00054 *  Further Details
00055 *  ===============
00056 *
00057 *  Any input parameter may be aliased with any output parameter.
00058 *
00059 *  Barring over/underflow and assuming a guard digit in subtraction, all
00060 *  output quantities are correct to within a few units in the last
00061 *  place (ulps).
00062 *
00063 *  In IEEE arithmetic, the code works correctly if one matrix element is
00064 *  infinite.
00065 *
00066 *  Overflow will not occur unless the largest singular value itself
00067 *  overflows or is within a few ulps of overflow. (On machines with
00068 *  partial overflow, like the Cray, overflow may occur if the largest
00069 *  singular value is within a factor of 2 of overflow.)
00070 *
00071 *  Underflow is harmless if underflow is gradual. Otherwise, results
00072 *  may correspond to a matrix modified by perturbations of size near
00073 *  the underflow threshold.
00074 *
00075 * =====================================================================
00076 *
00077 *     .. Parameters ..
00078       REAL               ZERO
00079       PARAMETER          ( ZERO = 0.0E0 )
00080       REAL               HALF
00081       PARAMETER          ( HALF = 0.5E0 )
00082       REAL               ONE
00083       PARAMETER          ( ONE = 1.0E0 )
00084       REAL               TWO
00085       PARAMETER          ( TWO = 2.0E0 )
00086       REAL               FOUR
00087       PARAMETER          ( FOUR = 4.0E0 )
00088 *     ..
00089 *     .. Local Scalars ..
00090       LOGICAL            GASMAL, SWAP
00091       INTEGER            PMAX
00092       REAL               A, CLT, CRT, D, FA, FT, GA, GT, HA, HT, L, M,
00093      $                   MM, R, S, SLT, SRT, T, TEMP, TSIGN, TT
00094 *     ..
00095 *     .. Intrinsic Functions ..
00096       INTRINSIC          ABS, SIGN, SQRT
00097 *     ..
00098 *     .. External Functions ..
00099       REAL               SLAMCH
00100       EXTERNAL           SLAMCH
00101 *     ..
00102 *     .. Executable Statements ..
00103 *
00104       FT = F
00105       FA = ABS( FT )
00106       HT = H
00107       HA = ABS( H )
00108 *
00109 *     PMAX points to the maximum absolute element of matrix
00110 *       PMAX = 1 if F largest in absolute values
00111 *       PMAX = 2 if G largest in absolute values
00112 *       PMAX = 3 if H largest in absolute values
00113 *
00114       PMAX = 1
00115       SWAP = ( HA.GT.FA )
00116       IF( SWAP ) THEN
00117          PMAX = 3
00118          TEMP = FT
00119          FT = HT
00120          HT = TEMP
00121          TEMP = FA
00122          FA = HA
00123          HA = TEMP
00124 *
00125 *        Now FA .ge. HA
00126 *
00127       END IF
00128       GT = G
00129       GA = ABS( GT )
00130       IF( GA.EQ.ZERO ) THEN
00131 *
00132 *        Diagonal matrix
00133 *
00134          SSMIN = HA
00135          SSMAX = FA
00136          CLT = ONE
00137          CRT = ONE
00138          SLT = ZERO
00139          SRT = ZERO
00140       ELSE
00141          GASMAL = .TRUE.
00142          IF( GA.GT.FA ) THEN
00143             PMAX = 2
00144             IF( ( FA / GA ).LT.SLAMCH( 'EPS' ) ) THEN
00145 *
00146 *              Case of very large GA
00147 *
00148                GASMAL = .FALSE.
00149                SSMAX = GA
00150                IF( HA.GT.ONE ) THEN
00151                   SSMIN = FA / ( GA / HA )
00152                ELSE
00153                   SSMIN = ( FA / GA )*HA
00154                END IF
00155                CLT = ONE
00156                SLT = HT / GT
00157                SRT = ONE
00158                CRT = FT / GT
00159             END IF
00160          END IF
00161          IF( GASMAL ) THEN
00162 *
00163 *           Normal case
00164 *
00165             D = FA - HA
00166             IF( D.EQ.FA ) THEN
00167 *
00168 *              Copes with infinite F or H
00169 *
00170                L = ONE
00171             ELSE
00172                L = D / FA
00173             END IF
00174 *
00175 *           Note that 0 .le. L .le. 1
00176 *
00177             M = GT / FT
00178 *
00179 *           Note that abs(M) .le. 1/macheps
00180 *
00181             T = TWO - L
00182 *
00183 *           Note that T .ge. 1
00184 *
00185             MM = M*M
00186             TT = T*T
00187             S = SQRT( TT+MM )
00188 *
00189 *           Note that 1 .le. S .le. 1 + 1/macheps
00190 *
00191             IF( L.EQ.ZERO ) THEN
00192                R = ABS( M )
00193             ELSE
00194                R = SQRT( L*L+MM )
00195             END IF
00196 *
00197 *           Note that 0 .le. R .le. 1 + 1/macheps
00198 *
00199             A = HALF*( S+R )
00200 *
00201 *           Note that 1 .le. A .le. 1 + abs(M)
00202 *
00203             SSMIN = HA / A
00204             SSMAX = FA*A
00205             IF( MM.EQ.ZERO ) THEN
00206 *
00207 *              Note that M is very tiny
00208 *
00209                IF( L.EQ.ZERO ) THEN
00210                   T = SIGN( TWO, FT )*SIGN( ONE, GT )
00211                ELSE
00212                   T = GT / SIGN( D, FT ) + M / T
00213                END IF
00214             ELSE
00215                T = ( M / ( S+T )+M / ( R+L ) )*( ONE+A )
00216             END IF
00217             L = SQRT( T*T+FOUR )
00218             CRT = TWO / L
00219             SRT = T / L
00220             CLT = ( CRT+SRT*M ) / A
00221             SLT = ( HT / FT )*SRT / A
00222          END IF
00223       END IF
00224       IF( SWAP ) THEN
00225          CSL = SRT
00226          SNL = CRT
00227          CSR = SLT
00228          SNR = CLT
00229       ELSE
00230          CSL = CLT
00231          SNL = SLT
00232          CSR = CRT
00233          SNR = SRT
00234       END IF
00235 *
00236 *     Correct signs of SSMAX and SSMIN
00237 *
00238       IF( PMAX.EQ.1 )
00239      $   TSIGN = SIGN( ONE, CSR )*SIGN( ONE, CSL )*SIGN( ONE, F )
00240       IF( PMAX.EQ.2 )
00241      $   TSIGN = SIGN( ONE, SNR )*SIGN( ONE, CSL )*SIGN( ONE, G )
00242       IF( PMAX.EQ.3 )
00243      $   TSIGN = SIGN( ONE, SNR )*SIGN( ONE, SNL )*SIGN( ONE, H )
00244       SSMAX = SIGN( SSMAX, TSIGN )
00245       SSMIN = SIGN( SSMIN, TSIGN*SIGN( ONE, F )*SIGN( ONE, H ) )
00246       RETURN
00247 *
00248 *     End of SLASV2
00249 *
00250       END