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