157
158
159
160
161
162
163 LOGICAL UPPER
164 DOUBLE PRECISION A1, A3, B1, B3, CSQ, CSU, CSV
165 COMPLEX*16 A2, B2, SNQ, SNU, SNV
166
167
168
169
170
171 DOUBLE PRECISION ZERO, ONE
172 parameter( zero = 0.0d+0, one = 1.0d+0 )
173
174
175 DOUBLE PRECISION A, AUA11, AUA12, AUA21, AUA22, AVB12, AVB11,
176 $ AVB21, AVB22, CSL, CSR, D, FB, FC, S1, S2,
177 $ SNL, SNR, UA11R, UA22R, VB11R, VB22R
178 COMPLEX*16 B, C, D1, R, T, UA11, UA12, UA21, UA22, VB11,
179 $ VB12, VB21, VB22
180
181
183
184
185 INTRINSIC abs, dble, dcmplx, dconjg, dimag
186
187
188 DOUBLE PRECISION ABS1
189
190
191 abs1( t ) = abs( dble( t ) ) + abs( dimag( 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 dlasv2( 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 zlartg( -dcmplx( vb11r ), dconjg( vb12 ), csq,
240 $ snq,
241 $ r )
242 ELSE IF( ( abs( vb11r )+abs1( vb12 ) ).EQ.zero ) THEN
243 CALL zlartg( -dcmplx( ua11r ), dconjg( ua12 ), csq,
244 $ snq,
245 $ r )
246 ELSE IF( aua12 / ( abs( ua11r )+abs1( ua12 ) ).LE.avb12 /
247 $ ( abs( vb11r )+abs1( vb12 ) ) ) THEN
248 CALL zlartg( -dcmplx( ua11r ), dconjg( ua12 ), csq,
249 $ snq,
250 $ r )
251 ELSE
252 CALL zlartg( -dcmplx( vb11r ), dconjg( vb12 ), csq,
253 $ snq,
254 $ r )
255 END IF
256
257 csu = csl
258 snu = -d1*snl
259 csv = csr
260 snv = -d1*snr
261
262 ELSE
263
264
265
266
267 ua21 = -dconjg( d1 )*snl*a1
268 ua22 = -dconjg( d1 )*snl*a2 + csl*a3
269
270 vb21 = -dconjg( d1 )*snr*b1
271 vb22 = -dconjg( d1 )*snr*b2 + csr*b3
272
273 aua22 = abs( snl )*abs1( a2 ) + abs( csl )*abs( a3 )
274 avb22 = abs( snr )*abs1( b2 ) + abs( csr )*abs( b3 )
275
276
277
278 IF( ( abs1( ua21 )+abs1( ua22 ) ).EQ.zero ) THEN
279 CALL zlartg( -dconjg( vb21 ), dconjg( vb22 ), csq,
280 $ snq,
281 $ r )
282 ELSE IF( ( abs1( vb21 )+abs( vb22 ) ).EQ.zero ) THEN
283 CALL zlartg( -dconjg( ua21 ), dconjg( ua22 ), csq,
284 $ snq,
285 $ r )
286 ELSE IF( aua22 / ( abs1( ua21 )+abs1( ua22 ) ).LE.avb22 /
287 $ ( abs1( vb21 )+abs1( vb22 ) ) ) THEN
288 CALL zlartg( -dconjg( ua21 ), dconjg( ua22 ), csq,
289 $ snq,
290 $ r )
291 ELSE
292 CALL zlartg( -dconjg( vb21 ), dconjg( vb22 ), csq,
293 $ snq,
294 $ r )
295 END IF
296
297 csu = snl
298 snu = d1*csl
299 csv = snr
300 snv = d1*csr
301
302 END IF
303
304 ELSE
305
306
307
308
309
310
311 a = a1*b3
312 d = a3*b1
313 c = a2*b3 - a3*b2
314 fc = abs( c )
315
316
317
318
319 d1 = one
320 IF( fc.NE.zero )
321 $ d1 = c / fc
322
323
324
325
326
327
328 CALL dlasv2( a, fc, d, s1, s2, snr, csr, snl, csl )
329
330 IF( abs( csr ).GE.abs( snr ) .OR. abs( csl ).GE.abs( snl ) )
331 $ THEN
332
333
334
335
336 ua21 = -d1*snr*a1 + csr*a2
337 ua22r = csr*a3
338
339 vb21 = -d1*snl*b1 + csl*b2
340 vb22r = csl*b3
341
342 aua21 = abs( snr )*abs( a1 ) + abs( csr )*abs1( a2 )
343 avb21 = abs( snl )*abs( b1 ) + abs( csl )*abs1( b2 )
344
345
346
347 IF( ( abs1( ua21 )+abs( ua22r ) ).EQ.zero ) THEN
348 CALL zlartg( dcmplx( vb22r ), vb21, csq, snq, r )
349 ELSE IF( ( abs1( vb21 )+abs( vb22r ) ).EQ.zero ) THEN
350 CALL zlartg( dcmplx( ua22r ), ua21, csq, snq, r )
351 ELSE IF( aua21 / ( abs1( ua21 )+abs( ua22r ) ).LE.avb21 /
352 $ ( abs1( vb21 )+abs( vb22r ) ) ) THEN
353 CALL zlartg( dcmplx( ua22r ), ua21, csq, snq, r )
354 ELSE
355 CALL zlartg( dcmplx( vb22r ), vb21, csq, snq, r )
356 END IF
357
358 csu = csr
359 snu = -dconjg( d1 )*snr
360 csv = csl
361 snv = -dconjg( d1 )*snl
362
363 ELSE
364
365
366
367
368 ua11 = csr*a1 + dconjg( d1 )*snr*a2
369 ua12 = dconjg( d1 )*snr*a3
370
371 vb11 = csl*b1 + dconjg( d1 )*snl*b2
372 vb12 = dconjg( d1 )*snl*b3
373
374 aua11 = abs( csr )*abs( a1 ) + abs( snr )*abs1( a2 )
375 avb11 = abs( csl )*abs( b1 ) + abs( snl )*abs1( b2 )
376
377
378
379 IF( ( abs1( ua11 )+abs1( ua12 ) ).EQ.zero ) THEN
380 CALL zlartg( vb12, vb11, csq, snq, r )
381 ELSE IF( ( abs1( vb11 )+abs1( vb12 ) ).EQ.zero ) THEN
382 CALL zlartg( ua12, ua11, csq, snq, r )
383 ELSE IF( aua11 / ( abs1( ua11 )+abs1( ua12 ) ).LE.avb11 /
384 $ ( abs1( vb11 )+abs1( vb12 ) ) ) THEN
385 CALL zlartg( ua12, ua11, csq, snq, r )
386 ELSE
387 CALL zlartg( vb12, vb11, csq, snq, r )
388 END IF
389
390 csu = snr
391 snu = dconjg( d1 )*csr
392 csv = snl
393 snv = dconjg( d1 )*csl
394
395 END IF
396
397 END IF
398
399 RETURN
400
401
402
subroutine zlartg(f, g, c, s, r)
ZLARTG generates a plane rotation with real cosine and complex sine.
subroutine dlasv2(f, g, h, ssmin, ssmax, snr, csr, snl, csl)
DLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix.