1*> \brief \b DORBDB
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DORBDB + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dorbdb.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dorbdb.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dorbdb.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE DORBDB( 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*       DOUBLE PRECISION   PHI( * ), THETA( * )
32*       DOUBLE PRECISION   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*> DORBDB 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 DORCSD
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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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*> \date November 2015
259*
260*> \ingroup doubleOTHERcomputational
261*
262*> \par Further Details:
263*  =====================
264*>
265*> \verbatim
266*>
267*>  The bidiagonal blocks B11, B12, B21, and B22 are represented
268*>  implicitly by angles THETA(1), ..., THETA(Q) and PHI(1), ...,
269*>  PHI(Q-1). B11 and B21 are upper bidiagonal, while B21 and B22 are
270*>  lower bidiagonal. Every entry in each bidiagonal band is a product
271*>  of a sine or cosine of a THETA with a sine or cosine of a PHI. See
272*>  [1] or DORCSD for details.
273*>
274*>  P1, P2, Q1, and Q2 are represented as products of elementary
275*>  reflectors. See DORCSD for details on generating P1, P2, Q1, and Q2
276*>  using DORGQR and DORGLQ.
277*> \endverbatim
278*
279*> \par References:
280*  ================
281*>
282*>  [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
283*>      Algorithms, 50(1):33-65, 2009.
284*>
285*  =====================================================================
286      SUBROUTINE DORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12,
287     $                   X21, LDX21, X22, LDX22, THETA, PHI, TAUP1,
288     $                   TAUP2, TAUQ1, TAUQ2, WORK, LWORK, INFO )
289*
290*  -- LAPACK computational routine (version 3.6.0) --
291*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
292*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
293*     November 2015
294*
295*     .. Scalar Arguments ..
296      CHARACTER          SIGNS, TRANS
297      INTEGER            INFO, LDX11, LDX12, LDX21, LDX22, LWORK, M, P,
298     $                   Q
299*     ..
300*     .. Array Arguments ..
301      DOUBLE PRECISION   PHI( * ), THETA( * )
302      DOUBLE PRECISION   TAUP1( * ), TAUP2( * ), TAUQ1( * ), TAUQ2( * ),
303     $                   WORK( * ), X11( LDX11, * ), X12( LDX12, * ),
304     $                   X21( LDX21, * ), X22( LDX22, * )
305*     ..
306*
307*  ====================================================================
308*
309*     .. Parameters ..
310      DOUBLE PRECISION   REALONE
311      PARAMETER          ( REALONE = 1.0D0 )
312      DOUBLE PRECISION   ONE
313      PARAMETER          ( ONE = 1.0D0 )
314*     ..
315*     .. Local Scalars ..
316      LOGICAL            COLMAJOR, LQUERY
317      INTEGER            I, LWORKMIN, LWORKOPT
318      DOUBLE PRECISION   Z1, Z2, Z3, Z4
319*     ..
320*     .. External Subroutines ..
321      EXTERNAL           DAXPY, DLARF, DLARFGP, DSCAL, XERBLA
322*     ..
323*     .. External Functions ..
324      DOUBLE PRECISION   DNRM2
325      LOGICAL            LSAME
326      EXTERNAL           DNRM2, LSAME
327*     ..
328*     .. Intrinsic Functions
329      INTRINSIC          ATAN2, COS, MAX, SIN
330*     ..
331*     .. Executable Statements ..
332*
333*     Test input arguments
334*
335      INFO = 0
336      COLMAJOR = .NOT. LSAME( TRANS, 'T' )
337      IF( .NOT. LSAME( SIGNS, 'O' ) ) THEN
338         Z1 = REALONE
339         Z2 = REALONE
340         Z3 = REALONE
341         Z4 = REALONE
342      ELSE
343         Z1 = REALONE
344         Z2 = -REALONE
345         Z3 = REALONE
346         Z4 = -REALONE
347      END IF
348      LQUERY = LWORK .EQ. -1
349*
350      IF( M .LT. 0 ) THEN
351         INFO = -3
352      ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN
353         INFO = -4
354      ELSE IF( Q .LT. 0 .OR. Q .GT. P .OR. Q .GT. M-P .OR.
355     $         Q .GT. M-Q ) THEN
356         INFO = -5
357      ELSE IF( COLMAJOR .AND. LDX11 .LT. MAX( 1, P ) ) THEN
358         INFO = -7
359      ELSE IF( .NOT.COLMAJOR .AND. LDX11 .LT. MAX( 1, Q ) ) THEN
360         INFO = -7
361      ELSE IF( COLMAJOR .AND. LDX12 .LT. MAX( 1, P ) ) THEN
362         INFO = -9
363      ELSE IF( .NOT.COLMAJOR .AND. LDX12 .LT. MAX( 1, M-Q ) ) THEN
364         INFO = -9
365      ELSE IF( COLMAJOR .AND. LDX21 .LT. MAX( 1, M-P ) ) THEN
366         INFO = -11
367      ELSE IF( .NOT.COLMAJOR .AND. LDX21 .LT. MAX( 1, Q ) ) THEN
368         INFO = -11
369      ELSE IF( COLMAJOR .AND. LDX22 .LT. MAX( 1, M-P ) ) THEN
370         INFO = -13
371      ELSE IF( .NOT.COLMAJOR .AND. LDX22 .LT. MAX( 1, M-Q ) ) THEN
372         INFO = -13
373      END IF
374*
375*     Compute workspace
376*
377      IF( INFO .EQ. 0 ) THEN
378         LWORKOPT = M - Q
379         LWORKMIN = M - Q
380         WORK(1) = LWORKOPT
381         IF( LWORK .LT. LWORKMIN .AND. .NOT. LQUERY ) THEN
382            INFO = -21
383         END IF
384      END IF
385      IF( INFO .NE. 0 ) THEN
386         CALL XERBLA( 'xORBDB', -INFO )
387         RETURN
388      ELSE IF( LQUERY ) THEN
389         RETURN
390      END IF
391*
392*     Handle column-major and row-major separately
393*
394      IF( COLMAJOR ) THEN
395*
396*        Reduce columns 1, ..., Q of X11, X12, X21, and X22
397*
398         DO I = 1, Q
399*
400            IF( I .EQ. 1 ) THEN
401               CALL DSCAL( P-I+1, Z1, X11(I,I), 1 )
402            ELSE
403               CALL DSCAL( P-I+1, Z1*COS(PHI(I-1)), X11(I,I), 1 )
404               CALL DAXPY( P-I+1, -Z1*Z3*Z4*SIN(PHI(I-1)), X12(I,I-1),
405     $                     1, X11(I,I), 1 )
406            END IF
407            IF( I .EQ. 1 ) THEN
408               CALL DSCAL( M-P-I+1, Z2, X21(I,I), 1 )
409            ELSE
410               CALL DSCAL( M-P-I+1, Z2*COS(PHI(I-1)), X21(I,I), 1 )
411               CALL DAXPY( M-P-I+1, -Z2*Z3*Z4*SIN(PHI(I-1)), X22(I,I-1),
412     $                     1, X21(I,I), 1 )
413            END IF
414*
415            THETA(I) = ATAN2( DNRM2( M-P-I+1, X21(I,I), 1 ),
416     $                 DNRM2( P-I+1, X11(I,I), 1 ) )
417*
418            IF( P .GT. I ) THEN
419               CALL DLARFGP( P-I+1, X11(I,I), X11(I+1,I), 1, TAUP1(I) )
420            ELSE IF( P .EQ. I ) THEN
421               CALL DLARFGP( P-I+1, X11(I,I), X11(I,I), 1, TAUP1(I) )
422            END IF
423            X11(I,I) = ONE
424            IF ( M-P .GT. I ) THEN
425               CALL DLARFGP( M-P-I+1, X21(I,I), X21(I+1,I), 1,
426     $                       TAUP2(I) )
427            ELSE IF ( M-P .EQ. I ) THEN
428               CALL DLARFGP( M-P-I+1, X21(I,I), X21(I,I), 1, TAUP2(I) )
429            END IF
430            X21(I,I) = ONE
431*
432            IF ( Q .GT. I ) THEN
433               CALL DLARF( 'L', P-I+1, Q-I, X11(I,I), 1, TAUP1(I),
434     $                     X11(I,I+1), LDX11, WORK )
435            END IF
436            IF ( M-Q+1 .GT. I ) THEN
437               CALL DLARF( 'L', P-I+1, M-Q-I+1, X11(I,I), 1, TAUP1(I),
438     $                     X12(I,I), LDX12, WORK )
439            END IF
440            IF ( Q .GT. I ) THEN
441               CALL DLARF( 'L', M-P-I+1, Q-I, X21(I,I), 1, TAUP2(I),
442     $                     X21(I,I+1), LDX21, WORK )
443            END IF
444            IF ( M-Q+1 .GT. I ) THEN
445               CALL DLARF( 'L', M-P-I+1, M-Q-I+1, X21(I,I), 1, TAUP2(I),
446     $                     X22(I,I), LDX22, WORK )
447            END IF
448*
449            IF( I .LT. Q ) THEN
450               CALL DSCAL( Q-I, -Z1*Z3*SIN(THETA(I)), X11(I,I+1),
451     $                     LDX11 )
452               CALL DAXPY( Q-I, Z2*Z3*COS(THETA(I)), X21(I,I+1), LDX21,
453     $                     X11(I,I+1), LDX11 )
454            END IF
455            CALL DSCAL( M-Q-I+1, -Z1*Z4*SIN(THETA(I)), X12(I,I), LDX12 )
456            CALL DAXPY( M-Q-I+1, Z2*Z4*COS(THETA(I)), X22(I,I), LDX22,
457     $                  X12(I,I), LDX12 )
458*
459            IF( I .LT. Q )
460     $         PHI(I) = ATAN2( DNRM2( Q-I, X11(I,I+1), LDX11 ),
461     $                  DNRM2( M-Q-I+1, X12(I,I), LDX12 ) )
462*
463            IF( I .LT. Q ) THEN
464               IF ( Q-I .EQ. 1 ) THEN
465                  CALL DLARFGP( Q-I, X11(I,I+1), X11(I,I+1), LDX11,
466     $                          TAUQ1(I) )
467               ELSE
468                  CALL DLARFGP( Q-I, X11(I,I+1), X11(I,I+2), LDX11,
469     $                          TAUQ1(I) )
470               END IF
471               X11(I,I+1) = ONE
472            END IF
473            IF ( Q+I-1 .LT. M ) THEN
474               IF ( M-Q .EQ. I ) THEN
475                  CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I,I), LDX12,
476     $                          TAUQ2(I) )
477               ELSE
478                  CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12,
479     $                          TAUQ2(I) )
480               END IF
481            END IF
482            X12(I,I) = ONE
483*
484            IF( I .LT. Q ) THEN
485               CALL DLARF( 'R', P-I, Q-I, X11(I,I+1), LDX11, TAUQ1(I),
486     $                     X11(I+1,I+1), LDX11, WORK )
487               CALL DLARF( 'R', M-P-I, Q-I, X11(I,I+1), LDX11, TAUQ1(I),
488     $                     X21(I+1,I+1), LDX21, WORK )
489            END IF
490            IF ( P .GT. I ) THEN
491               CALL DLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
492     $                     X12(I+1,I), LDX12, WORK )
493            END IF
494            IF ( M-P .GT. I ) THEN
495               CALL DLARF( 'R', M-P-I, M-Q-I+1, X12(I,I), LDX12,
496     $                     TAUQ2(I), X22(I+1,I), LDX22, WORK )
497            END IF
498*
499         END DO
500*
501*        Reduce columns Q + 1, ..., P of X12, X22
502*
503         DO I = Q + 1, P
504*
505            CALL DSCAL( M-Q-I+1, -Z1*Z4, X12(I,I), LDX12 )
506            IF ( I .GE. M-Q ) THEN
507               CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I,I), LDX12,
508     $                       TAUQ2(I) )
509            ELSE
510               CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12,
511     $                       TAUQ2(I) )
512            END IF
513            X12(I,I) = ONE
514*
515            IF ( P .GT. I ) THEN
516               CALL DLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, TAUQ2(I),
517     $                     X12(I+1,I), LDX12, WORK )
518            END IF
519            IF( M-P-Q .GE. 1 )
520     $         CALL DLARF( 'R', M-P-Q, M-Q-I+1, X12(I,I), LDX12,
521     $                     TAUQ2(I), X22(Q+1,I), LDX22, WORK )
522*
523         END DO
524*
525*        Reduce columns P + 1, ..., M - Q of X12, X22
526*
527         DO I = 1, M - P - Q
528*
529            CALL DSCAL( M-P-Q-I+1, Z2*Z4, X22(Q+I,P+I), LDX22 )
530            IF ( I .EQ. M-P-Q ) THEN
531               CALL DLARFGP( M-P-Q-I+1, X22(Q+I,P+I), X22(Q+I,P+I),
532     $                       LDX22, TAUQ2(P+I) )
533            ELSE
534               CALL DLARFGP( M-P-Q-I+1, X22(Q+I,P+I), X22(Q+I,P+I+1),
535     $                       LDX22, TAUQ2(P+I) )
536            END IF
537            X22(Q+I,P+I) = ONE
538            IF ( I .LT. M-P-Q ) THEN
539               CALL DLARF( 'R', M-P-Q-I, M-P-Q-I+1, X22(Q+I,P+I), LDX22,
540     $                     TAUQ2(P+I), X22(Q+I+1,P+I), LDX22, WORK )
541            END IF
542*
543         END DO
544*
545      ELSE
546*
547*        Reduce columns 1, ..., Q of X11, X12, X21, X22
548*
549         DO I = 1, Q
550*
551            IF( I .EQ. 1 ) THEN
552               CALL DSCAL( P-I+1, Z1, X11(I,I), LDX11 )
553            ELSE
554               CALL DSCAL( P-I+1, Z1*COS(PHI(I-1)), X11(I,I), LDX11 )
555               CALL DAXPY( P-I+1, -Z1*Z3*Z4*SIN(PHI(I-1)), X12(I-1,I),
556     $                     LDX12, X11(I,I), LDX11 )
557            END IF
558            IF( I .EQ. 1 ) THEN
559               CALL DSCAL( M-P-I+1, Z2, X21(I,I), LDX21 )
560            ELSE
561               CALL DSCAL( M-P-I+1, Z2*COS(PHI(I-1)), X21(I,I), LDX21 )
562               CALL DAXPY( M-P-I+1, -Z2*Z3*Z4*SIN(PHI(I-1)), X22(I-1,I),
563     $                     LDX22, X21(I,I), LDX21 )
564            END IF
565*
566            THETA(I) = ATAN2( DNRM2( M-P-I+1, X21(I,I), LDX21 ),
567     $                 DNRM2( P-I+1, X11(I,I), LDX11 ) )
568*
569            CALL DLARFGP( P-I+1, X11(I,I), X11(I,I+1), LDX11, TAUP1(I) )
570            X11(I,I) = ONE
571            IF ( I .EQ. M-P ) THEN
572               CALL DLARFGP( M-P-I+1, X21(I,I), X21(I,I), LDX21,
573     $                    TAUP2(I) )
574            ELSE
575               CALL DLARFGP( M-P-I+1, X21(I,I), X21(I,I+1), LDX21,
576     $                    TAUP2(I) )
577            END IF
578            X21(I,I) = ONE
579*
580            IF ( Q .GT. I ) THEN
581               CALL DLARF( 'R', Q-I, P-I+1, X11(I,I), LDX11, TAUP1(I),
582     $                     X11(I+1,I), LDX11, WORK )
583            END IF
584            IF ( M-Q+1 .GT. I ) THEN
585               CALL DLARF( 'R', M-Q-I+1, P-I+1, X11(I,I), LDX11,
586     $                     TAUP1(I), X12(I,I), LDX12, WORK )
587            END IF
588            IF ( Q .GT. I ) THEN
589               CALL DLARF( 'R', Q-I, M-P-I+1, X21(I,I), LDX21, TAUP2(I),
590     $                     X21(I+1,I), LDX21, WORK )
591            END IF
592            IF ( M-Q+1 .GT. I ) THEN
593               CALL DLARF( 'R', M-Q-I+1, M-P-I+1, X21(I,I), LDX21,
594     $                     TAUP2(I), X22(I,I), LDX22, WORK )
595            END IF
596*
597            IF( I .LT. Q ) THEN
598               CALL DSCAL( Q-I, -Z1*Z3*SIN(THETA(I)), X11(I+1,I), 1 )
599               CALL DAXPY( Q-I, Z2*Z3*COS(THETA(I)), X21(I+1,I), 1,
600     $                     X11(I+1,I), 1 )
601            END IF
602            CALL DSCAL( M-Q-I+1, -Z1*Z4*SIN(THETA(I)), X12(I,I), 1 )
603            CALL DAXPY( M-Q-I+1, Z2*Z4*COS(THETA(I)), X22(I,I), 1,
604     $                  X12(I,I), 1 )
605*
606            IF( I .LT. Q )
607     $         PHI(I) = ATAN2( DNRM2( Q-I, X11(I+1,I), 1 ),
608     $                  DNRM2( M-Q-I+1, X12(I,I), 1 ) )
609*
610            IF( I .LT. Q ) THEN
611               IF ( Q-I .EQ. 1) THEN
612                  CALL DLARFGP( Q-I, X11(I+1,I), X11(I+1,I), 1,
613     $                          TAUQ1(I) )
614               ELSE
615                  CALL DLARFGP( Q-I, X11(I+1,I), X11(I+2,I), 1,
616     $                          TAUQ1(I) )
617               END IF
618               X11(I+1,I) = ONE
619            END IF
620            IF ( M-Q .GT. I ) THEN
621               CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1,
622     $                       TAUQ2(I) )
623            ELSE
624               CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I,I), 1,
625     $                       TAUQ2(I) )
626            END IF
627            X12(I,I) = ONE
628*
629            IF( I .LT. Q ) THEN
630               CALL DLARF( 'L', Q-I, P-I, X11(I+1,I), 1, TAUQ1(I),
631     $                     X11(I+1,I+1), LDX11, WORK )
632               CALL DLARF( 'L', Q-I, M-P-I, X11(I+1,I), 1, TAUQ1(I),
633     $                     X21(I+1,I+1), LDX21, WORK )
634            END IF
635            CALL DLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1, TAUQ2(I),
636     $                  X12(I,I+1), LDX12, WORK )
637            IF ( M-P-I .GT. 0 ) THEN
638               CALL DLARF( 'L', M-Q-I+1, M-P-I, X12(I,I), 1, TAUQ2(I),
639     $                     X22(I,I+1), LDX22, WORK )
640            END IF
641*
642         END DO
643*
644*        Reduce columns Q + 1, ..., P of X12, X22
645*
646         DO I = Q + 1, P
647*
648            CALL DSCAL( M-Q-I+1, -Z1*Z4, X12(I,I), 1 )
649            CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1, TAUQ2(I) )
650            X12(I,I) = ONE
651*
652            IF ( P .GT. I ) THEN
653               CALL DLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1, TAUQ2(I),
654     $                  X12(I,I+1), LDX12, WORK )
655            END IF
656            IF( M-P-Q .GE. 1 )
657     $         CALL DLARF( 'L', M-Q-I+1, M-P-Q, X12(I,I), 1, TAUQ2(I),
658     $                     X22(I,Q+1), LDX22, WORK )
659*
660         END DO
661*
662*        Reduce columns P + 1, ..., M - Q of X12, X22
663*
664         DO I = 1, M - P - Q
665*
666            CALL DSCAL( M-P-Q-I+1, Z2*Z4, X22(P+I,Q+I), 1 )
667            IF ( M-P-Q .EQ. I ) THEN
668               CALL DLARFGP( M-P-Q-I+1, X22(P+I,Q+I), X22(P+I,Q+I), 1,
669     $                       TAUQ2(P+I) )
670            ELSE
671               CALL DLARFGP( M-P-Q-I+1, X22(P+I,Q+I), X22(P+I+1,Q+I), 1,
672     $                       TAUQ2(P+I) )
673               CALL DLARF( 'L', M-P-Q-I+1, M-P-Q-I, X22(P+I,Q+I), 1,
674     $                  TAUQ2(P+I), X22(P+I,Q+I+1), LDX22, WORK )
675            END IF
676            X22(P+I,Q+I) = ONE
677*
678         END DO
679*
680      END IF
681*
682      RETURN
683*
684*     End of DORBDB
685*
686      END
687
688