LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
cunbdb2.f
Go to the documentation of this file.
1*> \brief \b CUNBDB2
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download CUNBDB2 + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cunbdb2.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cunbdb2.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cunbdb2.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE CUNBDB2( M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI,
20* TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO )
21*
22* .. Scalar Arguments ..
23* INTEGER INFO, LWORK, M, P, Q, LDX11, LDX21
24* ..
25* .. Array Arguments ..
26* REAL PHI(*), THETA(*)
27* COMPLEX TAUP1(*), TAUP2(*), TAUQ1(*), WORK(*),
28* $ X11(LDX11,*), X21(LDX21,*)
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*>\verbatim
36*>
37*> CUNBDB2 simultaneously bidiagonalizes the blocks of a tall and skinny
38*> matrix X with orthonormal columns:
39*>
40*> [ B11 ]
41*> [ X11 ] [ P1 | ] [ 0 ]
42*> [-----] = [---------] [-----] Q1**T .
43*> [ X21 ] [ | P2 ] [ B21 ]
44*> [ 0 ]
45*>
46*> X11 is P-by-Q, and X21 is (M-P)-by-Q. P must be no larger than M-P,
47*> Q, or M-Q. Routines CUNBDB1, CUNBDB3, and CUNBDB4 handle cases in
48*> which P is not the minimum dimension.
49*>
50*> The unitary matrices P1, P2, and Q1 are P-by-P, (M-P)-by-(M-P),
51*> and (M-Q)-by-(M-Q), respectively. They are represented implicitly by
52*> Householder vectors.
53*>
54*> B11 and B12 are P-by-P bidiagonal matrices represented implicitly by
55*> angles THETA, PHI.
56*>
57*>\endverbatim
58*
59* Arguments:
60* ==========
61*
62*> \param[in] M
63*> \verbatim
64*> M is INTEGER
65*> The number of rows X11 plus the number of rows in X21.
66*> \endverbatim
67*>
68*> \param[in] P
69*> \verbatim
70*> P is INTEGER
71*> The number of rows in X11. 0 <= P <= min(M-P,Q,M-Q).
72*> \endverbatim
73*>
74*> \param[in] Q
75*> \verbatim
76*> Q is INTEGER
77*> The number of columns in X11 and X21. 0 <= Q <= M.
78*> \endverbatim
79*>
80*> \param[in,out] X11
81*> \verbatim
82*> X11 is COMPLEX array, dimension (LDX11,Q)
83*> On entry, the top block of the matrix X to be reduced. On
84*> exit, the columns of tril(X11) specify reflectors for P1 and
85*> the rows of triu(X11,1) specify reflectors for Q1.
86*> \endverbatim
87*>
88*> \param[in] LDX11
89*> \verbatim
90*> LDX11 is INTEGER
91*> The leading dimension of X11. LDX11 >= P.
92*> \endverbatim
93*>
94*> \param[in,out] X21
95*> \verbatim
96*> X21 is COMPLEX array, dimension (LDX21,Q)
97*> On entry, the bottom block of the matrix X to be reduced. On
98*> exit, the columns of tril(X21) specify reflectors for P2.
99*> \endverbatim
100*>
101*> \param[in] LDX21
102*> \verbatim
103*> LDX21 is INTEGER
104*> The leading dimension of X21. LDX21 >= M-P.
105*> \endverbatim
106*>
107*> \param[out] THETA
108*> \verbatim
109*> THETA is REAL array, dimension (Q)
110*> The entries of the bidiagonal blocks B11, B21 are defined by
111*> THETA and PHI. See Further Details.
112*> \endverbatim
113*>
114*> \param[out] PHI
115*> \verbatim
116*> PHI is REAL array, dimension (Q-1)
117*> The entries of the bidiagonal blocks B11, B21 are defined by
118*> THETA and PHI. See Further Details.
119*> \endverbatim
120*>
121*> \param[out] TAUP1
122*> \verbatim
123*> TAUP1 is COMPLEX array, dimension (P-1)
124*> The scalar factors of the elementary reflectors that define
125*> P1.
126*> \endverbatim
127*>
128*> \param[out] TAUP2
129*> \verbatim
130*> TAUP2 is COMPLEX array, dimension (Q)
131*> The scalar factors of the elementary reflectors that define
132*> P2.
133*> \endverbatim
134*>
135*> \param[out] TAUQ1
136*> \verbatim
137*> TAUQ1 is COMPLEX array, dimension (Q)
138*> The scalar factors of the elementary reflectors that define
139*> Q1.
140*> \endverbatim
141*>
142*> \param[out] WORK
143*> \verbatim
144*> WORK is COMPLEX array, dimension (LWORK)
145*> \endverbatim
146*>
147*> \param[in] LWORK
148*> \verbatim
149*> LWORK is INTEGER
150*> The dimension of the array WORK. LWORK >= M-Q.
151*>
152*> If LWORK = -1, then a workspace query is assumed; the routine
153*> only calculates the optimal size of the WORK array, returns
154*> this value as the first entry of the WORK array, and no error
155*> message related to LWORK is issued by XERBLA.
156*> \endverbatim
157*>
158*> \param[out] INFO
159*> \verbatim
160*> INFO is INTEGER
161*> = 0: successful exit.
162*> < 0: if INFO = -i, the i-th argument had an illegal value.
163*> \endverbatim
164*>
165*
166* Authors:
167* ========
168*
169*> \author Univ. of Tennessee
170*> \author Univ. of California Berkeley
171*> \author Univ. of Colorado Denver
172*> \author NAG Ltd.
173*
174*> \ingroup unbdb2
175*
176*> \par Further Details:
177* =====================
178*>
179*> \verbatim
180*>
181*> The upper-bidiagonal blocks B11, B21 are represented implicitly by
182*> angles THETA(1), ..., THETA(Q) and PHI(1), ..., PHI(Q-1). Every entry
183*> in each bidiagonal band is a product of a sine or cosine of a THETA
184*> with a sine or cosine of a PHI. See [1] or CUNCSD for details.
185*>
186*> P1, P2, and Q1 are represented as products of elementary reflectors.
187*> See CUNCSD2BY1 for details on generating P1, P2, and Q1 using CUNGQR
188*> and CUNGLQ.
189*> \endverbatim
190*
191*> \par References:
192* ================
193*>
194*> [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
195*> Algorithms, 50(1):33-65, 2009.
196*>
197* =====================================================================
198 SUBROUTINE cunbdb2( M, P, Q, X11, LDX11, X21, LDX21, THETA,
199 $ PHI,
200 $ TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO )
201*
202* -- LAPACK computational routine --
203* -- LAPACK is a software package provided by Univ. of Tennessee, --
204* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
205*
206* .. Scalar Arguments ..
207 INTEGER INFO, LWORK, M, P, Q, LDX11, LDX21
208* ..
209* .. Array Arguments ..
210 REAL PHI(*), THETA(*)
211 COMPLEX TAUP1(*), TAUP2(*), TAUQ1(*), WORK(*),
212 $ x11(ldx11,*), x21(ldx21,*)
213* ..
214*
215* ====================================================================
216*
217* .. Parameters ..
218 COMPLEX NEGONE
219 PARAMETER ( NEGONE = (-1.0e0,0.0e0) )
220* ..
221* .. Local Scalars ..
222 REAL C, S
223 INTEGER CHILDINFO, I, ILARF, IORBDB5, LLARF, LORBDB5,
224 $ lworkmin, lworkopt
225 LOGICAL LQUERY
226* ..
227* .. External Subroutines ..
228 EXTERNAL clarf1f, clarfgp, cunbdb5, csrot, cscal,
229 $ clacgv,
230 $ xerbla
231* ..
232* .. External Functions ..
233 REAL SCNRM2, SROUNDUP_LWORK
234 EXTERNAL SCNRM2, SROUNDUP_LWORK
235* ..
236* .. Intrinsic Function ..
237 INTRINSIC atan2, cos, max, sin, sqrt
238* ..
239* .. Executable Statements ..
240*
241* Test input arguments
242*
243 info = 0
244 lquery = lwork .EQ. -1
245*
246 IF( m .LT. 0 ) THEN
247 info = -1
248 ELSE IF( p .LT. 0 .OR. p .GT. m-p ) THEN
249 info = -2
250 ELSE IF( q .LT. 0 .OR. q .LT. p .OR. m-q .LT. p ) THEN
251 info = -3
252 ELSE IF( ldx11 .LT. max( 1, p ) ) THEN
253 info = -5
254 ELSE IF( ldx21 .LT. max( 1, m-p ) ) THEN
255 info = -7
256 END IF
257*
258* Compute workspace
259*
260 IF( info .EQ. 0 ) THEN
261 ilarf = 2
262 llarf = max( p-1, m-p, q-1 )
263 iorbdb5 = 2
264 lorbdb5 = q-1
265 lworkopt = max( ilarf+llarf-1, iorbdb5+lorbdb5-1 )
266 lworkmin = lworkopt
267 work(1) = sroundup_lwork(lworkopt)
268 IF( lwork .LT. lworkmin .AND. .NOT.lquery ) THEN
269 info = -14
270 END IF
271 END IF
272 IF( info .NE. 0 ) THEN
273 CALL xerbla( 'CUNBDB2', -info )
274 RETURN
275 ELSE IF( lquery ) THEN
276 RETURN
277 END IF
278*
279* Reduce rows 1, ..., P of X11 and X21
280*
281 DO i = 1, p
282*
283 IF( i .GT. 1 ) THEN
284 CALL csrot( q-i+1, x11(i,i), ldx11, x21(i-1,i), ldx21, c,
285 $ s )
286 END IF
287 CALL clacgv( q-i+1, x11(i,i), ldx11 )
288 CALL clarfgp( q-i+1, x11(i,i), x11(i,i+1), ldx11, tauq1(i) )
289 c = real( x11(i,i) )
290 CALL clarf1f( 'R', p-i, q-i+1, x11(i,i), ldx11, tauq1(i),
291 $ x11(i+1,i), ldx11, work(ilarf) )
292 CALL clarf1f( 'R', m-p-i+1, q-i+1, x11(i,i), ldx11,
293 $ tauq1(i), x21(i,i), ldx21, work(ilarf) )
294 CALL clacgv( q-i+1, x11(i,i), ldx11 )
295 s = sqrt( scnrm2( p-i, x11(i+1,i), 1 )**2
296 $ + scnrm2( m-p-i+1, x21(i,i), 1 )**2 )
297 theta(i) = atan2( s, c )
298*
299 CALL cunbdb5( p-i, m-p-i+1, q-i, x11(i+1,i), 1, x21(i,i), 1,
300 $ x11(i+1,i+1), ldx11, x21(i,i+1), ldx21,
301 $ work(iorbdb5), lorbdb5, childinfo )
302 CALL cscal( p-i, negone, x11(i+1,i), 1 )
303 CALL clarfgp( m-p-i+1, x21(i,i), x21(i+1,i), 1, taup2(i) )
304 IF( i .LT. p ) THEN
305 CALL clarfgp( p-i, x11(i+1,i), x11(i+2,i), 1, taup1(i) )
306 phi(i) = atan2( real( x11(i+1,i) ), real( x21(i,i) ) )
307 c = cos( phi(i) )
308 s = sin( phi(i) )
309 CALL clarf1f( 'L', p-i, q-i, x11(i+1,i), 1,
310 $ conjg(taup1(i)), x11(i+1,i+1), ldx11,
311 $ work(ilarf) )
312 END IF
313 CALL clarf1f( 'L', m-p-i+1, q-i, x21(i,i), 1,
314 $ conjg(taup2(i)), x21(i,i+1), ldx21,
315 $ work(ilarf) )
316*
317 END DO
318*
319* Reduce the bottom-right portion of X21 to the identity matrix
320*
321 DO i = p + 1, q
322 CALL clarfgp( m-p-i+1, x21(i,i), x21(i+1,i), 1, taup2(i) )
323 CALL clarf1f( 'L', m-p-i+1, q-i, x21(i,i), 1,
324 $ conjg(taup2(i)), x21(i,i+1), ldx21,
325 $ work(ilarf) )
326 END DO
327*
328 RETURN
329*
330* End of CUNBDB2
331*
332 END
333
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
clarf1f
subroutine clarf1f(side, m, n, v, incv, tau, c, ldc, work)
CLARF1F applies an elementary reflector to a general rectangular
Definition clarf1f.f:126
clacgv
subroutine clacgv(n, x, incx)
CLACGV conjugates a complex vector.
Definition clacgv.f:72
clarfgp
subroutine clarfgp(n, alpha, x, incx, tau)
CLARFGP generates an elementary reflector (Householder matrix) with non-negative beta.
Definition clarfgp.f:102
csrot
subroutine csrot(n, cx, incx, cy, incy, c, s)
CSROT
Definition csrot.f:98
cscal
subroutine cscal(n, ca, cx, incx)
CSCAL
Definition cscal.f:78
cunbdb2
subroutine cunbdb2(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, work, lwork, info)
CUNBDB2
Definition cunbdb2.f:201
cunbdb5
subroutine cunbdb5(m1, m2, n, x1, incx1, x2, incx2, q1, ldq1, q2, ldq2, work, lwork, info)
CUNBDB5
Definition cunbdb5.f:155