LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dorbdb4.f
Go to the documentation of this file.
1*> \brief \b DORBDB4
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download DORBDB4 + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dorbdb4.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dorbdb4.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dorbdb4.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE DORBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI,
20* TAUP1, TAUP2, TAUQ1, PHANTOM, WORK, LWORK,
21* INFO )
22*
23* .. Scalar Arguments ..
24* INTEGER INFO, LWORK, M, P, Q, LDX11, LDX21
25* ..
26* .. Array Arguments ..
27* DOUBLE PRECISION PHI(*), THETA(*)
28* DOUBLE PRECISION PHANTOM(*), TAUP1(*), TAUP2(*), TAUQ1(*),
29* $ WORK(*), X11(LDX11,*), X21(LDX21,*)
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*>\verbatim
37*>
38*> DORBDB4 simultaneously bidiagonalizes the blocks of a tall and skinny
39*> matrix X with orthonormal columns:
40*>
41*> [ B11 ]
42*> [ X11 ] [ P1 | ] [ 0 ]
43*> [-----] = [---------] [-----] Q1**T .
44*> [ X21 ] [ | P2 ] [ B21 ]
45*> [ 0 ]
46*>
47*> X11 is P-by-Q, and X21 is (M-P)-by-Q. M-Q must be no larger than P,
48*> M-P, or Q. Routines DORBDB1, DORBDB2, and DORBDB3 handle cases in
49*> which M-Q is not the minimum dimension.
50*>
51*> The orthogonal matrices P1, P2, and Q1 are P-by-P, (M-P)-by-(M-P),
52*> and (M-Q)-by-(M-Q), respectively. They are represented implicitly by
53*> Householder vectors.
54*>
55*> B11 and B12 are (M-Q)-by-(M-Q) bidiagonal matrices represented
56*> implicitly by angles THETA, PHI.
57*>
58*>\endverbatim
59*
60* Arguments:
61* ==========
62*
63*> \param[in] M
64*> \verbatim
65*> M is INTEGER
66*> The number of rows X11 plus the number of rows in X21.
67*> \endverbatim
68*>
69*> \param[in] P
70*> \verbatim
71*> P is INTEGER
72*> The number of rows in X11. 0 <= P <= M.
73*> \endverbatim
74*>
75*> \param[in] Q
76*> \verbatim
77*> Q is INTEGER
78*> The number of columns in X11 and X21. 0 <= Q <= M and
79*> M-Q <= min(P,M-P,Q).
80*> \endverbatim
81*>
82*> \param[in,out] X11
83*> \verbatim
84*> X11 is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (M-Q)
126*> The scalar factors of the elementary reflectors that define
127*> P1.
128*> \endverbatim
129*>
130*> \param[out] TAUP2
131*> \verbatim
132*> TAUP2 is DOUBLE PRECISION array, dimension (M-Q)
133*> The scalar factors of the elementary reflectors that define
134*> P2.
135*> \endverbatim
136*>
137*> \param[out] TAUQ1
138*> \verbatim
139*> TAUQ1 is DOUBLE PRECISION array, dimension (Q)
140*> The scalar factors of the elementary reflectors that define
141*> Q1.
142*> \endverbatim
143*>
144*> \param[out] PHANTOM
145*> \verbatim
146*> PHANTOM is DOUBLE PRECISION array, dimension (M)
147*> The routine computes an M-by-1 column vector Y that is
148*> orthogonal to the columns of [ X11; X21 ]. PHANTOM(1:P) and
149*> PHANTOM(P+1:M) contain Householder vectors for Y(1:P) and
150*> Y(P+1:M), respectively.
151*> \endverbatim
152*>
153*> \param[out] WORK
154*> \verbatim
155*> WORK is DOUBLE PRECISION array, dimension (LWORK)
156*> \endverbatim
157*>
158*> \param[in] LWORK
159*> \verbatim
160*> LWORK is INTEGER
161*> The dimension of the array WORK. LWORK >= M-Q.
162*>
163*> If LWORK = -1, then a workspace query is assumed; the routine
164*> only calculates the optimal size of the WORK array, returns
165*> this value as the first entry of the WORK array, and no error
166*> message related to LWORK is issued by XERBLA.
167*> \endverbatim
168*>
169*> \param[out] INFO
170*> \verbatim
171*> INFO is INTEGER
172*> = 0: successful exit.
173*> < 0: if INFO = -i, the i-th argument had an illegal value.
174*> \endverbatim
175*
176* Authors:
177* ========
178*
179*> \author Univ. of Tennessee
180*> \author Univ. of California Berkeley
181*> \author Univ. of Colorado Denver
182*> \author NAG Ltd.
183*
184*> \ingroup unbdb4
185*
186*> \par Further Details:
187* =====================
188*>
189*> \verbatim
190*>
191*> The upper-bidiagonal blocks B11, B21 are represented implicitly by
192*> angles THETA(1), ..., THETA(Q) and PHI(1), ..., PHI(Q-1). Every entry
193*> in each bidiagonal band is a product of a sine or cosine of a THETA
194*> with a sine or cosine of a PHI. See [1] or DORCSD for details.
195*>
196*> P1, P2, and Q1 are represented as products of elementary reflectors.
197*> See DORCSD2BY1 for details on generating P1, P2, and Q1 using DORGQR
198*> and DORGLQ.
199*> \endverbatim
200*
201*> \par References:
202* ================
203*>
204*> [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
205*> Algorithms, 50(1):33-65, 2009.
206*>
207* =====================================================================
208 SUBROUTINE dorbdb4( M, P, Q, X11, LDX11, X21, LDX21, THETA,
209 $ PHI,
210 $ TAUP1, TAUP2, TAUQ1, PHANTOM, WORK, LWORK,
211 $ INFO )
212*
213* -- LAPACK computational routine --
214* -- LAPACK is a software package provided by Univ. of Tennessee, --
215* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
216*
217* .. Scalar Arguments ..
218 INTEGER INFO, LWORK, M, P, Q, LDX11, LDX21
219* ..
220* .. Array Arguments ..
221 DOUBLE PRECISION PHI(*), THETA(*)
222 DOUBLE PRECISION PHANTOM(*), TAUP1(*), TAUP2(*), TAUQ1(*),
223 $ WORK(*), X11(LDX11,*), X21(LDX21,*)
224* ..
225*
226* ====================================================================
227*
228* .. Parameters ..
229 DOUBLE PRECISION NEGONE, ONE, ZERO
230 PARAMETER ( NEGONE = -1.0d0, one = 1.0d0, zero = 0.0d0 )
231* ..
232* .. Local Scalars ..
233 DOUBLE PRECISION C, S
234 INTEGER CHILDINFO, I, ILARF, IORBDB5, J, LLARF,
235 $ LORBDB5, LWORKMIN, LWORKOPT
236 LOGICAL LQUERY
237* ..
238* .. External Subroutines ..
239 EXTERNAL dlarf, dlarfgp, dorbdb5, drot, dscal,
240 $ xerbla
241* ..
242* .. External Functions ..
243 DOUBLE PRECISION DNRM2
244 EXTERNAL DNRM2
245* ..
246* .. Intrinsic Function ..
247 INTRINSIC atan2, cos, max, sin, sqrt
248* ..
249* .. Executable Statements ..
250*
251* Test input arguments
252*
253 info = 0
254 lquery = lwork .EQ. -1
255*
256 IF( m .LT. 0 ) THEN
257 info = -1
258 ELSE IF( p .LT. m-q .OR. m-p .LT. m-q ) THEN
259 info = -2
260 ELSE IF( q .LT. m-q .OR. q .GT. m ) THEN
261 info = -3
262 ELSE IF( ldx11 .LT. max( 1, p ) ) THEN
263 info = -5
264 ELSE IF( ldx21 .LT. max( 1, m-p ) ) THEN
265 info = -7
266 END IF
267*
268* Compute workspace
269*
270 IF( info .EQ. 0 ) THEN
271 ilarf = 2
272 llarf = max( q-1, p-1, m-p-1 )
273 iorbdb5 = 2
274 lorbdb5 = q
275 lworkopt = ilarf + llarf - 1
276 lworkopt = max( lworkopt, iorbdb5 + lorbdb5 - 1 )
277 lworkmin = lworkopt
278 work(1) = lworkopt
279 IF( lwork .LT. lworkmin .AND. .NOT.lquery ) THEN
280 info = -14
281 END IF
282 END IF
283 IF( info .NE. 0 ) THEN
284 CALL xerbla( 'DORBDB4', -info )
285 RETURN
286 ELSE IF( lquery ) THEN
287 RETURN
288 END IF
289*
290* Reduce columns 1, ..., M-Q of X11 and X21
291*
292 DO i = 1, m-q
293*
294 IF( i .EQ. 1 ) THEN
295 DO j = 1, m
296 phantom(j) = zero
297 END DO
298 CALL dorbdb5( p, m-p, q, phantom(1), 1, phantom(p+1), 1,
299 $ x11, ldx11, x21, ldx21, work(iorbdb5),
300 $ lorbdb5, childinfo )
301 CALL dscal( p, negone, phantom(1), 1 )
302 CALL dlarfgp( p, phantom(1), phantom(2), 1, taup1(1) )
303 CALL dlarfgp( m-p, phantom(p+1), phantom(p+2), 1,
304 $ taup2(1) )
305 theta(i) = atan2( phantom(1), phantom(p+1) )
306 c = cos( theta(i) )
307 s = sin( theta(i) )
308 CALL dlarf1f( 'L', p, q, phantom(1), 1, taup1(1), x11,
309 $ ldx11, work(ilarf) )
310 CALL dlarf1f( 'L', m-p, q, phantom(p+1), 1, taup2(1),
311 $ x21, ldx21, work(ilarf) )
312 ELSE
313 CALL dorbdb5( p-i+1, m-p-i+1, q-i+1, x11(i,i-1), 1,
314 $ x21(i,i-1), 1, x11(i,i), ldx11, x21(i,i),
315 $ ldx21, work(iorbdb5), lorbdb5, childinfo )
316 CALL dscal( p-i+1, negone, x11(i,i-1), 1 )
317 CALL dlarfgp( p-i+1, x11(i,i-1), x11(i+1,i-1), 1,
318 $ taup1(i) )
319 CALL dlarfgp( m-p-i+1, x21(i,i-1), x21(i+1,i-1), 1,
320 $ taup2(i) )
321 theta(i) = atan2( x11(i,i-1), x21(i,i-1) )
322 c = cos( theta(i) )
323 s = sin( theta(i) )
324 CALL dlarf1f( 'L', p-i+1, q-i+1, x11(i,i-1), 1, taup1(i),
325 $ x11(i,i), ldx11, work(ilarf) )
326 CALL dlarf1f( 'L', m-p-i+1, q-i+1, x21(i,i-1), 1,
327 $ taup2(i), x21(i,i), ldx21, work(ilarf) )
328 END IF
329*
330 CALL drot( q-i+1, x11(i,i), ldx11, x21(i,i), ldx21, s, -c )
331 CALL dlarfgp( q-i+1, x21(i,i), x21(i,i+1), ldx21, tauq1(i) )
332 c = x21(i,i)
333 CALL dlarf1f( 'R', p-i, q-i+1, x21(i,i), ldx21, tauq1(i),
334 $ x11(i+1,i), ldx11, work(ilarf) )
335 CALL dlarf1f( 'R', m-p-i, q-i+1, x21(i,i), ldx21, tauq1(i),
336 $ x21(i+1,i), ldx21, work(ilarf) )
337 IF( i .LT. m-q ) THEN
338 s = sqrt( dnrm2( p-i, x11(i+1,i), 1 )**2
339 $ + dnrm2( m-p-i, x21(i+1,i), 1 )**2 )
340 phi(i) = atan2( s, c )
341 END IF
342*
343 END DO
344*
345* Reduce the bottom-right portion of X11 to [ I 0 ]
346*
347 DO i = m - q + 1, p
348 CALL dlarfgp( q-i+1, x11(i,i), x11(i,i+1), ldx11, tauq1(i) )
349 CALL dlarf1f( 'R', p-i, q-i+1, x11(i,i), ldx11, tauq1(i),
350 $ x11(i+1,i), ldx11, work(ilarf) )
351 CALL dlarf1f( 'R', q-p, q-i+1, x11(i,i), ldx11, tauq1(i),
352 $ x21(m-q+1,i), ldx21, work(ilarf) )
353 END DO
354*
355* Reduce the bottom-right portion of X21 to [ 0 I ]
356*
357 DO i = p + 1, q
358 CALL dlarfgp( q-i+1, x21(m-q+i-p,i), x21(m-q+i-p,i+1),
359 $ ldx21,
360 $ tauq1(i) )
361 CALL dlarf1f( 'R', q-i, q-i+1, x21(m-q+i-p,i), ldx21,
362 $ tauq1(i),
363 $ x21(m-q+i-p+1,i), ldx21, work(ilarf) )
364 END DO
365*
366 RETURN
367*
368* End of DORBDB4
369*
370 END
371
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dlarf1f
subroutine dlarf1f(side, m, n, v, incv, tau, c, ldc, work)
DLARF1F applies an elementary reflector to a general rectangular
Definition dlarf1f.f:157
dlarf
subroutine dlarf(side, m, n, v, incv, tau, c, ldc, work)
DLARF applies an elementary reflector to a general rectangular matrix.
Definition dlarf.f:122
dlarfgp
subroutine dlarfgp(n, alpha, x, incx, tau)
DLARFGP generates an elementary reflector (Householder matrix) with non-negative beta.
Definition dlarfgp.f:102
drot
subroutine drot(n, dx, incx, dy, incy, c, s)
DROT
Definition drot.f:92
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
dorbdb4
subroutine dorbdb4(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, phantom, work, lwork, info)
DORBDB4
Definition dorbdb4.f:212
dorbdb5
subroutine dorbdb5(m1, m2, n, x1, incx1, x2, incx2, q1, ldq1, q2, ldq2, work, lwork, info)
DORBDB5
Definition dorbdb5.f:155