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