1*> \brief \b ZGGSVP
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download ZGGSVP + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zggsvp.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zggsvp.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zggsvp.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE ZGGSVP( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB,
22*                          TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ,
23*                          IWORK, RWORK, TAU, WORK, INFO )
24*
25*       .. Scalar Arguments ..
26*       CHARACTER          JOBQ, JOBU, JOBV
27*       INTEGER            INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P
28*       DOUBLE PRECISION   TOLA, TOLB
29*       ..
30*       .. Array Arguments ..
31*       INTEGER            IWORK( * )
32*       DOUBLE PRECISION   RWORK( * )
33*       COMPLEX*16         A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
34*      $                   TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * )
35*       ..
36*
37*
38*> \par Purpose:
39*  =============
40*>
41*> \verbatim
42*>
43*> This routine is deprecated and has been replaced by routine ZGGSVP3.
44*>
45*> ZGGSVP computes unitary matrices U, V and Q such that
46*>
47*>                    N-K-L  K    L
48*>  U**H*A*Q =     K ( 0    A12  A13 )  if M-K-L >= 0;
49*>                 L ( 0     0   A23 )
50*>             M-K-L ( 0     0    0  )
51*>
52*>                  N-K-L  K    L
53*>         =     K ( 0    A12  A13 )  if M-K-L < 0;
54*>             M-K ( 0     0   A23 )
55*>
56*>                  N-K-L  K    L
57*>  V**H*B*Q =   L ( 0     0   B13 )
58*>             P-L ( 0     0    0  )
59*>
60*> where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular
61*> upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0,
62*> otherwise A23 is (M-K)-by-L upper trapezoidal.  K+L = the effective
63*> numerical rank of the (M+P)-by-N matrix (A**H,B**H)**H.
64*>
65*> This decomposition is the preprocessing step for computing the
66*> Generalized Singular Value Decomposition (GSVD), see subroutine
67*> ZGGSVD.
68*> \endverbatim
69*
70*  Arguments:
71*  ==========
72*
73*> \param[in] JOBU
74*> \verbatim
75*>          JOBU is CHARACTER*1
76*>          = 'U':  Unitary matrix U is computed;
77*>          = 'N':  U is not computed.
78*> \endverbatim
79*>
80*> \param[in] JOBV
81*> \verbatim
82*>          JOBV is CHARACTER*1
83*>          = 'V':  Unitary matrix V is computed;
84*>          = 'N':  V is not computed.
85*> \endverbatim
86*>
87*> \param[in] JOBQ
88*> \verbatim
89*>          JOBQ is CHARACTER*1
90*>          = 'Q':  Unitary matrix Q is computed;
91*>          = 'N':  Q is not computed.
92*> \endverbatim
93*>
94*> \param[in] M
95*> \verbatim
96*>          M is INTEGER
97*>          The number of rows of the matrix A.  M >= 0.
98*> \endverbatim
99*>
100*> \param[in] P
101*> \verbatim
102*>          P is INTEGER
103*>          The number of rows of the matrix B.  P >= 0.
104*> \endverbatim
105*>
106*> \param[in] N
107*> \verbatim
108*>          N is INTEGER
109*>          The number of columns of the matrices A and B.  N >= 0.
110*> \endverbatim
111*>
112*> \param[in,out] A
113*> \verbatim
114*>          A is COMPLEX*16 array, dimension (LDA,N)
115*>          On entry, the M-by-N matrix A.
116*>          On exit, A contains the triangular (or trapezoidal) matrix
117*>          described in the Purpose section.
118*> \endverbatim
119*>
120*> \param[in] LDA
121*> \verbatim
122*>          LDA is INTEGER
123*>          The leading dimension of the array A. LDA >= max(1,M).
124*> \endverbatim
125*>
126*> \param[in,out] B
127*> \verbatim
128*>          B is COMPLEX*16 array, dimension (LDB,N)
129*>          On entry, the P-by-N matrix B.
130*>          On exit, B contains the triangular matrix described in
131*>          the Purpose section.
132*> \endverbatim
133*>
134*> \param[in] LDB
135*> \verbatim
136*>          LDB is INTEGER
137*>          The leading dimension of the array B. LDB >= max(1,P).
138*> \endverbatim
139*>
140*> \param[in] TOLA
141*> \verbatim
142*>          TOLA is DOUBLE PRECISION
143*> \endverbatim
144*>
145*> \param[in] TOLB
146*> \verbatim
147*>          TOLB is DOUBLE PRECISION
148*>
149*>          TOLA and TOLB are the thresholds to determine the effective
150*>          numerical rank of matrix B and a subblock of A. Generally,
151*>          they are set to
152*>             TOLA = MAX(M,N)*norm(A)*MAZHEPS,
153*>             TOLB = MAX(P,N)*norm(B)*MAZHEPS.
154*>          The size of TOLA and TOLB may affect the size of backward
155*>          errors of the decomposition.
156*> \endverbatim
157*>
158*> \param[out] K
159*> \verbatim
160*>          K is INTEGER
161*> \endverbatim
162*>
163*> \param[out] L
164*> \verbatim
165*>          L is INTEGER
166*>
167*>          On exit, K and L specify the dimension of the subblocks
168*>          described in Purpose section.
169*>          K + L = effective numerical rank of (A**H,B**H)**H.
170*> \endverbatim
171*>
172*> \param[out] U
173*> \verbatim
174*>          U is COMPLEX*16 array, dimension (LDU,M)
175*>          If JOBU = 'U', U contains the unitary matrix U.
176*>          If JOBU = 'N', U is not referenced.
177*> \endverbatim
178*>
179*> \param[in] LDU
180*> \verbatim
181*>          LDU is INTEGER
182*>          The leading dimension of the array U. LDU >= max(1,M) if
183*>          JOBU = 'U'; LDU >= 1 otherwise.
184*> \endverbatim
185*>
186*> \param[out] V
187*> \verbatim
188*>          V is COMPLEX*16 array, dimension (LDV,P)
189*>          If JOBV = 'V', V contains the unitary matrix V.
190*>          If JOBV = 'N', V is not referenced.
191*> \endverbatim
192*>
193*> \param[in] LDV
194*> \verbatim
195*>          LDV is INTEGER
196*>          The leading dimension of the array V. LDV >= max(1,P) if
197*>          JOBV = 'V'; LDV >= 1 otherwise.
198*> \endverbatim
199*>
200*> \param[out] Q
201*> \verbatim
202*>          Q is COMPLEX*16 array, dimension (LDQ,N)
203*>          If JOBQ = 'Q', Q contains the unitary matrix Q.
204*>          If JOBQ = 'N', Q is not referenced.
205*> \endverbatim
206*>
207*> \param[in] LDQ
208*> \verbatim
209*>          LDQ is INTEGER
210*>          The leading dimension of the array Q. LDQ >= max(1,N) if
211*>          JOBQ = 'Q'; LDQ >= 1 otherwise.
212*> \endverbatim
213*>
214*> \param[out] IWORK
215*> \verbatim
216*>          IWORK is INTEGER array, dimension (N)
217*> \endverbatim
218*>
219*> \param[out] RWORK
220*> \verbatim
221*>          RWORK is DOUBLE PRECISION array, dimension (2*N)
222*> \endverbatim
223*>
224*> \param[out] TAU
225*> \verbatim
226*>          TAU is COMPLEX*16 array, dimension (N)
227*> \endverbatim
228*>
229*> \param[out] WORK
230*> \verbatim
231*>          WORK is COMPLEX*16 array, dimension (max(3*N,M,P))
232*> \endverbatim
233*>
234*> \param[out] INFO
235*> \verbatim
236*>          INFO is INTEGER
237*>          = 0:  successful exit
238*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
239*> \endverbatim
240*
241*  Authors:
242*  ========
243*
244*> \author Univ. of Tennessee
245*> \author Univ. of California Berkeley
246*> \author Univ. of Colorado Denver
247*> \author NAG Ltd.
248*
249*> \date November 2011
250*
251*> \ingroup complex16OTHERcomputational
252*
253*> \par Further Details:
254*  =====================
255*>
256*> \verbatim
257*>
258*>  The subroutine uses LAPACK subroutine ZGEQPF for the QR factorization
259*>  with column pivoting to detect the effective numerical rank of the
260*>  a matrix. It may be replaced by a better rank determination strategy.
261*> \endverbatim
262*>
263*  =====================================================================
264      SUBROUTINE ZGGSVP( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB,
265     $                   TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ,
266     $                   IWORK, RWORK, TAU, WORK, INFO )
267*
268*  -- LAPACK computational routine (version 3.4.0) --
269*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
270*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
271*     November 2011
272*
273*     .. Scalar Arguments ..
274      CHARACTER          JOBQ, JOBU, JOBV
275      INTEGER            INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P
276      DOUBLE PRECISION   TOLA, TOLB
277*     ..
278*     .. Array Arguments ..
279      INTEGER            IWORK( * )
280      DOUBLE PRECISION   RWORK( * )
281      COMPLEX*16         A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
282     $                   TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * )
283*     ..
284*
285*  =====================================================================
286*
287*     .. Parameters ..
288      COMPLEX*16         CZERO, CONE
289      PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
290     $                   CONE = ( 1.0D+0, 0.0D+0 ) )
291*     ..
292*     .. Local Scalars ..
293      LOGICAL            FORWRD, WANTQ, WANTU, WANTV
294      INTEGER            I, J
295      COMPLEX*16         T
296*     ..
297*     .. External Functions ..
298      LOGICAL            LSAME
299      EXTERNAL           LSAME
300*     ..
301*     .. External Subroutines ..
302      EXTERNAL           XERBLA, ZGEQPF, ZGEQR2, ZGERQ2, ZLACPY, ZLAPMT,
303     $                   ZLASET, ZUNG2R, ZUNM2R, ZUNMR2
304*     ..
305*     .. Intrinsic Functions ..
306      INTRINSIC          ABS, DBLE, DIMAG, MAX, MIN
307*     ..
308*     .. Statement Functions ..
309      DOUBLE PRECISION   CABS1
310*     ..
311*     .. Statement Function definitions ..
312      CABS1( T ) = ABS( DBLE( T ) ) + ABS( DIMAG( T ) )
313*     ..
314*     .. Executable Statements ..
315*
316*     Test the input parameters
317*
318      WANTU = LSAME( JOBU, 'U' )
319      WANTV = LSAME( JOBV, 'V' )
320      WANTQ = LSAME( JOBQ, 'Q' )
321      FORWRD = .TRUE.
322*
323      INFO = 0
324      IF( .NOT.( WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN
325         INFO = -1
326      ELSE IF( .NOT.( WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN
327         INFO = -2
328      ELSE IF( .NOT.( WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN
329         INFO = -3
330      ELSE IF( M.LT.0 ) THEN
331         INFO = -4
332      ELSE IF( P.LT.0 ) THEN
333         INFO = -5
334      ELSE IF( N.LT.0 ) THEN
335         INFO = -6
336      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
337         INFO = -8
338      ELSE IF( LDB.LT.MAX( 1, P ) ) THEN
339         INFO = -10
340      ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN
341         INFO = -16
342      ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN
343         INFO = -18
344      ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
345         INFO = -20
346      END IF
347      IF( INFO.NE.0 ) THEN
348         CALL XERBLA( 'ZGGSVP', -INFO )
349         RETURN
350      END IF
351*
352*     QR with column pivoting of B: B*P = V*( S11 S12 )
353*                                           (  0   0  )
354*
355      DO 10 I = 1, N
356         IWORK( I ) = 0
357   10 CONTINUE
358      CALL ZGEQPF( P, N, B, LDB, IWORK, TAU, WORK, RWORK, INFO )
359*
360*     Update A := A*P
361*
362      CALL ZLAPMT( FORWRD, M, N, A, LDA, IWORK )
363*
364*     Determine the effective rank of matrix B.
365*
366      L = 0
367      DO 20 I = 1, MIN( P, N )
368         IF( CABS1( B( I, I ) ).GT.TOLB )
369     $      L = L + 1
370   20 CONTINUE
371*
372      IF( WANTV ) THEN
373*
374*        Copy the details of V, and form V.
375*
376         CALL ZLASET( 'Full', P, P, CZERO, CZERO, V, LDV )
377         IF( P.GT.1 )
378     $      CALL ZLACPY( 'Lower', P-1, N, B( 2, 1 ), LDB, V( 2, 1 ),
379     $                   LDV )
380         CALL ZUNG2R( P, P, MIN( P, N ), V, LDV, TAU, WORK, INFO )
381      END IF
382*
383*     Clean up B
384*
385      DO 40 J = 1, L - 1
386         DO 30 I = J + 1, L
387            B( I, J ) = CZERO
388   30    CONTINUE
389   40 CONTINUE
390      IF( P.GT.L )
391     $   CALL ZLASET( 'Full', P-L, N, CZERO, CZERO, B( L+1, 1 ), LDB )
392*
393      IF( WANTQ ) THEN
394*
395*        Set Q = I and Update Q := Q*P
396*
397         CALL ZLASET( 'Full', N, N, CZERO, CONE, Q, LDQ )
398         CALL ZLAPMT( FORWRD, N, N, Q, LDQ, IWORK )
399      END IF
400*
401      IF( P.GE.L .AND. N.NE.L ) THEN
402*
403*        RQ factorization of ( S11 S12 ) = ( 0 S12 )*Z
404*
405         CALL ZGERQ2( L, N, B, LDB, TAU, WORK, INFO )
406*
407*        Update A := A*Z**H
408*
409         CALL ZUNMR2( 'Right', 'Conjugate transpose', M, N, L, B, LDB,
410     $                TAU, A, LDA, WORK, INFO )
411         IF( WANTQ ) THEN
412*
413*           Update Q := Q*Z**H
414*
415            CALL ZUNMR2( 'Right', 'Conjugate transpose', N, N, L, B,
416     $                   LDB, TAU, Q, LDQ, WORK, INFO )
417         END IF
418*
419*        Clean up B
420*
421         CALL ZLASET( 'Full', L, N-L, CZERO, CZERO, B, LDB )
422         DO 60 J = N - L + 1, N
423            DO 50 I = J - N + L + 1, L
424               B( I, J ) = CZERO
425   50       CONTINUE
426   60    CONTINUE
427*
428      END IF
429*
430*     Let              N-L     L
431*                A = ( A11    A12 ) M,
432*
433*     then the following does the complete QR decomposition of A11:
434*
435*              A11 = U*(  0  T12 )*P1**H
436*                      (  0   0  )
437*
438      DO 70 I = 1, N - L
439         IWORK( I ) = 0
440   70 CONTINUE
441      CALL ZGEQPF( M, N-L, A, LDA, IWORK, TAU, WORK, RWORK, INFO )
442*
443*     Determine the effective rank of A11
444*
445      K = 0
446      DO 80 I = 1, MIN( M, N-L )
447         IF( CABS1( A( I, I ) ).GT.TOLA )
448     $      K = K + 1
449   80 CONTINUE
450*
451*     Update A12 := U**H*A12, where A12 = A( 1:M, N-L+1:N )
452*
453      CALL ZUNM2R( 'Left', 'Conjugate transpose', M, L, MIN( M, N-L ),
454     $             A, LDA, TAU, A( 1, N-L+1 ), LDA, WORK, INFO )
455*
456      IF( WANTU ) THEN
457*
458*        Copy the details of U, and form U
459*
460         CALL ZLASET( 'Full', M, M, CZERO, CZERO, U, LDU )
461         IF( M.GT.1 )
462     $      CALL ZLACPY( 'Lower', M-1, N-L, A( 2, 1 ), LDA, U( 2, 1 ),
463     $                   LDU )
464         CALL ZUNG2R( M, M, MIN( M, N-L ), U, LDU, TAU, WORK, INFO )
465      END IF
466*
467      IF( WANTQ ) THEN
468*
469*        Update Q( 1:N, 1:N-L )  = Q( 1:N, 1:N-L )*P1
470*
471         CALL ZLAPMT( FORWRD, N, N-L, Q, LDQ, IWORK )
472      END IF
473*
474*     Clean up A: set the strictly lower triangular part of
475*     A(1:K, 1:K) = 0, and A( K+1:M, 1:N-L ) = 0.
476*
477      DO 100 J = 1, K - 1
478         DO 90 I = J + 1, K
479            A( I, J ) = CZERO
480   90    CONTINUE
481  100 CONTINUE
482      IF( M.GT.K )
483     $   CALL ZLASET( 'Full', M-K, N-L, CZERO, CZERO, A( K+1, 1 ), LDA )
484*
485      IF( N-L.GT.K ) THEN
486*
487*        RQ factorization of ( T11 T12 ) = ( 0 T12 )*Z1
488*
489         CALL ZGERQ2( K, N-L, A, LDA, TAU, WORK, INFO )
490*
491         IF( WANTQ ) THEN
492*
493*           Update Q( 1:N,1:N-L ) = Q( 1:N,1:N-L )*Z1**H
494*
495            CALL ZUNMR2( 'Right', 'Conjugate transpose', N, N-L, K, A,
496     $                   LDA, TAU, Q, LDQ, WORK, INFO )
497         END IF
498*
499*        Clean up A
500*
501         CALL ZLASET( 'Full', K, N-L-K, CZERO, CZERO, A, LDA )
502         DO 120 J = N - L - K + 1, N - L
503            DO 110 I = J - N + L + K + 1, K
504               A( I, J ) = CZERO
505  110       CONTINUE
506  120    CONTINUE
507*
508      END IF
509*
510      IF( M.GT.K ) THEN
511*
512*        QR factorization of A( K+1:M,N-L+1:N )
513*
514         CALL ZGEQR2( M-K, L, A( K+1, N-L+1 ), LDA, TAU, WORK, INFO )
515*
516         IF( WANTU ) THEN
517*
518*           Update U(:,K+1:M) := U(:,K+1:M)*U1
519*
520            CALL ZUNM2R( 'Right', 'No transpose', M, M-K, MIN( M-K, L ),
521     $                   A( K+1, N-L+1 ), LDA, TAU, U( 1, K+1 ), LDU,
522     $                   WORK, INFO )
523         END IF
524*
525*        Clean up
526*
527         DO 140 J = N - L + 1, N
528            DO 130 I = J - N + K + L + 1, M
529               A( I, J ) = CZERO
530  130       CONTINUE
531  140    CONTINUE
532*
533      END IF
534*
535      RETURN
536*
537*     End of ZGGSVP
538*
539      END
540