158
159
160
161
162
163
164 LOGICAL UPPER
165 REAL A1, A3, B1, B3, CSQ, CSU, CSV
166 COMPLEX A2, B2, SNQ, SNU, SNV
167
168
169
170
171
172 REAL ZERO, ONE
173 parameter( zero = 0.0e+0, one = 1.0e+0 )
174
175
176 REAL A, AUA11, AUA12, AUA21, AUA22, AVB11, AVB12,
177 $ AVB21, AVB22, CSL, CSR, D, FB, FC, S1, S2, SNL,
178 $ SNR, UA11R, UA22R, VB11R, VB22R
179 COMPLEX B, C, D1, R, T, UA11, UA12, UA21, UA22, VB11,
180 $ VB12, VB21, VB22
181
182
184
185
186 INTRINSIC abs, aimag, cmplx, conjg, real
187
188
189 REAL ABS1
190
191
192 abs1( t ) = abs( real( t ) ) + abs( aimag( t ) )
193
194
195
196 IF( upper ) THEN
197
198
199
200
201
202
203 a = a1*b3
204 d = a3*b1
205 b = a2*b1 - a1*b2
206 fb = abs( b )
207
208
209
210
211 d1 = one
212 IF( fb.NE.zero )
213 $ d1 = b / fb
214
215
216
217
218
219
220 CALL slasv2( a, fb, d, s1, s2, snr, csr, snl, csl )
221
222 IF( abs( csl ).GE.abs( snl ) .OR. abs( csr ).GE.abs( snr ) )
223 $ THEN
224
225
226
227
228 ua11r = csl*a1
229 ua12 = csl*a2 + d1*snl*a3
230
231 vb11r = csr*b1
232 vb12 = csr*b2 + d1*snr*b3
233
234 aua12 = abs( csl )*abs1( a2 ) + abs( snl )*abs( a3 )
235 avb12 = abs( csr )*abs1( b2 ) + abs( snr )*abs( b3 )
236
237
238
239 IF( ( abs( ua11r )+abs1( ua12 ) ).EQ.zero ) THEN
240 CALL clartg( -cmplx( vb11r ), conjg( vb12 ), csq, snq,
241 $ r )
242 ELSE IF( ( abs( vb11r )+abs1( vb12 ) ).EQ.zero ) THEN
243 CALL clartg( -cmplx( ua11r ), conjg( ua12 ), csq, snq,
244 $ r )
245 ELSE IF( aua12 / ( abs( ua11r )+abs1( ua12 ) ).LE.avb12 /
246 $ ( abs( vb11r )+abs1( vb12 ) ) ) THEN
247 CALL clartg( -cmplx( ua11r ), conjg( ua12 ), csq, snq,
248 $ r )
249 ELSE
250 CALL clartg( -cmplx( vb11r ), conjg( vb12 ), csq, snq,
251 $ r )
252 END IF
253
254 csu = csl
255 snu = -d1*snl
256 csv = csr
257 snv = -d1*snr
258
259 ELSE
260
261
262
263
264 ua21 = -conjg( d1 )*snl*a1
265 ua22 = -conjg( d1 )*snl*a2 + csl*a3
266
267 vb21 = -conjg( d1 )*snr*b1
268 vb22 = -conjg( d1 )*snr*b2 + csr*b3
269
270 aua22 = abs( snl )*abs1( a2 ) + abs( csl )*abs( a3 )
271 avb22 = abs( snr )*abs1( b2 ) + abs( csr )*abs( b3 )
272
273
274
275 IF( ( abs1( ua21 )+abs1( ua22 ) ).EQ.zero ) THEN
276 CALL clartg( -conjg( vb21 ), conjg( vb22 ), csq, snq, r )
277 ELSE IF( ( abs1( vb21 )+abs( vb22 ) ).EQ.zero ) THEN
278 CALL clartg( -conjg( ua21 ), conjg( ua22 ), csq, snq, r )
279 ELSE IF( aua22 / ( abs1( ua21 )+abs1( ua22 ) ).LE.avb22 /
280 $ ( abs1( vb21 )+abs1( vb22 ) ) ) THEN
281 CALL clartg( -conjg( ua21 ), conjg( ua22 ), csq, snq, r )
282 ELSE
283 CALL clartg( -conjg( vb21 ), conjg( vb22 ), csq, snq, r )
284 END IF
285
286 csu = snl
287 snu = d1*csl
288 csv = snr
289 snv = d1*csr
290
291 END IF
292
293 ELSE
294
295
296
297
298
299
300 a = a1*b3
301 d = a3*b1
302 c = a2*b3 - a3*b2
303 fc = abs( c )
304
305
306
307
308 d1 = one
309 IF( fc.NE.zero )
310 $ d1 = c / fc
311
312
313
314
315
316
317 CALL slasv2( a, fc, d, s1, s2, snr, csr, snl, csl )
318
319 IF( abs( csr ).GE.abs( snr ) .OR. abs( csl ).GE.abs( snl ) )
320 $ THEN
321
322
323
324
325 ua21 = -d1*snr*a1 + csr*a2
326 ua22r = csr*a3
327
328 vb21 = -d1*snl*b1 + csl*b2
329 vb22r = csl*b3
330
331 aua21 = abs( snr )*abs( a1 ) + abs( csr )*abs1( a2 )
332 avb21 = abs( snl )*abs( b1 ) + abs( csl )*abs1( b2 )
333
334
335
336 IF( ( abs1( ua21 )+abs( ua22r ) ).EQ.zero ) THEN
337 CALL clartg( cmplx( vb22r ), vb21, csq, snq, r )
338 ELSE IF( ( abs1( vb21 )+abs( vb22r ) ).EQ.zero ) THEN
339 CALL clartg( cmplx( ua22r ), ua21, csq, snq, r )
340 ELSE IF( aua21 / ( abs1( ua21 )+abs( ua22r ) ).LE.avb21 /
341 $ ( abs1( vb21 )+abs( vb22r ) ) ) THEN
342 CALL clartg( cmplx( ua22r ), ua21, csq, snq, r )
343 ELSE
344 CALL clartg( cmplx( vb22r ), vb21, csq, snq, r )
345 END IF
346
347 csu = csr
348 snu = -conjg( d1 )*snr
349 csv = csl
350 snv = -conjg( d1 )*snl
351
352 ELSE
353
354
355
356
357 ua11 = csr*a1 + conjg( d1 )*snr*a2
358 ua12 = conjg( d1 )*snr*a3
359
360 vb11 = csl*b1 + conjg( d1 )*snl*b2
361 vb12 = conjg( d1 )*snl*b3
362
363 aua11 = abs( csr )*abs( a1 ) + abs( snr )*abs1( a2 )
364 avb11 = abs( csl )*abs( b1 ) + abs( snl )*abs1( b2 )
365
366
367
368 IF( ( abs1( ua11 )+abs1( ua12 ) ).EQ.zero ) THEN
369 CALL clartg( vb12, vb11, csq, snq, r )
370 ELSE IF( ( abs1( vb11 )+abs1( vb12 ) ).EQ.zero ) THEN
371 CALL clartg( ua12, ua11, csq, snq, r )
372 ELSE IF( aua11 / ( abs1( ua11 )+abs1( ua12 ) ).LE.avb11 /
373 $ ( abs1( vb11 )+abs1( vb12 ) ) ) THEN
374 CALL clartg( ua12, ua11, csq, snq, r )
375 ELSE
376 CALL clartg( vb12, vb11, csq, snq, r )
377 END IF
378
379 csu = snr
380 snu = conjg( d1 )*csr
381 csv = snl
382 snv = conjg( d1 )*csl
383
384 END IF
385
386 END IF
387
388 RETURN
389
390
391
subroutine clartg(f, g, c, s, r)
CLARTG generates a plane rotation with real cosine and complex sine.
subroutine slasv2(f, g, h, ssmin, ssmax, snr, csr, snl, csl)
SLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix.