LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ cgesvd()

subroutine cgesvd ( character  JOBU,
character  JOBVT,
integer  M,
integer  N,
complex, dimension( lda, * )  A,
integer  LDA,
real, dimension( * )  S,
complex, dimension( ldu, * )  U,
integer  LDU,
complex, dimension( ldvt, * )  VT,
integer  LDVT,
complex, dimension( * )  WORK,
integer  LWORK,
real, dimension( * )  RWORK,
integer  INFO 
)

CGESVD computes the singular value decomposition (SVD) for GE matrices

Download CGESVD + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 CGESVD computes the singular value decomposition (SVD) of a complex
 M-by-N matrix A, optionally computing the left and/or right singular
 vectors. The SVD is written

      A = U * SIGMA * conjugate-transpose(V)

 where SIGMA is an M-by-N matrix which is zero except for its
 min(m,n) diagonal elements, U is an M-by-M unitary matrix, and
 V is an N-by-N unitary matrix.  The diagonal elements of SIGMA
 are the singular values of A; they are real and non-negative, and
 are returned in descending order.  The first min(m,n) columns of
 U and V are the left and right singular vectors of A.

 Note that the routine returns V**H, not V.
Parameters
[in]JOBU
          JOBU is CHARACTER*1
          Specifies options for computing all or part of the matrix U:
          = 'A':  all M columns of U are returned in array U:
          = 'S':  the first min(m,n) columns of U (the left singular
                  vectors) are returned in the array U;
          = 'O':  the first min(m,n) columns of U (the left singular
                  vectors) are overwritten on the array A;
          = 'N':  no columns of U (no left singular vectors) are
                  computed.
[in]JOBVT
          JOBVT is CHARACTER*1
          Specifies options for computing all or part of the matrix
          V**H:
          = 'A':  all N rows of V**H are returned in the array VT;
          = 'S':  the first min(m,n) rows of V**H (the right singular
                  vectors) are returned in the array VT;
          = 'O':  the first min(m,n) rows of V**H (the right singular
                  vectors) are overwritten on the array A;
          = 'N':  no rows of V**H (no right singular vectors) are
                  computed.

          JOBVT and JOBU cannot both be 'O'.
[in]M
          M is INTEGER
          The number of rows of the input matrix A.  M >= 0.
[in]N
          N is INTEGER
          The number of columns of the input matrix A.  N >= 0.
[in,out]A
          A is COMPLEX array, dimension (LDA,N)
          On entry, the M-by-N matrix A.
          On exit,
          if JOBU = 'O',  A is overwritten with the first min(m,n)
                          columns of U (the left singular vectors,
                          stored columnwise);
          if JOBVT = 'O', A is overwritten with the first min(m,n)
                          rows of V**H (the right singular vectors,
                          stored rowwise);
          if JOBU .ne. 'O' and JOBVT .ne. 'O', the contents of A
                          are destroyed.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,M).
[out]S
          S is REAL array, dimension (min(M,N))
          The singular values of A, sorted so that S(i) >= S(i+1).
[out]U
          U is COMPLEX array, dimension (LDU,UCOL)
          (LDU,M) if JOBU = 'A' or (LDU,min(M,N)) if JOBU = 'S'.
          If JOBU = 'A', U contains the M-by-M unitary matrix U;
          if JOBU = 'S', U contains the first min(m,n) columns of U
          (the left singular vectors, stored columnwise);
          if JOBU = 'N' or 'O', U is not referenced.
[in]LDU
          LDU is INTEGER
          The leading dimension of the array U.  LDU >= 1; if
          JOBU = 'S' or 'A', LDU >= M.
[out]VT
          VT is COMPLEX array, dimension (LDVT,N)
          If JOBVT = 'A', VT contains the N-by-N unitary matrix
          V**H;
          if JOBVT = 'S', VT contains the first min(m,n) rows of
          V**H (the right singular vectors, stored rowwise);
          if JOBVT = 'N' or 'O', VT is not referenced.
[in]LDVT
          LDVT is INTEGER
          The leading dimension of the array VT.  LDVT >= 1; if
          JOBVT = 'A', LDVT >= N; if JOBVT = 'S', LDVT >= min(M,N).
[out]WORK
          WORK is COMPLEX array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.
          LWORK >=  MAX(1,2*MIN(M,N)+MAX(M,N)).
          For good performance, LWORK should generally be larger.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[out]RWORK
          RWORK is REAL array, dimension (5*min(M,N))
          On exit, if INFO > 0, RWORK(1:MIN(M,N)-1) contains the
          unconverged superdiagonal elements of an upper bidiagonal
          matrix B whose diagonal is in S (not necessarily sorted).
          B satisfies A = U * B * VT, so it has the same singular
          values as A, and singular vectors related by U and VT.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit.
          < 0:  if INFO = -i, the i-th argument had an illegal value.
          > 0:  if CBDSQR did not converge, INFO specifies how many
                superdiagonals of an intermediate bidiagonal form B
                did not converge to zero. See the description of RWORK
                above for details.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
April 2012

Definition at line 216 of file cgesvd.f.

216 *
217 * -- LAPACK driver routine (version 3.7.0) --
218 * -- LAPACK is a software package provided by Univ. of Tennessee, --
219 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
220 * April 2012
221 *
222 * .. Scalar Arguments ..
223  CHARACTER JOBU, JOBVT
224  INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N
225 * ..
226 * .. Array Arguments ..
227  REAL RWORK( * ), S( * )
228  COMPLEX A( LDA, * ), U( LDU, * ), VT( LDVT, * ),
229  $ WORK( * )
230 * ..
231 *
232 * =====================================================================
233 *
234 * .. Parameters ..
235  COMPLEX CZERO, CONE
236  parameter( czero = ( 0.0e0, 0.0e0 ),
237  $ cone = ( 1.0e0, 0.0e0 ) )
238  REAL ZERO, ONE
239  parameter( zero = 0.0e0, one = 1.0e0 )
240 * ..
241 * .. Local Scalars ..
242  LOGICAL LQUERY, WNTUA, WNTUAS, WNTUN, WNTUO, WNTUS,
243  $ WNTVA, WNTVAS, WNTVN, WNTVO, WNTVS
244  INTEGER BLK, CHUNK, I, IE, IERR, IR, IRWORK, ISCL,
245  $ ITAU, ITAUP, ITAUQ, IU, IWORK, LDWRKR, LDWRKU,
246  $ MAXWRK, MINMN, MINWRK, MNTHR, NCU, NCVT, NRU,
247  $ NRVT, WRKBL
248  INTEGER LWORK_CGEQRF, LWORK_CUNGQR_N, LWORK_CUNGQR_M,
249  $ LWORK_CGEBRD, LWORK_CUNGBR_P, LWORK_CUNGBR_Q,
250  $ LWORK_CGELQF, LWORK_CUNGLQ_N, LWORK_CUNGLQ_M
251  REAL ANRM, BIGNUM, EPS, SMLNUM
252 * ..
253 * .. Local Arrays ..
254  REAL DUM( 1 )
255  COMPLEX CDUM( 1 )
256 * ..
257 * .. External Subroutines ..
258  EXTERNAL cbdsqr, cgebrd, cgelqf, cgemm, cgeqrf, clacpy,
259  $ clascl, claset, cungbr, cunglq, cungqr, cunmbr,
260  $ slascl, xerbla
261 * ..
262 * .. External Functions ..
263  LOGICAL LSAME
264  INTEGER ILAENV
265  REAL CLANGE, SLAMCH
266  EXTERNAL lsame, ilaenv, clange, slamch
267 * ..
268 * .. Intrinsic Functions ..
269  INTRINSIC max, min, sqrt
270 * ..
271 * .. Executable Statements ..
272 *
273 * Test the input arguments
274 *
275  info = 0
276  minmn = min( m, n )
277  wntua = lsame( jobu, 'A' )
278  wntus = lsame( jobu, 'S' )
279  wntuas = wntua .OR. wntus
280  wntuo = lsame( jobu, 'O' )
281  wntun = lsame( jobu, 'N' )
282  wntva = lsame( jobvt, 'A' )
283  wntvs = lsame( jobvt, 'S' )
284  wntvas = wntva .OR. wntvs
285  wntvo = lsame( jobvt, 'O' )
286  wntvn = lsame( jobvt, 'N' )
287  lquery = ( lwork.EQ.-1 )
288 *
289  IF( .NOT.( wntua .OR. wntus .OR. wntuo .OR. wntun ) ) THEN
290  info = -1
291  ELSE IF( .NOT.( wntva .OR. wntvs .OR. wntvo .OR. wntvn ) .OR.
292  $ ( wntvo .AND. wntuo ) ) THEN
293  info = -2
294  ELSE IF( m.LT.0 ) THEN
295  info = -3
296  ELSE IF( n.LT.0 ) THEN
297  info = -4
298  ELSE IF( lda.LT.max( 1, m ) ) THEN
299  info = -6
300  ELSE IF( ldu.LT.1 .OR. ( wntuas .AND. ldu.LT.m ) ) THEN
301  info = -9
302  ELSE IF( ldvt.LT.1 .OR. ( wntva .AND. ldvt.LT.n ) .OR.
303  $ ( wntvs .AND. ldvt.LT.minmn ) ) THEN
304  info = -11
305  END IF
306 *
307 * Compute workspace
308 * (Note: Comments in the code beginning "Workspace:" describe the
309 * minimal amount of workspace needed at that point in the code,
310 * as well as the preferred amount for good performance.
311 * CWorkspace refers to complex workspace, and RWorkspace to
312 * real workspace. NB refers to the optimal block size for the
313 * immediately following subroutine, as returned by ILAENV.)
314 *
315  IF( info.EQ.0 ) THEN
316  minwrk = 1
317  maxwrk = 1
318  IF( m.GE.n .AND. minmn.GT.0 ) THEN
319 *
320 * Space needed for ZBDSQR is BDSPAC = 5*N
321 *
322  mnthr = ilaenv( 6, 'CGESVD', jobu // jobvt, m, n, 0, 0 )
323 * Compute space needed for CGEQRF
324  CALL cgeqrf( m, n, a, lda, cdum(1), cdum(1), -1, ierr )
325  lwork_cgeqrf = int( cdum(1) )
326 * Compute space needed for CUNGQR
327  CALL cungqr( m, n, n, a, lda, cdum(1), cdum(1), -1, ierr )
328  lwork_cungqr_n = int( cdum(1) )
329  CALL cungqr( m, m, n, a, lda, cdum(1), cdum(1), -1, ierr )
330  lwork_cungqr_m = int( cdum(1) )
331 * Compute space needed for CGEBRD
332  CALL cgebrd( n, n, a, lda, s, dum(1), cdum(1),
333  $ cdum(1), cdum(1), -1, ierr )
334  lwork_cgebrd = int( cdum(1) )
335 * Compute space needed for CUNGBR
336  CALL cungbr( 'P', n, n, n, a, lda, cdum(1),
337  $ cdum(1), -1, ierr )
338  lwork_cungbr_p = int( cdum(1) )
339  CALL cungbr( 'Q', n, n, n, a, lda, cdum(1),
340  $ cdum(1), -1, ierr )
341  lwork_cungbr_q = int( cdum(1) )
342 *
343  mnthr = ilaenv( 6, 'CGESVD', jobu // jobvt, m, n, 0, 0 )
344  IF( m.GE.mnthr ) THEN
345  IF( wntun ) THEN
346 *
347 * Path 1 (M much larger than N, JOBU='N')
348 *
349  maxwrk = n + lwork_cgeqrf
350  maxwrk = max( maxwrk, 2*n+lwork_cgebrd )
351  IF( wntvo .OR. wntvas )
352  $ maxwrk = max( maxwrk, 2*n+lwork_cungbr_p )
353  minwrk = 3*n
354  ELSE IF( wntuo .AND. wntvn ) THEN
355 *
356 * Path 2 (M much larger than N, JOBU='O', JOBVT='N')
357 *
358  wrkbl = n + lwork_cgeqrf
359  wrkbl = max( wrkbl, n+lwork_cungqr_n )
360  wrkbl = max( wrkbl, 2*n+lwork_cgebrd )
361  wrkbl = max( wrkbl, 2*n+lwork_cungbr_q )
362  maxwrk = max( n*n+wrkbl, n*n+m*n )
363  minwrk = 2*n + m
364  ELSE IF( wntuo .AND. wntvas ) THEN
365 *
366 * Path 3 (M much larger than N, JOBU='O', JOBVT='S' or
367 * 'A')
368 *
369  wrkbl = n + lwork_cgeqrf
370  wrkbl = max( wrkbl, n+lwork_cungqr_n )
371  wrkbl = max( wrkbl, 2*n+lwork_cgebrd )
372  wrkbl = max( wrkbl, 2*n+lwork_cungbr_q )
373  wrkbl = max( wrkbl, 2*n+lwork_cungbr_p )
374  maxwrk = max( n*n+wrkbl, n*n+m*n )
375  minwrk = 2*n + m
376  ELSE IF( wntus .AND. wntvn ) THEN
377 *
378 * Path 4 (M much larger than N, JOBU='S', JOBVT='N')
379 *
380  wrkbl = n + lwork_cgeqrf
381  wrkbl = max( wrkbl, n+lwork_cungqr_n )
382  wrkbl = max( wrkbl, 2*n+lwork_cgebrd )
383  wrkbl = max( wrkbl, 2*n+lwork_cungbr_q )
384  maxwrk = n*n + wrkbl
385  minwrk = 2*n + m
386  ELSE IF( wntus .AND. wntvo ) THEN
387 *
388 * Path 5 (M much larger than N, JOBU='S', JOBVT='O')
389 *
390  wrkbl = n + lwork_cgeqrf
391  wrkbl = max( wrkbl, n+lwork_cungqr_n )
392  wrkbl = max( wrkbl, 2*n+lwork_cgebrd )
393  wrkbl = max( wrkbl, 2*n+lwork_cungbr_q )
394  wrkbl = max( wrkbl, 2*n+lwork_cungbr_p )
395  maxwrk = 2*n*n + wrkbl
396  minwrk = 2*n + m
397  ELSE IF( wntus .AND. wntvas ) THEN
398 *
399 * Path 6 (M much larger than N, JOBU='S', JOBVT='S' or
400 * 'A')
401 *
402  wrkbl = n + lwork_cgeqrf
403  wrkbl = max( wrkbl, n+lwork_cungqr_n )
404  wrkbl = max( wrkbl, 2*n+lwork_cgebrd )
405  wrkbl = max( wrkbl, 2*n+lwork_cungbr_q )
406  wrkbl = max( wrkbl, 2*n+lwork_cungbr_p )
407  maxwrk = n*n + wrkbl
408  minwrk = 2*n + m
409  ELSE IF( wntua .AND. wntvn ) THEN
410 *
411 * Path 7 (M much larger than N, JOBU='A', JOBVT='N')
412 *
413  wrkbl = n + lwork_cgeqrf
414  wrkbl = max( wrkbl, n+lwork_cungqr_m )
415  wrkbl = max( wrkbl, 2*n+lwork_cgebrd )
416  wrkbl = max( wrkbl, 2*n+lwork_cungbr_q )
417  maxwrk = n*n + wrkbl
418  minwrk = 2*n + m
419  ELSE IF( wntua .AND. wntvo ) THEN
420 *
421 * Path 8 (M much larger than N, JOBU='A', JOBVT='O')
422 *
423  wrkbl = n + lwork_cgeqrf
424  wrkbl = max( wrkbl, n+lwork_cungqr_m )
425  wrkbl = max( wrkbl, 2*n+lwork_cgebrd )
426  wrkbl = max( wrkbl, 2*n+lwork_cungbr_q )
427  wrkbl = max( wrkbl, 2*n+lwork_cungbr_p )
428  maxwrk = 2*n*n + wrkbl
429  minwrk = 2*n + m
430  ELSE IF( wntua .AND. wntvas ) THEN
431 *
432 * Path 9 (M much larger than N, JOBU='A', JOBVT='S' or
433 * 'A')
434 *
435  wrkbl = n + lwork_cgeqrf
436  wrkbl = max( wrkbl, n+lwork_cungqr_m )
437  wrkbl = max( wrkbl, 2*n+lwork_cgebrd )
438  wrkbl = max( wrkbl, 2*n+lwork_cungbr_q )
439  wrkbl = max( wrkbl, 2*n+lwork_cungbr_p )
440  maxwrk = n*n + wrkbl
441  minwrk = 2*n + m
442  END IF
443  ELSE
444 *
445 * Path 10 (M at least N, but not much larger)
446 *
447  CALL cgebrd( m, n, a, lda, s, dum(1), cdum(1),
448  $ cdum(1), cdum(1), -1, ierr )
449  lwork_cgebrd = int( cdum(1) )
450  maxwrk = 2*n + lwork_cgebrd
451  IF( wntus .OR. wntuo ) THEN
452  CALL cungbr( 'Q', m, n, n, a, lda, cdum(1),
453  $ cdum(1), -1, ierr )
454  lwork_cungbr_q = int( cdum(1) )
455  maxwrk = max( maxwrk, 2*n+lwork_cungbr_q )
456  END IF
457  IF( wntua ) THEN
458  CALL cungbr( 'Q', m, m, n, a, lda, cdum(1),
459  $ cdum(1), -1, ierr )
460  lwork_cungbr_q = int( cdum(1) )
461  maxwrk = max( maxwrk, 2*n+lwork_cungbr_q )
462  END IF
463  IF( .NOT.wntvn ) THEN
464  maxwrk = max( maxwrk, 2*n+lwork_cungbr_p )
465  END IF
466  minwrk = 2*n + m
467  END IF
468  ELSE IF( minmn.GT.0 ) THEN
469 *
470 * Space needed for CBDSQR is BDSPAC = 5*M
471 *
472  mnthr = ilaenv( 6, 'CGESVD', jobu // jobvt, m, n, 0, 0 )
473 * Compute space needed for CGELQF
474  CALL cgelqf( m, n, a, lda, cdum(1), cdum(1), -1, ierr )
475  lwork_cgelqf = int( cdum(1) )
476 * Compute space needed for CUNGLQ
477  CALL cunglq( n, n, m, cdum(1), n, cdum(1), cdum(1), -1,
478  $ ierr )
479  lwork_cunglq_n = int( cdum(1) )
480  CALL cunglq( m, n, m, a, lda, cdum(1), cdum(1), -1, ierr )
481  lwork_cunglq_m = int( cdum(1) )
482 * Compute space needed for CGEBRD
483  CALL cgebrd( m, m, a, lda, s, dum(1), cdum(1),
484  $ cdum(1), cdum(1), -1, ierr )
485  lwork_cgebrd = int( cdum(1) )
486 * Compute space needed for CUNGBR P
487  CALL cungbr( 'P', m, m, m, a, n, cdum(1),
488  $ cdum(1), -1, ierr )
489  lwork_cungbr_p = int( cdum(1) )
490 * Compute space needed for CUNGBR Q
491  CALL cungbr( 'Q', m, m, m, a, n, cdum(1),
492  $ cdum(1), -1, ierr )
493  lwork_cungbr_q = int( cdum(1) )
494  IF( n.GE.mnthr ) THEN
495  IF( wntvn ) THEN
496 *
497 * Path 1t(N much larger than M, JOBVT='N')
498 *
499  maxwrk = m + lwork_cgelqf
500  maxwrk = max( maxwrk, 2*m+lwork_cgebrd )
501  IF( wntuo .OR. wntuas )
502  $ maxwrk = max( maxwrk, 2*m+lwork_cungbr_q )
503  minwrk = 3*m
504  ELSE IF( wntvo .AND. wntun ) THEN
505 *
506 * Path 2t(N much larger than M, JOBU='N', JOBVT='O')
507 *
508  wrkbl = m + lwork_cgelqf
509  wrkbl = max( wrkbl, m+lwork_cunglq_m )
510  wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
511  wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
512  maxwrk = max( m*m+wrkbl, m*m+m*n )
513  minwrk = 2*m + n
514  ELSE IF( wntvo .AND. wntuas ) THEN
515 *
516 * Path 3t(N much larger than M, JOBU='S' or 'A',
517 * JOBVT='O')
518 *
519  wrkbl = m + lwork_cgelqf
520  wrkbl = max( wrkbl, m+lwork_cunglq_m )
521  wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
522  wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
523  wrkbl = max( wrkbl, 2*m+lwork_cungbr_q )
524  maxwrk = max( m*m+wrkbl, m*m+m*n )
525  minwrk = 2*m + n
526  ELSE IF( wntvs .AND. wntun ) THEN
527 *
528 * Path 4t(N much larger than M, JOBU='N', JOBVT='S')
529 *
530  wrkbl = m + lwork_cgelqf
531  wrkbl = max( wrkbl, m+lwork_cunglq_m )
532  wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
533  wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
534  maxwrk = m*m + wrkbl
535  minwrk = 2*m + n
536  ELSE IF( wntvs .AND. wntuo ) THEN
537 *
538 * Path 5t(N much larger than M, JOBU='O', JOBVT='S')
539 *
540  wrkbl = m + lwork_cgelqf
541  wrkbl = max( wrkbl, m+lwork_cunglq_m )
542  wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
543  wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
544  wrkbl = max( wrkbl, 2*m+lwork_cungbr_q )
545  maxwrk = 2*m*m + wrkbl
546  minwrk = 2*m + n
547  ELSE IF( wntvs .AND. wntuas ) THEN
548 *
549 * Path 6t(N much larger than M, JOBU='S' or 'A',
550 * JOBVT='S')
551 *
552  wrkbl = m + lwork_cgelqf
553  wrkbl = max( wrkbl, m+lwork_cunglq_m )
554  wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
555  wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
556  wrkbl = max( wrkbl, 2*m+lwork_cungbr_q )
557  maxwrk = m*m + wrkbl
558  minwrk = 2*m + n
559  ELSE IF( wntva .AND. wntun ) THEN
560 *
561 * Path 7t(N much larger than M, JOBU='N', JOBVT='A')
562 *
563  wrkbl = m + lwork_cgelqf
564  wrkbl = max( wrkbl, m+lwork_cunglq_n )
565  wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
566  wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
567  maxwrk = m*m + wrkbl
568  minwrk = 2*m + n
569  ELSE IF( wntva .AND. wntuo ) THEN
570 *
571 * Path 8t(N much larger than M, JOBU='O', JOBVT='A')
572 *
573  wrkbl = m + lwork_cgelqf
574  wrkbl = max( wrkbl, m+lwork_cunglq_n )
575  wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
576  wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
577  wrkbl = max( wrkbl, 2*m+lwork_cungbr_q )
578  maxwrk = 2*m*m + wrkbl
579  minwrk = 2*m + n
580  ELSE IF( wntva .AND. wntuas ) THEN
581 *
582 * Path 9t(N much larger than M, JOBU='S' or 'A',
583 * JOBVT='A')
584 *
585  wrkbl = m + lwork_cgelqf
586  wrkbl = max( wrkbl, m+lwork_cunglq_n )
587  wrkbl = max( wrkbl, 2*m+lwork_cgebrd )
588  wrkbl = max( wrkbl, 2*m+lwork_cungbr_p )
589  wrkbl = max( wrkbl, 2*m+lwork_cungbr_q )
590  maxwrk = m*m + wrkbl
591  minwrk = 2*m + n
592  END IF
593  ELSE
594 *
595 * Path 10t(N greater than M, but not much larger)
596 *
597  CALL cgebrd( m, n, a, lda, s, dum(1), cdum(1),
598  $ cdum(1), cdum(1), -1, ierr )
599  lwork_cgebrd = int( cdum(1) )
600  maxwrk = 2*m + lwork_cgebrd
601  IF( wntvs .OR. wntvo ) THEN
602 * Compute space needed for CUNGBR P
603  CALL cungbr( 'P', m, n, m, a, n, cdum(1),
604  $ cdum(1), -1, ierr )
605  lwork_cungbr_p = int( cdum(1) )
606  maxwrk = max( maxwrk, 2*m+lwork_cungbr_p )
607  END IF
608  IF( wntva ) THEN
609  CALL cungbr( 'P', n, n, m, a, n, cdum(1),
610  $ cdum(1), -1, ierr )
611  lwork_cungbr_p = int( cdum(1) )
612  maxwrk = max( maxwrk, 2*m+lwork_cungbr_p )
613  END IF
614  IF( .NOT.wntun ) THEN
615  maxwrk = max( maxwrk, 2*m+lwork_cungbr_q )
616  END IF
617  minwrk = 2*m + n
618  END IF
619  END IF
620  maxwrk = max( minwrk, maxwrk )
621  work( 1 ) = maxwrk
622 *
623  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
624  info = -13
625  END IF
626  END IF
627 *
628  IF( info.NE.0 ) THEN
629  CALL xerbla( 'CGESVD', -info )
630  RETURN
631  ELSE IF( lquery ) THEN
632  RETURN
633  END IF
634 *
635 * Quick return if possible
636 *
637  IF( m.EQ.0 .OR. n.EQ.0 ) THEN
638  RETURN
639  END IF
640 *
641 * Get machine constants
642 *
643  eps = slamch( 'P' )
644  smlnum = sqrt( slamch( 'S' ) ) / eps
645  bignum = one / smlnum
646 *
647 * Scale A if max element outside range [SMLNUM,BIGNUM]
648 *
649  anrm = clange( 'M', m, n, a, lda, dum )
650  iscl = 0
651  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
652  iscl = 1
653  CALL clascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, ierr )
654  ELSE IF( anrm.GT.bignum ) THEN
655  iscl = 1
656  CALL clascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, ierr )
657  END IF
658 *
659  IF( m.GE.n ) THEN
660 *
661 * A has at least as many rows as columns. If A has sufficiently
662 * more rows than columns, first reduce using the QR
663 * decomposition (if sufficient workspace available)
664 *
665  IF( m.GE.mnthr ) THEN
666 *
667  IF( wntun ) THEN
668 *
669 * Path 1 (M much larger than N, JOBU='N')
670 * No left singular vectors to be computed
671 *
672  itau = 1
673  iwork = itau + n
674 *
675 * Compute A=Q*R
676 * (CWorkspace: need 2*N, prefer N+N*NB)
677 * (RWorkspace: need 0)
678 *
679  CALL cgeqrf( m, n, a, lda, work( itau ), work( iwork ),
680  $ lwork-iwork+1, ierr )
681 *
682 * Zero out below R
683 *
684  IF( n .GT. 1 ) THEN
685  CALL claset( 'L', n-1, n-1, czero, czero, a( 2, 1 ),
686  $ lda )
687  END IF
688  ie = 1
689  itauq = 1
690  itaup = itauq + n
691  iwork = itaup + n
692 *
693 * Bidiagonalize R in A
694 * (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
695 * (RWorkspace: need N)
696 *
697  CALL cgebrd( n, n, a, lda, s, rwork( ie ), work( itauq ),
698  $ work( itaup ), work( iwork ), lwork-iwork+1,
699  $ ierr )
700  ncvt = 0
701  IF( wntvo .OR. wntvas ) THEN
702 *
703 * If right singular vectors desired, generate P'.
704 * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
705 * (RWorkspace: 0)
706 *
707  CALL cungbr( 'P', n, n, n, a, lda, work( itaup ),
708  $ work( iwork ), lwork-iwork+1, ierr )
709  ncvt = n
710  END IF
711  irwork = ie + n
712 *
713 * Perform bidiagonal QR iteration, computing right
714 * singular vectors of A in A if desired
715 * (CWorkspace: 0)
716 * (RWorkspace: need BDSPAC)
717 *
718  CALL cbdsqr( 'U', n, ncvt, 0, 0, s, rwork( ie ), a, lda,
719  $ cdum, 1, cdum, 1, rwork( irwork ), info )
720 *
721 * If right singular vectors desired in VT, copy them there
722 *
723  IF( wntvas )
724  $ CALL clacpy( 'F', n, n, a, lda, vt, ldvt )
725 *
726  ELSE IF( wntuo .AND. wntvn ) THEN
727 *
728 * Path 2 (M much larger than N, JOBU='O', JOBVT='N')
729 * N left singular vectors to be overwritten on A and
730 * no right singular vectors to be computed
731 *
732  IF( lwork.GE.n*n+3*n ) THEN
733 *
734 * Sufficient workspace for a fast algorithm
735 *
736  ir = 1
737  IF( lwork.GE.max( wrkbl, lda*n )+lda*n ) THEN
738 *
739 * WORK(IU) is LDA by N, WORK(IR) is LDA by N
740 *
741  ldwrku = lda
742  ldwrkr = lda
743  ELSE IF( lwork.GE.max( wrkbl, lda*n )+n*n ) THEN
744 *
745 * WORK(IU) is LDA by N, WORK(IR) is N by N
746 *
747  ldwrku = lda
748  ldwrkr = n
749  ELSE
750 *
751 * WORK(IU) is LDWRKU by N, WORK(IR) is N by N
752 *
753  ldwrku = ( lwork-n*n ) / n
754  ldwrkr = n
755  END IF
756  itau = ir + ldwrkr*n
757  iwork = itau + n
758 *
759 * Compute A=Q*R
760 * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
761 * (RWorkspace: 0)
762 *
763  CALL cgeqrf( m, n, a, lda, work( itau ),
764  $ work( iwork ), lwork-iwork+1, ierr )
765 *
766 * Copy R to WORK(IR) and zero out below it
767 *
768  CALL clacpy( 'U', n, n, a, lda, work( ir ), ldwrkr )
769  CALL claset( 'L', n-1, n-1, czero, czero,
770  $ work( ir+1 ), ldwrkr )
771 *
772 * Generate Q in A
773 * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
774 * (RWorkspace: 0)
775 *
776  CALL cungqr( m, n, n, a, lda, work( itau ),
777  $ work( iwork ), lwork-iwork+1, ierr )
778  ie = 1
779  itauq = itau
780  itaup = itauq + n
781  iwork = itaup + n
782 *
783 * Bidiagonalize R in WORK(IR)
784 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
785 * (RWorkspace: need N)
786 *
787  CALL cgebrd( n, n, work( ir ), ldwrkr, s, rwork( ie ),
788  $ work( itauq ), work( itaup ),
789  $ work( iwork ), lwork-iwork+1, ierr )
790 *
791 * Generate left vectors bidiagonalizing R
792 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
793 * (RWorkspace: need 0)
794 *
795  CALL cungbr( 'Q', n, n, n, work( ir ), ldwrkr,
796  $ work( itauq ), work( iwork ),
797  $ lwork-iwork+1, ierr )
798  irwork = ie + n
799 *
800 * Perform bidiagonal QR iteration, computing left
801 * singular vectors of R in WORK(IR)
802 * (CWorkspace: need N*N)
803 * (RWorkspace: need BDSPAC)
804 *
805  CALL cbdsqr( 'U', n, 0, n, 0, s, rwork( ie ), cdum, 1,
806  $ work( ir ), ldwrkr, cdum, 1,
807  $ rwork( irwork ), info )
808  iu = itauq
809 *
810 * Multiply Q in A by left singular vectors of R in
811 * WORK(IR), storing result in WORK(IU) and copying to A
812 * (CWorkspace: need N*N+N, prefer N*N+M*N)
813 * (RWorkspace: 0)
814 *
815  DO 10 i = 1, m, ldwrku
816  chunk = min( m-i+1, ldwrku )
817  CALL cgemm( 'N', 'N', chunk, n, n, cone, a( i, 1 ),
818  $ lda, work( ir ), ldwrkr, czero,
819  $ work( iu ), ldwrku )
820  CALL clacpy( 'F', chunk, n, work( iu ), ldwrku,
821  $ a( i, 1 ), lda )
822  10 CONTINUE
823 *
824  ELSE
825 *
826 * Insufficient workspace for a fast algorithm
827 *
828  ie = 1
829  itauq = 1
830  itaup = itauq + n
831  iwork = itaup + n
832 *
833 * Bidiagonalize A
834 * (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB)
835 * (RWorkspace: N)
836 *
837  CALL cgebrd( m, n, a, lda, s, rwork( ie ),
838  $ work( itauq ), work( itaup ),
839  $ work( iwork ), lwork-iwork+1, ierr )
840 *
841 * Generate left vectors bidiagonalizing A
842 * (CWorkspace: need 3*N, prefer 2*N+N*NB)
843 * (RWorkspace: 0)
844 *
845  CALL cungbr( 'Q', m, n, n, a, lda, work( itauq ),
846  $ work( iwork ), lwork-iwork+1, ierr )
847  irwork = ie + n
848 *
849 * Perform bidiagonal QR iteration, computing left
850 * singular vectors of A in A
851 * (CWorkspace: need 0)
852 * (RWorkspace: need BDSPAC)
853 *
854  CALL cbdsqr( 'U', n, 0, m, 0, s, rwork( ie ), cdum, 1,
855  $ a, lda, cdum, 1, rwork( irwork ), info )
856 *
857  END IF
858 *
859  ELSE IF( wntuo .AND. wntvas ) THEN
860 *
861 * Path 3 (M much larger than N, JOBU='O', JOBVT='S' or 'A')
862 * N left singular vectors to be overwritten on A and
863 * N right singular vectors to be computed in VT
864 *
865  IF( lwork.GE.n*n+3*n ) THEN
866 *
867 * Sufficient workspace for a fast algorithm
868 *
869  ir = 1
870  IF( lwork.GE.max( wrkbl, lda*n )+lda*n ) THEN
871 *
872 * WORK(IU) is LDA by N and WORK(IR) is LDA by N
873 *
874  ldwrku = lda
875  ldwrkr = lda
876  ELSE IF( lwork.GE.max( wrkbl, lda*n )+n*n ) THEN
877 *
878 * WORK(IU) is LDA by N and WORK(IR) is N by N
879 *
880  ldwrku = lda
881  ldwrkr = n
882  ELSE
883 *
884 * WORK(IU) is LDWRKU by N and WORK(IR) is N by N
885 *
886  ldwrku = ( lwork-n*n ) / n
887  ldwrkr = n
888  END IF
889  itau = ir + ldwrkr*n
890  iwork = itau + n
891 *
892 * Compute A=Q*R
893 * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
894 * (RWorkspace: 0)
895 *
896  CALL cgeqrf( m, n, a, lda, work( itau ),
897  $ work( iwork ), lwork-iwork+1, ierr )
898 *
899 * Copy R to VT, zeroing out below it
900 *
901  CALL clacpy( 'U', n, n, a, lda, vt, ldvt )
902  IF( n.GT.1 )
903  $ CALL claset( 'L', n-1, n-1, czero, czero,
904  $ vt( 2, 1 ), ldvt )
905 *
906 * Generate Q in A
907 * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
908 * (RWorkspace: 0)
909 *
910  CALL cungqr( m, n, n, a, lda, work( itau ),
911  $ work( iwork ), lwork-iwork+1, ierr )
912  ie = 1
913  itauq = itau
914  itaup = itauq + n
915  iwork = itaup + n
916 *
917 * Bidiagonalize R in VT, copying result to WORK(IR)
918 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
919 * (RWorkspace: need N)
920 *
921  CALL cgebrd( n, n, vt, ldvt, s, rwork( ie ),
922  $ work( itauq ), work( itaup ),
923  $ work( iwork ), lwork-iwork+1, ierr )
924  CALL clacpy( 'L', n, n, vt, ldvt, work( ir ), ldwrkr )
925 *
926 * Generate left vectors bidiagonalizing R in WORK(IR)
927 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
928 * (RWorkspace: 0)
929 *
930  CALL cungbr( 'Q', n, n, n, work( ir ), ldwrkr,
931  $ work( itauq ), work( iwork ),
932  $ lwork-iwork+1, ierr )
933 *
934 * Generate right vectors bidiagonalizing R in VT
935 * (CWorkspace: need N*N+3*N-1, prefer N*N+2*N+(N-1)*NB)
936 * (RWorkspace: 0)
937 *
938  CALL cungbr( 'P', n, n, n, vt, ldvt, work( itaup ),
939  $ work( iwork ), lwork-iwork+1, ierr )
940  irwork = ie + n
941 *
942 * Perform bidiagonal QR iteration, computing left
943 * singular vectors of R in WORK(IR) and computing right
944 * singular vectors of R in VT
945 * (CWorkspace: need N*N)
946 * (RWorkspace: need BDSPAC)
947 *
948  CALL cbdsqr( 'U', n, n, n, 0, s, rwork( ie ), vt,
949  $ ldvt, work( ir ), ldwrkr, cdum, 1,
950  $ rwork( irwork ), info )
951  iu = itauq
952 *
953 * Multiply Q in A by left singular vectors of R in
954 * WORK(IR), storing result in WORK(IU) and copying to A
955 * (CWorkspace: need N*N+N, prefer N*N+M*N)
956 * (RWorkspace: 0)
957 *
958  DO 20 i = 1, m, ldwrku
959  chunk = min( m-i+1, ldwrku )
960  CALL cgemm( 'N', 'N', chunk, n, n, cone, a( i, 1 ),
961  $ lda, work( ir ), ldwrkr, czero,
962  $ work( iu ), ldwrku )
963  CALL clacpy( 'F', chunk, n, work( iu ), ldwrku,
964  $ a( i, 1 ), lda )
965  20 CONTINUE
966 *
967  ELSE
968 *
969 * Insufficient workspace for a fast algorithm
970 *
971  itau = 1
972  iwork = itau + n
973 *
974 * Compute A=Q*R
975 * (CWorkspace: need 2*N, prefer N+N*NB)
976 * (RWorkspace: 0)
977 *
978  CALL cgeqrf( m, n, a, lda, work( itau ),
979  $ work( iwork ), lwork-iwork+1, ierr )
980 *
981 * Copy R to VT, zeroing out below it
982 *
983  CALL clacpy( 'U', n, n, a, lda, vt, ldvt )
984  IF( n.GT.1 )
985  $ CALL claset( 'L', n-1, n-1, czero, czero,
986  $ vt( 2, 1 ), ldvt )
987 *
988 * Generate Q in A
989 * (CWorkspace: need 2*N, prefer N+N*NB)
990 * (RWorkspace: 0)
991 *
992  CALL cungqr( m, n, n, a, lda, work( itau ),
993  $ work( iwork ), lwork-iwork+1, ierr )
994  ie = 1
995  itauq = itau
996  itaup = itauq + n
997  iwork = itaup + n
998 *
999 * Bidiagonalize R in VT
1000 * (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
1001 * (RWorkspace: N)
1002 *
1003  CALL cgebrd( n, n, vt, ldvt, s, rwork( ie ),
1004  $ work( itauq ), work( itaup ),
1005  $ work( iwork ), lwork-iwork+1, ierr )
1006 *
1007 * Multiply Q in A by left vectors bidiagonalizing R
1008 * (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1009 * (RWorkspace: 0)
1010 *
1011  CALL cunmbr( 'Q', 'R', 'N', m, n, n, vt, ldvt,
1012  $ work( itauq ), a, lda, work( iwork ),
1013  $ lwork-iwork+1, ierr )
1014 *
1015 * Generate right vectors bidiagonalizing R in VT
1016 * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
1017 * (RWorkspace: 0)
1018 *
1019  CALL cungbr( 'P', n, n, n, vt, ldvt, work( itaup ),
1020  $ work( iwork ), lwork-iwork+1, ierr )
1021  irwork = ie + n
1022 *
1023 * Perform bidiagonal QR iteration, computing left
1024 * singular vectors of A in A and computing right
1025 * singular vectors of A in VT
1026 * (CWorkspace: 0)
1027 * (RWorkspace: need BDSPAC)
1028 *
1029  CALL cbdsqr( 'U', n, n, m, 0, s, rwork( ie ), vt,
1030  $ ldvt, a, lda, cdum, 1, rwork( irwork ),
1031  $ info )
1032 *
1033  END IF
1034 *
1035  ELSE IF( wntus ) THEN
1036 *
1037  IF( wntvn ) THEN
1038 *
1039 * Path 4 (M much larger than N, JOBU='S', JOBVT='N')
1040 * N left singular vectors to be computed in U and
1041 * no right singular vectors to be computed
1042 *
1043  IF( lwork.GE.n*n+3*n ) THEN
1044 *
1045 * Sufficient workspace for a fast algorithm
1046 *
1047  ir = 1
1048  IF( lwork.GE.wrkbl+lda*n ) THEN
1049 *
1050 * WORK(IR) is LDA by N
1051 *
1052  ldwrkr = lda
1053  ELSE
1054 *
1055 * WORK(IR) is N by N
1056 *
1057  ldwrkr = n
1058  END IF
1059  itau = ir + ldwrkr*n
1060  iwork = itau + n
1061 *
1062 * Compute A=Q*R
1063 * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1064 * (RWorkspace: 0)
1065 *
1066  CALL cgeqrf( m, n, a, lda, work( itau ),
1067  $ work( iwork ), lwork-iwork+1, ierr )
1068 *
1069 * Copy R to WORK(IR), zeroing out below it
1070 *
1071  CALL clacpy( 'U', n, n, a, lda, work( ir ),
1072  $ ldwrkr )
1073  CALL claset( 'L', n-1, n-1, czero, czero,
1074  $ work( ir+1 ), ldwrkr )
1075 *
1076 * Generate Q in A
1077 * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1078 * (RWorkspace: 0)
1079 *
1080  CALL cungqr( m, n, n, a, lda, work( itau ),
1081  $ work( iwork ), lwork-iwork+1, ierr )
1082  ie = 1
1083  itauq = itau
1084  itaup = itauq + n
1085  iwork = itaup + n
1086 *
1087 * Bidiagonalize R in WORK(IR)
1088 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
1089 * (RWorkspace: need N)
1090 *
1091  CALL cgebrd( n, n, work( ir ), ldwrkr, s,
1092  $ rwork( ie ), work( itauq ),
1093  $ work( itaup ), work( iwork ),
1094  $ lwork-iwork+1, ierr )
1095 *
1096 * Generate left vectors bidiagonalizing R in WORK(IR)
1097 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
1098 * (RWorkspace: 0)
1099 *
1100  CALL cungbr( 'Q', n, n, n, work( ir ), ldwrkr,
1101  $ work( itauq ), work( iwork ),
1102  $ lwork-iwork+1, ierr )
1103  irwork = ie + n
1104 *
1105 * Perform bidiagonal QR iteration, computing left
1106 * singular vectors of R in WORK(IR)
1107 * (CWorkspace: need N*N)
1108 * (RWorkspace: need BDSPAC)
1109 *
1110  CALL cbdsqr( 'U', n, 0, n, 0, s, rwork( ie ), cdum,
1111  $ 1, work( ir ), ldwrkr, cdum, 1,
1112  $ rwork( irwork ), info )
1113 *
1114 * Multiply Q in A by left singular vectors of R in
1115 * WORK(IR), storing result in U
1116 * (CWorkspace: need N*N)
1117 * (RWorkspace: 0)
1118 *
1119  CALL cgemm( 'N', 'N', m, n, n, cone, a, lda,
1120  $ work( ir ), ldwrkr, czero, u, ldu )
1121 *
1122  ELSE
1123 *
1124 * Insufficient workspace for a fast algorithm
1125 *
1126  itau = 1
1127  iwork = itau + n
1128 *
1129 * Compute A=Q*R, copying result to U
1130 * (CWorkspace: need 2*N, prefer N+N*NB)
1131 * (RWorkspace: 0)
1132 *
1133  CALL cgeqrf( m, n, a, lda, work( itau ),
1134  $ work( iwork ), lwork-iwork+1, ierr )
1135  CALL clacpy( 'L', m, n, a, lda, u, ldu )
1136 *
1137 * Generate Q in U
1138 * (CWorkspace: need 2*N, prefer N+N*NB)
1139 * (RWorkspace: 0)
1140 *
1141  CALL cungqr( m, n, n, u, ldu, work( itau ),
1142  $ work( iwork ), lwork-iwork+1, ierr )
1143  ie = 1
1144  itauq = itau
1145  itaup = itauq + n
1146  iwork = itaup + n
1147 *
1148 * Zero out below R in A
1149 *
1150  IF( n .GT. 1 ) THEN
1151  CALL claset( 'L', n-1, n-1, czero, czero,
1152  $ a( 2, 1 ), lda )
1153  END IF
1154 *
1155 * Bidiagonalize R in A
1156 * (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
1157 * (RWorkspace: need N)
1158 *
1159  CALL cgebrd( n, n, a, lda, s, rwork( ie ),
1160  $ work( itauq ), work( itaup ),
1161  $ work( iwork ), lwork-iwork+1, ierr )
1162 *
1163 * Multiply Q in U by left vectors bidiagonalizing R
1164 * (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1165 * (RWorkspace: 0)
1166 *
1167  CALL cunmbr( 'Q', 'R', 'N', m, n, n, a, lda,
1168  $ work( itauq ), u, ldu, work( iwork ),
1169  $ lwork-iwork+1, ierr )
1170  irwork = ie + n
1171 *
1172 * Perform bidiagonal QR iteration, computing left
1173 * singular vectors of A in U
1174 * (CWorkspace: 0)
1175 * (RWorkspace: need BDSPAC)
1176 *
1177  CALL cbdsqr( 'U', n, 0, m, 0, s, rwork( ie ), cdum,
1178  $ 1, u, ldu, cdum, 1, rwork( irwork ),
1179  $ info )
1180 *
1181  END IF
1182 *
1183  ELSE IF( wntvo ) THEN
1184 *
1185 * Path 5 (M much larger than N, JOBU='S', JOBVT='O')
1186 * N left singular vectors to be computed in U and
1187 * N right singular vectors to be overwritten on A
1188 *
1189  IF( lwork.GE.2*n*n+3*n ) THEN
1190 *
1191 * Sufficient workspace for a fast algorithm
1192 *
1193  iu = 1
1194  IF( lwork.GE.wrkbl+2*lda*n ) THEN
1195 *
1196 * WORK(IU) is LDA by N and WORK(IR) is LDA by N
1197 *
1198  ldwrku = lda
1199  ir = iu + ldwrku*n
1200  ldwrkr = lda
1201  ELSE IF( lwork.GE.wrkbl+( lda+n )*n ) THEN
1202 *
1203 * WORK(IU) is LDA by N and WORK(IR) is N by N
1204 *
1205  ldwrku = lda
1206  ir = iu + ldwrku*n
1207  ldwrkr = n
1208  ELSE
1209 *
1210 * WORK(IU) is N by N and WORK(IR) is N by N
1211 *
1212  ldwrku = n
1213  ir = iu + ldwrku*n
1214  ldwrkr = n
1215  END IF
1216  itau = ir + ldwrkr*n
1217  iwork = itau + n
1218 *
1219 * Compute A=Q*R
1220 * (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
1221 * (RWorkspace: 0)
1222 *
1223  CALL cgeqrf( m, n, a, lda, work( itau ),
1224  $ work( iwork ), lwork-iwork+1, ierr )
1225 *
1226 * Copy R to WORK(IU), zeroing out below it
1227 *
1228  CALL clacpy( 'U', n, n, a, lda, work( iu ),
1229  $ ldwrku )
1230  CALL claset( 'L', n-1, n-1, czero, czero,
1231  $ work( iu+1 ), ldwrku )
1232 *
1233 * Generate Q in A
1234 * (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
1235 * (RWorkspace: 0)
1236 *
1237  CALL cungqr( m, n, n, a, lda, work( itau ),
1238  $ work( iwork ), lwork-iwork+1, ierr )
1239  ie = 1
1240  itauq = itau
1241  itaup = itauq + n
1242  iwork = itaup + n
1243 *
1244 * Bidiagonalize R in WORK(IU), copying result to
1245 * WORK(IR)
1246 * (CWorkspace: need 2*N*N+3*N,
1247 * prefer 2*N*N+2*N+2*N*NB)
1248 * (RWorkspace: need N)
1249 *
1250  CALL cgebrd( n, n, work( iu ), ldwrku, s,
1251  $ rwork( ie ), work( itauq ),
1252  $ work( itaup ), work( iwork ),
1253  $ lwork-iwork+1, ierr )
1254  CALL clacpy( 'U', n, n, work( iu ), ldwrku,
1255  $ work( ir ), ldwrkr )
1256 *
1257 * Generate left bidiagonalizing vectors in WORK(IU)
1258 * (CWorkspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB)
1259 * (RWorkspace: 0)
1260 *
1261  CALL cungbr( 'Q', n, n, n, work( iu ), ldwrku,
1262  $ work( itauq ), work( iwork ),
1263  $ lwork-iwork+1, ierr )
1264 *
1265 * Generate right bidiagonalizing vectors in WORK(IR)
1266 * (CWorkspace: need 2*N*N+3*N-1,
1267 * prefer 2*N*N+2*N+(N-1)*NB)
1268 * (RWorkspace: 0)
1269 *
1270  CALL cungbr( 'P', n, n, n, work( ir ), ldwrkr,
1271  $ work( itaup ), work( iwork ),
1272  $ lwork-iwork+1, ierr )
1273  irwork = ie + n
1274 *
1275 * Perform bidiagonal QR iteration, computing left
1276 * singular vectors of R in WORK(IU) and computing
1277 * right singular vectors of R in WORK(IR)
1278 * (CWorkspace: need 2*N*N)
1279 * (RWorkspace: need BDSPAC)
1280 *
1281  CALL cbdsqr( 'U', n, n, n, 0, s, rwork( ie ),
1282  $ work( ir ), ldwrkr, work( iu ),
1283  $ ldwrku, cdum, 1, rwork( irwork ),
1284  $ info )
1285 *
1286 * Multiply Q in A by left singular vectors of R in
1287 * WORK(IU), storing result in U
1288 * (CWorkspace: need N*N)
1289 * (RWorkspace: 0)
1290 *
1291  CALL cgemm( 'N', 'N', m, n, n, cone, a, lda,
1292  $ work( iu ), ldwrku, czero, u, ldu )
1293 *
1294 * Copy right singular vectors of R to A
1295 * (CWorkspace: need N*N)
1296 * (RWorkspace: 0)
1297 *
1298  CALL clacpy( 'F', n, n, work( ir ), ldwrkr, a,
1299  $ lda )
1300 *
1301  ELSE
1302 *
1303 * Insufficient workspace for a fast algorithm
1304 *
1305  itau = 1
1306  iwork = itau + n
1307 *
1308 * Compute A=Q*R, copying result to U
1309 * (CWorkspace: need 2*N, prefer N+N*NB)
1310 * (RWorkspace: 0)
1311 *
1312  CALL cgeqrf( m, n, a, lda, work( itau ),
1313  $ work( iwork ), lwork-iwork+1, ierr )
1314  CALL clacpy( 'L', m, n, a, lda, u, ldu )
1315 *
1316 * Generate Q in U
1317 * (CWorkspace: need 2*N, prefer N+N*NB)
1318 * (RWorkspace: 0)
1319 *
1320  CALL cungqr( m, n, n, u, ldu, work( itau ),
1321  $ work( iwork ), lwork-iwork+1, ierr )
1322  ie = 1
1323  itauq = itau
1324  itaup = itauq + n
1325  iwork = itaup + n
1326 *
1327 * Zero out below R in A
1328 *
1329  IF( n .GT. 1 ) THEN
1330  CALL claset( 'L', n-1, n-1, czero, czero,
1331  $ a( 2, 1 ), lda )
1332  END IF
1333 *
1334 * Bidiagonalize R in A
1335 * (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
1336 * (RWorkspace: need N)
1337 *
1338  CALL cgebrd( n, n, a, lda, s, rwork( ie ),
1339  $ work( itauq ), work( itaup ),
1340  $ work( iwork ), lwork-iwork+1, ierr )
1341 *
1342 * Multiply Q in U by left vectors bidiagonalizing R
1343 * (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1344 * (RWorkspace: 0)
1345 *
1346  CALL cunmbr( 'Q', 'R', 'N', m, n, n, a, lda,
1347  $ work( itauq ), u, ldu, work( iwork ),
1348  $ lwork-iwork+1, ierr )
1349 *
1350 * Generate right vectors bidiagonalizing R in A
1351 * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
1352 * (RWorkspace: 0)
1353 *
1354  CALL cungbr( 'P', n, n, n, a, lda, work( itaup ),
1355  $ work( iwork ), lwork-iwork+1, ierr )
1356  irwork = ie + n
1357 *
1358 * Perform bidiagonal QR iteration, computing left
1359 * singular vectors of A in U and computing right
1360 * singular vectors of A in A
1361 * (CWorkspace: 0)
1362 * (RWorkspace: need BDSPAC)
1363 *
1364  CALL cbdsqr( 'U', n, n, m, 0, s, rwork( ie ), a,
1365  $ lda, u, ldu, cdum, 1, rwork( irwork ),
1366  $ info )
1367 *
1368  END IF
1369 *
1370  ELSE IF( wntvas ) THEN
1371 *
1372 * Path 6 (M much larger than N, JOBU='S', JOBVT='S'
1373 * or 'A')
1374 * N left singular vectors to be computed in U and
1375 * N right singular vectors to be computed in VT
1376 *
1377  IF( lwork.GE.n*n+3*n ) THEN
1378 *
1379 * Sufficient workspace for a fast algorithm
1380 *
1381  iu = 1
1382  IF( lwork.GE.wrkbl+lda*n ) THEN
1383 *
1384 * WORK(IU) is LDA by N
1385 *
1386  ldwrku = lda
1387  ELSE
1388 *
1389 * WORK(IU) is N by N
1390 *
1391  ldwrku = n
1392  END IF
1393  itau = iu + ldwrku*n
1394  iwork = itau + n
1395 *
1396 * Compute A=Q*R
1397 * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1398 * (RWorkspace: 0)
1399 *
1400  CALL cgeqrf( m, n, a, lda, work( itau ),
1401  $ work( iwork ), lwork-iwork+1, ierr )
1402 *
1403 * Copy R to WORK(IU), zeroing out below it
1404 *
1405  CALL clacpy( 'U', n, n, a, lda, work( iu ),
1406  $ ldwrku )
1407  CALL claset( 'L', n-1, n-1, czero, czero,
1408  $ work( iu+1 ), ldwrku )
1409 *
1410 * Generate Q in A
1411 * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1412 * (RWorkspace: 0)
1413 *
1414  CALL cungqr( m, n, n, a, lda, work( itau ),
1415  $ work( iwork ), lwork-iwork+1, ierr )
1416  ie = 1
1417  itauq = itau
1418  itaup = itauq + n
1419  iwork = itaup + n
1420 *
1421 * Bidiagonalize R in WORK(IU), copying result to VT
1422 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
1423 * (RWorkspace: need N)
1424 *
1425  CALL cgebrd( n, n, work( iu ), ldwrku, s,
1426  $ rwork( ie ), work( itauq ),
1427  $ work( itaup ), work( iwork ),
1428  $ lwork-iwork+1, ierr )
1429  CALL clacpy( 'U', n, n, work( iu ), ldwrku, vt,
1430  $ ldvt )
1431 *
1432 * Generate left bidiagonalizing vectors in WORK(IU)
1433 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
1434 * (RWorkspace: 0)
1435 *
1436  CALL cungbr( 'Q', n, n, n, work( iu ), ldwrku,
1437  $ work( itauq ), work( iwork ),
1438  $ lwork-iwork+1, ierr )
1439 *
1440 * Generate right bidiagonalizing vectors in VT
1441 * (CWorkspace: need N*N+3*N-1,
1442 * prefer N*N+2*N+(N-1)*NB)
1443 * (RWorkspace: 0)
1444 *
1445  CALL cungbr( 'P', n, n, n, vt, ldvt, work( itaup ),
1446  $ work( iwork ), lwork-iwork+1, ierr )
1447  irwork = ie + n
1448 *
1449 * Perform bidiagonal QR iteration, computing left
1450 * singular vectors of R in WORK(IU) and computing
1451 * right singular vectors of R in VT
1452 * (CWorkspace: need N*N)
1453 * (RWorkspace: need BDSPAC)
1454 *
1455  CALL cbdsqr( 'U', n, n, n, 0, s, rwork( ie ), vt,
1456  $ ldvt, work( iu ), ldwrku, cdum, 1,
1457  $ rwork( irwork ), info )
1458 *
1459 * Multiply Q in A by left singular vectors of R in
1460 * WORK(IU), storing result in U
1461 * (CWorkspace: need N*N)
1462 * (RWorkspace: 0)
1463 *
1464  CALL cgemm( 'N', 'N', m, n, n, cone, a, lda,
1465  $ work( iu ), ldwrku, czero, u, ldu )
1466 *
1467  ELSE
1468 *
1469 * Insufficient workspace for a fast algorithm
1470 *
1471  itau = 1
1472  iwork = itau + n
1473 *
1474 * Compute A=Q*R, copying result to U
1475 * (CWorkspace: need 2*N, prefer N+N*NB)
1476 * (RWorkspace: 0)
1477 *
1478  CALL cgeqrf( m, n, a, lda, work( itau ),
1479  $ work( iwork ), lwork-iwork+1, ierr )
1480  CALL clacpy( 'L', m, n, a, lda, u, ldu )
1481 *
1482 * Generate Q in U
1483 * (CWorkspace: need 2*N, prefer N+N*NB)
1484 * (RWorkspace: 0)
1485 *
1486  CALL cungqr( m, n, n, u, ldu, work( itau ),
1487  $ work( iwork ), lwork-iwork+1, ierr )
1488 *
1489 * Copy R to VT, zeroing out below it
1490 *
1491  CALL clacpy( 'U', n, n, a, lda, vt, ldvt )
1492  IF( n.GT.1 )
1493  $ CALL claset( 'L', n-1, n-1, czero, czero,
1494  $ vt( 2, 1 ), ldvt )
1495  ie = 1
1496  itauq = itau
1497  itaup = itauq + n
1498  iwork = itaup + n
1499 *
1500 * Bidiagonalize R in VT
1501 * (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
1502 * (RWorkspace: need N)
1503 *
1504  CALL cgebrd( n, n, vt, ldvt, s, rwork( ie ),
1505  $ work( itauq ), work( itaup ),
1506  $ work( iwork ), lwork-iwork+1, ierr )
1507 *
1508 * Multiply Q in U by left bidiagonalizing vectors
1509 * in VT
1510 * (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1511 * (RWorkspace: 0)
1512 *
1513  CALL cunmbr( 'Q', 'R', 'N', m, n, n, vt, ldvt,
1514  $ work( itauq ), u, ldu, work( iwork ),
1515  $ lwork-iwork+1, ierr )
1516 *
1517 * Generate right bidiagonalizing vectors in VT
1518 * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
1519 * (RWorkspace: 0)
1520 *
1521  CALL cungbr( 'P', n, n, n, vt, ldvt, work( itaup ),
1522  $ work( iwork ), lwork-iwork+1, ierr )
1523  irwork = ie + n
1524 *
1525 * Perform bidiagonal QR iteration, computing left
1526 * singular vectors of A in U and computing right
1527 * singular vectors of A in VT
1528 * (CWorkspace: 0)
1529 * (RWorkspace: need BDSPAC)
1530 *
1531  CALL cbdsqr( 'U', n, n, m, 0, s, rwork( ie ), vt,
1532  $ ldvt, u, ldu, cdum, 1,
1533  $ rwork( irwork ), info )
1534 *
1535  END IF
1536 *
1537  END IF
1538 *
1539  ELSE IF( wntua ) THEN
1540 *
1541  IF( wntvn ) THEN
1542 *
1543 * Path 7 (M much larger than N, JOBU='A', JOBVT='N')
1544 * M left singular vectors to be computed in U and
1545 * no right singular vectors to be computed
1546 *
1547  IF( lwork.GE.n*n+max( n+m, 3*n ) ) THEN
1548 *
1549 * Sufficient workspace for a fast algorithm
1550 *
1551  ir = 1
1552  IF( lwork.GE.wrkbl+lda*n ) THEN
1553 *
1554 * WORK(IR) is LDA by N
1555 *
1556  ldwrkr = lda
1557  ELSE
1558 *
1559 * WORK(IR) is N by N
1560 *
1561  ldwrkr = n
1562  END IF
1563  itau = ir + ldwrkr*n
1564  iwork = itau + n
1565 *
1566 * Compute A=Q*R, copying result to U
1567 * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1568 * (RWorkspace: 0)
1569 *
1570  CALL cgeqrf( m, n, a, lda, work( itau ),
1571  $ work( iwork ), lwork-iwork+1, ierr )
1572  CALL clacpy( 'L', m, n, a, lda, u, ldu )
1573 *
1574 * Copy R to WORK(IR), zeroing out below it
1575 *
1576  CALL clacpy( 'U', n, n, a, lda, work( ir ),
1577  $ ldwrkr )
1578  CALL claset( 'L', n-1, n-1, czero, czero,
1579  $ work( ir+1 ), ldwrkr )
1580 *
1581 * Generate Q in U
1582 * (CWorkspace: need N*N+N+M, prefer N*N+N+M*NB)
1583 * (RWorkspace: 0)
1584 *
1585  CALL cungqr( m, m, n, u, ldu, work( itau ),
1586  $ work( iwork ), lwork-iwork+1, ierr )
1587  ie = 1
1588  itauq = itau
1589  itaup = itauq + n
1590  iwork = itaup + n
1591 *
1592 * Bidiagonalize R in WORK(IR)
1593 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
1594 * (RWorkspace: need N)
1595 *
1596  CALL cgebrd( n, n, work( ir ), ldwrkr, s,
1597  $ rwork( ie ), work( itauq ),
1598  $ work( itaup ), work( iwork ),
1599  $ lwork-iwork+1, ierr )
1600 *
1601 * Generate left bidiagonalizing vectors in WORK(IR)
1602 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
1603 * (RWorkspace: 0)
1604 *
1605  CALL cungbr( 'Q', n, n, n, work( ir ), ldwrkr,
1606  $ work( itauq ), work( iwork ),
1607  $ lwork-iwork+1, ierr )
1608  irwork = ie + n
1609 *
1610 * Perform bidiagonal QR iteration, computing left
1611 * singular vectors of R in WORK(IR)
1612 * (CWorkspace: need N*N)
1613 * (RWorkspace: need BDSPAC)
1614 *
1615  CALL cbdsqr( 'U', n, 0, n, 0, s, rwork( ie ), cdum,
1616  $ 1, work( ir ), ldwrkr, cdum, 1,
1617  $ rwork( irwork ), info )
1618 *
1619 * Multiply Q in U by left singular vectors of R in
1620 * WORK(IR), storing result in A
1621 * (CWorkspace: need N*N)
1622 * (RWorkspace: 0)
1623 *
1624  CALL cgemm( 'N', 'N', m, n, n, cone, u, ldu,
1625  $ work( ir ), ldwrkr, czero, a, lda )
1626 *
1627 * Copy left singular vectors of A from A to U
1628 *
1629  CALL clacpy( 'F', m, n, a, lda, u, ldu )
1630 *
1631  ELSE
1632 *
1633 * Insufficient workspace for a fast algorithm
1634 *
1635  itau = 1
1636  iwork = itau + n
1637 *
1638 * Compute A=Q*R, copying result to U
1639 * (CWorkspace: need 2*N, prefer N+N*NB)
1640 * (RWorkspace: 0)
1641 *
1642  CALL cgeqrf( m, n, a, lda, work( itau ),
1643  $ work( iwork ), lwork-iwork+1, ierr )
1644  CALL clacpy( 'L', m, n, a, lda, u, ldu )
1645 *
1646 * Generate Q in U
1647 * (CWorkspace: need N+M, prefer N+M*NB)
1648 * (RWorkspace: 0)
1649 *
1650  CALL cungqr( m, m, n, u, ldu, work( itau ),
1651  $ work( iwork ), lwork-iwork+1, ierr )
1652  ie = 1
1653  itauq = itau
1654  itaup = itauq + n
1655  iwork = itaup + n
1656 *
1657 * Zero out below R in A
1658 *
1659  IF( n .GT. 1 ) THEN
1660  CALL claset( 'L', n-1, n-1, czero, czero,
1661  $ a( 2, 1 ), lda )
1662  END IF
1663 *
1664 * Bidiagonalize R in A
1665 * (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
1666 * (RWorkspace: need N)
1667 *
1668  CALL cgebrd( n, n, a, lda, s, rwork( ie ),
1669  $ work( itauq ), work( itaup ),
1670  $ work( iwork ), lwork-iwork+1, ierr )
1671 *
1672 * Multiply Q in U by left bidiagonalizing vectors
1673 * in A
1674 * (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1675 * (RWorkspace: 0)
1676 *
1677  CALL cunmbr( 'Q', 'R', 'N', m, n, n, a, lda,
1678  $ work( itauq ), u, ldu, work( iwork ),
1679  $ lwork-iwork+1, ierr )
1680  irwork = ie + n
1681 *
1682 * Perform bidiagonal QR iteration, computing left
1683 * singular vectors of A in U
1684 * (CWorkspace: 0)
1685 * (RWorkspace: need BDSPAC)
1686 *
1687  CALL cbdsqr( 'U', n, 0, m, 0, s, rwork( ie ), cdum,
1688  $ 1, u, ldu, cdum, 1, rwork( irwork ),
1689  $ info )
1690 *
1691  END IF
1692 *
1693  ELSE IF( wntvo ) THEN
1694 *
1695 * Path 8 (M much larger than N, JOBU='A', JOBVT='O')
1696 * M left singular vectors to be computed in U and
1697 * N right singular vectors to be overwritten on A
1698 *
1699  IF( lwork.GE.2*n*n+max( n+m, 3*n ) ) THEN
1700 *
1701 * Sufficient workspace for a fast algorithm
1702 *
1703  iu = 1
1704  IF( lwork.GE.wrkbl+2*lda*n ) THEN
1705 *
1706 * WORK(IU) is LDA by N and WORK(IR) is LDA by N
1707 *
1708  ldwrku = lda
1709  ir = iu + ldwrku*n
1710  ldwrkr = lda
1711  ELSE IF( lwork.GE.wrkbl+( lda+n )*n ) THEN
1712 *
1713 * WORK(IU) is LDA by N and WORK(IR) is N by N
1714 *
1715  ldwrku = lda
1716  ir = iu + ldwrku*n
1717  ldwrkr = n
1718  ELSE
1719 *
1720 * WORK(IU) is N by N and WORK(IR) is N by N
1721 *
1722  ldwrku = n
1723  ir = iu + ldwrku*n
1724  ldwrkr = n
1725  END IF
1726  itau = ir + ldwrkr*n
1727  iwork = itau + n
1728 *
1729 * Compute A=Q*R, copying result to U
1730 * (CWorkspace: need 2*N*N+2*N, prefer 2*N*N+N+N*NB)
1731 * (RWorkspace: 0)
1732 *
1733  CALL cgeqrf( m, n, a, lda, work( itau ),
1734  $ work( iwork ), lwork-iwork+1, ierr )
1735  CALL clacpy( 'L', m, n, a, lda, u, ldu )
1736 *
1737 * Generate Q in U
1738 * (CWorkspace: need 2*N*N+N+M, prefer 2*N*N+N+M*NB)
1739 * (RWorkspace: 0)
1740 *
1741  CALL cungqr( m, m, n, u, ldu, work( itau ),
1742  $ work( iwork ), lwork-iwork+1, ierr )
1743 *
1744 * Copy R to WORK(IU), zeroing out below it
1745 *
1746  CALL clacpy( 'U', n, n, a, lda, work( iu ),
1747  $ ldwrku )
1748  CALL claset( 'L', n-1, n-1, czero, czero,
1749  $ work( iu+1 ), ldwrku )
1750  ie = 1
1751  itauq = itau
1752  itaup = itauq + n
1753  iwork = itaup + n
1754 *
1755 * Bidiagonalize R in WORK(IU), copying result to
1756 * WORK(IR)
1757 * (CWorkspace: need 2*N*N+3*N,
1758 * prefer 2*N*N+2*N+2*N*NB)
1759 * (RWorkspace: need N)
1760 *
1761  CALL cgebrd( n, n, work( iu ), ldwrku, s,
1762  $ rwork( ie ), work( itauq ),
1763  $ work( itaup ), work( iwork ),
1764  $ lwork-iwork+1, ierr )
1765  CALL clacpy( 'U', n, n, work( iu ), ldwrku,
1766  $ work( ir ), ldwrkr )
1767 *
1768 * Generate left bidiagonalizing vectors in WORK(IU)
1769 * (CWorkspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB)
1770 * (RWorkspace: 0)
1771 *
1772  CALL cungbr( 'Q', n, n, n, work( iu ), ldwrku,
1773  $ work( itauq ), work( iwork ),
1774  $ lwork-iwork+1, ierr )
1775 *
1776 * Generate right bidiagonalizing vectors in WORK(IR)
1777 * (CWorkspace: need 2*N*N+3*N-1,
1778 * prefer 2*N*N+2*N+(N-1)*NB)
1779 * (RWorkspace: 0)
1780 *
1781  CALL cungbr( 'P', n, n, n, work( ir ), ldwrkr,
1782  $ work( itaup ), work( iwork ),
1783  $ lwork-iwork+1, ierr )
1784  irwork = ie + n
1785 *
1786 * Perform bidiagonal QR iteration, computing left
1787 * singular vectors of R in WORK(IU) and computing
1788 * right singular vectors of R in WORK(IR)
1789 * (CWorkspace: need 2*N*N)
1790 * (RWorkspace: need BDSPAC)
1791 *
1792  CALL cbdsqr( 'U', n, n, n, 0, s, rwork( ie ),
1793  $ work( ir ), ldwrkr, work( iu ),
1794  $ ldwrku, cdum, 1, rwork( irwork ),
1795  $ info )
1796 *
1797 * Multiply Q in U by left singular vectors of R in
1798 * WORK(IU), storing result in A
1799 * (CWorkspace: need N*N)
1800 * (RWorkspace: 0)
1801 *
1802  CALL cgemm( 'N', 'N', m, n, n, cone, u, ldu,
1803  $ work( iu ), ldwrku, czero, a, lda )
1804 *
1805 * Copy left singular vectors of A from A to U
1806 *
1807  CALL clacpy( 'F', m, n, a, lda, u, ldu )
1808 *
1809 * Copy right singular vectors of R from WORK(IR) to A
1810 *
1811  CALL clacpy( 'F', n, n, work( ir ), ldwrkr, a,
1812  $ lda )
1813 *
1814  ELSE
1815 *
1816 * Insufficient workspace for a fast algorithm
1817 *
1818  itau = 1
1819  iwork = itau + n
1820 *
1821 * Compute A=Q*R, copying result to U
1822 * (CWorkspace: need 2*N, prefer N+N*NB)
1823 * (RWorkspace: 0)
1824 *
1825  CALL cgeqrf( m, n, a, lda, work( itau ),
1826  $ work( iwork ), lwork-iwork+1, ierr )
1827  CALL clacpy( 'L', m, n, a, lda, u, ldu )
1828 *
1829 * Generate Q in U
1830 * (CWorkspace: need N+M, prefer N+M*NB)
1831 * (RWorkspace: 0)
1832 *
1833  CALL cungqr( m, m, n, u, ldu, work( itau ),
1834  $ work( iwork ), lwork-iwork+1, ierr )
1835  ie = 1
1836  itauq = itau
1837  itaup = itauq + n
1838  iwork = itaup + n
1839 *
1840 * Zero out below R in A
1841 *
1842  IF( n .GT. 1 ) THEN
1843  CALL claset( 'L', n-1, n-1, czero, czero,
1844  $ a( 2, 1 ), lda )
1845  END IF
1846 *
1847 * Bidiagonalize R in A
1848 * (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
1849 * (RWorkspace: need N)
1850 *
1851  CALL cgebrd( n, n, a, lda, s, rwork( ie ),
1852  $ work( itauq ), work( itaup ),
1853  $ work( iwork ), lwork-iwork+1, ierr )
1854 *
1855 * Multiply Q in U by left bidiagonalizing vectors
1856 * in A
1857 * (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1858 * (RWorkspace: 0)
1859 *
1860  CALL cunmbr( 'Q', 'R', 'N', m, n, n, a, lda,
1861  $ work( itauq ), u, ldu, work( iwork ),
1862  $ lwork-iwork+1, ierr )
1863 *
1864 * Generate right bidiagonalizing vectors in A
1865 * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
1866 * (RWorkspace: 0)
1867 *
1868  CALL cungbr( 'P', n, n, n, a, lda, work( itaup ),
1869  $ work( iwork ), lwork-iwork+1, ierr )
1870  irwork = ie + n
1871 *
1872 * Perform bidiagonal QR iteration, computing left
1873 * singular vectors of A in U and computing right
1874 * singular vectors of A in A
1875 * (CWorkspace: 0)
1876 * (RWorkspace: need BDSPAC)
1877 *
1878  CALL cbdsqr( 'U', n, n, m, 0, s, rwork( ie ), a,
1879  $ lda, u, ldu, cdum, 1, rwork( irwork ),
1880  $ info )
1881 *
1882  END IF
1883 *
1884  ELSE IF( wntvas ) THEN
1885 *
1886 * Path 9 (M much larger than N, JOBU='A', JOBVT='S'
1887 * or 'A')
1888 * M left singular vectors to be computed in U and
1889 * N right singular vectors to be computed in VT
1890 *
1891  IF( lwork.GE.n*n+max( n+m, 3*n ) ) THEN
1892 *
1893 * Sufficient workspace for a fast algorithm
1894 *
1895  iu = 1
1896  IF( lwork.GE.wrkbl+lda*n ) THEN
1897 *
1898 * WORK(IU) is LDA by N
1899 *
1900  ldwrku = lda
1901  ELSE
1902 *
1903 * WORK(IU) is N by N
1904 *
1905  ldwrku = n
1906  END IF
1907  itau = iu + ldwrku*n
1908  iwork = itau + n
1909 *
1910 * Compute A=Q*R, copying result to U
1911 * (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
1912 * (RWorkspace: 0)
1913 *
1914  CALL cgeqrf( m, n, a, lda, work( itau ),
1915  $ work( iwork ), lwork-iwork+1, ierr )
1916  CALL clacpy( 'L', m, n, a, lda, u, ldu )
1917 *
1918 * Generate Q in U
1919 * (CWorkspace: need N*N+N+M, prefer N*N+N+M*NB)
1920 * (RWorkspace: 0)
1921 *
1922  CALL cungqr( m, m, n, u, ldu, work( itau ),
1923  $ work( iwork ), lwork-iwork+1, ierr )
1924 *
1925 * Copy R to WORK(IU), zeroing out below it
1926 *
1927  CALL clacpy( 'U', n, n, a, lda, work( iu ),
1928  $ ldwrku )
1929  CALL claset( 'L', n-1, n-1, czero, czero,
1930  $ work( iu+1 ), ldwrku )
1931  ie = 1
1932  itauq = itau
1933  itaup = itauq + n
1934  iwork = itaup + n
1935 *
1936 * Bidiagonalize R in WORK(IU), copying result to VT
1937 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
1938 * (RWorkspace: need N)
1939 *
1940  CALL cgebrd( n, n, work( iu ), ldwrku, s,
1941  $ rwork( ie ), work( itauq ),
1942  $ work( itaup ), work( iwork ),
1943  $ lwork-iwork+1, ierr )
1944  CALL clacpy( 'U', n, n, work( iu ), ldwrku, vt,
1945  $ ldvt )
1946 *
1947 * Generate left bidiagonalizing vectors in WORK(IU)
1948 * (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
1949 * (RWorkspace: 0)
1950 *
1951  CALL cungbr( 'Q', n, n, n, work( iu ), ldwrku,
1952  $ work( itauq ), work( iwork ),
1953  $ lwork-iwork+1, ierr )
1954 *
1955 * Generate right bidiagonalizing vectors in VT
1956 * (CWorkspace: need N*N+3*N-1,
1957 * prefer N*N+2*N+(N-1)*NB)
1958 * (RWorkspace: need 0)
1959 *
1960  CALL cungbr( 'P', n, n, n, vt, ldvt, work( itaup ),
1961  $ work( iwork ), lwork-iwork+1, ierr )
1962  irwork = ie + n
1963 *
1964 * Perform bidiagonal QR iteration, computing left
1965 * singular vectors of R in WORK(IU) and computing
1966 * right singular vectors of R in VT
1967 * (CWorkspace: need N*N)
1968 * (RWorkspace: need BDSPAC)
1969 *
1970  CALL cbdsqr( 'U', n, n, n, 0, s, rwork( ie ), vt,
1971  $ ldvt, work( iu ), ldwrku, cdum, 1,
1972  $ rwork( irwork ), info )
1973 *
1974 * Multiply Q in U by left singular vectors of R in
1975 * WORK(IU), storing result in A
1976 * (CWorkspace: need N*N)
1977 * (RWorkspace: 0)
1978 *
1979  CALL cgemm( 'N', 'N', m, n, n, cone, u, ldu,
1980  $ work( iu ), ldwrku, czero, a, lda )
1981 *
1982 * Copy left singular vectors of A from A to U
1983 *
1984  CALL clacpy( 'F', m, n, a, lda, u, ldu )
1985 *
1986  ELSE
1987 *
1988 * Insufficient workspace for a fast algorithm
1989 *
1990  itau = 1
1991  iwork = itau + n
1992 *
1993 * Compute A=Q*R, copying result to U
1994 * (CWorkspace: need 2*N, prefer N+N*NB)
1995 * (RWorkspace: 0)
1996 *
1997  CALL cgeqrf( m, n, a, lda, work( itau ),
1998  $ work( iwork ), lwork-iwork+1, ierr )
1999  CALL clacpy( 'L', m, n, a, lda, u, ldu )
2000 *
2001 * Generate Q in U
2002 * (CWorkspace: need N+M, prefer N+M*NB)
2003 * (RWorkspace: 0)
2004 *
2005  CALL cungqr( m, m, n, u, ldu, work( itau ),
2006  $ work( iwork ), lwork-iwork+1, ierr )
2007 *
2008 * Copy R from A to VT, zeroing out below it
2009 *
2010  CALL clacpy( 'U', n, n, a, lda, vt, ldvt )
2011  IF( n.GT.1 )
2012  $ CALL claset( 'L', n-1, n-1, czero, czero,
2013  $ vt( 2, 1 ), ldvt )
2014  ie = 1
2015  itauq = itau
2016  itaup = itauq + n
2017  iwork = itaup + n
2018 *
2019 * Bidiagonalize R in VT
2020 * (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
2021 * (RWorkspace: need N)
2022 *
2023  CALL cgebrd( n, n, vt, ldvt, s, rwork( ie ),
2024  $ work( itauq ), work( itaup ),
2025  $ work( iwork ), lwork-iwork+1, ierr )
2026 *
2027 * Multiply Q in U by left bidiagonalizing vectors
2028 * in VT
2029 * (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
2030 * (RWorkspace: 0)
2031 *
2032  CALL cunmbr( 'Q', 'R', 'N', m, n, n, vt, ldvt,
2033  $ work( itauq ), u, ldu, work( iwork ),
2034  $ lwork-iwork+1, ierr )
2035 *
2036 * Generate right bidiagonalizing vectors in VT
2037 * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
2038 * (RWorkspace: 0)
2039 *
2040  CALL cungbr( 'P', n, n, n, vt, ldvt, work( itaup ),
2041  $ work( iwork ), lwork-iwork+1, ierr )
2042  irwork = ie + n
2043 *
2044 * Perform bidiagonal QR iteration, computing left
2045 * singular vectors of A in U and computing right
2046 * singular vectors of A in VT
2047 * (CWorkspace: 0)
2048 * (RWorkspace: need BDSPAC)
2049 *
2050  CALL cbdsqr( 'U', n, n, m, 0, s, rwork( ie ), vt,
2051  $ ldvt, u, ldu, cdum, 1,
2052  $ rwork( irwork ), info )
2053 *
2054  END IF
2055 *
2056  END IF
2057 *
2058  END IF
2059 *
2060  ELSE
2061 *
2062 * M .LT. MNTHR
2063 *
2064 * Path 10 (M at least N, but not much larger)
2065 * Reduce to bidiagonal form without QR decomposition
2066 *
2067  ie = 1
2068  itauq = 1
2069  itaup = itauq + n
2070  iwork = itaup + n
2071 *
2072 * Bidiagonalize A
2073 * (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB)
2074 * (RWorkspace: need N)
2075 *
2076  CALL cgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
2077  $ work( itaup ), work( iwork ), lwork-iwork+1,
2078  $ ierr )
2079  IF( wntuas ) THEN
2080 *
2081 * If left singular vectors desired in U, copy result to U
2082 * and generate left bidiagonalizing vectors in U
2083 * (CWorkspace: need 2*N+NCU, prefer 2*N+NCU*NB)
2084 * (RWorkspace: 0)
2085 *
2086  CALL clacpy( 'L', m, n, a, lda, u, ldu )
2087  IF( wntus )
2088  $ ncu = n
2089  IF( wntua )
2090  $ ncu = m
2091  CALL cungbr( 'Q', m, ncu, n, u, ldu, work( itauq ),
2092  $ work( iwork ), lwork-iwork+1, ierr )
2093  END IF
2094  IF( wntvas ) THEN
2095 *
2096 * If right singular vectors desired in VT, copy result to
2097 * VT and generate right bidiagonalizing vectors in VT
2098 * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
2099 * (RWorkspace: 0)
2100 *
2101  CALL clacpy( 'U', n, n, a, lda, vt, ldvt )
2102  CALL cungbr( 'P', n, n, n, vt, ldvt, work( itaup ),
2103  $ work( iwork ), lwork-iwork+1, ierr )
2104  END IF
2105  IF( wntuo ) THEN
2106 *
2107 * If left singular vectors desired in A, generate left
2108 * bidiagonalizing vectors in A
2109 * (CWorkspace: need 3*N, prefer 2*N+N*NB)
2110 * (RWorkspace: 0)
2111 *
2112  CALL cungbr( 'Q', m, n, n, a, lda, work( itauq ),
2113  $ work( iwork ), lwork-iwork+1, ierr )
2114  END IF
2115  IF( wntvo ) THEN
2116 *
2117 * If right singular vectors desired in A, generate right
2118 * bidiagonalizing vectors in A
2119 * (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
2120 * (RWorkspace: 0)
2121 *
2122  CALL cungbr( 'P', n, n, n, a, lda, work( itaup ),
2123  $ work( iwork ), lwork-iwork+1, ierr )
2124  END IF
2125  irwork = ie + n
2126  IF( wntuas .OR. wntuo )
2127  $ nru = m
2128  IF( wntun )
2129  $ nru = 0
2130  IF( wntvas .OR. wntvo )
2131  $ ncvt = n
2132  IF( wntvn )
2133  $ ncvt = 0
2134  IF( ( .NOT.wntuo ) .AND. ( .NOT.wntvo ) ) THEN
2135 *
2136 * Perform bidiagonal QR iteration, if desired, computing
2137 * left singular vectors in U and computing right singular
2138 * vectors in VT
2139 * (CWorkspace: 0)
2140 * (RWorkspace: need BDSPAC)
2141 *
2142  CALL cbdsqr( 'U', n, ncvt, nru, 0, s, rwork( ie ), vt,
2143  $ ldvt, u, ldu, cdum, 1, rwork( irwork ),
2144  $ info )
2145  ELSE IF( ( .NOT.wntuo ) .AND. wntvo ) THEN
2146 *
2147 * Perform bidiagonal QR iteration, if desired, computing
2148 * left singular vectors in U and computing right singular
2149 * vectors in A
2150 * (CWorkspace: 0)
2151 * (RWorkspace: need BDSPAC)
2152 *
2153  CALL cbdsqr( 'U', n, ncvt, nru, 0, s, rwork( ie ), a,
2154  $ lda, u, ldu, cdum, 1, rwork( irwork ),
2155  $ info )
2156  ELSE
2157 *
2158 * Perform bidiagonal QR iteration, if desired, computing
2159 * left singular vectors in A and computing right singular
2160 * vectors in VT
2161 * (CWorkspace: 0)
2162 * (RWorkspace: need BDSPAC)
2163 *
2164  CALL cbdsqr( 'U', n, ncvt, nru, 0, s, rwork( ie ), vt,
2165  $ ldvt, a, lda, cdum, 1, rwork( irwork ),
2166  $ info )
2167  END IF
2168 *
2169  END IF
2170 *
2171  ELSE
2172 *
2173 * A has more columns than rows. If A has sufficiently more
2174 * columns than rows, first reduce using the LQ decomposition (if
2175 * sufficient workspace available)
2176 *
2177  IF( n.GE.mnthr ) THEN
2178 *
2179  IF( wntvn ) THEN
2180 *
2181 * Path 1t(N much larger than M, JOBVT='N')
2182 * No right singular vectors to be computed
2183 *
2184  itau = 1
2185  iwork = itau + m
2186 *
2187 * Compute A=L*Q
2188 * (CWorkspace: need 2*M, prefer M+M*NB)
2189 * (RWorkspace: 0)
2190 *
2191  CALL cgelqf( m, n, a, lda, work( itau ), work( iwork ),
2192  $ lwork-iwork+1, ierr )
2193 *
2194 * Zero out above L
2195 *
2196  CALL claset( 'U', m-1, m-1, czero, czero, a( 1, 2 ),
2197  $ lda )
2198  ie = 1
2199  itauq = 1
2200  itaup = itauq + m
2201  iwork = itaup + m
2202 *
2203 * Bidiagonalize L in A
2204 * (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
2205 * (RWorkspace: need M)
2206 *
2207  CALL cgebrd( m, m, a, lda, s, rwork( ie ), work( itauq ),
2208  $ work( itaup ), work( iwork ), lwork-iwork+1,
2209  $ ierr )
2210  IF( wntuo .OR. wntuas ) THEN
2211 *
2212 * If left singular vectors desired, generate Q
2213 * (CWorkspace: need 3*M, prefer 2*M+M*NB)
2214 * (RWorkspace: 0)
2215 *
2216  CALL cungbr( 'Q', m, m, m, a, lda, work( itauq ),
2217  $ work( iwork ), lwork-iwork+1, ierr )
2218  END IF
2219  irwork = ie + m
2220  nru = 0
2221  IF( wntuo .OR. wntuas )
2222  $ nru = m
2223 *
2224 * Perform bidiagonal QR iteration, computing left singular
2225 * vectors of A in A if desired
2226 * (CWorkspace: 0)
2227 * (RWorkspace: need BDSPAC)
2228 *
2229  CALL cbdsqr( 'U', m, 0, nru, 0, s, rwork( ie ), cdum, 1,
2230  $ a, lda, cdum, 1, rwork( irwork ), info )
2231 *
2232 * If left singular vectors desired in U, copy them there
2233 *
2234  IF( wntuas )
2235  $ CALL clacpy( 'F', m, m, a, lda, u, ldu )
2236 *
2237  ELSE IF( wntvo .AND. wntun ) THEN
2238 *
2239 * Path 2t(N much larger than M, JOBU='N', JOBVT='O')
2240 * M right singular vectors to be overwritten on A and
2241 * no left singular vectors to be computed
2242 *
2243  IF( lwork.GE.m*m+3*m ) THEN
2244 *
2245 * Sufficient workspace for a fast algorithm
2246 *
2247  ir = 1
2248  IF( lwork.GE.max( wrkbl, lda*n )+lda*m ) THEN
2249 *
2250 * WORK(IU) is LDA by N and WORK(IR) is LDA by M
2251 *
2252  ldwrku = lda
2253  chunk = n
2254  ldwrkr = lda
2255  ELSE IF( lwork.GE.max( wrkbl, lda*n )+m*m ) THEN
2256 *
2257 * WORK(IU) is LDA by N and WORK(IR) is M by M
2258 *
2259  ldwrku = lda
2260  chunk = n
2261  ldwrkr = m
2262  ELSE
2263 *
2264 * WORK(IU) is M by CHUNK and WORK(IR) is M by M
2265 *
2266  ldwrku = m
2267  chunk = ( lwork-m*m ) / m
2268  ldwrkr = m
2269  END IF
2270  itau = ir + ldwrkr*m
2271  iwork = itau + m
2272 *
2273 * Compute A=L*Q
2274 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2275 * (RWorkspace: 0)
2276 *
2277  CALL cgelqf( m, n, a, lda, work( itau ),
2278  $ work( iwork ), lwork-iwork+1, ierr )
2279 *
2280 * Copy L to WORK(IR) and zero out above it
2281 *
2282  CALL clacpy( 'L', m, m, a, lda, work( ir ), ldwrkr )
2283  CALL claset( 'U', m-1, m-1, czero, czero,
2284  $ work( ir+ldwrkr ), ldwrkr )
2285 *
2286 * Generate Q in A
2287 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2288 * (RWorkspace: 0)
2289 *
2290  CALL cunglq( m, n, m, a, lda, work( itau ),
2291  $ work( iwork ), lwork-iwork+1, ierr )
2292  ie = 1
2293  itauq = itau
2294  itaup = itauq + m
2295  iwork = itaup + m
2296 *
2297 * Bidiagonalize L in WORK(IR)
2298 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
2299 * (RWorkspace: need M)
2300 *
2301  CALL cgebrd( m, m, work( ir ), ldwrkr, s, rwork( ie ),
2302  $ work( itauq ), work( itaup ),
2303  $ work( iwork ), lwork-iwork+1, ierr )
2304 *
2305 * Generate right vectors bidiagonalizing L
2306 * (CWorkspace: need M*M+3*M-1, prefer M*M+2*M+(M-1)*NB)
2307 * (RWorkspace: 0)
2308 *
2309  CALL cungbr( 'P', m, m, m, work( ir ), ldwrkr,
2310  $ work( itaup ), work( iwork ),
2311  $ lwork-iwork+1, ierr )
2312  irwork = ie + m
2313 *
2314 * Perform bidiagonal QR iteration, computing right
2315 * singular vectors of L in WORK(IR)
2316 * (CWorkspace: need M*M)
2317 * (RWorkspace: need BDSPAC)
2318 *
2319  CALL cbdsqr( 'U', m, m, 0, 0, s, rwork( ie ),
2320  $ work( ir ), ldwrkr, cdum, 1, cdum, 1,
2321  $ rwork( irwork ), info )
2322  iu = itauq
2323 *
2324 * Multiply right singular vectors of L in WORK(IR) by Q
2325 * in A, storing result in WORK(IU) and copying to A
2326 * (CWorkspace: need M*M+M, prefer M*M+M*N)
2327 * (RWorkspace: 0)
2328 *
2329  DO 30 i = 1, n, chunk
2330  blk = min( n-i+1, chunk )
2331  CALL cgemm( 'N', 'N', m, blk, m, cone, work( ir ),
2332  $ ldwrkr, a( 1, i ), lda, czero,
2333  $ work( iu ), ldwrku )
2334  CALL clacpy( 'F', m, blk, work( iu ), ldwrku,
2335  $ a( 1, i ), lda )
2336  30 CONTINUE
2337 *
2338  ELSE
2339 *
2340 * Insufficient workspace for a fast algorithm
2341 *
2342  ie = 1
2343  itauq = 1
2344  itaup = itauq + m
2345  iwork = itaup + m
2346 *
2347 * Bidiagonalize A
2348 * (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)
2349 * (RWorkspace: need M)
2350 *
2351  CALL cgebrd( m, n, a, lda, s, rwork( ie ),
2352  $ work( itauq ), work( itaup ),
2353  $ work( iwork ), lwork-iwork+1, ierr )
2354 *
2355 * Generate right vectors bidiagonalizing A
2356 * (CWorkspace: need 3*M, prefer 2*M+M*NB)
2357 * (RWorkspace: 0)
2358 *
2359  CALL cungbr( 'P', m, n, m, a, lda, work( itaup ),
2360  $ work( iwork ), lwork-iwork+1, ierr )
2361  irwork = ie + m
2362 *
2363 * Perform bidiagonal QR iteration, computing right
2364 * singular vectors of A in A
2365 * (CWorkspace: 0)
2366 * (RWorkspace: need BDSPAC)
2367 *
2368  CALL cbdsqr( 'L', m, n, 0, 0, s, rwork( ie ), a, lda,
2369  $ cdum, 1, cdum, 1, rwork( irwork ), info )
2370 *
2371  END IF
2372 *
2373  ELSE IF( wntvo .AND. wntuas ) THEN
2374 *
2375 * Path 3t(N much larger than M, JOBU='S' or 'A', JOBVT='O')
2376 * M right singular vectors to be overwritten on A and
2377 * M left singular vectors to be computed in U
2378 *
2379  IF( lwork.GE.m*m+3*m ) THEN
2380 *
2381 * Sufficient workspace for a fast algorithm
2382 *
2383  ir = 1
2384  IF( lwork.GE.max( wrkbl, lda*n )+lda*m ) THEN
2385 *
2386 * WORK(IU) is LDA by N and WORK(IR) is LDA by M
2387 *
2388  ldwrku = lda
2389  chunk = n
2390  ldwrkr = lda
2391  ELSE IF( lwork.GE.max( wrkbl, lda*n )+m*m ) THEN
2392 *
2393 * WORK(IU) is LDA by N and WORK(IR) is M by M
2394 *
2395  ldwrku = lda
2396  chunk = n
2397  ldwrkr = m
2398  ELSE
2399 *
2400 * WORK(IU) is M by CHUNK and WORK(IR) is M by M
2401 *
2402  ldwrku = m
2403  chunk = ( lwork-m*m ) / m
2404  ldwrkr = m
2405  END IF
2406  itau = ir + ldwrkr*m
2407  iwork = itau + m
2408 *
2409 * Compute A=L*Q
2410 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2411 * (RWorkspace: 0)
2412 *
2413  CALL cgelqf( m, n, a, lda, work( itau ),
2414  $ work( iwork ), lwork-iwork+1, ierr )
2415 *
2416 * Copy L to U, zeroing about above it
2417 *
2418  CALL clacpy( 'L', m, m, a, lda, u, ldu )
2419  CALL claset( 'U', m-1, m-1, czero, czero, u( 1, 2 ),
2420  $ ldu )
2421 *
2422 * Generate Q in A
2423 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2424 * (RWorkspace: 0)
2425 *
2426  CALL cunglq( m, n, m, a, lda, work( itau ),
2427  $ work( iwork ), lwork-iwork+1, ierr )
2428  ie = 1
2429  itauq = itau
2430  itaup = itauq + m
2431  iwork = itaup + m
2432 *
2433 * Bidiagonalize L in U, copying result to WORK(IR)
2434 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
2435 * (RWorkspace: need M)
2436 *
2437  CALL cgebrd( m, m, u, ldu, s, rwork( ie ),
2438  $ work( itauq ), work( itaup ),
2439  $ work( iwork ), lwork-iwork+1, ierr )
2440  CALL clacpy( 'U', m, m, u, ldu, work( ir ), ldwrkr )
2441 *
2442 * Generate right vectors bidiagonalizing L in WORK(IR)
2443 * (CWorkspace: need M*M+3*M-1, prefer M*M+2*M+(M-1)*NB)
2444 * (RWorkspace: 0)
2445 *
2446  CALL cungbr( 'P', m, m, m, work( ir ), ldwrkr,
2447  $ work( itaup ), work( iwork ),
2448  $ lwork-iwork+1, ierr )
2449 *
2450 * Generate left vectors bidiagonalizing L in U
2451 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
2452 * (RWorkspace: 0)
2453 *
2454  CALL cungbr( 'Q', m, m, m, u, ldu, work( itauq ),
2455  $ work( iwork ), lwork-iwork+1, ierr )
2456  irwork = ie + m
2457 *
2458 * Perform bidiagonal QR iteration, computing left
2459 * singular vectors of L in U, and computing right
2460 * singular vectors of L in WORK(IR)
2461 * (CWorkspace: need M*M)
2462 * (RWorkspace: need BDSPAC)
2463 *
2464  CALL cbdsqr( 'U', m, m, m, 0, s, rwork( ie ),
2465  $ work( ir ), ldwrkr, u, ldu, cdum, 1,
2466  $ rwork( irwork ), info )
2467  iu = itauq
2468 *
2469 * Multiply right singular vectors of L in WORK(IR) by Q
2470 * in A, storing result in WORK(IU) and copying to A
2471 * (CWorkspace: need M*M+M, prefer M*M+M*N))
2472 * (RWorkspace: 0)
2473 *
2474  DO 40 i = 1, n, chunk
2475  blk = min( n-i+1, chunk )
2476  CALL cgemm( 'N', 'N', m, blk, m, cone, work( ir ),
2477  $ ldwrkr, a( 1, i ), lda, czero,
2478  $ work( iu ), ldwrku )
2479  CALL clacpy( 'F', m, blk, work( iu ), ldwrku,
2480  $ a( 1, i ), lda )
2481  40 CONTINUE
2482 *
2483  ELSE
2484 *
2485 * Insufficient workspace for a fast algorithm
2486 *
2487  itau = 1
2488  iwork = itau + m
2489 *
2490 * Compute A=L*Q
2491 * (CWorkspace: need 2*M, prefer M+M*NB)
2492 * (RWorkspace: 0)
2493 *
2494  CALL cgelqf( m, n, a, lda, work( itau ),
2495  $ work( iwork ), lwork-iwork+1, ierr )
2496 *
2497 * Copy L to U, zeroing out above it
2498 *
2499  CALL clacpy( 'L', m, m, a, lda, u, ldu )
2500  CALL claset( 'U', m-1, m-1, czero, czero, u( 1, 2 ),
2501  $ ldu )
2502 *
2503 * Generate Q in A
2504 * (CWorkspace: need 2*M, prefer M+M*NB)
2505 * (RWorkspace: 0)
2506 *
2507  CALL cunglq( m, n, m, a, lda, work( itau ),
2508  $ work( iwork ), lwork-iwork+1, ierr )
2509  ie = 1
2510  itauq = itau
2511  itaup = itauq + m
2512  iwork = itaup + m
2513 *
2514 * Bidiagonalize L in U
2515 * (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
2516 * (RWorkspace: need M)
2517 *
2518  CALL cgebrd( m, m, u, ldu, s, rwork( ie ),
2519  $ work( itauq ), work( itaup ),
2520  $ work( iwork ), lwork-iwork+1, ierr )
2521 *
2522 * Multiply right vectors bidiagonalizing L by Q in A
2523 * (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
2524 * (RWorkspace: 0)
2525 *
2526  CALL cunmbr( 'P', 'L', 'C', m, n, m, u, ldu,
2527  $ work( itaup ), a, lda, work( iwork ),
2528  $ lwork-iwork+1, ierr )
2529 *
2530 * Generate left vectors bidiagonalizing L in U
2531 * (CWorkspace: need 3*M, prefer 2*M+M*NB)
2532 * (RWorkspace: 0)
2533 *
2534  CALL cungbr( 'Q', m, m, m, u, ldu, work( itauq ),
2535  $ work( iwork ), lwork-iwork+1, ierr )
2536  irwork = ie + m
2537 *
2538 * Perform bidiagonal QR iteration, computing left
2539 * singular vectors of A in U and computing right
2540 * singular vectors of A in A
2541 * (CWorkspace: 0)
2542 * (RWorkspace: need BDSPAC)
2543 *
2544  CALL cbdsqr( 'U', m, n, m, 0, s, rwork( ie ), a, lda,
2545  $ u, ldu, cdum, 1, rwork( irwork ), info )
2546 *
2547  END IF
2548 *
2549  ELSE IF( wntvs ) THEN
2550 *
2551  IF( wntun ) THEN
2552 *
2553 * Path 4t(N much larger than M, JOBU='N', JOBVT='S')
2554 * M right singular vectors to be computed in VT and
2555 * no left singular vectors to be computed
2556 *
2557  IF( lwork.GE.m*m+3*m ) THEN
2558 *
2559 * Sufficient workspace for a fast algorithm
2560 *
2561  ir = 1
2562  IF( lwork.GE.wrkbl+lda*m ) THEN
2563 *
2564 * WORK(IR) is LDA by M
2565 *
2566  ldwrkr = lda
2567  ELSE
2568 *
2569 * WORK(IR) is M by M
2570 *
2571  ldwrkr = m
2572  END IF
2573  itau = ir + ldwrkr*m
2574  iwork = itau + m
2575 *
2576 * Compute A=L*Q
2577 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2578 * (RWorkspace: 0)
2579 *
2580  CALL cgelqf( m, n, a, lda, work( itau ),
2581  $ work( iwork ), lwork-iwork+1, ierr )
2582 *
2583 * Copy L to WORK(IR), zeroing out above it
2584 *
2585  CALL clacpy( 'L', m, m, a, lda, work( ir ),
2586  $ ldwrkr )
2587  CALL claset( 'U', m-1, m-1, czero, czero,
2588  $ work( ir+ldwrkr ), ldwrkr )
2589 *
2590 * Generate Q in A
2591 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2592 * (RWorkspace: 0)
2593 *
2594  CALL cunglq( m, n, m, a, lda, work( itau ),
2595  $ work( iwork ), lwork-iwork+1, ierr )
2596  ie = 1
2597  itauq = itau
2598  itaup = itauq + m
2599  iwork = itaup + m
2600 *
2601 * Bidiagonalize L in WORK(IR)
2602 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
2603 * (RWorkspace: need M)
2604 *
2605  CALL cgebrd( m, m, work( ir ), ldwrkr, s,
2606  $ rwork( ie ), work( itauq ),
2607  $ work( itaup ), work( iwork ),
2608  $ lwork-iwork+1, ierr )
2609 *
2610 * Generate right vectors bidiagonalizing L in
2611 * WORK(IR)
2612 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+(M-1)*NB)
2613 * (RWorkspace: 0)
2614 *
2615  CALL cungbr( 'P', m, m, m, work( ir ), ldwrkr,
2616  $ work( itaup ), work( iwork ),
2617  $ lwork-iwork+1, ierr )
2618  irwork = ie + m
2619 *
2620 * Perform bidiagonal QR iteration, computing right
2621 * singular vectors of L in WORK(IR)
2622 * (CWorkspace: need M*M)
2623 * (RWorkspace: need BDSPAC)
2624 *
2625  CALL cbdsqr( 'U', m, m, 0, 0, s, rwork( ie ),
2626  $ work( ir ), ldwrkr, cdum, 1, cdum, 1,
2627  $ rwork( irwork ), info )
2628 *
2629 * Multiply right singular vectors of L in WORK(IR) by
2630 * Q in A, storing result in VT
2631 * (CWorkspace: need M*M)
2632 * (RWorkspace: 0)
2633 *
2634  CALL cgemm( 'N', 'N', m, n, m, cone, work( ir ),
2635  $ ldwrkr, a, lda, czero, vt, ldvt )
2636 *
2637  ELSE
2638 *
2639 * Insufficient workspace for a fast algorithm
2640 *
2641  itau = 1
2642  iwork = itau + m
2643 *
2644 * Compute A=L*Q
2645 * (CWorkspace: need 2*M, prefer M+M*NB)
2646 * (RWorkspace: 0)
2647 *
2648  CALL cgelqf( m, n, a, lda, work( itau ),
2649  $ work( iwork ), lwork-iwork+1, ierr )
2650 *
2651 * Copy result to VT
2652 *
2653  CALL clacpy( 'U', m, n, a, lda, vt, ldvt )
2654 *
2655 * Generate Q in VT
2656 * (CWorkspace: need 2*M, prefer M+M*NB)
2657 * (RWorkspace: 0)
2658 *
2659  CALL cunglq( m, n, m, vt, ldvt, work( itau ),
2660  $ work( iwork ), lwork-iwork+1, ierr )
2661  ie = 1
2662  itauq = itau
2663  itaup = itauq + m
2664  iwork = itaup + m
2665 *
2666 * Zero out above L in A
2667 *
2668  CALL claset( 'U', m-1, m-1, czero, czero,
2669  $ a( 1, 2 ), lda )
2670 *
2671 * Bidiagonalize L in A
2672 * (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
2673 * (RWorkspace: need M)
2674 *
2675  CALL cgebrd( m, m, a, lda, s, rwork( ie ),
2676  $ work( itauq ), work( itaup ),
2677  $ work( iwork ), lwork-iwork+1, ierr )
2678 *
2679 * Multiply right vectors bidiagonalizing L by Q in VT
2680 * (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
2681 * (RWorkspace: 0)
2682 *
2683  CALL cunmbr( 'P', 'L', 'C', m, n, m, a, lda,
2684  $ work( itaup ), vt, ldvt,
2685  $ work( iwork ), lwork-iwork+1, ierr )
2686  irwork = ie + m
2687 *
2688 * Perform bidiagonal QR iteration, computing right
2689 * singular vectors of A in VT
2690 * (CWorkspace: 0)
2691 * (RWorkspace: need BDSPAC)
2692 *
2693  CALL cbdsqr( 'U', m, n, 0, 0, s, rwork( ie ), vt,
2694  $ ldvt, cdum, 1, cdum, 1,
2695  $ rwork( irwork ), info )
2696 *
2697  END IF
2698 *
2699  ELSE IF( wntuo ) THEN
2700 *
2701 * Path 5t(N much larger than M, JOBU='O', JOBVT='S')
2702 * M right singular vectors to be computed in VT and
2703 * M left singular vectors to be overwritten on A
2704 *
2705  IF( lwork.GE.2*m*m+3*m ) THEN
2706 *
2707 * Sufficient workspace for a fast algorithm
2708 *
2709  iu = 1
2710  IF( lwork.GE.wrkbl+2*lda*m ) THEN
2711 *
2712 * WORK(IU) is LDA by M and WORK(IR) is LDA by M
2713 *
2714  ldwrku = lda
2715  ir = iu + ldwrku*m
2716  ldwrkr = lda
2717  ELSE IF( lwork.GE.wrkbl+( lda+m )*m ) THEN
2718 *
2719 * WORK(IU) is LDA by M and WORK(IR) is M by M
2720 *
2721  ldwrku = lda
2722  ir = iu + ldwrku*m
2723  ldwrkr = m
2724  ELSE
2725 *
2726 * WORK(IU) is M by M and WORK(IR) is M by M
2727 *
2728  ldwrku = m
2729  ir = iu + ldwrku*m
2730  ldwrkr = m
2731  END IF
2732  itau = ir + ldwrkr*m
2733  iwork = itau + m
2734 *
2735 * Compute A=L*Q
2736 * (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
2737 * (RWorkspace: 0)
2738 *
2739  CALL cgelqf( m, n, a, lda, work( itau ),
2740  $ work( iwork ), lwork-iwork+1, ierr )
2741 *
2742 * Copy L to WORK(IU), zeroing out below it
2743 *
2744  CALL clacpy( 'L', m, m, a, lda, work( iu ),
2745  $ ldwrku )
2746  CALL claset( 'U', m-1, m-1, czero, czero,
2747  $ work( iu+ldwrku ), ldwrku )
2748 *
2749 * Generate Q in A
2750 * (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
2751 * (RWorkspace: 0)
2752 *
2753  CALL cunglq( m, n, m, a, lda, work( itau ),
2754  $ work( iwork ), lwork-iwork+1, ierr )
2755  ie = 1
2756  itauq = itau
2757  itaup = itauq + m
2758  iwork = itaup + m
2759 *
2760 * Bidiagonalize L in WORK(IU), copying result to
2761 * WORK(IR)
2762 * (CWorkspace: need 2*M*M+3*M,
2763 * prefer 2*M*M+2*M+2*M*NB)
2764 * (RWorkspace: need M)
2765 *
2766  CALL cgebrd( m, m, work( iu ), ldwrku, s,
2767  $ rwork( ie ), work( itauq ),
2768  $ work( itaup ), work( iwork ),
2769  $ lwork-iwork+1, ierr )
2770  CALL clacpy( 'L', m, m, work( iu ), ldwrku,
2771  $ work( ir ), ldwrkr )
2772 *
2773 * Generate right bidiagonalizing vectors in WORK(IU)
2774 * (CWorkspace: need 2*M*M+3*M-1,
2775 * prefer 2*M*M+2*M+(M-1)*NB)
2776 * (RWorkspace: 0)
2777 *
2778  CALL cungbr( 'P', m, m, m, work( iu ), ldwrku,
2779  $ work( itaup ), work( iwork ),
2780  $ lwork-iwork+1, ierr )
2781 *
2782 * Generate left bidiagonalizing vectors in WORK(IR)
2783 * (CWorkspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB)
2784 * (RWorkspace: 0)
2785 *
2786  CALL cungbr( 'Q', m, m, m, work( ir ), ldwrkr,
2787  $ work( itauq ), work( iwork ),
2788  $ lwork-iwork+1, ierr )
2789  irwork = ie + m
2790 *
2791 * Perform bidiagonal QR iteration, computing left
2792 * singular vectors of L in WORK(IR) and computing
2793 * right singular vectors of L in WORK(IU)
2794 * (CWorkspace: need 2*M*M)
2795 * (RWorkspace: need BDSPAC)
2796 *
2797  CALL cbdsqr( 'U', m, m, m, 0, s, rwork( ie ),
2798  $ work( iu ), ldwrku, work( ir ),
2799  $ ldwrkr, cdum, 1, rwork( irwork ),
2800  $ info )
2801 *
2802 * Multiply right singular vectors of L in WORK(IU) by
2803 * Q in A, storing result in VT
2804 * (CWorkspace: need M*M)
2805 * (RWorkspace: 0)
2806 *
2807  CALL cgemm( 'N', 'N', m, n, m, cone, work( iu ),
2808  $ ldwrku, a, lda, czero, vt, ldvt )
2809 *
2810 * Copy left singular vectors of L to A
2811 * (CWorkspace: need M*M)
2812 * (RWorkspace: 0)
2813 *
2814  CALL clacpy( 'F', m, m, work( ir ), ldwrkr, a,
2815  $ lda )
2816 *
2817  ELSE
2818 *
2819 * Insufficient workspace for a fast algorithm
2820 *
2821  itau = 1
2822  iwork = itau + m
2823 *
2824 * Compute A=L*Q, copying result to VT
2825 * (CWorkspace: need 2*M, prefer M+M*NB)
2826 * (RWorkspace: 0)
2827 *
2828  CALL cgelqf( m, n, a, lda, work( itau ),
2829  $ work( iwork ), lwork-iwork+1, ierr )
2830  CALL clacpy( 'U', m, n, a, lda, vt, ldvt )
2831 *
2832 * Generate Q in VT
2833 * (CWorkspace: need 2*M, prefer M+M*NB)
2834 * (RWorkspace: 0)
2835 *
2836  CALL cunglq( m, n, m, vt, ldvt, work( itau ),
2837  $ work( iwork ), lwork-iwork+1, ierr )
2838  ie = 1
2839  itauq = itau
2840  itaup = itauq + m
2841  iwork = itaup + m
2842 *
2843 * Zero out above L in A
2844 *
2845  CALL claset( 'U', m-1, m-1, czero, czero,
2846  $ a( 1, 2 ), lda )
2847 *
2848 * Bidiagonalize L in A
2849 * (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
2850 * (RWorkspace: need M)
2851 *
2852  CALL cgebrd( m, m, a, lda, s, rwork( ie ),
2853  $ work( itauq ), work( itaup ),
2854  $ work( iwork ), lwork-iwork+1, ierr )
2855 *
2856 * Multiply right vectors bidiagonalizing L by Q in VT
2857 * (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
2858 * (RWorkspace: 0)
2859 *
2860  CALL cunmbr( 'P', 'L', 'C', m, n, m, a, lda,
2861  $ work( itaup ), vt, ldvt,
2862  $ work( iwork ), lwork-iwork+1, ierr )
2863 *
2864 * Generate left bidiagonalizing vectors of L in A
2865 * (CWorkspace: need 3*M, prefer 2*M+M*NB)
2866 * (RWorkspace: 0)
2867 *
2868  CALL cungbr( 'Q', m, m, m, a, lda, work( itauq ),
2869  $ work( iwork ), lwork-iwork+1, ierr )
2870  irwork = ie + m
2871 *
2872 * Perform bidiagonal QR iteration, computing left
2873 * singular vectors of A in A and computing right
2874 * singular vectors of A in VT
2875 * (CWorkspace: 0)
2876 * (RWorkspace: need BDSPAC)
2877 *
2878  CALL cbdsqr( 'U', m, n, m, 0, s, rwork( ie ), vt,
2879  $ ldvt, a, lda, cdum, 1,
2880  $ rwork( irwork ), info )
2881 *
2882  END IF
2883 *
2884  ELSE IF( wntuas ) THEN
2885 *
2886 * Path 6t(N much larger than M, JOBU='S' or 'A',
2887 * JOBVT='S')
2888 * M right singular vectors to be computed in VT and
2889 * M left singular vectors to be computed in U
2890 *
2891  IF( lwork.GE.m*m+3*m ) THEN
2892 *
2893 * Sufficient workspace for a fast algorithm
2894 *
2895  iu = 1
2896  IF( lwork.GE.wrkbl+lda*m ) THEN
2897 *
2898 * WORK(IU) is LDA by N
2899 *
2900  ldwrku = lda
2901  ELSE
2902 *
2903 * WORK(IU) is LDA by M
2904 *
2905  ldwrku = m
2906  END IF
2907  itau = iu + ldwrku*m
2908  iwork = itau + m
2909 *
2910 * Compute A=L*Q
2911 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2912 * (RWorkspace: 0)
2913 *
2914  CALL cgelqf( m, n, a, lda, work( itau ),
2915  $ work( iwork ), lwork-iwork+1, ierr )
2916 *
2917 * Copy L to WORK(IU), zeroing out above it
2918 *
2919  CALL clacpy( 'L', m, m, a, lda, work( iu ),
2920  $ ldwrku )
2921  CALL claset( 'U', m-1, m-1, czero, czero,
2922  $ work( iu+ldwrku ), ldwrku )
2923 *
2924 * Generate Q in A
2925 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
2926 * (RWorkspace: 0)
2927 *
2928  CALL cunglq( m, n, m, a, lda, work( itau ),
2929  $ work( iwork ), lwork-iwork+1, ierr )
2930  ie = 1
2931  itauq = itau
2932  itaup = itauq + m
2933  iwork = itaup + m
2934 *
2935 * Bidiagonalize L in WORK(IU), copying result to U
2936 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
2937 * (RWorkspace: need M)
2938 *
2939  CALL cgebrd( m, m, work( iu ), ldwrku, s,
2940  $ rwork( ie ), work( itauq ),
2941  $ work( itaup ), work( iwork ),
2942  $ lwork-iwork+1, ierr )
2943  CALL clacpy( 'L', m, m, work( iu ), ldwrku, u,
2944  $ ldu )
2945 *
2946 * Generate right bidiagonalizing vectors in WORK(IU)
2947 * (CWorkspace: need M*M+3*M-1,
2948 * prefer M*M+2*M+(M-1)*NB)
2949 * (RWorkspace: 0)
2950 *
2951  CALL cungbr( 'P', m, m, m, work( iu ), ldwrku,
2952  $ work( itaup ), work( iwork ),
2953  $ lwork-iwork+1, ierr )
2954 *
2955 * Generate left bidiagonalizing vectors in U
2956 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
2957 * (RWorkspace: 0)
2958 *
2959  CALL cungbr( 'Q', m, m, m, u, ldu, work( itauq ),
2960  $ work( iwork ), lwork-iwork+1, ierr )
2961  irwork = ie + m
2962 *
2963 * Perform bidiagonal QR iteration, computing left
2964 * singular vectors of L in U and computing right
2965 * singular vectors of L in WORK(IU)
2966 * (CWorkspace: need M*M)
2967 * (RWorkspace: need BDSPAC)
2968 *
2969  CALL cbdsqr( 'U', m, m, m, 0, s, rwork( ie ),
2970  $ work( iu ), ldwrku, u, ldu, cdum, 1,
2971  $ rwork( irwork ), info )
2972 *
2973 * Multiply right singular vectors of L in WORK(IU) by
2974 * Q in A, storing result in VT
2975 * (CWorkspace: need M*M)
2976 * (RWorkspace: 0)
2977 *
2978  CALL cgemm( 'N', 'N', m, n, m, cone, work( iu ),
2979  $ ldwrku, a, lda, czero, vt, ldvt )
2980 *
2981  ELSE
2982 *
2983 * Insufficient workspace for a fast algorithm
2984 *
2985  itau = 1
2986  iwork = itau + m
2987 *
2988 * Compute A=L*Q, copying result to VT
2989 * (CWorkspace: need 2*M, prefer M+M*NB)
2990 * (RWorkspace: 0)
2991 *
2992  CALL cgelqf( m, n, a, lda, work( itau ),
2993  $ work( iwork ), lwork-iwork+1, ierr )
2994  CALL clacpy( 'U', m, n, a, lda, vt, ldvt )
2995 *
2996 * Generate Q in VT
2997 * (CWorkspace: need 2*M, prefer M+M*NB)
2998 * (RWorkspace: 0)
2999 *
3000  CALL cunglq( m, n, m, vt, ldvt, work( itau ),
3001  $ work( iwork ), lwork-iwork+1, ierr )
3002 *
3003 * Copy L to U, zeroing out above it
3004 *
3005  CALL clacpy( 'L', m, m, a, lda, u, ldu )
3006  CALL claset( 'U', m-1, m-1, czero, czero,
3007  $ u( 1, 2 ), ldu )
3008  ie = 1
3009  itauq = itau
3010  itaup = itauq + m
3011  iwork = itaup + m
3012 *
3013 * Bidiagonalize L in U
3014 * (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
3015 * (RWorkspace: need M)
3016 *
3017  CALL cgebrd( m, m, u, ldu, s, rwork( ie ),
3018  $ work( itauq ), work( itaup ),
3019  $ work( iwork ), lwork-iwork+1, ierr )
3020 *
3021 * Multiply right bidiagonalizing vectors in U by Q
3022 * in VT
3023 * (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
3024 * (RWorkspace: 0)
3025 *
3026  CALL cunmbr( 'P', 'L', 'C', m, n, m, u, ldu,
3027  $ work( itaup ), vt, ldvt,
3028  $ work( iwork ), lwork-iwork+1, ierr )
3029 *
3030 * Generate left bidiagonalizing vectors in U
3031 * (CWorkspace: need 3*M, prefer 2*M+M*NB)
3032 * (RWorkspace: 0)
3033 *
3034  CALL cungbr( 'Q', m, m, m, u, ldu, work( itauq ),
3035  $ work( iwork ), lwork-iwork+1, ierr )
3036  irwork = ie + m
3037 *
3038 * Perform bidiagonal QR iteration, computing left
3039 * singular vectors of A in U and computing right
3040 * singular vectors of A in VT
3041 * (CWorkspace: 0)
3042 * (RWorkspace: need BDSPAC)
3043 *
3044  CALL cbdsqr( 'U', m, n, m, 0, s, rwork( ie ), vt,
3045  $ ldvt, u, ldu, cdum, 1,
3046  $ rwork( irwork ), info )
3047 *
3048  END IF
3049 *
3050  END IF
3051 *
3052  ELSE IF( wntva ) THEN
3053 *
3054  IF( wntun ) THEN
3055 *
3056 * Path 7t(N much larger than M, JOBU='N', JOBVT='A')
3057 * N right singular vectors to be computed in VT and
3058 * no left singular vectors to be computed
3059 *
3060  IF( lwork.GE.m*m+max( n+m, 3*m ) ) THEN
3061 *
3062 * Sufficient workspace for a fast algorithm
3063 *
3064  ir = 1
3065  IF( lwork.GE.wrkbl+lda*m ) THEN
3066 *
3067 * WORK(IR) is LDA by M
3068 *
3069  ldwrkr = lda
3070  ELSE
3071 *
3072 * WORK(IR) is M by M
3073 *
3074  ldwrkr = m
3075  END IF
3076  itau = ir + ldwrkr*m
3077  iwork = itau + m
3078 *
3079 * Compute A=L*Q, copying result to VT
3080 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
3081 * (RWorkspace: 0)
3082 *
3083  CALL cgelqf( m, n, a, lda, work( itau ),
3084  $ work( iwork ), lwork-iwork+1, ierr )
3085  CALL clacpy( 'U', m, n, a, lda, vt, ldvt )
3086 *
3087 * Copy L to WORK(IR), zeroing out above it
3088 *
3089  CALL clacpy( 'L', m, m, a, lda, work( ir ),
3090  $ ldwrkr )
3091  CALL claset( 'U', m-1, m-1, czero, czero,
3092  $ work( ir+ldwrkr ), ldwrkr )
3093 *
3094 * Generate Q in VT
3095 * (CWorkspace: need M*M+M+N, prefer M*M+M+N*NB)
3096 * (RWorkspace: 0)
3097 *
3098  CALL cunglq( n, n, m, vt, ldvt, work( itau ),
3099  $ work( iwork ), lwork-iwork+1, ierr )
3100  ie = 1
3101  itauq = itau
3102  itaup = itauq + m
3103  iwork = itaup + m
3104 *
3105 * Bidiagonalize L in WORK(IR)
3106 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
3107 * (RWorkspace: need M)
3108 *
3109  CALL cgebrd( m, m, work( ir ), ldwrkr, s,
3110  $ rwork( ie ), work( itauq ),
3111  $ work( itaup ), work( iwork ),
3112  $ lwork-iwork+1, ierr )
3113 *
3114 * Generate right bidiagonalizing vectors in WORK(IR)
3115 * (CWorkspace: need M*M+3*M-1,
3116 * prefer M*M+2*M+(M-1)*NB)
3117 * (RWorkspace: 0)
3118 *
3119  CALL cungbr( 'P', m, m, m, work( ir ), ldwrkr,
3120  $ work( itaup ), work( iwork ),
3121  $ lwork-iwork+1, ierr )
3122  irwork = ie + m
3123 *
3124 * Perform bidiagonal QR iteration, computing right
3125 * singular vectors of L in WORK(IR)
3126 * (CWorkspace: need M*M)
3127 * (RWorkspace: need BDSPAC)
3128 *
3129  CALL cbdsqr( 'U', m, m, 0, 0, s, rwork( ie ),
3130  $ work( ir ), ldwrkr, cdum, 1, cdum, 1,
3131  $ rwork( irwork ), info )
3132 *
3133 * Multiply right singular vectors of L in WORK(IR) by
3134 * Q in VT, storing result in A
3135 * (CWorkspace: need M*M)
3136 * (RWorkspace: 0)
3137 *
3138  CALL cgemm( 'N', 'N', m, n, m, cone, work( ir ),
3139  $ ldwrkr, vt, ldvt, czero, a, lda )
3140 *
3141 * Copy right singular vectors of A from A to VT
3142 *
3143  CALL clacpy( 'F', m, n, a, lda, vt, ldvt )
3144 *
3145  ELSE
3146 *
3147 * Insufficient workspace for a fast algorithm
3148 *
3149  itau = 1
3150  iwork = itau + m
3151 *
3152 * Compute A=L*Q, copying result to VT
3153 * (CWorkspace: need 2*M, prefer M+M*NB)
3154 * (RWorkspace: 0)
3155 *
3156  CALL cgelqf( m, n, a, lda, work( itau ),
3157  $ work( iwork ), lwork-iwork+1, ierr )
3158  CALL clacpy( 'U', m, n, a, lda, vt, ldvt )
3159 *
3160 * Generate Q in VT
3161 * (CWorkspace: need M+N, prefer M+N*NB)
3162 * (RWorkspace: 0)
3163 *
3164  CALL cunglq( n, n, m, vt, ldvt, work( itau ),
3165  $ work( iwork ), lwork-iwork+1, ierr )
3166  ie = 1
3167  itauq = itau
3168  itaup = itauq + m
3169  iwork = itaup + m
3170 *
3171 * Zero out above L in A
3172 *
3173  CALL claset( 'U', m-1, m-1, czero, czero,
3174  $ a( 1, 2 ), lda )
3175 *
3176 * Bidiagonalize L in A
3177 * (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
3178 * (RWorkspace: need M)
3179 *
3180  CALL cgebrd( m, m, a, lda, s, rwork( ie ),
3181  $ work( itauq ), work( itaup ),
3182  $ work( iwork ), lwork-iwork+1, ierr )
3183 *
3184 * Multiply right bidiagonalizing vectors in A by Q
3185 * in VT
3186 * (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
3187 * (RWorkspace: 0)
3188 *
3189  CALL cunmbr( 'P', 'L', 'C', m, n, m, a, lda,
3190  $ work( itaup ), vt, ldvt,
3191  $ work( iwork ), lwork-iwork+1, ierr )
3192  irwork = ie + m
3193 *
3194 * Perform bidiagonal QR iteration, computing right
3195 * singular vectors of A in VT
3196 * (CWorkspace: 0)
3197 * (RWorkspace: need BDSPAC)
3198 *
3199  CALL cbdsqr( 'U', m, n, 0, 0, s, rwork( ie ), vt,
3200  $ ldvt, cdum, 1, cdum, 1,
3201  $ rwork( irwork ), info )
3202 *
3203  END IF
3204 *
3205  ELSE IF( wntuo ) THEN
3206 *
3207 * Path 8t(N much larger than M, JOBU='O', JOBVT='A')
3208 * N right singular vectors to be computed in VT and
3209 * M left singular vectors to be overwritten on A
3210 *
3211  IF( lwork.GE.2*m*m+max( n+m, 3*m ) ) THEN
3212 *
3213 * Sufficient workspace for a fast algorithm
3214 *
3215  iu = 1
3216  IF( lwork.GE.wrkbl+2*lda*m ) THEN
3217 *
3218 * WORK(IU) is LDA by M and WORK(IR) is LDA by M
3219 *
3220  ldwrku = lda
3221  ir = iu + ldwrku*m
3222  ldwrkr = lda
3223  ELSE IF( lwork.GE.wrkbl+( lda+m )*m ) THEN
3224 *
3225 * WORK(IU) is LDA by M and WORK(IR) is M by M
3226 *
3227  ldwrku = lda
3228  ir = iu + ldwrku*m
3229  ldwrkr = m
3230  ELSE
3231 *
3232 * WORK(IU) is M by M and WORK(IR) is M by M
3233 *
3234  ldwrku = m
3235  ir = iu + ldwrku*m
3236  ldwrkr = m
3237  END IF
3238  itau = ir + ldwrkr*m
3239  iwork = itau + m
3240 *
3241 * Compute A=L*Q, copying result to VT
3242 * (CWorkspace: need 2*M*M+2*M, prefer 2*M*M+M+M*NB)
3243 * (RWorkspace: 0)
3244 *
3245  CALL cgelqf( m, n, a, lda, work( itau ),
3246  $ work( iwork ), lwork-iwork+1, ierr )
3247  CALL clacpy( 'U', m, n, a, lda, vt, ldvt )
3248 *
3249 * Generate Q in VT
3250 * (CWorkspace: need 2*M*M+M+N, prefer 2*M*M+M+N*NB)
3251 * (RWorkspace: 0)
3252 *
3253  CALL cunglq( n, n, m, vt, ldvt, work( itau ),
3254  $ work( iwork ), lwork-iwork+1, ierr )
3255 *
3256 * Copy L to WORK(IU), zeroing out above it
3257 *
3258  CALL clacpy( 'L', m, m, a, lda, work( iu ),
3259  $ ldwrku )
3260  CALL claset( 'U', m-1, m-1, czero, czero,
3261  $ work( iu+ldwrku ), ldwrku )
3262  ie = 1
3263  itauq = itau
3264  itaup = itauq + m
3265  iwork = itaup + m
3266 *
3267 * Bidiagonalize L in WORK(IU), copying result to
3268 * WORK(IR)
3269 * (CWorkspace: need 2*M*M+3*M,
3270 * prefer 2*M*M+2*M+2*M*NB)
3271 * (RWorkspace: need M)
3272 *
3273  CALL cgebrd( m, m, work( iu ), ldwrku, s,
3274  $ rwork( ie ), work( itauq ),
3275  $ work( itaup ), work( iwork ),
3276  $ lwork-iwork+1, ierr )
3277  CALL clacpy( 'L', m, m, work( iu ), ldwrku,
3278  $ work( ir ), ldwrkr )
3279 *
3280 * Generate right bidiagonalizing vectors in WORK(IU)
3281 * (CWorkspace: need 2*M*M+3*M-1,
3282 * prefer 2*M*M+2*M+(M-1)*NB)
3283 * (RWorkspace: 0)
3284 *
3285  CALL cungbr( 'P', m, m, m, work( iu ), ldwrku,
3286  $ work( itaup ), work( iwork ),
3287  $ lwork-iwork+1, ierr )
3288 *
3289 * Generate left bidiagonalizing vectors in WORK(IR)
3290 * (CWorkspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB)
3291 * (RWorkspace: 0)
3292 *
3293  CALL cungbr( 'Q', m, m, m, work( ir ), ldwrkr,
3294  $ work( itauq ), work( iwork ),
3295  $ lwork-iwork+1, ierr )
3296  irwork = ie + m
3297 *
3298 * Perform bidiagonal QR iteration, computing left
3299 * singular vectors of L in WORK(IR) and computing
3300 * right singular vectors of L in WORK(IU)
3301 * (CWorkspace: need 2*M*M)
3302 * (RWorkspace: need BDSPAC)
3303 *
3304  CALL cbdsqr( 'U', m, m, m, 0, s, rwork( ie ),
3305  $ work( iu ), ldwrku, work( ir ),
3306  $ ldwrkr, cdum, 1, rwork( irwork ),
3307  $ info )
3308 *
3309 * Multiply right singular vectors of L in WORK(IU) by
3310 * Q in VT, storing result in A
3311 * (CWorkspace: need M*M)
3312 * (RWorkspace: 0)
3313 *
3314  CALL cgemm( 'N', 'N', m, n, m, cone, work( iu ),
3315  $ ldwrku, vt, ldvt, czero, a, lda )
3316 *
3317 * Copy right singular vectors of A from A to VT
3318 *
3319  CALL clacpy( 'F', m, n, a, lda, vt, ldvt )
3320 *
3321 * Copy left singular vectors of A from WORK(IR) to A
3322 *
3323  CALL clacpy( 'F', m, m, work( ir ), ldwrkr, a,
3324  $ lda )
3325 *
3326  ELSE
3327 *
3328 * Insufficient workspace for a fast algorithm
3329 *
3330  itau = 1
3331  iwork = itau + m
3332 *
3333 * Compute A=L*Q, copying result to VT
3334 * (CWorkspace: need 2*M, prefer M+M*NB)
3335 * (RWorkspace: 0)
3336 *
3337  CALL cgelqf( m, n, a, lda, work( itau ),
3338  $ work( iwork ), lwork-iwork+1, ierr )
3339  CALL clacpy( 'U', m, n, a, lda, vt, ldvt )
3340 *
3341 * Generate Q in VT
3342 * (CWorkspace: need M+N, prefer M+N*NB)
3343 * (RWorkspace: 0)
3344 *
3345  CALL cunglq( n, n, m, vt, ldvt, work( itau ),
3346  $ work( iwork ), lwork-iwork+1, ierr )
3347  ie = 1
3348  itauq = itau
3349  itaup = itauq + m
3350  iwork = itaup + m
3351 *
3352 * Zero out above L in A
3353 *
3354  CALL claset( 'U', m-1, m-1, czero, czero,
3355  $ a( 1, 2 ), lda )
3356 *
3357 * Bidiagonalize L in A
3358 * (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
3359 * (RWorkspace: need M)
3360 *
3361  CALL cgebrd( m, m, a, lda, s, rwork( ie ),
3362  $ work( itauq ), work( itaup ),
3363  $ work( iwork ), lwork-iwork+1, ierr )
3364 *
3365 * Multiply right bidiagonalizing vectors in A by Q
3366 * in VT
3367 * (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
3368 * (RWorkspace: 0)
3369 *
3370  CALL cunmbr( 'P', 'L', 'C', m, n, m, a, lda,
3371  $ work( itaup ), vt, ldvt,
3372  $ work( iwork ), lwork-iwork+1, ierr )
3373 *
3374 * Generate left bidiagonalizing vectors in A
3375 * (CWorkspace: need 3*M, prefer 2*M+M*NB)
3376 * (RWorkspace: 0)
3377 *
3378  CALL cungbr( 'Q', m, m, m, a, lda, work( itauq ),
3379  $ work( iwork ), lwork-iwork+1, ierr )
3380  irwork = ie + m
3381 *
3382 * Perform bidiagonal QR iteration, computing left
3383 * singular vectors of A in A and computing right
3384 * singular vectors of A in VT
3385 * (CWorkspace: 0)
3386 * (RWorkspace: need BDSPAC)
3387 *
3388  CALL cbdsqr( 'U', m, n, m, 0, s, rwork( ie ), vt,
3389  $ ldvt, a, lda, cdum, 1,
3390  $ rwork( irwork ), info )
3391 *
3392  END IF
3393 *
3394  ELSE IF( wntuas ) THEN
3395 *
3396 * Path 9t(N much larger than M, JOBU='S' or 'A',
3397 * JOBVT='A')
3398 * N right singular vectors to be computed in VT and
3399 * M left singular vectors to be computed in U
3400 *
3401  IF( lwork.GE.m*m+max( n+m, 3*m ) ) THEN
3402 *
3403 * Sufficient workspace for a fast algorithm
3404 *
3405  iu = 1
3406  IF( lwork.GE.wrkbl+lda*m ) THEN
3407 *
3408 * WORK(IU) is LDA by M
3409 *
3410  ldwrku = lda
3411  ELSE
3412 *
3413 * WORK(IU) is M by M
3414 *
3415  ldwrku = m
3416  END IF
3417  itau = iu + ldwrku*m
3418  iwork = itau + m
3419 *
3420 * Compute A=L*Q, copying result to VT
3421 * (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
3422 * (RWorkspace: 0)
3423 *
3424  CALL cgelqf( m, n, a, lda, work( itau ),
3425  $ work( iwork ), lwork-iwork+1, ierr )
3426  CALL clacpy( 'U', m, n, a, lda, vt, ldvt )
3427 *
3428 * Generate Q in VT
3429 * (CWorkspace: need M*M+M+N, prefer M*M+M+N*NB)
3430 * (RWorkspace: 0)
3431 *
3432  CALL cunglq( n, n, m, vt, ldvt, work( itau ),
3433  $ work( iwork ), lwork-iwork+1, ierr )
3434 *
3435 * Copy L to WORK(IU), zeroing out above it
3436 *
3437  CALL clacpy( 'L', m, m, a, lda, work( iu ),
3438  $ ldwrku )
3439  CALL claset( 'U', m-1, m-1, czero, czero,
3440  $ work( iu+ldwrku ), ldwrku )
3441  ie = 1
3442  itauq = itau
3443  itaup = itauq + m
3444  iwork = itaup + m
3445 *
3446 * Bidiagonalize L in WORK(IU), copying result to U
3447 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
3448 * (RWorkspace: need M)
3449 *
3450  CALL cgebrd( m, m, work( iu ), ldwrku, s,
3451  $ rwork( ie ), work( itauq ),
3452  $ work( itaup ), work( iwork ),
3453  $ lwork-iwork+1, ierr )
3454  CALL clacpy( 'L', m, m, work( iu ), ldwrku, u,
3455  $ ldu )
3456 *
3457 * Generate right bidiagonalizing vectors in WORK(IU)
3458 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+(M-1)*NB)
3459 * (RWorkspace: 0)
3460 *
3461  CALL cungbr( 'P', m, m, m, work( iu ), ldwrku,
3462  $ work( itaup ), work( iwork ),
3463  $ lwork-iwork+1, ierr )
3464 *
3465 * Generate left bidiagonalizing vectors in U
3466 * (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
3467 * (RWorkspace: 0)
3468 *
3469  CALL cungbr( 'Q', m, m, m, u, ldu, work( itauq ),
3470  $ work( iwork ), lwork-iwork+1, ierr )
3471  irwork = ie + m
3472 *
3473 * Perform bidiagonal QR iteration, computing left
3474 * singular vectors of L in U and computing right
3475 * singular vectors of L in WORK(IU)
3476 * (CWorkspace: need M*M)
3477 * (RWorkspace: need BDSPAC)
3478 *
3479  CALL cbdsqr( 'U', m, m, m, 0, s, rwork( ie ),
3480  $ work( iu ), ldwrku, u, ldu, cdum, 1,
3481  $ rwork( irwork ), info )
3482 *
3483 * Multiply right singular vectors of L in WORK(IU) by
3484 * Q in VT, storing result in A
3485 * (CWorkspace: need M*M)
3486 * (RWorkspace: 0)
3487 *
3488  CALL cgemm( 'N', 'N', m, n, m, cone, work( iu ),
3489  $ ldwrku, vt, ldvt, czero, a, lda )
3490 *
3491 * Copy right singular vectors of A from A to VT
3492 *
3493  CALL clacpy( 'F', m, n, a, lda, vt, ldvt )
3494 *
3495  ELSE
3496 *
3497 * Insufficient workspace for a fast algorithm
3498 *
3499  itau = 1
3500  iwork = itau + m
3501 *
3502 * Compute A=L*Q, copying result to VT
3503 * (CWorkspace: need 2*M, prefer M+M*NB)
3504 * (RWorkspace: 0)
3505 *
3506  CALL cgelqf( m, n, a, lda, work( itau ),
3507  $ work( iwork ), lwork-iwork+1, ierr )
3508  CALL clacpy( 'U', m, n, a, lda, vt, ldvt )
3509 *
3510 * Generate Q in VT
3511 * (CWorkspace: need M+N, prefer M+N*NB)
3512 * (RWorkspace: 0)
3513 *
3514  CALL cunglq( n, n, m, vt, ldvt, work( itau ),
3515  $ work( iwork ), lwork-iwork+1, ierr )
3516 *
3517 * Copy L to U, zeroing out above it
3518 *
3519  CALL clacpy( 'L', m, m, a, lda, u, ldu )
3520  CALL claset( 'U', m-1, m-1, czero, czero,
3521  $ u( 1, 2 ), ldu )
3522  ie = 1
3523  itauq = itau
3524  itaup = itauq + m
3525  iwork = itaup + m
3526 *
3527 * Bidiagonalize L in U
3528 * (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
3529 * (RWorkspace: need M)
3530 *
3531  CALL cgebrd( m, m, u, ldu, s, rwork( ie ),
3532  $ work( itauq ), work( itaup ),
3533  $ work( iwork ), lwork-iwork+1, ierr )
3534 *
3535 * Multiply right bidiagonalizing vectors in U by Q
3536 * in VT
3537 * (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
3538 * (RWorkspace: 0)
3539 *
3540  CALL cunmbr( 'P', 'L', 'C', m, n, m, u, ldu,
3541  $ work( itaup ), vt, ldvt,
3542  $ work( iwork ), lwork-iwork+1, ierr )
3543 *
3544 * Generate left bidiagonalizing vectors in U
3545 * (CWorkspace: need 3*M, prefer 2*M+M*NB)
3546 * (RWorkspace: 0)
3547 *
3548  CALL cungbr( 'Q', m, m, m, u, ldu, work( itauq ),
3549  $ work( iwork ), lwork-iwork+1, ierr )
3550  irwork = ie + m
3551 *
3552 * Perform bidiagonal QR iteration, computing left
3553 * singular vectors of A in U and computing right
3554 * singular vectors of A in VT
3555 * (CWorkspace: 0)
3556 * (RWorkspace: need BDSPAC)
3557 *
3558  CALL cbdsqr( 'U', m, n, m, 0, s, rwork( ie ), vt,
3559  $ ldvt, u, ldu, cdum, 1,
3560  $ rwork( irwork ), info )
3561 *
3562  END IF
3563 *
3564  END IF
3565 *
3566  END IF
3567 *
3568  ELSE
3569 *
3570 * N .LT. MNTHR
3571 *
3572 * Path 10t(N greater than M, but not much larger)
3573 * Reduce to bidiagonal form without LQ decomposition
3574 *
3575  ie = 1
3576  itauq = 1
3577  itaup = itauq + m
3578  iwork = itaup + m
3579 *
3580 * Bidiagonalize A
3581 * (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)
3582 * (RWorkspace: M)
3583 *
3584  CALL cgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
3585  $ work( itaup ), work( iwork ), lwork-iwork+1,
3586  $ ierr )
3587  IF( wntuas ) THEN
3588 *
3589 * If left singular vectors desired in U, copy result to U
3590 * and generate left bidiagonalizing vectors in U
3591 * (CWorkspace: need 3*M-1, prefer 2*M+(M-1)*NB)
3592 * (RWorkspace: 0)
3593 *
3594  CALL clacpy( 'L', m, m, a, lda, u, ldu )
3595  CALL cungbr( 'Q', m, m, n, u, ldu, work( itauq ),
3596  $ work( iwork ), lwork-iwork+1, ierr )
3597  END IF
3598  IF( wntvas ) THEN
3599 *
3600 * If right singular vectors desired in VT, copy result to
3601 * VT and generate right bidiagonalizing vectors in VT
3602 * (CWorkspace: need 2*M+NRVT, prefer 2*M+NRVT*NB)
3603 * (RWorkspace: 0)
3604 *
3605  CALL clacpy( 'U', m, n, a, lda, vt, ldvt )
3606  IF( wntva )
3607  $ nrvt = n
3608  IF( wntvs )
3609  $ nrvt = m
3610  CALL cungbr( 'P', nrvt, n, m, vt, ldvt, work( itaup ),
3611  $ work( iwork ), lwork-iwork+1, ierr )
3612  END IF
3613  IF( wntuo ) THEN
3614 *
3615 * If left singular vectors desired in A, generate left
3616 * bidiagonalizing vectors in A
3617 * (CWorkspace: need 3*M-1, prefer 2*M+(M-1)*NB)
3618 * (RWorkspace: 0)
3619 *
3620  CALL cungbr( 'Q', m, m, n, a, lda, work( itauq ),
3621  $ work( iwork ), lwork-iwork+1, ierr )
3622  END IF
3623  IF( wntvo ) THEN
3624 *
3625 * If right singular vectors desired in A, generate right
3626 * bidiagonalizing vectors in A
3627 * (CWorkspace: need 3*M, prefer 2*M+M*NB)
3628 * (RWorkspace: 0)
3629 *
3630  CALL cungbr( 'P', m, n, m, a, lda, work( itaup ),
3631  $ work( iwork ), lwork-iwork+1, ierr )
3632  END IF
3633  irwork = ie + m
3634  IF( wntuas .OR. wntuo )
3635  $ nru = m
3636  IF( wntun )
3637  $ nru = 0
3638  IF( wntvas .OR. wntvo )
3639  $ ncvt = n
3640  IF( wntvn )
3641  $ ncvt = 0
3642  IF( ( .NOT.wntuo ) .AND. ( .NOT.wntvo ) ) THEN
3643 *
3644 * Perform bidiagonal QR iteration, if desired, computing
3645 * left singular vectors in U and computing right singular
3646 * vectors in VT
3647 * (CWorkspace: 0)
3648 * (RWorkspace: need BDSPAC)
3649 *
3650  CALL cbdsqr( 'L', m, ncvt, nru, 0, s, rwork( ie ), vt,
3651  $ ldvt, u, ldu, cdum, 1, rwork( irwork ),
3652  $ info )
3653  ELSE IF( ( .NOT.wntuo ) .AND. wntvo ) THEN
3654 *
3655 * Perform bidiagonal QR iteration, if desired, computing
3656 * left singular vectors in U and computing right singular
3657 * vectors in A
3658 * (CWorkspace: 0)
3659 * (RWorkspace: need BDSPAC)
3660 *
3661  CALL cbdsqr( 'L', m, ncvt, nru, 0, s, rwork( ie ), a,
3662  $ lda, u, ldu, cdum, 1, rwork( irwork ),
3663  $ info )
3664  ELSE
3665 *
3666 * Perform bidiagonal QR iteration, if desired, computing
3667 * left singular vectors in A and computing right singular
3668 * vectors in VT
3669 * (CWorkspace: 0)
3670 * (RWorkspace: need BDSPAC)
3671 *
3672  CALL cbdsqr( 'L', m, ncvt, nru, 0, s, rwork( ie ), vt,
3673  $ ldvt, a, lda, cdum, 1, rwork( irwork ),
3674  $ info )
3675  END IF
3676 *
3677  END IF
3678 *
3679  END IF
3680 *
3681 * Undo scaling if necessary
3682 *
3683  IF( iscl.EQ.1 ) THEN
3684  IF( anrm.GT.bignum )
3685  $ CALL slascl( 'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
3686  $ ierr )
3687  IF( info.NE.0 .AND. anrm.GT.bignum )
3688  $ CALL slascl( 'G', 0, 0, bignum, anrm, minmn-1, 1,
3689  $ rwork( ie ), minmn, ierr )
3690  IF( anrm.LT.smlnum )
3691  $ CALL slascl( 'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
3692  $ ierr )
3693  IF( info.NE.0 .AND. anrm.LT.smlnum )
3694  $ CALL slascl( 'G', 0, 0, smlnum, anrm, minmn-1, 1,
3695  $ rwork( ie ), minmn, ierr )
3696  END IF
3697 *
3698 * Return optimal workspace in WORK(1)
3699 *
3700  work( 1 ) = maxwrk
3701 *
3702  RETURN
3703 *
3704 * End of CGESVD
3705 *
Here is the call graph for this function:
Here is the caller graph for this function:
cunmbr
subroutine cunmbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMBR
Definition: cunmbr.f:199
clange
real function clange(NORM, M, N, A, LDA, WORK)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: clange.f:117
cgemm
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:189
cunglq
subroutine cunglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGLQ
Definition: cunglq.f:129
clacpy
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:105
slascl
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: slascl.f:145
cungqr
subroutine cungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGQR
Definition: cungqr.f:130
cgebrd
subroutine cgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
CGEBRD
Definition: cgebrd.f:208
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
cgelqf
subroutine cgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGELQF
Definition: cgelqf.f:145
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
claset
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: claset.f:108
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:70
cungbr
subroutine cungbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGBR
Definition: cungbr.f:159
cgeqrf
subroutine cgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGEQRF
Definition: cgeqrf.f:147
cbdsqr
subroutine cbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, RWORK, INFO)
CBDSQR
Definition: cbdsqr.f:224
ilaenv
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: tstiee.f:83
clascl
subroutine clascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: clascl.f:145