001:       SUBROUTINE ZGGSVP( JOBU, JOBV, JOBQ, M, P, N, A, LDA, B, LDB,
002:      $                   TOLA, TOLB, K, L, U, LDU, V, LDV, Q, LDQ,
003:      $                   IWORK, RWORK, TAU, WORK, INFO )
004: *
005: *  -- LAPACK routine (version 3.2) --
006: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
007: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
008: *     November 2006
009: *
010: *     .. Scalar Arguments ..
011:       CHARACTER          JOBQ, JOBU, JOBV
012:       INTEGER            INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N, P
013:       DOUBLE PRECISION   TOLA, TOLB
014: *     ..
015: *     .. Array Arguments ..
016:       INTEGER            IWORK( * )
017:       DOUBLE PRECISION   RWORK( * )
018:       COMPLEX*16         A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
019:      $                   TAU( * ), U( LDU, * ), V( LDV, * ), WORK( * )
020: *     ..
021: *
022: *  Purpose
023: *  =======
024: *
025: *  ZGGSVP computes unitary matrices U, V and Q such that
026: *
027: *                   N-K-L  K    L
028: *   U'*A*Q =     K ( 0    A12  A13 )  if M-K-L >= 0;
029: *                L ( 0     0   A23 )
030: *            M-K-L ( 0     0    0  )
031: *
032: *                   N-K-L  K    L
033: *          =     K ( 0    A12  A13 )  if M-K-L < 0;
034: *              M-K ( 0     0   A23 )
035: *
036: *                 N-K-L  K    L
037: *   V'*B*Q =   L ( 0     0   B13 )
038: *            P-L ( 0     0    0  )
039: *
040: *  where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular
041: *  upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0,
042: *  otherwise A23 is (M-K)-by-L upper trapezoidal.  K+L = the effective
043: *  numerical rank of the (M+P)-by-N matrix (A',B')'.  Z' denotes the
044: *  conjugate transpose of Z.
045: *
046: *  This decomposition is the preprocessing step for computing the
047: *  Generalized Singular Value Decomposition (GSVD), see subroutine
048: *  ZGGSVD.
049: *
050: *  Arguments
051: *  =========
052: *
053: *  JOBU    (input) CHARACTER*1
054: *          = 'U':  Unitary matrix U is computed;
055: *          = 'N':  U is not computed.
056: *
057: *  JOBV    (input) CHARACTER*1
058: *          = 'V':  Unitary matrix V is computed;
059: *          = 'N':  V is not computed.
060: *
061: *  JOBQ    (input) CHARACTER*1
062: *          = 'Q':  Unitary matrix Q is computed;
063: *          = 'N':  Q is not computed.
064: *
065: *  M       (input) INTEGER
066: *          The number of rows of the matrix A.  M >= 0.
067: *
068: *  P       (input) INTEGER
069: *          The number of rows of the matrix B.  P >= 0.
070: *
071: *  N       (input) INTEGER
072: *          The number of columns of the matrices A and B.  N >= 0.
073: *
074: *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
075: *          On entry, the M-by-N matrix A.
076: *          On exit, A contains the triangular (or trapezoidal) matrix
077: *          described in the Purpose section.
078: *
079: *  LDA     (input) INTEGER
080: *          The leading dimension of the array A. LDA >= max(1,M).
081: *
082: *  B       (input/output) COMPLEX*16 array, dimension (LDB,N)
083: *          On entry, the P-by-N matrix B.
084: *          On exit, B contains the triangular matrix described in
085: *          the Purpose section.
086: *
087: *  LDB     (input) INTEGER
088: *          The leading dimension of the array B. LDB >= max(1,P).
089: *
090: *  TOLA    (input) DOUBLE PRECISION
091: *  TOLB    (input) DOUBLE PRECISION
092: *          TOLA and TOLB are the thresholds to determine the effective
093: *          numerical rank of matrix B and a subblock of A. Generally,
094: *          they are set to
095: *             TOLA = MAX(M,N)*norm(A)*MAZHEPS,
096: *             TOLB = MAX(P,N)*norm(B)*MAZHEPS.
097: *          The size of TOLA and TOLB may affect the size of backward
098: *          errors of the decomposition.
099: *
100: *  K       (output) INTEGER
101: *  L       (output) INTEGER
102: *          On exit, K and L specify the dimension of the subblocks
103: *          described in Purpose section.
104: *          K + L = effective numerical rank of (A',B')'.
105: *
106: *  U       (output) COMPLEX*16 array, dimension (LDU,M)
107: *          If JOBU = 'U', U contains the unitary matrix U.
108: *          If JOBU = 'N', U is not referenced.
109: *
110: *  LDU     (input) INTEGER
111: *          The leading dimension of the array U. LDU >= max(1,M) if
112: *          JOBU = 'U'; LDU >= 1 otherwise.
113: *
114: *  V       (output) COMPLEX*16 array, dimension (LDV,P)
115: *          If JOBV = 'V', V contains the unitary matrix V.
116: *          If JOBV = 'N', V is not referenced.
117: *
118: *  LDV     (input) INTEGER
119: *          The leading dimension of the array V. LDV >= max(1,P) if
120: *          JOBV = 'V'; LDV >= 1 otherwise.
121: *
122: *  Q       (output) COMPLEX*16 array, dimension (LDQ,N)
123: *          If JOBQ = 'Q', Q contains the unitary matrix Q.
124: *          If JOBQ = 'N', Q is not referenced.
125: *
126: *  LDQ     (input) INTEGER
127: *          The leading dimension of the array Q. LDQ >= max(1,N) if
128: *          JOBQ = 'Q'; LDQ >= 1 otherwise.
129: *
130: *  IWORK   (workspace) INTEGER array, dimension (N)
131: *
132: *  RWORK   (workspace) DOUBLE PRECISION array, dimension (2*N)
133: *
134: *  TAU     (workspace) COMPLEX*16 array, dimension (N)
135: *
136: *  WORK    (workspace) COMPLEX*16 array, dimension (max(3*N,M,P))
137: *
138: *  INFO    (output) INTEGER
139: *          = 0:  successful exit
140: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
141: *
142: *  Further Details
143: *  ===============
144: *
145: *  The subroutine uses LAPACK subroutine ZGEQPF for the QR factorization
146: *  with column pivoting to detect the effective numerical rank of the
147: *  a matrix. It may be replaced by a better rank determination strategy.
148: *
149: *  =====================================================================
150: *
151: *     .. Parameters ..
152:       COMPLEX*16         CZERO, CONE
153:       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
154:      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
155: *     ..
156: *     .. Local Scalars ..
157:       LOGICAL            FORWRD, WANTQ, WANTU, WANTV
158:       INTEGER            I, J
159:       COMPLEX*16         T
160: *     ..
161: *     .. External Functions ..
162:       LOGICAL            LSAME
163:       EXTERNAL           LSAME
164: *     ..
165: *     .. External Subroutines ..
166:       EXTERNAL           XERBLA, ZGEQPF, ZGEQR2, ZGERQ2, ZLACPY, ZLAPMT,
167:      $                   ZLASET, ZUNG2R, ZUNM2R, ZUNMR2
168: *     ..
169: *     .. Intrinsic Functions ..
170:       INTRINSIC          ABS, DBLE, DIMAG, MAX, MIN
171: *     ..
172: *     .. Statement Functions ..
173:       DOUBLE PRECISION   CABS1
174: *     ..
175: *     .. Statement Function definitions ..
176:       CABS1( T ) = ABS( DBLE( T ) ) + ABS( DIMAG( T ) )
177: *     ..
178: *     .. Executable Statements ..
179: *
180: *     Test the input parameters
181: *
182:       WANTU = LSAME( JOBU, 'U' )
183:       WANTV = LSAME( JOBV, 'V' )
184:       WANTQ = LSAME( JOBQ, 'Q' )
185:       FORWRD = .TRUE.
186: *
187:       INFO = 0
188:       IF( .NOT.( WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN
189:          INFO = -1
190:       ELSE IF( .NOT.( WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN
191:          INFO = -2
192:       ELSE IF( .NOT.( WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN
193:          INFO = -3
194:       ELSE IF( M.LT.0 ) THEN
195:          INFO = -4
196:       ELSE IF( P.LT.0 ) THEN
197:          INFO = -5
198:       ELSE IF( N.LT.0 ) THEN
199:          INFO = -6
200:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
201:          INFO = -8
202:       ELSE IF( LDB.LT.MAX( 1, P ) ) THEN
203:          INFO = -10
204:       ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN
205:          INFO = -16
206:       ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN
207:          INFO = -18
208:       ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
209:          INFO = -20
210:       END IF
211:       IF( INFO.NE.0 ) THEN
212:          CALL XERBLA( 'ZGGSVP', -INFO )
213:          RETURN
214:       END IF
215: *
216: *     QR with column pivoting of B: B*P = V*( S11 S12 )
217: *                                           (  0   0  )
218: *
219:       DO 10 I = 1, N
220:          IWORK( I ) = 0
221:    10 CONTINUE
222:       CALL ZGEQPF( P, N, B, LDB, IWORK, TAU, WORK, RWORK, INFO )
223: *
224: *     Update A := A*P
225: *
226:       CALL ZLAPMT( FORWRD, M, N, A, LDA, IWORK )
227: *
228: *     Determine the effective rank of matrix B.
229: *
230:       L = 0
231:       DO 20 I = 1, MIN( P, N )
232:          IF( CABS1( B( I, I ) ).GT.TOLB )
233:      $      L = L + 1
234:    20 CONTINUE
235: *
236:       IF( WANTV ) THEN
237: *
238: *        Copy the details of V, and form V.
239: *
240:          CALL ZLASET( 'Full', P, P, CZERO, CZERO, V, LDV )
241:          IF( P.GT.1 )
242:      $      CALL ZLACPY( 'Lower', P-1, N, B( 2, 1 ), LDB, V( 2, 1 ),
243:      $                   LDV )
244:          CALL ZUNG2R( P, P, MIN( P, N ), V, LDV, TAU, WORK, INFO )
245:       END IF
246: *
247: *     Clean up B
248: *
249:       DO 40 J = 1, L - 1
250:          DO 30 I = J + 1, L
251:             B( I, J ) = CZERO
252:    30    CONTINUE
253:    40 CONTINUE
254:       IF( P.GT.L )
255:      $   CALL ZLASET( 'Full', P-L, N, CZERO, CZERO, B( L+1, 1 ), LDB )
256: *
257:       IF( WANTQ ) THEN
258: *
259: *        Set Q = I and Update Q := Q*P
260: *
261:          CALL ZLASET( 'Full', N, N, CZERO, CONE, Q, LDQ )
262:          CALL ZLAPMT( FORWRD, N, N, Q, LDQ, IWORK )
263:       END IF
264: *
265:       IF( P.GE.L .AND. N.NE.L ) THEN
266: *
267: *        RQ factorization of ( S11 S12 ) = ( 0 S12 )*Z
268: *
269:          CALL ZGERQ2( L, N, B, LDB, TAU, WORK, INFO )
270: *
271: *        Update A := A*Z'
272: *
273:          CALL ZUNMR2( 'Right', 'Conjugate transpose', M, N, L, B, LDB,
274:      $                TAU, A, LDA, WORK, INFO )
275:          IF( WANTQ ) THEN
276: *
277: *           Update Q := Q*Z'
278: *
279:             CALL ZUNMR2( 'Right', 'Conjugate transpose', N, N, L, B,
280:      $                   LDB, TAU, Q, LDQ, WORK, INFO )
281:          END IF
282: *
283: *        Clean up B
284: *
285:          CALL ZLASET( 'Full', L, N-L, CZERO, CZERO, B, LDB )
286:          DO 60 J = N - L + 1, N
287:             DO 50 I = J - N + L + 1, L
288:                B( I, J ) = CZERO
289:    50       CONTINUE
290:    60    CONTINUE
291: *
292:       END IF
293: *
294: *     Let              N-L     L
295: *                A = ( A11    A12 ) M,
296: *
297: *     then the following does the complete QR decomposition of A11:
298: *
299: *              A11 = U*(  0  T12 )*P1'
300: *                      (  0   0  )
301: *
302:       DO 70 I = 1, N - L
303:          IWORK( I ) = 0
304:    70 CONTINUE
305:       CALL ZGEQPF( M, N-L, A, LDA, IWORK, TAU, WORK, RWORK, INFO )
306: *
307: *     Determine the effective rank of A11
308: *
309:       K = 0
310:       DO 80 I = 1, MIN( M, N-L )
311:          IF( CABS1( A( I, I ) ).GT.TOLA )
312:      $      K = K + 1
313:    80 CONTINUE
314: *
315: *     Update A12 := U'*A12, where A12 = A( 1:M, N-L+1:N )
316: *
317:       CALL ZUNM2R( 'Left', 'Conjugate transpose', M, L, MIN( M, N-L ),
318:      $             A, LDA, TAU, A( 1, N-L+1 ), LDA, WORK, INFO )
319: *
320:       IF( WANTU ) THEN
321: *
322: *        Copy the details of U, and form U
323: *
324:          CALL ZLASET( 'Full', M, M, CZERO, CZERO, U, LDU )
325:          IF( M.GT.1 )
326:      $      CALL ZLACPY( 'Lower', M-1, N-L, A( 2, 1 ), LDA, U( 2, 1 ),
327:      $                   LDU )
328:          CALL ZUNG2R( M, M, MIN( M, N-L ), U, LDU, TAU, WORK, INFO )
329:       END IF
330: *
331:       IF( WANTQ ) THEN
332: *
333: *        Update Q( 1:N, 1:N-L )  = Q( 1:N, 1:N-L )*P1
334: *
335:          CALL ZLAPMT( FORWRD, N, N-L, Q, LDQ, IWORK )
336:       END IF
337: *
338: *     Clean up A: set the strictly lower triangular part of
339: *     A(1:K, 1:K) = 0, and A( K+1:M, 1:N-L ) = 0.
340: *
341:       DO 100 J = 1, K - 1
342:          DO 90 I = J + 1, K
343:             A( I, J ) = CZERO
344:    90    CONTINUE
345:   100 CONTINUE
346:       IF( M.GT.K )
347:      $   CALL ZLASET( 'Full', M-K, N-L, CZERO, CZERO, A( K+1, 1 ), LDA )
348: *
349:       IF( N-L.GT.K ) THEN
350: *
351: *        RQ factorization of ( T11 T12 ) = ( 0 T12 )*Z1
352: *
353:          CALL ZGERQ2( K, N-L, A, LDA, TAU, WORK, INFO )
354: *
355:          IF( WANTQ ) THEN
356: *
357: *           Update Q( 1:N,1:N-L ) = Q( 1:N,1:N-L )*Z1'
358: *
359:             CALL ZUNMR2( 'Right', 'Conjugate transpose', N, N-L, K, A,
360:      $                   LDA, TAU, Q, LDQ, WORK, INFO )
361:          END IF
362: *
363: *        Clean up A
364: *
365:          CALL ZLASET( 'Full', K, N-L-K, CZERO, CZERO, A, LDA )
366:          DO 120 J = N - L - K + 1, N - L
367:             DO 110 I = J - N + L + K + 1, K
368:                A( I, J ) = CZERO
369:   110       CONTINUE
370:   120    CONTINUE
371: *
372:       END IF
373: *
374:       IF( M.GT.K ) THEN
375: *
376: *        QR factorization of A( K+1:M,N-L+1:N )
377: *
378:          CALL ZGEQR2( M-K, L, A( K+1, N-L+1 ), LDA, TAU, WORK, INFO )
379: *
380:          IF( WANTU ) THEN
381: *
382: *           Update U(:,K+1:M) := U(:,K+1:M)*U1
383: *
384:             CALL ZUNM2R( 'Right', 'No transpose', M, M-K, MIN( M-K, L ),
385:      $                   A( K+1, N-L+1 ), LDA, TAU, U( 1, K+1 ), LDU,
386:      $                   WORK, INFO )
387:          END IF
388: *
389: *        Clean up
390: *
391:          DO 140 J = N - L + 1, N
392:             DO 130 I = J - N + K + L + 1, M
393:                A( I, J ) = CZERO
394:   130       CONTINUE
395:   140    CONTINUE
396: *
397:       END IF
398: *
399:       RETURN
400: *
401: *     End of ZGGSVP
402: *
403:       END
404: