![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
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