1 SUBROUTINE DORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12,
2 $ X21, LDX21, X22, LDX22, THETA, PHI, TAUP1,
3 $ TAUP2, TAUQ1, TAUQ2, WORK, LWORK, INFO )
4 IMPLICIT NONE
5 *
6 * -- LAPACK routine (version 3.3.0) --
7 *
8 * -- Contributed by Brian Sutton of the Randolph-Macon College --
9 * -- November 2010
10 *
11 * -- LAPACK is a software package provided by Univ. of Tennessee, --
12 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
13 *
14 * .. Scalar Arguments ..
15 CHARACTER SIGNS, TRANS
16 INTEGER INFO, LDX11, LDX12, LDX21, LDX22, LWORK, M, P,
17 $ Q
18 * ..
19 * .. Array Arguments ..
20 DOUBLE PRECISION PHI( * ), THETA( * )
21 DOUBLE PRECISION TAUP1( * ), TAUP2( * ), TAUQ1( * ), TAUQ2( * ),
22 $ WORK( * ), X11( LDX11, * ), X12( LDX12, * ),
23 $ X21( LDX21, * ), X22( LDX22, * )
24 * ..
25 *
26 * Purpose
27 * =======
28 *
29 * DORBDB simultaneously bidiagonalizes the blocks of an M-by-M
30 * partitioned orthogonal matrix X:
31 *
32 * [ B11 | B12 0 0 ]
33 * [ X11 | X12 ] [ P1 | ] [ 0 | 0 -I 0 ] [ Q1 | ]**T
34 * X = [-----------] = [---------] [----------------] [---------] .
35 * [ X21 | X22 ] [ | P2 ] [ B21 | B22 0 0 ] [ | Q2 ]
36 * [ 0 | 0 0 I ]
37 *
38 * X11 is P-by-Q. Q must be no larger than P, M-P, or M-Q. (If this is
39 * not the case, then X must be transposed and/or permuted. This can be
40 * done in constant time using the TRANS and SIGNS options. See DORCSD
41 * for details.)
42 *
43 * The orthogonal matrices P1, P2, Q1, and Q2 are P-by-P, (M-P)-by-
44 * (M-P), Q-by-Q, and (M-Q)-by-(M-Q), respectively. They are
45 * represented implicitly by Householder vectors.
46 *
47 * B11, B12, B21, and B22 are Q-by-Q bidiagonal matrices represented
48 * implicitly by angles THETA, PHI.
49 *
50 * Arguments
51 * =========
52 *
53 * TRANS (input) CHARACTER
54 * = 'T': X, U1, U2, V1T, and V2T are stored in row-major
55 * order;
56 * otherwise: X, U1, U2, V1T, and V2T are stored in column-
57 * major order.
58 *
59 * SIGNS (input) CHARACTER
60 * = 'O': The lower-left block is made nonpositive (the
61 * "other" convention);
62 * otherwise: The upper-right block is made nonpositive (the
63 * "default" convention).
64 *
65 * M (input) INTEGER
66 * The number of rows and columns in X.
67 *
68 * P (input) INTEGER
69 * The number of rows in X11 and X12. 0 <= P <= M.
70 *
71 * Q (input) INTEGER
72 * The number of columns in X11 and X21. 0 <= Q <=
73 * MIN(P,M-P,M-Q).
74 *
75 * X11 (input/output) DOUBLE PRECISION array, dimension (LDX11,Q)
76 * On entry, the top-left block of the orthogonal matrix to be
77 * reduced. On exit, the form depends on TRANS:
78 * If TRANS = 'N', then
79 * the columns of tril(X11) specify reflectors for P1,
80 * the rows of triu(X11,1) specify reflectors for Q1;
81 * else TRANS = 'T', and
82 * the rows of triu(X11) specify reflectors for P1,
83 * the columns of tril(X11,-1) specify reflectors for Q1.
84 *
85 * LDX11 (input) INTEGER
86 * The leading dimension of X11. If TRANS = 'N', then LDX11 >=
87 * P; else LDX11 >= Q.
88 *
89 * X12 (input/output) DOUBLE PRECISION array, dimension (LDX12,M-Q)
90 * On entry, the top-right block of the orthogonal matrix to
91 * be reduced. On exit, the form depends on TRANS:
92 * If TRANS = 'N', then
93 * the rows of triu(X12) specify the first P reflectors for
94 * Q2;
95 * else TRANS = 'T', and
96 * the columns of tril(X12) specify the first P reflectors
97 * for Q2.
98 *
99 * LDX12 (input) INTEGER
100 * The leading dimension of X12. If TRANS = 'N', then LDX12 >=
101 * P; else LDX11 >= M-Q.
102 *
103 * X21 (input/output) DOUBLE PRECISION array, dimension (LDX21,Q)
104 * On entry, the bottom-left block of the orthogonal matrix to
105 * be reduced. On exit, the form depends on TRANS:
106 * If TRANS = 'N', then
107 * the columns of tril(X21) specify reflectors for P2;
108 * else TRANS = 'T', and
109 * the rows of triu(X21) specify reflectors for P2.
110 *
111 * LDX21 (input) INTEGER
112 * The leading dimension of X21. If TRANS = 'N', then LDX21 >=
113 * M-P; else LDX21 >= Q.
114 *
115 * X22 (input/output) DOUBLE PRECISION array, dimension (LDX22,M-Q)
116 * On entry, the bottom-right block of the orthogonal matrix to
117 * be reduced. On exit, the form depends on TRANS:
118 * If TRANS = 'N', then
119 * the rows of triu(X22(Q+1:M-P,P+1:M-Q)) specify the last
120 * M-P-Q reflectors for Q2,
121 * else TRANS = 'T', and
122 * the columns of tril(X22(P+1:M-Q,Q+1:M-P)) specify the last
123 * M-P-Q reflectors for P2.
124 *
125 * LDX22 (input) INTEGER
126 * The leading dimension of X22. If TRANS = 'N', then LDX22 >=
127 * M-P; else LDX22 >= M-Q.
128 *
129 * THETA (output) DOUBLE PRECISION array, dimension (Q)
130 * The entries of the bidiagonal blocks B11, B12, B21, B22 can
131 * be computed from the angles THETA and PHI. See Further
132 * Details.
133 *
134 * PHI (output) DOUBLE PRECISION array, dimension (Q-1)
135 * The entries of the bidiagonal blocks B11, B12, B21, B22 can
136 * be computed from the angles THETA and PHI. See Further
137 * Details.
138 *
139 * TAUP1 (output) DOUBLE PRECISION array, dimension (P)
140 * The scalar factors of the elementary reflectors that define
141 * P1.
142 *
143 * TAUP2 (output) DOUBLE PRECISION array, dimension (M-P)
144 * The scalar factors of the elementary reflectors that define
145 * P2.
146 *
147 * TAUQ1 (output) DOUBLE PRECISION array, dimension (Q)
148 * The scalar factors of the elementary reflectors that define
149 * Q1.
150 *
151 * TAUQ2 (output) DOUBLE PRECISION array, dimension (M-Q)
152 * The scalar factors of the elementary reflectors that define
153 * Q2.
154 *
155 * WORK (workspace) DOUBLE PRECISION array, dimension (LWORK)
156 *
157 * LWORK (input) INTEGER
158 * The dimension of the array WORK. LWORK >= M-Q.
159 *
160 * If LWORK = -1, then a workspace query is assumed; the routine
161 * only calculates the optimal size of the WORK array, returns
162 * this value as the first entry of the WORK array, and no error
163 * message related to LWORK is issued by XERBLA.
164 *
165 * INFO (output) INTEGER
166 * = 0: successful exit.
167 * < 0: if INFO = -i, the i-th argument had an illegal value.
168 *
169 * Further Details
170 * ===============
171 *
172 * The bidiagonal blocks B11, B12, B21, and B22 are represented
173 * implicitly by angles THETA(1), ..., THETA(Q) and PHI(1), ...,
174 * PHI(Q-1). B11 and B21 are upper bidiagonal, while B21 and B22 are
175 * lower bidiagonal. Every entry in each bidiagonal band is a product
176 * of a sine or cosine of a THETA with a sine or cosine of a PHI. See
177 * [1] or DORCSD for details.
178 *
179 * P1, P2, Q1, and Q2 are represented as products of elementary
180 * reflectors. See DORCSD for details on generating P1, P2, Q1, and Q2
181 * using DORGQR and DORGLQ.
182 *
183 * Reference
184 * =========
185 *
186 * [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
187 * Algorithms, 50(1):33-65, 2009.
188 *
189 * ====================================================================
190 *
191 * .. Parameters ..
192 DOUBLE PRECISION REALONE
193 PARAMETER ( REALONE = 1.0D0 )
194 DOUBLE PRECISION NEGONE, ONE
195 PARAMETER ( NEGONE = -1.0D0, ONE = 1.0D0 )
196 * ..
197 * .. Local Scalars ..
198 LOGICAL COLMAJOR, LQUERY
199 INTEGER I, LWORKMIN, LWORKOPT
200 DOUBLE PRECISION Z1, Z2, Z3, Z4
201 * ..
202 * .. External Subroutines ..
203 EXTERNAL DAXPY, DLARF, DLARFGP, DSCAL, XERBLA
204 * ..
205 * .. External Functions ..
206 DOUBLE PRECISION DNRM2
207 LOGICAL LSAME
208 EXTERNAL DNRM2, LSAME
209 * ..
210 * .. Intrinsic Functions
211 INTRINSIC ATAN2, COS, MAX, MIN, SIN
212 * ..
213 * .. Executable Statements ..
214 *
215 * Test input arguments
216 *
217 INFO = 0
218 COLMAJOR = .NOT. LSAME( TRANS, 'T' )
219 IF( .NOT. LSAME( SIGNS, 'O' ) ) THEN
220 Z1 = REALONE
221 Z2 = REALONE
222 Z3 = REALONE
223 Z4 = REALONE
224 ELSE
225 Z1 = REALONE
226 Z2 = -REALONE
227 Z3 = REALONE
228 Z4 = -REALONE
229 END IF
230 LQUERY = LWORK .EQ. -1
231 *
232 IF( M .LT. 0 ) THEN
233 INFO = -3
234 ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN
235 INFO = -4
236 ELSE IF( Q .LT. 0 .OR. Q .GT. P .OR. Q .GT. M-P .OR.
237 $ Q .GT. M-Q ) THEN
238 INFO = -5
239 ELSE IF( COLMAJOR .AND. LDX11 .LT. MAX( 1, P ) ) THEN
240 INFO = -7
241 ELSE IF( .NOT.COLMAJOR .AND. LDX11 .LT. MAX( 1, Q ) ) THEN
242 INFO = -7
243 ELSE IF( COLMAJOR .AND. LDX12 .LT. MAX( 1, P ) ) THEN
244 INFO = -9
245 ELSE IF( .NOT.COLMAJOR .AND. LDX12 .LT. MAX( 1, M-Q ) ) THEN
246 INFO = -9
247 ELSE IF( COLMAJOR .AND. LDX21 .LT. MAX( 1, M-P ) ) THEN
248 INFO = -11
249 ELSE IF( .NOT.COLMAJOR .AND. LDX21 .LT. MAX( 1, Q ) ) THEN
250 INFO = -11
251 ELSE IF( COLMAJOR .AND. LDX22 .LT. MAX( 1, M-P ) ) THEN
252 INFO = -13
253 ELSE IF( .NOT.COLMAJOR .AND. LDX22 .LT. MAX( 1, M-Q ) ) THEN
254 INFO = -13
255 END IF
256 *
257 * Compute workspace
258 *
259 IF( INFO .EQ. 0 ) THEN
260 LWORKOPT = M - Q
261 LWORKMIN = M - Q
262 WORK(1) = LWORKOPT
263 IF( LWORK .LT. LWORKMIN .AND. .NOT. LQUERY ) THEN
264 INFO = -21
265 END IF
266 END IF
267 IF( INFO .NE. 0 ) THEN
268 CALL XERBLA( 'xORBDB', -INFO )
269 RETURN
270 ELSE IF( LQUERY ) THEN
271 RETURN
272 END IF
273 *
274 * Handle column-major and row-major separately
275 *
276 IF( COLMAJOR ) THEN
277 *
278 * Reduce columns 1, ..., Q of X11, X12, X21, and X22
279 *
280 DO I = 1, Q
281 *
282 IF( I .EQ. 1 ) THEN
283 CALL DSCAL( P-I+1, Z1, X11(I,I), 1 )
284 ELSE
285 CALL DSCAL( P-I+1, Z1*COS(PHI(I-1)), X11(I,I), 1 )
286 CALL DAXPY( P-I+1, -Z1*Z3*Z4*SIN(PHI(I-1)), X12(I,I-1),
287 $ 1, X11(I,I), 1 )
288 END IF
289 IF( I .EQ. 1 ) THEN
290 CALL DSCAL( M-P-I+1, Z2, X21(I,I), 1 )
291 ELSE
292 CALL DSCAL( M-P-I+1, Z2*COS(PHI(I-1)), X21(I,I), 1 )
293 CALL DAXPY( M-P-I+1, -Z2*Z3*Z4*SIN(PHI(I-1)), X22(I,I-1),
294 $ 1, X21(I,I), 1 )
295 END IF
296 *
297 THETA(I) = ATAN2( DNRM2( M-P-I+1, X21(I,I), 1 ),
298 $ DNRM2( P-I+1, X11(I,I), 1 ) )
299 *
300 CALL DLARFGP( P-I+1, X11(I,I), X11(I+1,I), 1, TAUP1(I) )
301 X11(I,I) = ONE
302 CALL DLARFGP( M-P-I+1, X21(I,I), X21(I+1,I), 1, TAUP2(I) )
303 X21(I,I) = ONE
304 *
305 CALL DLARF( 'L', P-I+1, Q-I, X11(I,I), 1, TAUP1(I),
306 $ X11(I,I+1), LDX11, WORK )
307 CALL DLARF( 'L', P-I+1, M-Q-I+1, X11(I,I), 1, TAUP1(I),
308 $ X12(I,I), LDX12, WORK )
309 CALL DLARF( 'L', M-P-I+1, Q-I, X21(I,I), 1, TAUP2(I),
310 $ X21(I,I+1), LDX21, WORK )
311 CALL DLARF( 'L', M-P-I+1, M-Q-I+1, X21(I,I), 1, TAUP2(I),
312 $ X22(I,I), LDX22, WORK )
313 *
314 IF( I .LT. Q ) THEN
315 CALL DSCAL( Q-I, -Z1*Z3*SIN(THETA(I)), X11(I,I+1),
316 $ LDX11 )
317 CALL DAXPY( Q-I, Z2*Z3*COS(THETA(I)), X21(I,I+1), LDX21,
318 $ X11(I,I+1), LDX11 )
319 END IF
320 CALL DSCAL( M-Q-I+1, -Z1*Z4*SIN(THETA(I)), X12(I,I), LDX12 )
321 CALL DAXPY( M-Q-I+1, Z2*Z4*COS(THETA(I)), X22(I,I), LDX22,
322 $ X12(I,I), LDX12 )
323 *
324 IF( I .LT. Q )
325 $ PHI(I) = ATAN2( DNRM2( Q-I, X11(I,I+1), LDX11 ),
326 $ DNRM2( M-Q-I+1, X12(I,I), LDX12 ) )
327 *
328 IF( I .LT. Q ) THEN
329 CALL DLARFGP( Q-I, X11(I,I+1), X11(I,I+2), LDX11,
330 $ TAUQ1(I) )
331 X11(I,I+1) = ONE
332 END IF
333 CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12,
334 $ TAUQ2(I) )
335 X12(I,I) = ONE
336 *
337 IF( I .LT. Q ) THEN
338 CALL DLARF( 'R', P-I, Q-I, X11(I,I+1), LDX11, TAUQ1(I),
339 $ X11(I+1,I+1), LDX11, WORK )
340 CALL DLARF( 'R', M-P-I, Q-I, X11(I,I+1), LDX11, TAUQ1(I),
341 $ X21(I+1,I+1), LDX21, WORK )
342 END IF
343 CALL DLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
344 $ X12(I+1,I), LDX12, WORK )
345 CALL DLARF( 'R', M-P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
346 $ X22(I+1,I), LDX22, WORK )
347 *
348 END DO
349 *
350 * Reduce columns Q + 1, ..., P of X12, X22
351 *
352 DO I = Q + 1, P
353 *
354 CALL DSCAL( M-Q-I+1, -Z1*Z4, X12(I,I), LDX12 )
355 CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12,
356 $ TAUQ2(I) )
357 X12(I,I) = ONE
358 *
359 CALL DLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
360 $ X12(I+1,I), LDX12, WORK )
361 IF( M-P-Q .GE. 1 )
362 $ CALL DLARF( 'R', M-P-Q, M-Q-I+1, X12(I,I), LDX12,
363 $ TAUQ2(I), X22(Q+1,I), LDX22, WORK )
364 *
365 END DO
366 *
367 * Reduce columns P + 1, ..., M - Q of X12, X22
368 *
369 DO I = 1, M - P - Q
370 *
371 CALL DSCAL( M-P-Q-I+1, Z2*Z4, X22(Q+I,P+I), LDX22 )
372 CALL DLARFGP( M-P-Q-I+1, X22(Q+I,P+I), X22(Q+I,P+I+1),
373 $ LDX22, TAUQ2(P+I) )
374 X22(Q+I,P+I) = ONE
375 CALL DLARF( 'R', M-P-Q-I, M-P-Q-I+1, X22(Q+I,P+I), LDX22,
376 $ TAUQ2(P+I), X22(Q+I+1,P+I), LDX22, WORK )
377 *
378 END DO
379 *
380 ELSE
381 *
382 * Reduce columns 1, ..., Q of X11, X12, X21, X22
383 *
384 DO I = 1, Q
385 *
386 IF( I .EQ. 1 ) THEN
387 CALL DSCAL( P-I+1, Z1, X11(I,I), LDX11 )
388 ELSE
389 CALL DSCAL( P-I+1, Z1*COS(PHI(I-1)), X11(I,I), LDX11 )
390 CALL DAXPY( P-I+1, -Z1*Z3*Z4*SIN(PHI(I-1)), X12(I-1,I),
391 $ LDX12, X11(I,I), LDX11 )
392 END IF
393 IF( I .EQ. 1 ) THEN
394 CALL DSCAL( M-P-I+1, Z2, X21(I,I), LDX21 )
395 ELSE
396 CALL DSCAL( M-P-I+1, Z2*COS(PHI(I-1)), X21(I,I), LDX21 )
397 CALL DAXPY( M-P-I+1, -Z2*Z3*Z4*SIN(PHI(I-1)), X22(I-1,I),
398 $ LDX22, X21(I,I), LDX21 )
399 END IF
400 *
401 THETA(I) = ATAN2( DNRM2( M-P-I+1, X21(I,I), LDX21 ),
402 $ DNRM2( P-I+1, X11(I,I), LDX11 ) )
403 *
404 CALL DLARFGP( P-I+1, X11(I,I), X11(I,I+1), LDX11, TAUP1(I) )
405 X11(I,I) = ONE
406 CALL DLARFGP( M-P-I+1, X21(I,I), X21(I,I+1), LDX21,
407 $ TAUP2(I) )
408 X21(I,I) = ONE
409 *
410 CALL DLARF( 'R', Q-I, P-I+1, X11(I,I), LDX11, TAUP1(I),
411 $ X11(I+1,I), LDX11, WORK )
412 CALL DLARF( 'R', M-Q-I+1, P-I+1, X11(I,I), LDX11, TAUP1(I),
413 $ X12(I,I), LDX12, WORK )
414 CALL DLARF( 'R', Q-I, M-P-I+1, X21(I,I), LDX21, TAUP2(I),
415 $ X21(I+1,I), LDX21, WORK )
416 CALL DLARF( 'R', M-Q-I+1, M-P-I+1, X21(I,I), LDX21,
417 $ TAUP2(I), X22(I,I), LDX22, WORK )
418 *
419 IF( I .LT. Q ) THEN
420 CALL DSCAL( Q-I, -Z1*Z3*SIN(THETA(I)), X11(I+1,I), 1 )
421 CALL DAXPY( Q-I, Z2*Z3*COS(THETA(I)), X21(I+1,I), 1,
422 $ X11(I+1,I), 1 )
423 END IF
424 CALL DSCAL( M-Q-I+1, -Z1*Z4*SIN(THETA(I)), X12(I,I), 1 )
425 CALL DAXPY( M-Q-I+1, Z2*Z4*COS(THETA(I)), X22(I,I), 1,
426 $ X12(I,I), 1 )
427 *
428 IF( I .LT. Q )
429 $ PHI(I) = ATAN2( DNRM2( Q-I, X11(I+1,I), 1 ),
430 $ DNRM2( M-Q-I+1, X12(I,I), 1 ) )
431 *
432 IF( I .LT. Q ) THEN
433 CALL DLARFGP( Q-I, X11(I+1,I), X11(I+2,I), 1, TAUQ1(I) )
434 X11(I+1,I) = ONE
435 END IF
436 CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1, TAUQ2(I) )
437 X12(I,I) = ONE
438 *
439 IF( I .LT. Q ) THEN
440 CALL DLARF( 'L', Q-I, P-I, X11(I+1,I), 1, TAUQ1(I),
441 $ X11(I+1,I+1), LDX11, WORK )
442 CALL DLARF( 'L', Q-I, M-P-I, X11(I+1,I), 1, TAUQ1(I),
443 $ X21(I+1,I+1), LDX21, WORK )
444 END IF
445 CALL DLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1, TAUQ2(I),
446 $ X12(I,I+1), LDX12, WORK )
447 CALL DLARF( 'L', M-Q-I+1, M-P-I, X12(I,I), 1, TAUQ2(I),
448 $ X22(I,I+1), LDX22, WORK )
449 *
450 END DO
451 *
452 * Reduce columns Q + 1, ..., P of X12, X22
453 *
454 DO I = Q + 1, P
455 *
456 CALL DSCAL( M-Q-I+1, -Z1*Z4, X12(I,I), 1 )
457 CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1, TAUQ2(I) )
458 X12(I,I) = ONE
459 *
460 CALL DLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1, TAUQ2(I),
461 $ X12(I,I+1), LDX12, WORK )
462 IF( M-P-Q .GE. 1 )
463 $ CALL DLARF( 'L', M-Q-I+1, M-P-Q, X12(I,I), 1, TAUQ2(I),
464 $ X22(I,Q+1), LDX22, WORK )
465 *
466 END DO
467 *
468 * Reduce columns P + 1, ..., M - Q of X12, X22
469 *
470 DO I = 1, M - P - Q
471 *
472 CALL DSCAL( M-P-Q-I+1, Z2*Z4, X22(P+I,Q+I), 1 )
473 CALL DLARFGP( M-P-Q-I+1, X22(P+I,Q+I), X22(P+I+1,Q+I), 1,
474 $ TAUQ2(P+I) )
475 X22(P+I,Q+I) = ONE
476 *
477 CALL DLARF( 'L', M-P-Q-I+1, M-P-Q-I, X22(P+I,Q+I), 1,
478 $ TAUQ2(P+I), X22(P+I,Q+I+1), LDX22, WORK )
479 *
480 END DO
481 *
482 END IF
483 *
484 RETURN
485 *
486 * End of DORBDB
487 *
488 END
489
2 $ X21, LDX21, X22, LDX22, THETA, PHI, TAUP1,
3 $ TAUP2, TAUQ1, TAUQ2, WORK, LWORK, INFO )
4 IMPLICIT NONE
5 *
6 * -- LAPACK routine (version 3.3.0) --
7 *
8 * -- Contributed by Brian Sutton of the Randolph-Macon College --
9 * -- November 2010
10 *
11 * -- LAPACK is a software package provided by Univ. of Tennessee, --
12 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
13 *
14 * .. Scalar Arguments ..
15 CHARACTER SIGNS, TRANS
16 INTEGER INFO, LDX11, LDX12, LDX21, LDX22, LWORK, M, P,
17 $ Q
18 * ..
19 * .. Array Arguments ..
20 DOUBLE PRECISION PHI( * ), THETA( * )
21 DOUBLE PRECISION TAUP1( * ), TAUP2( * ), TAUQ1( * ), TAUQ2( * ),
22 $ WORK( * ), X11( LDX11, * ), X12( LDX12, * ),
23 $ X21( LDX21, * ), X22( LDX22, * )
24 * ..
25 *
26 * Purpose
27 * =======
28 *
29 * DORBDB simultaneously bidiagonalizes the blocks of an M-by-M
30 * partitioned orthogonal matrix X:
31 *
32 * [ B11 | B12 0 0 ]
33 * [ X11 | X12 ] [ P1 | ] [ 0 | 0 -I 0 ] [ Q1 | ]**T
34 * X = [-----------] = [---------] [----------------] [---------] .
35 * [ X21 | X22 ] [ | P2 ] [ B21 | B22 0 0 ] [ | Q2 ]
36 * [ 0 | 0 0 I ]
37 *
38 * X11 is P-by-Q. Q must be no larger than P, M-P, or M-Q. (If this is
39 * not the case, then X must be transposed and/or permuted. This can be
40 * done in constant time using the TRANS and SIGNS options. See DORCSD
41 * for details.)
42 *
43 * The orthogonal matrices P1, P2, Q1, and Q2 are P-by-P, (M-P)-by-
44 * (M-P), Q-by-Q, and (M-Q)-by-(M-Q), respectively. They are
45 * represented implicitly by Householder vectors.
46 *
47 * B11, B12, B21, and B22 are Q-by-Q bidiagonal matrices represented
48 * implicitly by angles THETA, PHI.
49 *
50 * Arguments
51 * =========
52 *
53 * TRANS (input) CHARACTER
54 * = 'T': X, U1, U2, V1T, and V2T are stored in row-major
55 * order;
56 * otherwise: X, U1, U2, V1T, and V2T are stored in column-
57 * major order.
58 *
59 * SIGNS (input) CHARACTER
60 * = 'O': The lower-left block is made nonpositive (the
61 * "other" convention);
62 * otherwise: The upper-right block is made nonpositive (the
63 * "default" convention).
64 *
65 * M (input) INTEGER
66 * The number of rows and columns in X.
67 *
68 * P (input) INTEGER
69 * The number of rows in X11 and X12. 0 <= P <= M.
70 *
71 * Q (input) INTEGER
72 * The number of columns in X11 and X21. 0 <= Q <=
73 * MIN(P,M-P,M-Q).
74 *
75 * X11 (input/output) DOUBLE PRECISION array, dimension (LDX11,Q)
76 * On entry, the top-left block of the orthogonal matrix to be
77 * reduced. On exit, the form depends on TRANS:
78 * If TRANS = 'N', then
79 * the columns of tril(X11) specify reflectors for P1,
80 * the rows of triu(X11,1) specify reflectors for Q1;
81 * else TRANS = 'T', and
82 * the rows of triu(X11) specify reflectors for P1,
83 * the columns of tril(X11,-1) specify reflectors for Q1.
84 *
85 * LDX11 (input) INTEGER
86 * The leading dimension of X11. If TRANS = 'N', then LDX11 >=
87 * P; else LDX11 >= Q.
88 *
89 * X12 (input/output) DOUBLE PRECISION array, dimension (LDX12,M-Q)
90 * On entry, the top-right block of the orthogonal matrix to
91 * be reduced. On exit, the form depends on TRANS:
92 * If TRANS = 'N', then
93 * the rows of triu(X12) specify the first P reflectors for
94 * Q2;
95 * else TRANS = 'T', and
96 * the columns of tril(X12) specify the first P reflectors
97 * for Q2.
98 *
99 * LDX12 (input) INTEGER
100 * The leading dimension of X12. If TRANS = 'N', then LDX12 >=
101 * P; else LDX11 >= M-Q.
102 *
103 * X21 (input/output) DOUBLE PRECISION array, dimension (LDX21,Q)
104 * On entry, the bottom-left block of the orthogonal matrix to
105 * be reduced. On exit, the form depends on TRANS:
106 * If TRANS = 'N', then
107 * the columns of tril(X21) specify reflectors for P2;
108 * else TRANS = 'T', and
109 * the rows of triu(X21) specify reflectors for P2.
110 *
111 * LDX21 (input) INTEGER
112 * The leading dimension of X21. If TRANS = 'N', then LDX21 >=
113 * M-P; else LDX21 >= Q.
114 *
115 * X22 (input/output) DOUBLE PRECISION array, dimension (LDX22,M-Q)
116 * On entry, the bottom-right block of the orthogonal matrix to
117 * be reduced. On exit, the form depends on TRANS:
118 * If TRANS = 'N', then
119 * the rows of triu(X22(Q+1:M-P,P+1:M-Q)) specify the last
120 * M-P-Q reflectors for Q2,
121 * else TRANS = 'T', and
122 * the columns of tril(X22(P+1:M-Q,Q+1:M-P)) specify the last
123 * M-P-Q reflectors for P2.
124 *
125 * LDX22 (input) INTEGER
126 * The leading dimension of X22. If TRANS = 'N', then LDX22 >=
127 * M-P; else LDX22 >= M-Q.
128 *
129 * THETA (output) DOUBLE PRECISION array, dimension (Q)
130 * The entries of the bidiagonal blocks B11, B12, B21, B22 can
131 * be computed from the angles THETA and PHI. See Further
132 * Details.
133 *
134 * PHI (output) DOUBLE PRECISION array, dimension (Q-1)
135 * The entries of the bidiagonal blocks B11, B12, B21, B22 can
136 * be computed from the angles THETA and PHI. See Further
137 * Details.
138 *
139 * TAUP1 (output) DOUBLE PRECISION array, dimension (P)
140 * The scalar factors of the elementary reflectors that define
141 * P1.
142 *
143 * TAUP2 (output) DOUBLE PRECISION array, dimension (M-P)
144 * The scalar factors of the elementary reflectors that define
145 * P2.
146 *
147 * TAUQ1 (output) DOUBLE PRECISION array, dimension (Q)
148 * The scalar factors of the elementary reflectors that define
149 * Q1.
150 *
151 * TAUQ2 (output) DOUBLE PRECISION array, dimension (M-Q)
152 * The scalar factors of the elementary reflectors that define
153 * Q2.
154 *
155 * WORK (workspace) DOUBLE PRECISION array, dimension (LWORK)
156 *
157 * LWORK (input) INTEGER
158 * The dimension of the array WORK. LWORK >= M-Q.
159 *
160 * If LWORK = -1, then a workspace query is assumed; the routine
161 * only calculates the optimal size of the WORK array, returns
162 * this value as the first entry of the WORK array, and no error
163 * message related to LWORK is issued by XERBLA.
164 *
165 * INFO (output) INTEGER
166 * = 0: successful exit.
167 * < 0: if INFO = -i, the i-th argument had an illegal value.
168 *
169 * Further Details
170 * ===============
171 *
172 * The bidiagonal blocks B11, B12, B21, and B22 are represented
173 * implicitly by angles THETA(1), ..., THETA(Q) and PHI(1), ...,
174 * PHI(Q-1). B11 and B21 are upper bidiagonal, while B21 and B22 are
175 * lower bidiagonal. Every entry in each bidiagonal band is a product
176 * of a sine or cosine of a THETA with a sine or cosine of a PHI. See
177 * [1] or DORCSD for details.
178 *
179 * P1, P2, Q1, and Q2 are represented as products of elementary
180 * reflectors. See DORCSD for details on generating P1, P2, Q1, and Q2
181 * using DORGQR and DORGLQ.
182 *
183 * Reference
184 * =========
185 *
186 * [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
187 * Algorithms, 50(1):33-65, 2009.
188 *
189 * ====================================================================
190 *
191 * .. Parameters ..
192 DOUBLE PRECISION REALONE
193 PARAMETER ( REALONE = 1.0D0 )
194 DOUBLE PRECISION NEGONE, ONE
195 PARAMETER ( NEGONE = -1.0D0, ONE = 1.0D0 )
196 * ..
197 * .. Local Scalars ..
198 LOGICAL COLMAJOR, LQUERY
199 INTEGER I, LWORKMIN, LWORKOPT
200 DOUBLE PRECISION Z1, Z2, Z3, Z4
201 * ..
202 * .. External Subroutines ..
203 EXTERNAL DAXPY, DLARF, DLARFGP, DSCAL, XERBLA
204 * ..
205 * .. External Functions ..
206 DOUBLE PRECISION DNRM2
207 LOGICAL LSAME
208 EXTERNAL DNRM2, LSAME
209 * ..
210 * .. Intrinsic Functions
211 INTRINSIC ATAN2, COS, MAX, MIN, SIN
212 * ..
213 * .. Executable Statements ..
214 *
215 * Test input arguments
216 *
217 INFO = 0
218 COLMAJOR = .NOT. LSAME( TRANS, 'T' )
219 IF( .NOT. LSAME( SIGNS, 'O' ) ) THEN
220 Z1 = REALONE
221 Z2 = REALONE
222 Z3 = REALONE
223 Z4 = REALONE
224 ELSE
225 Z1 = REALONE
226 Z2 = -REALONE
227 Z3 = REALONE
228 Z4 = -REALONE
229 END IF
230 LQUERY = LWORK .EQ. -1
231 *
232 IF( M .LT. 0 ) THEN
233 INFO = -3
234 ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN
235 INFO = -4
236 ELSE IF( Q .LT. 0 .OR. Q .GT. P .OR. Q .GT. M-P .OR.
237 $ Q .GT. M-Q ) THEN
238 INFO = -5
239 ELSE IF( COLMAJOR .AND. LDX11 .LT. MAX( 1, P ) ) THEN
240 INFO = -7
241 ELSE IF( .NOT.COLMAJOR .AND. LDX11 .LT. MAX( 1, Q ) ) THEN
242 INFO = -7
243 ELSE IF( COLMAJOR .AND. LDX12 .LT. MAX( 1, P ) ) THEN
244 INFO = -9
245 ELSE IF( .NOT.COLMAJOR .AND. LDX12 .LT. MAX( 1, M-Q ) ) THEN
246 INFO = -9
247 ELSE IF( COLMAJOR .AND. LDX21 .LT. MAX( 1, M-P ) ) THEN
248 INFO = -11
249 ELSE IF( .NOT.COLMAJOR .AND. LDX21 .LT. MAX( 1, Q ) ) THEN
250 INFO = -11
251 ELSE IF( COLMAJOR .AND. LDX22 .LT. MAX( 1, M-P ) ) THEN
252 INFO = -13
253 ELSE IF( .NOT.COLMAJOR .AND. LDX22 .LT. MAX( 1, M-Q ) ) THEN
254 INFO = -13
255 END IF
256 *
257 * Compute workspace
258 *
259 IF( INFO .EQ. 0 ) THEN
260 LWORKOPT = M - Q
261 LWORKMIN = M - Q
262 WORK(1) = LWORKOPT
263 IF( LWORK .LT. LWORKMIN .AND. .NOT. LQUERY ) THEN
264 INFO = -21
265 END IF
266 END IF
267 IF( INFO .NE. 0 ) THEN
268 CALL XERBLA( 'xORBDB', -INFO )
269 RETURN
270 ELSE IF( LQUERY ) THEN
271 RETURN
272 END IF
273 *
274 * Handle column-major and row-major separately
275 *
276 IF( COLMAJOR ) THEN
277 *
278 * Reduce columns 1, ..., Q of X11, X12, X21, and X22
279 *
280 DO I = 1, Q
281 *
282 IF( I .EQ. 1 ) THEN
283 CALL DSCAL( P-I+1, Z1, X11(I,I), 1 )
284 ELSE
285 CALL DSCAL( P-I+1, Z1*COS(PHI(I-1)), X11(I,I), 1 )
286 CALL DAXPY( P-I+1, -Z1*Z3*Z4*SIN(PHI(I-1)), X12(I,I-1),
287 $ 1, X11(I,I), 1 )
288 END IF
289 IF( I .EQ. 1 ) THEN
290 CALL DSCAL( M-P-I+1, Z2, X21(I,I), 1 )
291 ELSE
292 CALL DSCAL( M-P-I+1, Z2*COS(PHI(I-1)), X21(I,I), 1 )
293 CALL DAXPY( M-P-I+1, -Z2*Z3*Z4*SIN(PHI(I-1)), X22(I,I-1),
294 $ 1, X21(I,I), 1 )
295 END IF
296 *
297 THETA(I) = ATAN2( DNRM2( M-P-I+1, X21(I,I), 1 ),
298 $ DNRM2( P-I+1, X11(I,I), 1 ) )
299 *
300 CALL DLARFGP( P-I+1, X11(I,I), X11(I+1,I), 1, TAUP1(I) )
301 X11(I,I) = ONE
302 CALL DLARFGP( M-P-I+1, X21(I,I), X21(I+1,I), 1, TAUP2(I) )
303 X21(I,I) = ONE
304 *
305 CALL DLARF( 'L', P-I+1, Q-I, X11(I,I), 1, TAUP1(I),
306 $ X11(I,I+1), LDX11, WORK )
307 CALL DLARF( 'L', P-I+1, M-Q-I+1, X11(I,I), 1, TAUP1(I),
308 $ X12(I,I), LDX12, WORK )
309 CALL DLARF( 'L', M-P-I+1, Q-I, X21(I,I), 1, TAUP2(I),
310 $ X21(I,I+1), LDX21, WORK )
311 CALL DLARF( 'L', M-P-I+1, M-Q-I+1, X21(I,I), 1, TAUP2(I),
312 $ X22(I,I), LDX22, WORK )
313 *
314 IF( I .LT. Q ) THEN
315 CALL DSCAL( Q-I, -Z1*Z3*SIN(THETA(I)), X11(I,I+1),
316 $ LDX11 )
317 CALL DAXPY( Q-I, Z2*Z3*COS(THETA(I)), X21(I,I+1), LDX21,
318 $ X11(I,I+1), LDX11 )
319 END IF
320 CALL DSCAL( M-Q-I+1, -Z1*Z4*SIN(THETA(I)), X12(I,I), LDX12 )
321 CALL DAXPY( M-Q-I+1, Z2*Z4*COS(THETA(I)), X22(I,I), LDX22,
322 $ X12(I,I), LDX12 )
323 *
324 IF( I .LT. Q )
325 $ PHI(I) = ATAN2( DNRM2( Q-I, X11(I,I+1), LDX11 ),
326 $ DNRM2( M-Q-I+1, X12(I,I), LDX12 ) )
327 *
328 IF( I .LT. Q ) THEN
329 CALL DLARFGP( Q-I, X11(I,I+1), X11(I,I+2), LDX11,
330 $ TAUQ1(I) )
331 X11(I,I+1) = ONE
332 END IF
333 CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12,
334 $ TAUQ2(I) )
335 X12(I,I) = ONE
336 *
337 IF( I .LT. Q ) THEN
338 CALL DLARF( 'R', P-I, Q-I, X11(I,I+1), LDX11, TAUQ1(I),
339 $ X11(I+1,I+1), LDX11, WORK )
340 CALL DLARF( 'R', M-P-I, Q-I, X11(I,I+1), LDX11, TAUQ1(I),
341 $ X21(I+1,I+1), LDX21, WORK )
342 END IF
343 CALL DLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
344 $ X12(I+1,I), LDX12, WORK )
345 CALL DLARF( 'R', M-P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
346 $ X22(I+1,I), LDX22, WORK )
347 *
348 END DO
349 *
350 * Reduce columns Q + 1, ..., P of X12, X22
351 *
352 DO I = Q + 1, P
353 *
354 CALL DSCAL( M-Q-I+1, -Z1*Z4, X12(I,I), LDX12 )
355 CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12,
356 $ TAUQ2(I) )
357 X12(I,I) = ONE
358 *
359 CALL DLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
360 $ X12(I+1,I), LDX12, WORK )
361 IF( M-P-Q .GE. 1 )
362 $ CALL DLARF( 'R', M-P-Q, M-Q-I+1, X12(I,I), LDX12,
363 $ TAUQ2(I), X22(Q+1,I), LDX22, WORK )
364 *
365 END DO
366 *
367 * Reduce columns P + 1, ..., M - Q of X12, X22
368 *
369 DO I = 1, M - P - Q
370 *
371 CALL DSCAL( M-P-Q-I+1, Z2*Z4, X22(Q+I,P+I), LDX22 )
372 CALL DLARFGP( M-P-Q-I+1, X22(Q+I,P+I), X22(Q+I,P+I+1),
373 $ LDX22, TAUQ2(P+I) )
374 X22(Q+I,P+I) = ONE
375 CALL DLARF( 'R', M-P-Q-I, M-P-Q-I+1, X22(Q+I,P+I), LDX22,
376 $ TAUQ2(P+I), X22(Q+I+1,P+I), LDX22, WORK )
377 *
378 END DO
379 *
380 ELSE
381 *
382 * Reduce columns 1, ..., Q of X11, X12, X21, X22
383 *
384 DO I = 1, Q
385 *
386 IF( I .EQ. 1 ) THEN
387 CALL DSCAL( P-I+1, Z1, X11(I,I), LDX11 )
388 ELSE
389 CALL DSCAL( P-I+1, Z1*COS(PHI(I-1)), X11(I,I), LDX11 )
390 CALL DAXPY( P-I+1, -Z1*Z3*Z4*SIN(PHI(I-1)), X12(I-1,I),
391 $ LDX12, X11(I,I), LDX11 )
392 END IF
393 IF( I .EQ. 1 ) THEN
394 CALL DSCAL( M-P-I+1, Z2, X21(I,I), LDX21 )
395 ELSE
396 CALL DSCAL( M-P-I+1, Z2*COS(PHI(I-1)), X21(I,I), LDX21 )
397 CALL DAXPY( M-P-I+1, -Z2*Z3*Z4*SIN(PHI(I-1)), X22(I-1,I),
398 $ LDX22, X21(I,I), LDX21 )
399 END IF
400 *
401 THETA(I) = ATAN2( DNRM2( M-P-I+1, X21(I,I), LDX21 ),
402 $ DNRM2( P-I+1, X11(I,I), LDX11 ) )
403 *
404 CALL DLARFGP( P-I+1, X11(I,I), X11(I,I+1), LDX11, TAUP1(I) )
405 X11(I,I) = ONE
406 CALL DLARFGP( M-P-I+1, X21(I,I), X21(I,I+1), LDX21,
407 $ TAUP2(I) )
408 X21(I,I) = ONE
409 *
410 CALL DLARF( 'R', Q-I, P-I+1, X11(I,I), LDX11, TAUP1(I),
411 $ X11(I+1,I), LDX11, WORK )
412 CALL DLARF( 'R', M-Q-I+1, P-I+1, X11(I,I), LDX11, TAUP1(I),
413 $ X12(I,I), LDX12, WORK )
414 CALL DLARF( 'R', Q-I, M-P-I+1, X21(I,I), LDX21, TAUP2(I),
415 $ X21(I+1,I), LDX21, WORK )
416 CALL DLARF( 'R', M-Q-I+1, M-P-I+1, X21(I,I), LDX21,
417 $ TAUP2(I), X22(I,I), LDX22, WORK )
418 *
419 IF( I .LT. Q ) THEN
420 CALL DSCAL( Q-I, -Z1*Z3*SIN(THETA(I)), X11(I+1,I), 1 )
421 CALL DAXPY( Q-I, Z2*Z3*COS(THETA(I)), X21(I+1,I), 1,
422 $ X11(I+1,I), 1 )
423 END IF
424 CALL DSCAL( M-Q-I+1, -Z1*Z4*SIN(THETA(I)), X12(I,I), 1 )
425 CALL DAXPY( M-Q-I+1, Z2*Z4*COS(THETA(I)), X22(I,I), 1,
426 $ X12(I,I), 1 )
427 *
428 IF( I .LT. Q )
429 $ PHI(I) = ATAN2( DNRM2( Q-I, X11(I+1,I), 1 ),
430 $ DNRM2( M-Q-I+1, X12(I,I), 1 ) )
431 *
432 IF( I .LT. Q ) THEN
433 CALL DLARFGP( Q-I, X11(I+1,I), X11(I+2,I), 1, TAUQ1(I) )
434 X11(I+1,I) = ONE
435 END IF
436 CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1, TAUQ2(I) )
437 X12(I,I) = ONE
438 *
439 IF( I .LT. Q ) THEN
440 CALL DLARF( 'L', Q-I, P-I, X11(I+1,I), 1, TAUQ1(I),
441 $ X11(I+1,I+1), LDX11, WORK )
442 CALL DLARF( 'L', Q-I, M-P-I, X11(I+1,I), 1, TAUQ1(I),
443 $ X21(I+1,I+1), LDX21, WORK )
444 END IF
445 CALL DLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1, TAUQ2(I),
446 $ X12(I,I+1), LDX12, WORK )
447 CALL DLARF( 'L', M-Q-I+1, M-P-I, X12(I,I), 1, TAUQ2(I),
448 $ X22(I,I+1), LDX22, WORK )
449 *
450 END DO
451 *
452 * Reduce columns Q + 1, ..., P of X12, X22
453 *
454 DO I = Q + 1, P
455 *
456 CALL DSCAL( M-Q-I+1, -Z1*Z4, X12(I,I), 1 )
457 CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1, TAUQ2(I) )
458 X12(I,I) = ONE
459 *
460 CALL DLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1, TAUQ2(I),
461 $ X12(I,I+1), LDX12, WORK )
462 IF( M-P-Q .GE. 1 )
463 $ CALL DLARF( 'L', M-Q-I+1, M-P-Q, X12(I,I), 1, TAUQ2(I),
464 $ X22(I,Q+1), LDX22, WORK )
465 *
466 END DO
467 *
468 * Reduce columns P + 1, ..., M - Q of X12, X22
469 *
470 DO I = 1, M - P - Q
471 *
472 CALL DSCAL( M-P-Q-I+1, Z2*Z4, X22(P+I,Q+I), 1 )
473 CALL DLARFGP( M-P-Q-I+1, X22(P+I,Q+I), X22(P+I+1,Q+I), 1,
474 $ TAUQ2(P+I) )
475 X22(P+I,Q+I) = ONE
476 *
477 CALL DLARF( 'L', M-P-Q-I+1, M-P-Q-I, X22(P+I,Q+I), 1,
478 $ TAUQ2(P+I), X22(P+I,Q+I+1), LDX22, WORK )
479 *
480 END DO
481 *
482 END IF
483 *
484 RETURN
485 *
486 * End of DORBDB
487 *
488 END
489