LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
sorbdb.f
Go to the documentation of this file.
1*> \brief \b SORBDB
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download SORBDB + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sorbdb.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sorbdb.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sorbdb.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12,
22* X21, LDX21, X22, LDX22, THETA, PHI, TAUP1,
23* TAUP2, TAUQ1, TAUQ2, WORK, LWORK, INFO )
24*
25* .. Scalar Arguments ..
26* CHARACTER SIGNS, TRANS
27* INTEGER INFO, LDX11, LDX12, LDX21, LDX22, LWORK, M, P,
28* $ Q
29* ..
30* .. Array Arguments ..
31* REAL PHI( * ), THETA( * )
32* REAL TAUP1( * ), TAUP2( * ), TAUQ1( * ), TAUQ2( * ),
33* $ WORK( * ), X11( LDX11, * ), X12( LDX12, * ),
34* $ X21( LDX21, * ), X22( LDX22, * )
35* ..
36*
37*
38*> \par Purpose:
39* =============
40*>
41*> \verbatim
42*>
43*> SORBDB simultaneously bidiagonalizes the blocks of an M-by-M
44*> partitioned orthogonal matrix X:
45*>
46*> [ B11 | B12 0 0 ]
47*> [ X11 | X12 ] [ P1 | ] [ 0 | 0 -I 0 ] [ Q1 | ]**T
48*> X = [-----------] = [---------] [----------------] [---------] .
49*> [ X21 | X22 ] [ | P2 ] [ B21 | B22 0 0 ] [ | Q2 ]
50*> [ 0 | 0 0 I ]
51*>
52*> X11 is P-by-Q. Q must be no larger than P, M-P, or M-Q. (If this is
53*> not the case, then X must be transposed and/or permuted. This can be
54*> done in constant time using the TRANS and SIGNS options. See SORCSD
55*> for details.)
56*>
57*> The orthogonal matrices P1, P2, Q1, and Q2 are P-by-P, (M-P)-by-
58*> (M-P), Q-by-Q, and (M-Q)-by-(M-Q), respectively. They are
59*> represented implicitly by Householder vectors.
60*>
61*> B11, B12, B21, and B22 are Q-by-Q bidiagonal matrices represented
62*> implicitly by angles THETA, PHI.
63*> \endverbatim
64*
65* Arguments:
66* ==========
67*
68*> \param[in] TRANS
69*> \verbatim
70*> TRANS is CHARACTER
71*> = 'T': X, U1, U2, V1T, and V2T are stored in row-major
72*> order;
73*> otherwise: X, U1, U2, V1T, and V2T are stored in column-
74*> major order.
75*> \endverbatim
76*>
77*> \param[in] SIGNS
78*> \verbatim
79*> SIGNS is CHARACTER
80*> = 'O': The lower-left block is made nonpositive (the
81*> "other" convention);
82*> otherwise: The upper-right block is made nonpositive (the
83*> "default" convention).
84*> \endverbatim
85*>
86*> \param[in] M
87*> \verbatim
88*> M is INTEGER
89*> The number of rows and columns in X.
90*> \endverbatim
91*>
92*> \param[in] P
93*> \verbatim
94*> P is INTEGER
95*> The number of rows in X11 and X12. 0 <= P <= M.
96*> \endverbatim
97*>
98*> \param[in] Q
99*> \verbatim
100*> Q is INTEGER
101*> The number of columns in X11 and X21. 0 <= Q <=
102*> MIN(P,M-P,M-Q).
103*> \endverbatim
104*>
105*> \param[in,out] X11
106*> \verbatim
107*> X11 is REAL array, dimension (LDX11,Q)
108*> On entry, the top-left block of the orthogonal matrix to be
109*> reduced. On exit, the form depends on TRANS:
110*> If TRANS = 'N', then
111*> the columns of tril(X11) specify reflectors for P1,
112*> the rows of triu(X11,1) specify reflectors for Q1;
113*> else TRANS = 'T', and
114*> the rows of triu(X11) specify reflectors for P1,
115*> the columns of tril(X11,-1) specify reflectors for Q1.
116*> \endverbatim
117*>
118*> \param[in] LDX11
119*> \verbatim
120*> LDX11 is INTEGER
121*> The leading dimension of X11. If TRANS = 'N', then LDX11 >=
122*> P; else LDX11 >= Q.
123*> \endverbatim
124*>
125*> \param[in,out] X12
126*> \verbatim
127*> X12 is REAL array, dimension (LDX12,M-Q)
128*> On entry, the top-right block of the orthogonal matrix to
129*> be reduced. On exit, the form depends on TRANS:
130*> If TRANS = 'N', then
131*> the rows of triu(X12) specify the first P reflectors for
132*> Q2;
133*> else TRANS = 'T', and
134*> the columns of tril(X12) specify the first P reflectors
135*> for Q2.
136*> \endverbatim
137*>
138*> \param[in] LDX12
139*> \verbatim
140*> LDX12 is INTEGER
141*> The leading dimension of X12. If TRANS = 'N', then LDX12 >=
142*> P; else LDX11 >= M-Q.
143*> \endverbatim
144*>
145*> \param[in,out] X21
146*> \verbatim
147*> X21 is REAL array, dimension (LDX21,Q)
148*> On entry, the bottom-left block of the orthogonal matrix to
149*> be reduced. On exit, the form depends on TRANS:
150*> If TRANS = 'N', then
151*> the columns of tril(X21) specify reflectors for P2;
152*> else TRANS = 'T', and
153*> the rows of triu(X21) specify reflectors for P2.
154*> \endverbatim
155*>
156*> \param[in] LDX21
157*> \verbatim
158*> LDX21 is INTEGER
159*> The leading dimension of X21. If TRANS = 'N', then LDX21 >=
160*> M-P; else LDX21 >= Q.
161*> \endverbatim
162*>
163*> \param[in,out] X22
164*> \verbatim
165*> X22 is REAL array, dimension (LDX22,M-Q)
166*> On entry, the bottom-right block of the orthogonal matrix to
167*> be reduced. On exit, the form depends on TRANS:
168*> If TRANS = 'N', then
169*> the rows of triu(X22(Q+1:M-P,P+1:M-Q)) specify the last
170*> M-P-Q reflectors for Q2,
171*> else TRANS = 'T', and
172*> the columns of tril(X22(P+1:M-Q,Q+1:M-P)) specify the last
173*> M-P-Q reflectors for P2.
174*> \endverbatim
175*>
176*> \param[in] LDX22
177*> \verbatim
178*> LDX22 is INTEGER
179*> The leading dimension of X22. If TRANS = 'N', then LDX22 >=
180*> M-P; else LDX22 >= M-Q.
181*> \endverbatim
182*>
183*> \param[out] THETA
184*> \verbatim
185*> THETA is REAL array, dimension (Q)
186*> The entries of the bidiagonal blocks B11, B12, B21, B22 can
187*> be computed from the angles THETA and PHI. See Further
188*> Details.
189*> \endverbatim
190*>
191*> \param[out] PHI
192*> \verbatim
193*> PHI is REAL array, dimension (Q-1)
194*> The entries of the bidiagonal blocks B11, B12, B21, B22 can
195*> be computed from the angles THETA and PHI. See Further
196*> Details.
197*> \endverbatim
198*>
199*> \param[out] TAUP1
200*> \verbatim
201*> TAUP1 is REAL array, dimension (P)
202*> The scalar factors of the elementary reflectors that define
203*> P1.
204*> \endverbatim
205*>
206*> \param[out] TAUP2
207*> \verbatim
208*> TAUP2 is REAL array, dimension (M-P)
209*> The scalar factors of the elementary reflectors that define
210*> P2.
211*> \endverbatim
212*>
213*> \param[out] TAUQ1
214*> \verbatim
215*> TAUQ1 is REAL array, dimension (Q)
216*> The scalar factors of the elementary reflectors that define
217*> Q1.
218*> \endverbatim
219*>
220*> \param[out] TAUQ2
221*> \verbatim
222*> TAUQ2 is REAL array, dimension (M-Q)
223*> The scalar factors of the elementary reflectors that define
224*> Q2.
225*> \endverbatim
226*>
227*> \param[out] WORK
228*> \verbatim
229*> WORK is REAL array, dimension (LWORK)
230*> \endverbatim
231*>
232*> \param[in] LWORK
233*> \verbatim
234*> LWORK is INTEGER
235*> The dimension of the array WORK. LWORK >= M-Q.
236*>
237*> If LWORK = -1, then a workspace query is assumed; the routine
238*> only calculates the optimal size of the WORK array, returns
239*> this value as the first entry of the WORK array, and no error
240*> message related to LWORK is issued by XERBLA.
241*> \endverbatim
242*>
243*> \param[out] INFO
244*> \verbatim
245*> INFO is INTEGER
246*> = 0: successful exit.
247*> < 0: if INFO = -i, the i-th argument had an illegal value.
248*> \endverbatim
249*
250* Authors:
251* ========
252*
253*> \author Univ. of Tennessee
254*> \author Univ. of California Berkeley
255*> \author Univ. of Colorado Denver
256*> \author NAG Ltd.
257*
258*> \ingroup unbdb
259*
260*> \par Further Details:
261* =====================
262*>
263*> \verbatim
264*>
265*> The bidiagonal blocks B11, B12, B21, and B22 are represented
266*> implicitly by angles THETA(1), ..., THETA(Q) and PHI(1), ...,
267*> PHI(Q-1). B11 and B21 are upper bidiagonal, while B21 and B22 are
268*> lower bidiagonal. Every entry in each bidiagonal band is a product
269*> of a sine or cosine of a THETA with a sine or cosine of a PHI. See
270*> [1] or SORCSD for details.
271*>
272*> P1, P2, Q1, and Q2 are represented as products of elementary
273*> reflectors. See SORCSD for details on generating P1, P2, Q1, and Q2
274*> using SORGQR and SORGLQ.
275*> \endverbatim
276*
277*> \par References:
278* ================
279*>
280*> [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
281*> Algorithms, 50(1):33-65, 2009.
282*>
283* =====================================================================
284 SUBROUTINE sorbdb( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12,
285 $ X21, LDX21, X22, LDX22, THETA, PHI, TAUP1,
286 $ TAUP2, TAUQ1, TAUQ2, WORK, LWORK, INFO )
287*
288* -- LAPACK computational routine --
289* -- LAPACK is a software package provided by Univ. of Tennessee, --
290* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
291*
292* .. Scalar Arguments ..
293 CHARACTER SIGNS, TRANS
294 INTEGER INFO, LDX11, LDX12, LDX21, LDX22, LWORK, M, P,
295 $ q
296* ..
297* .. Array Arguments ..
298 REAL PHI( * ), THETA( * )
299 REAL TAUP1( * ), TAUP2( * ), TAUQ1( * ), TAUQ2( * ),
300 $ work( * ), x11( ldx11, * ), x12( ldx12, * ),
301 $ x21( ldx21, * ), x22( ldx22, * )
302* ..
303*
304* ====================================================================
305*
306* .. Parameters ..
307 REAL REALONE
308 PARAMETER ( REALONE = 1.0e0 )
309 REAL ONE
310 parameter( one = 1.0e0 )
311* ..
312* .. Local Scalars ..
313 LOGICAL COLMAJOR, LQUERY
314 INTEGER I, LWORKMIN, LWORKOPT
315 REAL Z1, Z2, Z3, Z4
316* ..
317* .. External Subroutines ..
318 EXTERNAL saxpy, slarf, slarfgp, sscal, xerbla
319* ..
320* .. External Functions ..
321 REAL SNRM2
322 LOGICAL LSAME
323 EXTERNAL snrm2, lsame
324* ..
325* .. Intrinsic Functions
326 INTRINSIC atan2, cos, max, sin
327* ..
328* .. Executable Statements ..
329*
330* Test input arguments
331*
332 info = 0
333 colmajor = .NOT. lsame( trans, 'T' )
334 IF( .NOT. lsame( signs, 'O' ) ) THEN
335 z1 = realone
336 z2 = realone
337 z3 = realone
338 z4 = realone
339 ELSE
340 z1 = realone
341 z2 = -realone
342 z3 = realone
343 z4 = -realone
344 END IF
345 lquery = lwork .EQ. -1
346*
347 IF( m .LT. 0 ) THEN
348 info = -3
349 ELSE IF( p .LT. 0 .OR. p .GT. m ) THEN
350 info = -4
351 ELSE IF( q .LT. 0 .OR. q .GT. p .OR. q .GT. m-p .OR.
352 $ q .GT. m-q ) THEN
353 info = -5
354 ELSE IF( colmajor .AND. ldx11 .LT. max( 1, p ) ) THEN
355 info = -7
356 ELSE IF( .NOT.colmajor .AND. ldx11 .LT. max( 1, q ) ) THEN
357 info = -7
358 ELSE IF( colmajor .AND. ldx12 .LT. max( 1, p ) ) THEN
359 info = -9
360 ELSE IF( .NOT.colmajor .AND. ldx12 .LT. max( 1, m-q ) ) THEN
361 info = -9
362 ELSE IF( colmajor .AND. ldx21 .LT. max( 1, m-p ) ) THEN
363 info = -11
364 ELSE IF( .NOT.colmajor .AND. ldx21 .LT. max( 1, q ) ) THEN
365 info = -11
366 ELSE IF( colmajor .AND. ldx22 .LT. max( 1, m-p ) ) THEN
367 info = -13
368 ELSE IF( .NOT.colmajor .AND. ldx22 .LT. max( 1, m-q ) ) THEN
369 info = -13
370 END IF
371*
372* Compute workspace
373*
374 IF( info .EQ. 0 ) THEN
375 lworkopt = m - q
376 lworkmin = m - q
377 work(1) = lworkopt
378 IF( lwork .LT. lworkmin .AND. .NOT. lquery ) THEN
379 info = -21
380 END IF
381 END IF
382 IF( info .NE. 0 ) THEN
383 CALL xerbla( 'xORBDB', -info )
384 RETURN
385 ELSE IF( lquery ) THEN
386 RETURN
387 END IF
388*
389* Handle column-major and row-major separately
390*
391 IF( colmajor ) THEN
392*
393* Reduce columns 1, ..., Q of X11, X12, X21, and X22
394*
395 DO i = 1, q
396*
397 IF( i .EQ. 1 ) THEN
398 CALL sscal( p-i+1, z1, x11(i,i), 1 )
399 ELSE
400 CALL sscal( p-i+1, z1*cos(phi(i-1)), x11(i,i), 1 )
401 CALL saxpy( p-i+1, -z1*z3*z4*sin(phi(i-1)), x12(i,i-1),
402 $ 1, x11(i,i), 1 )
403 END IF
404 IF( i .EQ. 1 ) THEN
405 CALL sscal( m-p-i+1, z2, x21(i,i), 1 )
406 ELSE
407 CALL sscal( m-p-i+1, z2*cos(phi(i-1)), x21(i,i), 1 )
408 CALL saxpy( m-p-i+1, -z2*z3*z4*sin(phi(i-1)), x22(i,i-1),
409 $ 1, x21(i,i), 1 )
410 END IF
411*
412 theta(i) = atan2( snrm2( m-p-i+1, x21(i,i), 1 ),
413 $ snrm2( p-i+1, x11(i,i), 1 ) )
414*
415 IF( p .GT. i ) THEN
416 CALL slarfgp( p-i+1, x11(i,i), x11(i+1,i), 1, taup1(i) )
417 ELSE IF( p .EQ. i ) THEN
418 CALL slarfgp( p-i+1, x11(i,i), x11(i,i), 1, taup1(i) )
419 END IF
420 x11(i,i) = one
421 IF ( m-p .GT. i ) THEN
422 CALL slarfgp( m-p-i+1, x21(i,i), x21(i+1,i), 1,
423 $ taup2(i) )
424 ELSE IF ( m-p .EQ. i ) THEN
425 CALL slarfgp( m-p-i+1, x21(i,i), x21(i,i), 1, taup2(i) )
426 END IF
427 x21(i,i) = one
428*
429 IF ( q .GT. i ) THEN
430 CALL slarf( 'L', p-i+1, q-i, x11(i,i), 1, taup1(i),
431 $ x11(i,i+1), ldx11, work )
432 END IF
433 IF ( m-q+1 .GT. i ) THEN
434 CALL slarf( 'L', p-i+1, m-q-i+1, x11(i,i), 1, taup1(i),
435 $ x12(i,i), ldx12, work )
436 END IF
437 IF ( q .GT. i ) THEN
438 CALL slarf( 'L', m-p-i+1, q-i, x21(i,i), 1, taup2(i),
439 $ x21(i,i+1), ldx21, work )
440 END IF
441 IF ( m-q+1 .GT. i ) THEN
442 CALL slarf( 'L', m-p-i+1, m-q-i+1, x21(i,i), 1, taup2(i),
443 $ x22(i,i), ldx22, work )
444 END IF
445*
446 IF( i .LT. q ) THEN
447 CALL sscal( q-i, -z1*z3*sin(theta(i)), x11(i,i+1),
448 $ ldx11 )
449 CALL saxpy( q-i, z2*z3*cos(theta(i)), x21(i,i+1), ldx21,
450 $ x11(i,i+1), ldx11 )
451 END IF
452 CALL sscal( m-q-i+1, -z1*z4*sin(theta(i)), x12(i,i), ldx12 )
453 CALL saxpy( m-q-i+1, z2*z4*cos(theta(i)), x22(i,i), ldx22,
454 $ x12(i,i), ldx12 )
455*
456 IF( i .LT. q )
457 $ phi(i) = atan2( snrm2( q-i, x11(i,i+1), ldx11 ),
458 $ snrm2( m-q-i+1, x12(i,i), ldx12 ) )
459*
460 IF( i .LT. q ) THEN
461 IF ( q-i .EQ. 1 ) THEN
462 CALL slarfgp( q-i, x11(i,i+1), x11(i,i+1), ldx11,
463 $ tauq1(i) )
464 ELSE
465 CALL slarfgp( q-i, x11(i,i+1), x11(i,i+2), ldx11,
466 $ tauq1(i) )
467 END IF
468 x11(i,i+1) = one
469 END IF
470 IF ( q+i-1 .LT. m ) THEN
471 IF ( m-q .EQ. i ) THEN
472 CALL slarfgp( m-q-i+1, x12(i,i), x12(i,i), ldx12,
473 $ tauq2(i) )
474 ELSE
475 CALL slarfgp( m-q-i+1, x12(i,i), x12(i,i+1), ldx12,
476 $ tauq2(i) )
477 END IF
478 END IF
479 x12(i,i) = one
480*
481 IF( i .LT. q ) THEN
482 CALL slarf( 'R', p-i, q-i, x11(i,i+1), ldx11, tauq1(i),
483 $ x11(i+1,i+1), ldx11, work )
484 CALL slarf( 'R', m-p-i, q-i, x11(i,i+1), ldx11, tauq1(i),
485 $ x21(i+1,i+1), ldx21, work )
486 END IF
487 IF ( p .GT. i ) THEN
488 CALL slarf( 'R', p-i, m-q-i+1, x12(i,i), ldx12, tauq2(i),
489 $ x12(i+1,i), ldx12, work )
490 END IF
491 IF ( m-p .GT. i ) THEN
492 CALL slarf( 'R', m-p-i, m-q-i+1, x12(i,i), ldx12,
493 $ tauq2(i), x22(i+1,i), ldx22, work )
494 END IF
495*
496 END DO
497*
498* Reduce columns Q + 1, ..., P of X12, X22
499*
500 DO i = q + 1, p
501*
502 CALL sscal( m-q-i+1, -z1*z4, x12(i,i), ldx12 )
503 IF ( i .GE. m-q ) THEN
504 CALL slarfgp( m-q-i+1, x12(i,i), x12(i,i), ldx12,
505 $ tauq2(i) )
506 ELSE
507 CALL slarfgp( m-q-i+1, x12(i,i), x12(i,i+1), ldx12,
508 $ tauq2(i) )
509 END IF
510 x12(i,i) = one
511*
512 IF ( p .GT. i ) THEN
513 CALL slarf( 'R', p-i, m-q-i+1, x12(i,i), ldx12, tauq2(i),
514 $ x12(i+1,i), ldx12, work )
515 END IF
516 IF( m-p-q .GE. 1 )
517 $ CALL slarf( 'R', m-p-q, m-q-i+1, x12(i,i), ldx12,
518 $ tauq2(i), x22(q+1,i), ldx22, work )
519*
520 END DO
521*
522* Reduce columns P + 1, ..., M - Q of X12, X22
523*
524 DO i = 1, m - p - q
525*
526 CALL sscal( m-p-q-i+1, z2*z4, x22(q+i,p+i), ldx22 )
527 IF ( i .EQ. m-p-q ) THEN
528 CALL slarfgp( m-p-q-i+1, x22(q+i,p+i), x22(q+i,p+i),
529 $ ldx22, tauq2(p+i) )
530 ELSE
531 CALL slarfgp( m-p-q-i+1, x22(q+i,p+i), x22(q+i,p+i+1),
532 $ ldx22, tauq2(p+i) )
533 END IF
534 x22(q+i,p+i) = one
535 IF ( i .LT. m-p-q ) THEN
536 CALL slarf( 'R', m-p-q-i, m-p-q-i+1, x22(q+i,p+i), ldx22,
537 $ tauq2(p+i), x22(q+i+1,p+i), ldx22, work )
538 END IF
539*
540 END DO
541*
542 ELSE
543*
544* Reduce columns 1, ..., Q of X11, X12, X21, X22
545*
546 DO i = 1, q
547*
548 IF( i .EQ. 1 ) THEN
549 CALL sscal( p-i+1, z1, x11(i,i), ldx11 )
550 ELSE
551 CALL sscal( p-i+1, z1*cos(phi(i-1)), x11(i,i), ldx11 )
552 CALL saxpy( p-i+1, -z1*z3*z4*sin(phi(i-1)), x12(i-1,i),
553 $ ldx12, x11(i,i), ldx11 )
554 END IF
555 IF( i .EQ. 1 ) THEN
556 CALL sscal( m-p-i+1, z2, x21(i,i), ldx21 )
557 ELSE
558 CALL sscal( m-p-i+1, z2*cos(phi(i-1)), x21(i,i), ldx21 )
559 CALL saxpy( m-p-i+1, -z2*z3*z4*sin(phi(i-1)), x22(i-1,i),
560 $ ldx22, x21(i,i), ldx21 )
561 END IF
562*
563 theta(i) = atan2( snrm2( m-p-i+1, x21(i,i), ldx21 ),
564 $ snrm2( p-i+1, x11(i,i), ldx11 ) )
565*
566 CALL slarfgp( p-i+1, x11(i,i), x11(i,i+1), ldx11, taup1(i) )
567 x11(i,i) = one
568 IF ( i .EQ. m-p ) THEN
569 CALL slarfgp( m-p-i+1, x21(i,i), x21(i,i), ldx21,
570 $ taup2(i) )
571 ELSE
572 CALL slarfgp( m-p-i+1, x21(i,i), x21(i,i+1), ldx21,
573 $ taup2(i) )
574 END IF
575 x21(i,i) = one
576*
577 IF ( q .GT. i ) THEN
578 CALL slarf( 'R', q-i, p-i+1, x11(i,i), ldx11, taup1(i),
579 $ x11(i+1,i), ldx11, work )
580 END IF
581 IF ( m-q+1 .GT. i ) THEN
582 CALL slarf( 'R', m-q-i+1, p-i+1, x11(i,i), ldx11,
583 $ taup1(i), x12(i,i), ldx12, work )
584 END IF
585 IF ( q .GT. i ) THEN
586 CALL slarf( 'R', q-i, m-p-i+1, x21(i,i), ldx21, taup2(i),
587 $ x21(i+1,i), ldx21, work )
588 END IF
589 IF ( m-q+1 .GT. i ) THEN
590 CALL slarf( 'R', m-q-i+1, m-p-i+1, x21(i,i), ldx21,
591 $ taup2(i), x22(i,i), ldx22, work )
592 END IF
593*
594 IF( i .LT. q ) THEN
595 CALL sscal( q-i, -z1*z3*sin(theta(i)), x11(i+1,i), 1 )
596 CALL saxpy( q-i, z2*z3*cos(theta(i)), x21(i+1,i), 1,
597 $ x11(i+1,i), 1 )
598 END IF
599 CALL sscal( m-q-i+1, -z1*z4*sin(theta(i)), x12(i,i), 1 )
600 CALL saxpy( m-q-i+1, z2*z4*cos(theta(i)), x22(i,i), 1,
601 $ x12(i,i), 1 )
602*
603 IF( i .LT. q )
604 $ phi(i) = atan2( snrm2( q-i, x11(i+1,i), 1 ),
605 $ snrm2( m-q-i+1, x12(i,i), 1 ) )
606*
607 IF( i .LT. q ) THEN
608 IF ( q-i .EQ. 1) THEN
609 CALL slarfgp( q-i, x11(i+1,i), x11(i+1,i), 1,
610 $ tauq1(i) )
611 ELSE
612 CALL slarfgp( q-i, x11(i+1,i), x11(i+2,i), 1,
613 $ tauq1(i) )
614 END IF
615 x11(i+1,i) = one
616 END IF
617 IF ( m-q .GT. i ) THEN
618 CALL slarfgp( m-q-i+1, x12(i,i), x12(i+1,i), 1,
619 $ tauq2(i) )
620 ELSE
621 CALL slarfgp( m-q-i+1, x12(i,i), x12(i,i), 1,
622 $ tauq2(i) )
623 END IF
624 x12(i,i) = one
625*
626 IF( i .LT. q ) THEN
627 CALL slarf( 'L', q-i, p-i, x11(i+1,i), 1, tauq1(i),
628 $ x11(i+1,i+1), ldx11, work )
629 CALL slarf( 'L', q-i, m-p-i, x11(i+1,i), 1, tauq1(i),
630 $ x21(i+1,i+1), ldx21, work )
631 END IF
632 CALL slarf( 'L', m-q-i+1, p-i, x12(i,i), 1, tauq2(i),
633 $ x12(i,i+1), ldx12, work )
634 IF ( m-p-i .GT. 0 ) THEN
635 CALL slarf( 'L', m-q-i+1, m-p-i, x12(i,i), 1, tauq2(i),
636 $ x22(i,i+1), ldx22, work )
637 END IF
638*
639 END DO
640*
641* Reduce columns Q + 1, ..., P of X12, X22
642*
643 DO i = q + 1, p
644*
645 CALL sscal( m-q-i+1, -z1*z4, x12(i,i), 1 )
646 CALL slarfgp( m-q-i+1, x12(i,i), x12(i+1,i), 1, tauq2(i) )
647 x12(i,i) = one
648*
649 IF ( p .GT. i ) THEN
650 CALL slarf( 'L', m-q-i+1, p-i, x12(i,i), 1, tauq2(i),
651 $ x12(i,i+1), ldx12, work )
652 END IF
653 IF( m-p-q .GE. 1 )
654 $ CALL slarf( 'L', m-q-i+1, m-p-q, x12(i,i), 1, tauq2(i),
655 $ x22(i,q+1), ldx22, work )
656*
657 END DO
658*
659* Reduce columns P + 1, ..., M - Q of X12, X22
660*
661 DO i = 1, m - p - q
662*
663 CALL sscal( m-p-q-i+1, z2*z4, x22(p+i,q+i), 1 )
664 IF ( m-p-q .EQ. i ) THEN
665 CALL slarfgp( m-p-q-i+1, x22(p+i,q+i), x22(p+i,q+i), 1,
666 $ tauq2(p+i) )
667 x22(p+i,q+i) = one
668 ELSE
669 CALL slarfgp( m-p-q-i+1, x22(p+i,q+i), x22(p+i+1,q+i), 1,
670 $ tauq2(p+i) )
671 x22(p+i,q+i) = one
672 CALL slarf( 'L', m-p-q-i+1, m-p-q-i, x22(p+i,q+i), 1,
673 $ tauq2(p+i), x22(p+i,q+i+1), ldx22, work )
674 END IF
675*
676*
677 END DO
678*
679 END IF
680*
681 RETURN
682*
683* End of SORBDB
684*
685 END
686
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine saxpy(n, sa, sx, incx, sy, incy)
SAXPY
Definition saxpy.f:89
subroutine slarf(side, m, n, v, incv, tau, c, ldc, work)
SLARF applies an elementary reflector to a general rectangular matrix.
Definition slarf.f:124
subroutine slarfgp(n, alpha, x, incx, tau)
SLARFGP generates an elementary reflector (Householder matrix) with non-negative beta.
Definition slarfgp.f:104
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
subroutine sorbdb(trans, signs, m, p, q, x11, ldx11, x12, ldx12, x21, ldx21, x22, ldx22, theta, phi, taup1, taup2, tauq1, tauq2, work, lwork, info)
SORBDB
Definition sorbdb.f:287