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