LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dlasv2.f
Go to the documentation of this file.
00001 *> \brief \b DLASV2
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DLASV2 + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasv2.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasv2.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasv2.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DLASV2( F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       DOUBLE PRECISION   CSL, CSR, F, G, H, SNL, SNR, SSMAX, SSMIN
00025 *       ..
00026 *  
00027 *
00028 *> \par Purpose:
00029 *  =============
00030 *>
00031 *> \verbatim
00032 *>
00033 *> DLASV2 computes the singular value decomposition of a 2-by-2
00034 *> triangular matrix
00035 *>    [  F   G  ]
00036 *>    [  0   H  ].
00037 *> On return, abs(SSMAX) is the larger singular value, abs(SSMIN) is the
00038 *> smaller singular value, and (CSL,SNL) and (CSR,SNR) are the left and
00039 *> right singular vectors for abs(SSMAX), giving the decomposition
00040 *>
00041 *>    [ CSL  SNL ] [  F   G  ] [ CSR -SNR ]  =  [ SSMAX   0   ]
00042 *>    [-SNL  CSL ] [  0   H  ] [ SNR  CSR ]     [  0    SSMIN ].
00043 *> \endverbatim
00044 *
00045 *  Arguments:
00046 *  ==========
00047 *
00048 *> \param[in] F
00049 *> \verbatim
00050 *>          F is DOUBLE PRECISION
00051 *>          The (1,1) element of the 2-by-2 matrix.
00052 *> \endverbatim
00053 *>
00054 *> \param[in] G
00055 *> \verbatim
00056 *>          G is DOUBLE PRECISION
00057 *>          The (1,2) element of the 2-by-2 matrix.
00058 *> \endverbatim
00059 *>
00060 *> \param[in] H
00061 *> \verbatim
00062 *>          H is DOUBLE PRECISION
00063 *>          The (2,2) element of the 2-by-2 matrix.
00064 *> \endverbatim
00065 *>
00066 *> \param[out] SSMIN
00067 *> \verbatim
00068 *>          SSMIN is DOUBLE PRECISION
00069 *>          abs(SSMIN) is the smaller singular value.
00070 *> \endverbatim
00071 *>
00072 *> \param[out] SSMAX
00073 *> \verbatim
00074 *>          SSMAX is DOUBLE PRECISION
00075 *>          abs(SSMAX) is the larger singular value.
00076 *> \endverbatim
00077 *>
00078 *> \param[out] SNL
00079 *> \verbatim
00080 *>          SNL is DOUBLE PRECISION
00081 *> \endverbatim
00082 *>
00083 *> \param[out] CSL
00084 *> \verbatim
00085 *>          CSL is DOUBLE PRECISION
00086 *>          The vector (CSL, SNL) is a unit left singular vector for the
00087 *>          singular value abs(SSMAX).
00088 *> \endverbatim
00089 *>
00090 *> \param[out] SNR
00091 *> \verbatim
00092 *>          SNR is DOUBLE PRECISION
00093 *> \endverbatim
00094 *>
00095 *> \param[out] CSR
00096 *> \verbatim
00097 *>          CSR is DOUBLE PRECISION
00098 *>          The vector (CSR, SNR) is a unit right singular vector for the
00099 *>          singular value abs(SSMAX).
00100 *> \endverbatim
00101 *
00102 *  Authors:
00103 *  ========
00104 *
00105 *> \author Univ. of Tennessee 
00106 *> \author Univ. of California Berkeley 
00107 *> \author Univ. of Colorado Denver 
00108 *> \author NAG Ltd. 
00109 *
00110 *> \date November 2011
00111 *
00112 *> \ingroup auxOTHERauxiliary
00113 *
00114 *> \par Further Details:
00115 *  =====================
00116 *>
00117 *> \verbatim
00118 *>
00119 *>  Any input parameter may be aliased with any output parameter.
00120 *>
00121 *>  Barring over/underflow and assuming a guard digit in subtraction, all
00122 *>  output quantities are correct to within a few units in the last
00123 *>  place (ulps).
00124 *>
00125 *>  In IEEE arithmetic, the code works correctly if one matrix element is
00126 *>  infinite.
00127 *>
00128 *>  Overflow will not occur unless the largest singular value itself
00129 *>  overflows or is within a few ulps of overflow. (On machines with
00130 *>  partial overflow, like the Cray, overflow may occur if the largest
00131 *>  singular value is within a factor of 2 of overflow.)
00132 *>
00133 *>  Underflow is harmless if underflow is gradual. Otherwise, results
00134 *>  may correspond to a matrix modified by perturbations of size near
00135 *>  the underflow threshold.
00136 *> \endverbatim
00137 *>
00138 *  =====================================================================
00139       SUBROUTINE DLASV2( F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL )
00140 *
00141 *  -- LAPACK auxiliary routine (version 3.4.0) --
00142 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00143 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00144 *     November 2011
00145 *
00146 *     .. Scalar Arguments ..
00147       DOUBLE PRECISION   CSL, CSR, F, G, H, SNL, SNR, SSMAX, SSMIN
00148 *     ..
00149 *
00150 * =====================================================================
00151 *
00152 *     .. Parameters ..
00153       DOUBLE PRECISION   ZERO
00154       PARAMETER          ( ZERO = 0.0D0 )
00155       DOUBLE PRECISION   HALF
00156       PARAMETER          ( HALF = 0.5D0 )
00157       DOUBLE PRECISION   ONE
00158       PARAMETER          ( ONE = 1.0D0 )
00159       DOUBLE PRECISION   TWO
00160       PARAMETER          ( TWO = 2.0D0 )
00161       DOUBLE PRECISION   FOUR
00162       PARAMETER          ( FOUR = 4.0D0 )
00163 *     ..
00164 *     .. Local Scalars ..
00165       LOGICAL            GASMAL, SWAP
00166       INTEGER            PMAX
00167       DOUBLE PRECISION   A, CLT, CRT, D, FA, FT, GA, GT, HA, HT, L, M,
00168      $                   MM, R, S, SLT, SRT, T, TEMP, TSIGN, TT
00169 *     ..
00170 *     .. Intrinsic Functions ..
00171       INTRINSIC          ABS, SIGN, SQRT
00172 *     ..
00173 *     .. External Functions ..
00174       DOUBLE PRECISION   DLAMCH
00175       EXTERNAL           DLAMCH
00176 *     ..
00177 *     .. Executable Statements ..
00178 *
00179       FT = F
00180       FA = ABS( FT )
00181       HT = H
00182       HA = ABS( H )
00183 *
00184 *     PMAX points to the maximum absolute element of matrix
00185 *       PMAX = 1 if F largest in absolute values
00186 *       PMAX = 2 if G largest in absolute values
00187 *       PMAX = 3 if H largest in absolute values
00188 *
00189       PMAX = 1
00190       SWAP = ( HA.GT.FA )
00191       IF( SWAP ) THEN
00192          PMAX = 3
00193          TEMP = FT
00194          FT = HT
00195          HT = TEMP
00196          TEMP = FA
00197          FA = HA
00198          HA = TEMP
00199 *
00200 *        Now FA .ge. HA
00201 *
00202       END IF
00203       GT = G
00204       GA = ABS( GT )
00205       IF( GA.EQ.ZERO ) THEN
00206 *
00207 *        Diagonal matrix
00208 *
00209          SSMIN = HA
00210          SSMAX = FA
00211          CLT = ONE
00212          CRT = ONE
00213          SLT = ZERO
00214          SRT = ZERO
00215       ELSE
00216          GASMAL = .TRUE.
00217          IF( GA.GT.FA ) THEN
00218             PMAX = 2
00219             IF( ( FA / GA ).LT.DLAMCH( 'EPS' ) ) THEN
00220 *
00221 *              Case of very large GA
00222 *
00223                GASMAL = .FALSE.
00224                SSMAX = GA
00225                IF( HA.GT.ONE ) THEN
00226                   SSMIN = FA / ( GA / HA )
00227                ELSE
00228                   SSMIN = ( FA / GA )*HA
00229                END IF
00230                CLT = ONE
00231                SLT = HT / GT
00232                SRT = ONE
00233                CRT = FT / GT
00234             END IF
00235          END IF
00236          IF( GASMAL ) THEN
00237 *
00238 *           Normal case
00239 *
00240             D = FA - HA
00241             IF( D.EQ.FA ) THEN
00242 *
00243 *              Copes with infinite F or H
00244 *
00245                L = ONE
00246             ELSE
00247                L = D / FA
00248             END IF
00249 *
00250 *           Note that 0 .le. L .le. 1
00251 *
00252             M = GT / FT
00253 *
00254 *           Note that abs(M) .le. 1/macheps
00255 *
00256             T = TWO - L
00257 *
00258 *           Note that T .ge. 1
00259 *
00260             MM = M*M
00261             TT = T*T
00262             S = SQRT( TT+MM )
00263 *
00264 *           Note that 1 .le. S .le. 1 + 1/macheps
00265 *
00266             IF( L.EQ.ZERO ) THEN
00267                R = ABS( M )
00268             ELSE
00269                R = SQRT( L*L+MM )
00270             END IF
00271 *
00272 *           Note that 0 .le. R .le. 1 + 1/macheps
00273 *
00274             A = HALF*( S+R )
00275 *
00276 *           Note that 1 .le. A .le. 1 + abs(M)
00277 *
00278             SSMIN = HA / A
00279             SSMAX = FA*A
00280             IF( MM.EQ.ZERO ) THEN
00281 *
00282 *              Note that M is very tiny
00283 *
00284                IF( L.EQ.ZERO ) THEN
00285                   T = SIGN( TWO, FT )*SIGN( ONE, GT )
00286                ELSE
00287                   T = GT / SIGN( D, FT ) + M / T
00288                END IF
00289             ELSE
00290                T = ( M / ( S+T )+M / ( R+L ) )*( ONE+A )
00291             END IF
00292             L = SQRT( T*T+FOUR )
00293             CRT = TWO / L
00294             SRT = T / L
00295             CLT = ( CRT+SRT*M ) / A
00296             SLT = ( HT / FT )*SRT / A
00297          END IF
00298       END IF
00299       IF( SWAP ) THEN
00300          CSL = SRT
00301          SNL = CRT
00302          CSR = SLT
00303          SNR = CLT
00304       ELSE
00305          CSL = CLT
00306          SNL = SLT
00307          CSR = CRT
00308          SNR = SRT
00309       END IF
00310 *
00311 *     Correct signs of SSMAX and SSMIN
00312 *
00313       IF( PMAX.EQ.1 )
00314      $   TSIGN = SIGN( ONE, CSR )*SIGN( ONE, CSL )*SIGN( ONE, F )
00315       IF( PMAX.EQ.2 )
00316      $   TSIGN = SIGN( ONE, SNR )*SIGN( ONE, CSL )*SIGN( ONE, G )
00317       IF( PMAX.EQ.3 )
00318      $   TSIGN = SIGN( ONE, SNR )*SIGN( ONE, SNL )*SIGN( ONE, H )
00319       SSMAX = SIGN( SSMAX, TSIGN )
00320       SSMIN = SIGN( SSMIN, TSIGN*SIGN( ONE, F )*SIGN( ONE, H ) )
00321       RETURN
00322 *
00323 *     End of DLASV2
00324 *
00325       END
 All Files Functions