LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zlags2.f
Go to the documentation of this file.
1*> \brief \b ZLAGS2
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download ZLAGS2 + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlags2.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlags2.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlags2.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE ZLAGS2( UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU, CSV,
20* SNV, CSQ, SNQ )
21*
22* .. Scalar Arguments ..
23* LOGICAL UPPER
24* DOUBLE PRECISION A1, A3, B1, B3, CSQ, CSU, CSV
25* COMPLEX*16 A2, B2, SNQ, SNU, SNV
26* ..
27*
28*
29*> \par Purpose:
30* =============
31*>
32*> \verbatim
33*>
34*> ZLAGS2 computes 2-by-2 unitary matrices U, V and Q, such
35*> that if ( UPPER ) then
36*>
37*> U**H *A*Q = U**H *( A1 A2 )*Q = ( x 0 )
38*> ( 0 A3 ) ( x x )
39*> and
40*> V**H*B*Q = V**H *( B1 B2 )*Q = ( x 0 )
41*> ( 0 B3 ) ( x x )
42*>
43*> or if ( .NOT.UPPER ) then
44*>
45*> U**H *A*Q = U**H *( A1 0 )*Q = ( x x )
46*> ( A2 A3 ) ( 0 x )
47*> and
48*> V**H *B*Q = V**H *( B1 0 )*Q = ( x x )
49*> ( B2 B3 ) ( 0 x )
50*> where
51*>
52*> U = ( CSU SNU ), V = ( CSV SNV ),
53*> ( -SNU**H CSU ) ( -SNV**H CSV )
54*>
55*> Q = ( CSQ SNQ )
56*> ( -SNQ**H CSQ )
57*>
58*> The rows of the transformed A and B are parallel. Moreover, if the
59*> input 2-by-2 matrix A is not zero, then the transformed (1,1) entry
60*> of A is not zero. If the input matrices A and B are both not zero,
61*> then the transformed (2,2) element of B is not zero, except when the
62*> first rows of input A and B are parallel and the second rows are
63*> zero.
64*> \endverbatim
65*
66* Arguments:
67* ==========
68*
69*> \param[in] UPPER
70*> \verbatim
71*> UPPER is LOGICAL
72*> = .TRUE.: the input matrices A and B are upper triangular.
73*> = .FALSE.: the input matrices A and B are lower triangular.
74*> \endverbatim
75*>
76*> \param[in] A1
77*> \verbatim
78*> A1 is DOUBLE PRECISION
79*> \endverbatim
80*>
81*> \param[in] A2
82*> \verbatim
83*> A2 is COMPLEX*16
84*> \endverbatim
85*>
86*> \param[in] A3
87*> \verbatim
88*> A3 is DOUBLE PRECISION
89*> On entry, A1, A2 and A3 are elements of the input 2-by-2
90*> upper (lower) triangular matrix A.
91*> \endverbatim
92*>
93*> \param[in] B1
94*> \verbatim
95*> B1 is DOUBLE PRECISION
96*> \endverbatim
97*>
98*> \param[in] B2
99*> \verbatim
100*> B2 is COMPLEX*16
101*> \endverbatim
102*>
103*> \param[in] B3
104*> \verbatim
105*> B3 is DOUBLE PRECISION
106*> On entry, B1, B2 and B3 are elements of the input 2-by-2
107*> upper (lower) triangular matrix B.
108*> \endverbatim
109*>
110*> \param[out] CSU
111*> \verbatim
112*> CSU is DOUBLE PRECISION
113*> \endverbatim
114*>
115*> \param[out] SNU
116*> \verbatim
117*> SNU is COMPLEX*16
118*> The desired unitary matrix U.
119*> \endverbatim
120*>
121*> \param[out] CSV
122*> \verbatim
123*> CSV is DOUBLE PRECISION
124*> \endverbatim
125*>
126*> \param[out] SNV
127*> \verbatim
128*> SNV is COMPLEX*16
129*> The desired unitary matrix V.
130*> \endverbatim
131*>
132*> \param[out] CSQ
133*> \verbatim
134*> CSQ is DOUBLE PRECISION
135*> \endverbatim
136*>
137*> \param[out] SNQ
138*> \verbatim
139*> SNQ is COMPLEX*16
140*> The desired unitary matrix Q.
141*> \endverbatim
142*
143* Authors:
144* ========
145*
146*> \author Univ. of Tennessee
147*> \author Univ. of California Berkeley
148*> \author Univ. of Colorado Denver
149*> \author NAG Ltd.
150*
151*> \ingroup lags2
152*
153* =====================================================================
154 SUBROUTINE zlags2( UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU,
155 $ CSV,
156 $ SNV, CSQ, SNQ )
157*
158* -- LAPACK auxiliary routine --
159* -- LAPACK is a software package provided by Univ. of Tennessee, --
160* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
161*
162* .. Scalar Arguments ..
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* .. Parameters ..
171 DOUBLE PRECISION ZERO, ONE
172 PARAMETER ( ZERO = 0.0d+0, one = 1.0d+0 )
173* ..
174* .. Local Scalars ..
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* .. External Subroutines ..
182 EXTERNAL dlasv2, zlartg
183* ..
184* .. Intrinsic Functions ..
185 INTRINSIC abs, dble, dcmplx, dconjg, dimag
186* ..
187* .. Statement Functions ..
188 DOUBLE PRECISION ABS1
189* ..
190* .. Statement Function definitions ..
191 abs1( t ) = abs( dble( t ) ) + abs( dimag( t ) )
192* ..
193* .. Executable Statements ..
194*
195 IF( upper ) THEN
196*
197* Input matrices A and B are upper triangular matrices
198*
199* Form matrix C = A*adj(B) = ( a b )
200* ( 0 d )
201*
202 a = a1*b3
203 d = a3*b1
204 b = a2*b1 - a1*b2
205 fb = abs( b )
206*
207* Transform complex 2-by-2 matrix C to real matrix by unitary
208* diagonal matrix diag(1,D1).
209*
210 d1 = one
211 IF( fb.NE.zero )
212 $ d1 = b / fb
213*
214* The SVD of real 2 by 2 triangular C
215*
216* ( CSL -SNL )*( A B )*( CSR SNR ) = ( R 0 )
217* ( SNL CSL ) ( 0 D ) ( -SNR CSR ) ( 0 T )
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* Compute the (1,1) and (1,2) elements of U**H *A and V**H *B,
225* and (1,2) element of |U|**H *|A| and |V|**H *|B|.
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* zero (1,2) elements of U**H *A and V**H *B
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* Compute the (2,1) and (2,2) elements of U**H *A and V**H *B,
265* and (2,2) element of |U|**H *|A| and |V|**H *|B|.
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* zero (2,2) elements of U**H *A and V**H *B, and then swap.
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* Input matrices A and B are lower triangular matrices
307*
308* Form matrix C = A*adj(B) = ( a 0 )
309* ( c d )
310*
311 a = a1*b3
312 d = a3*b1
313 c = a2*b3 - a3*b2
314 fc = abs( c )
315*
316* Transform complex 2-by-2 matrix C to real matrix by unitary
317* diagonal matrix diag(d1,1).
318*
319 d1 = one
320 IF( fc.NE.zero )
321 $ d1 = c / fc
322*
323* The SVD of real 2 by 2 triangular C
324*
325* ( CSL -SNL )*( A 0 )*( CSR SNR ) = ( R 0 )
326* ( SNL CSL ) ( C D ) ( -SNR CSR ) ( 0 T )
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* Compute the (2,1) and (2,2) elements of U**H *A and V**H *B,
334* and (2,1) element of |U|**H *|A| and |V|**H *|B|.
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* zero (2,1) elements of U**H *A and V**H *B.
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* Compute the (1,1) and (1,2) elements of U**H *A and V**H *B,
366* and (1,1) element of |U|**H *|A| and |V|**H *|B|.
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* zero (1,1) elements of U**H *A and V**H *B, and then swap.
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* End of ZLAGS2
402*
403 END
subroutine zlags2(upper, a1, a2, a3, b1, b2, b3, csu, snu, csv, snv, csq, snq)
ZLAGS2
Definition zlags2.f:157
subroutine zlartg(f, g, c, s, r)
ZLARTG generates a plane rotation with real cosine and complex sine.
Definition zlartg.f90:116
subroutine dlasv2(f, g, h, ssmin, ssmax, snr, csr, snl, csl)
DLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix.
Definition dlasv2.f:134