001:       SUBROUTINE SLASV2( F, G, H, SSMIN, SSMAX, SNR, CSR, SNL, CSL )
002: *
003: *  -- LAPACK auxiliary routine (version 3.2) --
004: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
005: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       REAL               CSL, CSR, F, G, H, SNL, SNR, SSMAX, SSMIN
010: *     ..
011: *
012: *  Purpose
013: *  =======
014: *
015: *  SLASV2 computes the singular value decomposition of a 2-by-2
016: *  triangular matrix
017: *     [  F   G  ]
018: *     [  0   H  ].
019: *  On return, abs(SSMAX) is the larger singular value, abs(SSMIN) is the
020: *  smaller singular value, and (CSL,SNL) and (CSR,SNR) are the left and
021: *  right singular vectors for abs(SSMAX), giving the decomposition
022: *
023: *     [ CSL  SNL ] [  F   G  ] [ CSR -SNR ]  =  [ SSMAX   0   ]
024: *     [-SNL  CSL ] [  0   H  ] [ SNR  CSR ]     [  0    SSMIN ].
025: *
026: *  Arguments
027: *  =========
028: *
029: *  F       (input) REAL
030: *          The (1,1) element of the 2-by-2 matrix.
031: *
032: *  G       (input) REAL
033: *          The (1,2) element of the 2-by-2 matrix.
034: *
035: *  H       (input) REAL
036: *          The (2,2) element of the 2-by-2 matrix.
037: *
038: *  SSMIN   (output) REAL
039: *          abs(SSMIN) is the smaller singular value.
040: *
041: *  SSMAX   (output) REAL
042: *          abs(SSMAX) is the larger singular value.
043: *
044: *  SNL     (output) REAL
045: *  CSL     (output) REAL
046: *          The vector (CSL, SNL) is a unit left singular vector for the
047: *          singular value abs(SSMAX).
048: *
049: *  SNR     (output) REAL
050: *  CSR     (output) REAL
051: *          The vector (CSR, SNR) is a unit right singular vector for the
052: *          singular value abs(SSMAX).
053: *
054: *  Further Details
055: *  ===============
056: *
057: *  Any input parameter may be aliased with any output parameter.
058: *
059: *  Barring over/underflow and assuming a guard digit in subtraction, all
060: *  output quantities are correct to within a few units in the last
061: *  place (ulps).
062: *
063: *  In IEEE arithmetic, the code works correctly if one matrix element is
064: *  infinite.
065: *
066: *  Overflow will not occur unless the largest singular value itself
067: *  overflows or is within a few ulps of overflow. (On machines with
068: *  partial overflow, like the Cray, overflow may occur if the largest
069: *  singular value is within a factor of 2 of overflow.)
070: *
071: *  Underflow is harmless if underflow is gradual. Otherwise, results
072: *  may correspond to a matrix modified by perturbations of size near
073: *  the underflow threshold.
074: *
075: * =====================================================================
076: *
077: *     .. Parameters ..
078:       REAL               ZERO
079:       PARAMETER          ( ZERO = 0.0E0 )
080:       REAL               HALF
081:       PARAMETER          ( HALF = 0.5E0 )
082:       REAL               ONE
083:       PARAMETER          ( ONE = 1.0E0 )
084:       REAL               TWO
085:       PARAMETER          ( TWO = 2.0E0 )
086:       REAL               FOUR
087:       PARAMETER          ( FOUR = 4.0E0 )
088: *     ..
089: *     .. Local Scalars ..
090:       LOGICAL            GASMAL, SWAP
091:       INTEGER            PMAX
092:       REAL               A, CLT, CRT, D, FA, FT, GA, GT, HA, HT, L, M,
093:      $                   MM, R, S, SLT, SRT, T, TEMP, TSIGN, TT
094: *     ..
095: *     .. Intrinsic Functions ..
096:       INTRINSIC          ABS, SIGN, SQRT
097: *     ..
098: *     .. External Functions ..
099:       REAL               SLAMCH
100:       EXTERNAL           SLAMCH
101: *     ..
102: *     .. Executable Statements ..
103: *
104:       FT = F
105:       FA = ABS( FT )
106:       HT = H
107:       HA = ABS( H )
108: *
109: *     PMAX points to the maximum absolute element of matrix
110: *       PMAX = 1 if F largest in absolute values
111: *       PMAX = 2 if G largest in absolute values
112: *       PMAX = 3 if H largest in absolute values
113: *
114:       PMAX = 1
115:       SWAP = ( HA.GT.FA )
116:       IF( SWAP ) THEN
117:          PMAX = 3
118:          TEMP = FT
119:          FT = HT
120:          HT = TEMP
121:          TEMP = FA
122:          FA = HA
123:          HA = TEMP
124: *
125: *        Now FA .ge. HA
126: *
127:       END IF
128:       GT = G
129:       GA = ABS( GT )
130:       IF( GA.EQ.ZERO ) THEN
131: *
132: *        Diagonal matrix
133: *
134:          SSMIN = HA
135:          SSMAX = FA
136:          CLT = ONE
137:          CRT = ONE
138:          SLT = ZERO
139:          SRT = ZERO
140:       ELSE
141:          GASMAL = .TRUE.
142:          IF( GA.GT.FA ) THEN
143:             PMAX = 2
144:             IF( ( FA / GA ).LT.SLAMCH( 'EPS' ) ) THEN
145: *
146: *              Case of very large GA
147: *
148:                GASMAL = .FALSE.
149:                SSMAX = GA
150:                IF( HA.GT.ONE ) THEN
151:                   SSMIN = FA / ( GA / HA )
152:                ELSE
153:                   SSMIN = ( FA / GA )*HA
154:                END IF
155:                CLT = ONE
156:                SLT = HT / GT
157:                SRT = ONE
158:                CRT = FT / GT
159:             END IF
160:          END IF
161:          IF( GASMAL ) THEN
162: *
163: *           Normal case
164: *
165:             D = FA - HA
166:             IF( D.EQ.FA ) THEN
167: *
168: *              Copes with infinite F or H
169: *
170:                L = ONE
171:             ELSE
172:                L = D / FA
173:             END IF
174: *
175: *           Note that 0 .le. L .le. 1
176: *
177:             M = GT / FT
178: *
179: *           Note that abs(M) .le. 1/macheps
180: *
181:             T = TWO - L
182: *
183: *           Note that T .ge. 1
184: *
185:             MM = M*M
186:             TT = T*T
187:             S = SQRT( TT+MM )
188: *
189: *           Note that 1 .le. S .le. 1 + 1/macheps
190: *
191:             IF( L.EQ.ZERO ) THEN
192:                R = ABS( M )
193:             ELSE
194:                R = SQRT( L*L+MM )
195:             END IF
196: *
197: *           Note that 0 .le. R .le. 1 + 1/macheps
198: *
199:             A = HALF*( S+R )
200: *
201: *           Note that 1 .le. A .le. 1 + abs(M)
202: *
203:             SSMIN = HA / A
204:             SSMAX = FA*A
205:             IF( MM.EQ.ZERO ) THEN
206: *
207: *              Note that M is very tiny
208: *
209:                IF( L.EQ.ZERO ) THEN
210:                   T = SIGN( TWO, FT )*SIGN( ONE, GT )
211:                ELSE
212:                   T = GT / SIGN( D, FT ) + M / T
213:                END IF
214:             ELSE
215:                T = ( M / ( S+T )+M / ( R+L ) )*( ONE+A )
216:             END IF
217:             L = SQRT( T*T+FOUR )
218:             CRT = TWO / L
219:             SRT = T / L
220:             CLT = ( CRT+SRT*M ) / A
221:             SLT = ( HT / FT )*SRT / A
222:          END IF
223:       END IF
224:       IF( SWAP ) THEN
225:          CSL = SRT
226:          SNL = CRT
227:          CSR = SLT
228:          SNR = CLT
229:       ELSE
230:          CSL = CLT
231:          SNL = SLT
232:          CSR = CRT
233:          SNR = SRT
234:       END IF
235: *
236: *     Correct signs of SSMAX and SSMIN
237: *
238:       IF( PMAX.EQ.1 )
239:      $   TSIGN = SIGN( ONE, CSR )*SIGN( ONE, CSL )*SIGN( ONE, F )
240:       IF( PMAX.EQ.2 )
241:      $   TSIGN = SIGN( ONE, SNR )*SIGN( ONE, CSL )*SIGN( ONE, G )
242:       IF( PMAX.EQ.3 )
243:      $   TSIGN = SIGN( ONE, SNR )*SIGN( ONE, SNL )*SIGN( ONE, H )
244:       SSMAX = SIGN( SSMAX, TSIGN )
245:       SSMIN = SIGN( SSMIN, TSIGN*SIGN( ONE, F )*SIGN( ONE, H ) )
246:       RETURN
247: *
248: *     End of SLASV2
249: *
250:       END
251: