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