LAPACK  3.6.1 LAPACK: Linear Algebra PACKage
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 *> \date November 2011
154 *
155 *> \ingroup complexOTHERauxiliary
156 *
157 * =====================================================================
158  SUBROUTINE clags2( UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU, CSV,
159  \$ snv, csq, snq )
160 *
161 * -- LAPACK auxiliary routine (version 3.4.0) --
162 * -- LAPACK is a software package provided by Univ. of Tennessee, --
163 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
164 * November 2011
165 *
166 * .. Scalar Arguments ..
167  LOGICAL UPPER
168  REAL A1, A3, B1, B3, CSQ, CSU, CSV
169  COMPLEX A2, B2, SNQ, SNU, SNV
170 * ..
171 *
172 * =====================================================================
173 *
174 * .. Parameters ..
175  REAL ZERO, ONE
176  parameter ( zero = 0.0e+0, one = 1.0e+0 )
177 * ..
178 * .. Local Scalars ..
179  REAL A, AUA11, AUA12, AUA21, AUA22, AVB11, AVB12,
180  \$ avb21, avb22, csl, csr, d, fb, fc, s1, s2, snl,
181  \$ snr, ua11r, ua22r, vb11r, vb22r
182  COMPLEX B, C, D1, R, T, UA11, UA12, UA21, UA22, VB11,
183  \$ vb12, vb21, vb22
184 * ..
185 * .. External Subroutines ..
186  EXTERNAL clartg, slasv2
187 * ..
188 * .. Intrinsic Functions ..
189  INTRINSIC abs, aimag, cmplx, conjg, real
190 * ..
191 * .. Statement Functions ..
192  REAL ABS1
193 * ..
194 * .. Statement Function definitions ..
195  abs1( t ) = abs( REAL( T ) ) + abs( AIMAG( t ) )
196 * ..
197 * .. Executable Statements ..
198 *
199  IF( upper ) THEN
200 *
201 * Input matrices A and B are upper triangular matrices
202 *
203 * Form matrix C = A*adj(B) = ( a b )
204 * ( 0 d )
205 *
206  a = a1*b3
207  d = a3*b1
208  b = a2*b1 - a1*b2
209  fb = abs( b )
210 *
211 * Transform complex 2-by-2 matrix C to real matrix by unitary
212 * diagonal matrix diag(1,D1).
213 *
214  d1 = one
215  IF( fb.NE.zero )
216  \$ d1 = b / fb
217 *
218 * The SVD of real 2 by 2 triangular C
219 *
220 * ( CSL -SNL )*( A B )*( CSR SNR ) = ( R 0 )
221 * ( SNL CSL ) ( 0 D ) ( -SNR CSR ) ( 0 T )
222 *
223  CALL slasv2( a, fb, d, s1, s2, snr, csr, snl, csl )
224 *
225  IF( abs( csl ).GE.abs( snl ) .OR. abs( csr ).GE.abs( snr ) )
226  \$ THEN
227 *
228 * Compute the (1,1) and (1,2) elements of U**H *A and V**H *B,
229 * and (1,2) element of |U|**H *|A| and |V|**H *|B|.
230 *
231  ua11r = csl*a1
232  ua12 = csl*a2 + d1*snl*a3
233 *
234  vb11r = csr*b1
235  vb12 = csr*b2 + d1*snr*b3
236 *
237  aua12 = abs( csl )*abs1( a2 ) + abs( snl )*abs( a3 )
238  avb12 = abs( csr )*abs1( b2 ) + abs( snr )*abs( b3 )
239 *
240 * zero (1,2) elements of U**H *A and V**H *B
241 *
242  IF( ( abs( ua11r )+abs1( ua12 ) ).EQ.zero ) THEN
243  CALL clartg( -cmplx( vb11r ), conjg( vb12 ), csq, snq,
244  \$ r )
245  ELSE IF( ( abs( vb11r )+abs1( vb12 ) ).EQ.zero ) THEN
246  CALL clartg( -cmplx( ua11r ), conjg( ua12 ), csq, snq,
247  \$ r )
248  ELSE IF( aua12 / ( abs( ua11r )+abs1( ua12 ) ).LE.avb12 /
249  \$ ( abs( vb11r )+abs1( vb12 ) ) ) THEN
250  CALL clartg( -cmplx( ua11r ), conjg( ua12 ), csq, snq,
251  \$ r )
252  ELSE
253  CALL clartg( -cmplx( vb11r ), conjg( vb12 ), csq, 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 = -conjg( d1 )*snl*a1
268  ua22 = -conjg( d1 )*snl*a2 + csl*a3
269 *
270  vb21 = -conjg( d1 )*snr*b1
271  vb22 = -conjg( 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 clartg( -conjg( vb21 ), conjg( vb22 ), csq, snq, r )
280  ELSE IF( ( abs1( vb21 )+abs( vb22 ) ).EQ.zero ) THEN
281  CALL clartg( -conjg( ua21 ), conjg( ua22 ), csq, snq, r )
282  ELSE IF( aua22 / ( abs1( ua21 )+abs1( ua22 ) ).LE.avb22 /
283  \$ ( abs1( vb21 )+abs1( vb22 ) ) ) THEN
284  CALL clartg( -conjg( ua21 ), conjg( ua22 ), csq, snq, r )
285  ELSE
286  CALL clartg( -conjg( vb21 ), conjg( vb22 ), csq, snq, 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 * Input matrices A and B are lower triangular matrices
299 *
300 * Form matrix C = A*adj(B) = ( a 0 )
301 * ( c d )
302 *
303  a = a1*b3
304  d = a3*b1
305  c = a2*b3 - a3*b2
306  fc = abs( c )
307 *
308 * Transform complex 2-by-2 matrix C to real matrix by unitary
309 * diagonal matrix diag(d1,1).
310 *
311  d1 = one
312  IF( fc.NE.zero )
313  \$ d1 = c / fc
314 *
315 * The SVD of real 2 by 2 triangular C
316 *
317 * ( CSL -SNL )*( A 0 )*( CSR SNR ) = ( R 0 )
318 * ( SNL CSL ) ( C D ) ( -SNR CSR ) ( 0 T )
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 * Compute the (2,1) and (2,2) elements of U**H *A and V**H *B,
326 * and (2,1) element of |U|**H *|A| and |V|**H *|B|.
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 * zero (2,1) elements of U**H *A and V**H *B.
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 * Compute the (1,1) and (1,2) elements of U**H *A and V**H *B,
358 * and (1,1) element of |U|**H *|A| and |V|**H *|B|.
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 * zero (1,1) elements of U**H *A and V**H *B, and then swap.
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 * End of CLAGS2
394 *
395  END
subroutine clartg(F, G, CS, SN, R)
CLARTG generates a plane rotation with real cosine and complex sine.
Definition: clartg.f:105
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:140
subroutine clags2(UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU, CSV, SNV, CSQ, SNQ)
CLAGS2
Definition: clags2.f:160