LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
dorbdb1.f
Go to the documentation of this file.
1 *> \brief \b DORBDB1
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DORBDB1 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dorbdb1.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dorbdb1.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dorbdb1.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DORBDB1( 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 * DOUBLE PRECISION TAUP1(*), TAUP2(*), TAUQ1(*), WORK(*),
30 * $ X11(LDX11,*), X21(LDX21,*)
31 * ..
32 *
33 *
34 *> \par Purpose:
35 *> =============
36 *>
37 *>\verbatim
38 *>
39 *> DORBDB1 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. Q must be no larger than P,
49 *> M-P, or M-Q. Routines DORBDB2, DORBDB3, and DORBDB4 handle cases in
50 *> which Q is not the minimum dimension.
51 *>
52 *> The orthogonal 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 Q-by-Q bidiagonal matrices represented implicitly by
57 *> 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.
74 *> \endverbatim
75 *>
76 *> \param[in] Q
77 *> \verbatim
78 *> Q is INTEGER
79 *> The number of columns in X11 and X21. 0 <= Q <=
80 *> MIN(P,M-P,M-Q).
81 *> \endverbatim
82 *>
83 *> \param[in,out] X11
84 *> \verbatim
85 *> X11 is DOUBLE PRECISION array, dimension (LDX11,Q)
86 *> On entry, the top block of the matrix X to be reduced. On
87 *> exit, the columns of tril(X11) specify reflectors for P1 and
88 *> the rows of triu(X11,1) specify reflectors for Q1.
89 *> \endverbatim
90 *>
91 *> \param[in] LDX11
92 *> \verbatim
93 *> LDX11 is INTEGER
94 *> The leading dimension of X11. LDX11 >= P.
95 *> \endverbatim
96 *>
97 *> \param[in,out] X21
98 *> \verbatim
99 *> X21 is DOUBLE PRECISION array, dimension (LDX21,Q)
100 *> On entry, the bottom block of the matrix X to be reduced. On
101 *> exit, the columns of tril(X21) specify reflectors for P2.
102 *> \endverbatim
103 *>
104 *> \param[in] LDX21
105 *> \verbatim
106 *> LDX21 is INTEGER
107 *> The leading dimension of X21. LDX21 >= M-P.
108 *> \endverbatim
109 *>
110 *> \param[out] THETA
111 *> \verbatim
112 *> THETA is DOUBLE PRECISION array, dimension (Q)
113 *> The entries of the bidiagonal blocks B11, B21 are defined by
114 *> THETA and PHI. See Further Details.
115 *> \endverbatim
116 *>
117 *> \param[out] PHI
118 *> \verbatim
119 *> PHI is DOUBLE PRECISION array, dimension (Q-1)
120 *> The entries of the bidiagonal blocks B11, B21 are defined by
121 *> THETA and PHI. See Further Details.
122 *> \endverbatim
123 *>
124 *> \param[out] TAUP1
125 *> \verbatim
126 *> TAUP1 is DOUBLE PRECISION array, dimension (P)
127 *> The scalar factors of the elementary reflectors that define
128 *> P1.
129 *> \endverbatim
130 *>
131 *> \param[out] TAUP2
132 *> \verbatim
133 *> TAUP2 is DOUBLE PRECISION array, dimension (M-P)
134 *> The scalar factors of the elementary reflectors that define
135 *> P2.
136 *> \endverbatim
137 *>
138 *> \param[out] TAUQ1
139 *> \verbatim
140 *> TAUQ1 is DOUBLE PRECISION array, dimension (Q)
141 *> The scalar factors of the elementary reflectors that define
142 *> Q1.
143 *> \endverbatim
144 *>
145 *> \param[out] WORK
146 *> \verbatim
147 *> WORK is DOUBLE PRECISION array, dimension (LWORK)
148 *> \endverbatim
149 *>
150 *> \param[in] LWORK
151 *> \verbatim
152 *> LWORK is INTEGER
153 *> The dimension of the array WORK. LWORK >= M-Q.
154 *>
155 *> If LWORK = -1, then a workspace query is assumed; the routine
156 *> only calculates the optimal size of the WORK array, returns
157 *> this value as the first entry of the WORK array, and no error
158 *> message related to LWORK is issued by XERBLA.
159 *> \endverbatim
160 *>
161 *> \param[out] INFO
162 *> \verbatim
163 *> INFO is INTEGER
164 *> = 0: successful exit.
165 *> < 0: if INFO = -i, the i-th argument had an illegal value.
166 *> \endverbatim
167 *>
168 *
169 * Authors:
170 * ========
171 *
172 *> \author Univ. of Tennessee
173 *> \author Univ. of California Berkeley
174 *> \author Univ. of Colorado Denver
175 *> \author NAG Ltd.
176 *
177 *> \date July 2012
178 *
179 *> \ingroup doubleOTHERcomputational
180 *
181 *> \par Further Details:
182 * =====================
183 *>
184 *> \verbatim
185 *>
186 *> The upper-bidiagonal blocks B11, B21 are represented implicitly by
187 *> angles THETA(1), ..., THETA(Q) and PHI(1), ..., PHI(Q-1). Every entry
188 *> in each bidiagonal band is a product of a sine or cosine of a THETA
189 *> with a sine or cosine of a PHI. See [1] or DORCSD for details.
190 *>
191 *> P1, P2, and Q1 are represented as products of elementary reflectors.
192 *> See DORCSD2BY1 for details on generating P1, P2, and Q1 using DORGQR
193 *> and DORGLQ.
194 *> \endverbatim
195 *
196 *> \par References:
197 * ================
198 *>
199 *> [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
200 *> Algorithms, 50(1):33-65, 2009.
201 *>
202 * =====================================================================
203  SUBROUTINE dorbdb1( M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI,
204  $ taup1, taup2, tauq1, work, lwork, info )
205 *
206 * -- LAPACK computational routine (version 3.6.1) --
207 * -- LAPACK is a software package provided by Univ. of Tennessee, --
208 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
209 * July 2012
210 *
211 * .. Scalar Arguments ..
212  INTEGER INFO, LWORK, M, P, Q, LDX11, LDX21
213 * ..
214 * .. Array Arguments ..
215  DOUBLE PRECISION PHI(*), THETA(*)
216  DOUBLE PRECISION TAUP1(*), TAUP2(*), TAUQ1(*), WORK(*),
217  $ x11(ldx11,*), x21(ldx21,*)
218 * ..
219 *
220 * ====================================================================
221 *
222 * .. Parameters ..
223  DOUBLE PRECISION ONE
224  parameter ( one = 1.0d0 )
225 * ..
226 * .. Local Scalars ..
227  DOUBLE PRECISION C, S
228  INTEGER CHILDINFO, I, ILARF, IORBDB5, LLARF, LORBDB5,
229  $ lworkmin, lworkopt
230  LOGICAL LQUERY
231 * ..
232 * .. External Subroutines ..
233  EXTERNAL dlarf, dlarfgp, dorbdb5, drot, xerbla
234 * ..
235 * .. External Functions ..
236  DOUBLE PRECISION DNRM2
237  EXTERNAL dnrm2
238 * ..
239 * .. Intrinsic Function ..
240  INTRINSIC atan2, cos, max, sin, sqrt
241 * ..
242 * .. Executable Statements ..
243 *
244 * Test input arguments
245 *
246  info = 0
247  lquery = lwork .EQ. -1
248 *
249  IF( m .LT. 0 ) THEN
250  info = -1
251  ELSE IF( p .LT. q .OR. m-p .LT. q ) THEN
252  info = -2
253  ELSE IF( q .LT. 0 .OR. m-q .LT. q ) THEN
254  info = -3
255  ELSE IF( ldx11 .LT. max( 1, p ) ) THEN
256  info = -5
257  ELSE IF( ldx21 .LT. max( 1, m-p ) ) THEN
258  info = -7
259  END IF
260 *
261 * Compute workspace
262 *
263  IF( info .EQ. 0 ) THEN
264  ilarf = 2
265  llarf = max( p-1, m-p-1, q-1 )
266  iorbdb5 = 2
267  lorbdb5 = q-2
268  lworkopt = max( ilarf+llarf-1, iorbdb5+lorbdb5-1 )
269  lworkmin = lworkopt
270  work(1) = lworkopt
271  IF( lwork .LT. lworkmin .AND. .NOT.lquery ) THEN
272  info = -14
273  END IF
274  END IF
275  IF( info .NE. 0 ) THEN
276  CALL xerbla( 'DORBDB1', -info )
277  RETURN
278  ELSE IF( lquery ) THEN
279  RETURN
280  END IF
281 *
282 * Reduce columns 1, ..., Q of X11 and X21
283 *
284  DO i = 1, q
285 *
286  CALL dlarfgp( p-i+1, x11(i,i), x11(i+1,i), 1, taup1(i) )
287  CALL dlarfgp( m-p-i+1, x21(i,i), x21(i+1,i), 1, taup2(i) )
288  theta(i) = atan2( x21(i,i), x11(i,i) )
289  c = cos( theta(i) )
290  s = sin( theta(i) )
291  x11(i,i) = one
292  x21(i,i) = one
293  CALL dlarf( 'L', p-i+1, q-i, x11(i,i), 1, taup1(i), x11(i,i+1),
294  $ ldx11, work(ilarf) )
295  CALL dlarf( 'L', m-p-i+1, q-i, x21(i,i), 1, taup2(i),
296  $ x21(i,i+1), ldx21, work(ilarf) )
297 *
298  IF( i .LT. q ) THEN
299  CALL drot( q-i, x11(i,i+1), ldx11, x21(i,i+1), ldx21, c, s )
300  CALL dlarfgp( q-i, x21(i,i+1), x21(i,i+2), ldx21, tauq1(i) )
301  s = x21(i,i+1)
302  x21(i,i+1) = one
303  CALL dlarf( 'R', p-i, q-i, x21(i,i+1), ldx21, tauq1(i),
304  $ x11(i+1,i+1), ldx11, work(ilarf) )
305  CALL dlarf( 'R', m-p-i, q-i, x21(i,i+1), ldx21, tauq1(i),
306  $ x21(i+1,i+1), ldx21, work(ilarf) )
307  c = sqrt( dnrm2( p-i, x11(i+1,i+1), 1 )**2
308  $ + dnrm2( m-p-i, x21(i+1,i+1), 1 )**2 )
309  phi(i) = atan2( s, c )
310  CALL dorbdb5( p-i, m-p-i, q-i-1, x11(i+1,i+1), 1,
311  $ x21(i+1,i+1), 1, x11(i+1,i+2), ldx11,
312  $ x21(i+1,i+2), ldx21, work(iorbdb5), lorbdb5,
313  $ childinfo )
314  END IF
315 *
316  END DO
317 *
318  RETURN
319 *
320 * End of DORBDB1
321 *
322  END
323 
subroutine dlarfgp(N, ALPHA, X, INCX, TAU)
DLARFGP generates an elementary reflector (Householder matrix) with non-negative beta.
Definition: dlarfgp.f:106
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
Definition: drot.f:53
subroutine dlarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
DLARF applies an elementary reflector to a general rectangular matrix.
Definition: dlarf.f:126
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dorbdb1(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
DORBDB1
Definition: dorbdb1.f:205
subroutine dorbdb5(M1, M2, N, X1, INCX1, X2, INCX2, Q1, LDQ1, Q2, LDQ2, WORK, LWORK, INFO)
DORBDB5
Definition: dorbdb5.f:158