LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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 orthonomal 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 complex16OTHERcomputational
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)
XERBLA
Definition: xerbla.f:60
subroutine zdrot(N, ZX, INCX, ZY, INCY, C, S)
ZDROT
Definition: zdrot.f:98
subroutine zlarfgp(N, ALPHA, X, INCX, TAU)
ZLARFGP generates an elementary reflector (Householder matrix) with non-negative beta.
Definition: zlarfgp.f:104
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 zunbdb5(M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, LDQ2, WORK, LWORK, INFO)
ZUNBDB5
Definition: zunbdb5.f:156
subroutine zunbdb3(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
ZUNBDB3
Definition: zunbdb3.f:201