LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ zgesdd()

subroutine zgesdd ( character  JOBZ,
integer  M,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  S,
complex*16, dimension( ldu, * )  U,
integer  LDU,
complex*16, dimension( ldvt, * )  VT,
integer  LDVT,
complex*16, dimension( * )  WORK,
integer  LWORK,
double precision, dimension( * )  RWORK,
integer, dimension( * )  IWORK,
integer  INFO 
)

ZGESDD

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

Purpose:
 ZGESDD computes the singular value decomposition (SVD) of a complex
 M-by-N matrix A, optionally computing the left and/or right singular
 vectors, by using divide-and-conquer method. 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 VT = V**H, not V.

 The divide and conquer algorithm makes very mild assumptions about
 floating point arithmetic. It will work on machines with a guard
 digit in add/subtract, or on those binary machines without guard
 digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
 Cray-2. It could conceivably fail on hexadecimal or decimal machines
 without guard digits, but we know of none.
Parameters
[in]JOBZ
          JOBZ is CHARACTER*1
          Specifies options for computing all or part of the matrix U:
          = 'A':  all M columns of U and all N rows of V**H are
                  returned in the arrays U and VT;
          = 'S':  the first min(M,N) columns of U and the first
                  min(M,N) rows of V**H are returned in the arrays U
                  and VT;
          = 'O':  If M >= N, the first N columns of U are overwritten
                  in the array A and all rows of V**H are returned in
                  the array VT;
                  otherwise, all columns of U are returned in the
                  array U and the first M rows of V**H are overwritten
                  in the array A;
          = 'N':  no columns of U or rows of V**H are computed.
[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*16 array, dimension (LDA,N)
          On entry, the M-by-N matrix A.
          On exit,
          if JOBZ = 'O',  A is overwritten with the first N columns
                          of U (the left singular vectors, stored
                          columnwise) if M >= N;
                          A is overwritten with the first M rows
                          of V**H (the right singular vectors, stored
                          rowwise) otherwise.
          if JOBZ .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 DOUBLE PRECISION array, dimension (min(M,N))
          The singular values of A, sorted so that S(i) >= S(i+1).
[out]U
          U is COMPLEX*16 array, dimension (LDU,UCOL)
          UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N;
          UCOL = min(M,N) if JOBZ = 'S'.
          If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M
          unitary matrix U;
          if JOBZ = 'S', U contains the first min(M,N) columns of U
          (the left singular vectors, stored columnwise);
          if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced.
[in]LDU
          LDU is INTEGER
          The leading dimension of the array U.  LDU >= 1;
          if JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M.
[out]VT
          VT is COMPLEX*16 array, dimension (LDVT,N)
          If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the
          N-by-N unitary matrix V**H;
          if JOBZ = 'S', VT contains the first min(M,N) rows of
          V**H (the right singular vectors, stored rowwise);
          if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced.
[in]LDVT
          LDVT is INTEGER
          The leading dimension of the array VT.  LDVT >= 1;
          if JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N;
          if JOBZ = 'S', LDVT >= min(M,N).
[out]WORK
          WORK is COMPLEX*16 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 >= 1.
          If LWORK = -1, a workspace query is assumed.  The optimal
          size for the WORK array is calculated and stored in WORK(1),
          and no other work except argument checking is performed.

          Let mx = max(M,N) and mn = min(M,N).
          If JOBZ = 'N', LWORK >= 2*mn + mx.
          If JOBZ = 'O', LWORK >= 2*mn*mn + 2*mn + mx.
          If JOBZ = 'S', LWORK >=   mn*mn + 3*mn.
          If JOBZ = 'A', LWORK >=   mn*mn + 2*mn + mx.
          These are not tight minimums in all cases; see comments inside code.
          For good performance, LWORK should generally be larger;
          a query is recommended.
[out]RWORK
          RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
          Let mx = max(M,N) and mn = min(M,N).
          If JOBZ = 'N',    LRWORK >= 5*mn (LAPACK <= 3.6 needs 7*mn);
          else if mx >> mn, LRWORK >= 5*mn*mn + 5*mn;
          else              LRWORK >= max( 5*mn*mn + 5*mn,
                                           2*mx*mn + 2*mn*mn + mn ).
[out]IWORK
          IWORK is INTEGER array, dimension (8*min(M,N))
[out]INFO
          INFO is INTEGER
          <  0:  if INFO = -i, the i-th argument had an illegal value.
          = -4:  if A had a NAN entry.
          >  0:  The updating process of DBDSDC did not converge.
          =  0:  successful exit.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Ming Gu and Huan Ren, Computer Science Division, University of California at Berkeley, USA

Definition at line 225 of file zgesdd.f.

227  implicit none
228 *
229 * -- LAPACK driver routine --
230 * -- LAPACK is a software package provided by Univ. of Tennessee, --
231 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
232 *
233 * .. Scalar Arguments ..
234  CHARACTER JOBZ
235  INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N
236 * ..
237 * .. Array Arguments ..
238  INTEGER IWORK( * )
239  DOUBLE PRECISION RWORK( * ), S( * )
240  COMPLEX*16 A( LDA, * ), U( LDU, * ), VT( LDVT, * ),
241  $ WORK( * )
242 * ..
243 *
244 * =====================================================================
245 *
246 * .. Parameters ..
247  COMPLEX*16 CZERO, CONE
248  parameter( czero = ( 0.0d+0, 0.0d+0 ),
249  $ cone = ( 1.0d+0, 0.0d+0 ) )
250  DOUBLE PRECISION ZERO, ONE
251  parameter( zero = 0.0d+0, one = 1.0d+0 )
252 * ..
253 * .. Local Scalars ..
254  LOGICAL LQUERY, WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS
255  INTEGER BLK, CHUNK, I, IE, IERR, IL, IR, IRU, IRVT,
256  $ ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT,
257  $ LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK,
258  $ MNTHR1, MNTHR2, NRWORK, NWORK, WRKBL
259  INTEGER LWORK_ZGEBRD_MN, LWORK_ZGEBRD_MM,
260  $ LWORK_ZGEBRD_NN, LWORK_ZGELQF_MN,
261  $ LWORK_ZGEQRF_MN,
262  $ LWORK_ZUNGBR_P_MN, LWORK_ZUNGBR_P_NN,
263  $ LWORK_ZUNGBR_Q_MN, LWORK_ZUNGBR_Q_MM,
264  $ LWORK_ZUNGLQ_MN, LWORK_ZUNGLQ_NN,
265  $ LWORK_ZUNGQR_MM, LWORK_ZUNGQR_MN,
266  $ LWORK_ZUNMBR_PRC_MM, LWORK_ZUNMBR_QLN_MM,
267  $ LWORK_ZUNMBR_PRC_MN, LWORK_ZUNMBR_QLN_MN,
268  $ LWORK_ZUNMBR_PRC_NN, LWORK_ZUNMBR_QLN_NN
269  DOUBLE PRECISION ANRM, BIGNUM, EPS, SMLNUM
270 * ..
271 * .. Local Arrays ..
272  INTEGER IDUM( 1 )
273  DOUBLE PRECISION DUM( 1 )
274  COMPLEX*16 CDUM( 1 )
275 * ..
276 * .. External Subroutines ..
277  EXTERNAL dbdsdc, dlascl, xerbla, zgebrd, zgelqf, zgemm,
280 * ..
281 * .. External Functions ..
282  LOGICAL LSAME, DISNAN
283  DOUBLE PRECISION DLAMCH, ZLANGE, DROUNDUP_LWORK
284  EXTERNAL lsame, dlamch, zlange, disnan,
286 * ..
287 * .. Intrinsic Functions ..
288  INTRINSIC int, max, min, sqrt
289 * ..
290 * .. Executable Statements ..
291 *
292 * Test the input arguments
293 *
294  info = 0
295  minmn = min( m, n )
296  mnthr1 = int( minmn*17.0d0 / 9.0d0 )
297  mnthr2 = int( minmn*5.0d0 / 3.0d0 )
298  wntqa = lsame( jobz, 'A' )
299  wntqs = lsame( jobz, 'S' )
300  wntqas = wntqa .OR. wntqs
301  wntqo = lsame( jobz, 'O' )
302  wntqn = lsame( jobz, 'N' )
303  lquery = ( lwork.EQ.-1 )
304  minwrk = 1
305  maxwrk = 1
306 *
307  IF( .NOT.( wntqa .OR. wntqs .OR. wntqo .OR. wntqn ) ) THEN
308  info = -1
309  ELSE IF( m.LT.0 ) THEN
310  info = -2
311  ELSE IF( n.LT.0 ) THEN
312  info = -3
313  ELSE IF( lda.LT.max( 1, m ) ) THEN
314  info = -5
315  ELSE IF( ldu.LT.1 .OR. ( wntqas .AND. ldu.LT.m ) .OR.
316  $ ( wntqo .AND. m.LT.n .AND. ldu.LT.m ) ) THEN
317  info = -8
318  ELSE IF( ldvt.LT.1 .OR. ( wntqa .AND. ldvt.LT.n ) .OR.
319  $ ( wntqs .AND. ldvt.LT.minmn ) .OR.
320  $ ( wntqo .AND. m.GE.n .AND. ldvt.LT.n ) ) THEN
321  info = -10
322  END IF
323 *
324 * Compute workspace
325 * Note: Comments in the code beginning "Workspace:" describe the
326 * minimal amount of workspace allocated at that point in the code,
327 * as well as the preferred amount for good performance.
328 * CWorkspace refers to complex workspace, and RWorkspace to
329 * real workspace. NB refers to the optimal block size for the
330 * immediately following subroutine, as returned by ILAENV.)
331 *
332  IF( info.EQ.0 ) THEN
333  minwrk = 1
334  maxwrk = 1
335  IF( m.GE.n .AND. minmn.GT.0 ) THEN
336 *
337 * There is no complex work space needed for bidiagonal SVD
338 * The real work space needed for bidiagonal SVD (dbdsdc) is
339 * BDSPAC = 3*N*N + 4*N for singular values and vectors;
340 * BDSPAC = 4*N for singular values only;
341 * not including e, RU, and RVT matrices.
342 *
343 * Compute space preferred for each routine
344  CALL zgebrd( m, n, cdum(1), m, dum(1), dum(1), cdum(1),
345  $ cdum(1), cdum(1), -1, ierr )
346  lwork_zgebrd_mn = int( cdum(1) )
347 *
348  CALL zgebrd( n, n, cdum(1), n, dum(1), dum(1), cdum(1),
349  $ cdum(1), cdum(1), -1, ierr )
350  lwork_zgebrd_nn = int( cdum(1) )
351 *
352  CALL zgeqrf( m, n, cdum(1), m, cdum(1), cdum(1), -1, ierr )
353  lwork_zgeqrf_mn = int( cdum(1) )
354 *
355  CALL zungbr( 'P', n, n, n, cdum(1), n, cdum(1), cdum(1),
356  $ -1, ierr )
357  lwork_zungbr_p_nn = int( cdum(1) )
358 *
359  CALL zungbr( 'Q', m, m, n, cdum(1), m, cdum(1), cdum(1),
360  $ -1, ierr )
361  lwork_zungbr_q_mm = int( cdum(1) )
362 *
363  CALL zungbr( 'Q', m, n, n, cdum(1), m, cdum(1), cdum(1),
364  $ -1, ierr )
365  lwork_zungbr_q_mn = int( cdum(1) )
366 *
367  CALL zungqr( m, m, n, cdum(1), m, cdum(1), cdum(1),
368  $ -1, ierr )
369  lwork_zungqr_mm = int( cdum(1) )
370 *
371  CALL zungqr( m, n, n, cdum(1), m, cdum(1), cdum(1),
372  $ -1, ierr )
373  lwork_zungqr_mn = int( cdum(1) )
374 *
375  CALL zunmbr( 'P', 'R', 'C', n, n, n, cdum(1), n, cdum(1),
376  $ cdum(1), n, cdum(1), -1, ierr )
377  lwork_zunmbr_prc_nn = int( cdum(1) )
378 *
379  CALL zunmbr( 'Q', 'L', 'N', m, m, n, cdum(1), m, cdum(1),
380  $ cdum(1), m, cdum(1), -1, ierr )
381  lwork_zunmbr_qln_mm = int( cdum(1) )
382 *
383  CALL zunmbr( 'Q', 'L', 'N', m, n, n, cdum(1), m, cdum(1),
384  $ cdum(1), m, cdum(1), -1, ierr )
385  lwork_zunmbr_qln_mn = int( cdum(1) )
386 *
387  CALL zunmbr( 'Q', 'L', 'N', n, n, n, cdum(1), n, cdum(1),
388  $ cdum(1), n, cdum(1), -1, ierr )
389  lwork_zunmbr_qln_nn = int( cdum(1) )
390 *
391  IF( m.GE.mnthr1 ) THEN
392  IF( wntqn ) THEN
393 *
394 * Path 1 (M >> N, JOBZ='N')
395 *
396  maxwrk = n + lwork_zgeqrf_mn
397  maxwrk = max( maxwrk, 2*n + lwork_zgebrd_nn )
398  minwrk = 3*n
399  ELSE IF( wntqo ) THEN
400 *
401 * Path 2 (M >> N, JOBZ='O')
402 *
403  wrkbl = n + lwork_zgeqrf_mn
404  wrkbl = max( wrkbl, n + lwork_zungqr_mn )
405  wrkbl = max( wrkbl, 2*n + lwork_zgebrd_nn )
406  wrkbl = max( wrkbl, 2*n + lwork_zunmbr_qln_nn )
407  wrkbl = max( wrkbl, 2*n + lwork_zunmbr_prc_nn )
408  maxwrk = m*n + n*n + wrkbl
409  minwrk = 2*n*n + 3*n
410  ELSE IF( wntqs ) THEN
411 *
412 * Path 3 (M >> N, JOBZ='S')
413 *
414  wrkbl = n + lwork_zgeqrf_mn
415  wrkbl = max( wrkbl, n + lwork_zungqr_mn )
416  wrkbl = max( wrkbl, 2*n + lwork_zgebrd_nn )
417  wrkbl = max( wrkbl, 2*n + lwork_zunmbr_qln_nn )
418  wrkbl = max( wrkbl, 2*n + lwork_zunmbr_prc_nn )
419  maxwrk = n*n + wrkbl
420  minwrk = n*n + 3*n
421  ELSE IF( wntqa ) THEN
422 *
423 * Path 4 (M >> N, JOBZ='A')
424 *
425  wrkbl = n + lwork_zgeqrf_mn
426  wrkbl = max( wrkbl, n + lwork_zungqr_mm )
427  wrkbl = max( wrkbl, 2*n + lwork_zgebrd_nn )
428  wrkbl = max( wrkbl, 2*n + lwork_zunmbr_qln_nn )
429  wrkbl = max( wrkbl, 2*n + lwork_zunmbr_prc_nn )
430  maxwrk = n*n + wrkbl
431  minwrk = n*n + max( 3*n, n + m )
432  END IF
433  ELSE IF( m.GE.mnthr2 ) THEN
434 *
435 * Path 5 (M >> N, but not as much as MNTHR1)
436 *
437  maxwrk = 2*n + lwork_zgebrd_mn
438  minwrk = 2*n + m
439  IF( wntqo ) THEN
440 * Path 5o (M >> N, JOBZ='O')
441  maxwrk = max( maxwrk, 2*n + lwork_zungbr_p_nn )
442  maxwrk = max( maxwrk, 2*n + lwork_zungbr_q_mn )
443  maxwrk = maxwrk + m*n
444  minwrk = minwrk + n*n
445  ELSE IF( wntqs ) THEN
446 * Path 5s (M >> N, JOBZ='S')
447  maxwrk = max( maxwrk, 2*n + lwork_zungbr_p_nn )
448  maxwrk = max( maxwrk, 2*n + lwork_zungbr_q_mn )
449  ELSE IF( wntqa ) THEN
450 * Path 5a (M >> N, JOBZ='A')
451  maxwrk = max( maxwrk, 2*n + lwork_zungbr_p_nn )
452  maxwrk = max( maxwrk, 2*n + lwork_zungbr_q_mm )
453  END IF
454  ELSE
455 *
456 * Path 6 (M >= N, but not much larger)
457 *
458  maxwrk = 2*n + lwork_zgebrd_mn
459  minwrk = 2*n + m
460  IF( wntqo ) THEN
461 * Path 6o (M >= N, JOBZ='O')
462  maxwrk = max( maxwrk, 2*n + lwork_zunmbr_prc_nn )
463  maxwrk = max( maxwrk, 2*n + lwork_zunmbr_qln_mn )
464  maxwrk = maxwrk + m*n
465  minwrk = minwrk + n*n
466  ELSE IF( wntqs ) THEN
467 * Path 6s (M >= N, JOBZ='S')
468  maxwrk = max( maxwrk, 2*n + lwork_zunmbr_qln_mn )
469  maxwrk = max( maxwrk, 2*n + lwork_zunmbr_prc_nn )
470  ELSE IF( wntqa ) THEN
471 * Path 6a (M >= N, JOBZ='A')
472  maxwrk = max( maxwrk, 2*n + lwork_zunmbr_qln_mm )
473  maxwrk = max( maxwrk, 2*n + lwork_zunmbr_prc_nn )
474  END IF
475  END IF
476  ELSE IF( minmn.GT.0 ) THEN
477 *
478 * There is no complex work space needed for bidiagonal SVD
479 * The real work space needed for bidiagonal SVD (dbdsdc) is
480 * BDSPAC = 3*M*M + 4*M for singular values and vectors;
481 * BDSPAC = 4*M for singular values only;
482 * not including e, RU, and RVT matrices.
483 *
484 * Compute space preferred for each routine
485  CALL zgebrd( m, n, cdum(1), m, dum(1), dum(1), cdum(1),
486  $ cdum(1), cdum(1), -1, ierr )
487  lwork_zgebrd_mn = int( cdum(1) )
488 *
489  CALL zgebrd( m, m, cdum(1), m, dum(1), dum(1), cdum(1),
490  $ cdum(1), cdum(1), -1, ierr )
491  lwork_zgebrd_mm = int( cdum(1) )
492 *
493  CALL zgelqf( m, n, cdum(1), m, cdum(1), cdum(1), -1, ierr )
494  lwork_zgelqf_mn = int( cdum(1) )
495 *
496  CALL zungbr( 'P', m, n, m, cdum(1), m, cdum(1), cdum(1),
497  $ -1, ierr )
498  lwork_zungbr_p_mn = int( cdum(1) )
499 *
500  CALL zungbr( 'P', n, n, m, cdum(1), n, cdum(1), cdum(1),
501  $ -1, ierr )
502  lwork_zungbr_p_nn = int( cdum(1) )
503 *
504  CALL zungbr( 'Q', m, m, n, cdum(1), m, cdum(1), cdum(1),
505  $ -1, ierr )
506  lwork_zungbr_q_mm = int( cdum(1) )
507 *
508  CALL zunglq( m, n, m, cdum(1), m, cdum(1), cdum(1),
509  $ -1, ierr )
510  lwork_zunglq_mn = int( cdum(1) )
511 *
512  CALL zunglq( n, n, m, cdum(1), n, cdum(1), cdum(1),
513  $ -1, ierr )
514  lwork_zunglq_nn = int( cdum(1) )
515 *
516  CALL zunmbr( 'P', 'R', 'C', m, m, m, cdum(1), m, cdum(1),
517  $ cdum(1), m, cdum(1), -1, ierr )
518  lwork_zunmbr_prc_mm = int( cdum(1) )
519 *
520  CALL zunmbr( 'P', 'R', 'C', m, n, m, cdum(1), m, cdum(1),
521  $ cdum(1), m, cdum(1), -1, ierr )
522  lwork_zunmbr_prc_mn = int( cdum(1) )
523 *
524  CALL zunmbr( 'P', 'R', 'C', n, n, m, cdum(1), n, cdum(1),
525  $ cdum(1), n, cdum(1), -1, ierr )
526  lwork_zunmbr_prc_nn = int( cdum(1) )
527 *
528  CALL zunmbr( 'Q', 'L', 'N', m, m, m, cdum(1), m, cdum(1),
529  $ cdum(1), m, cdum(1), -1, ierr )
530  lwork_zunmbr_qln_mm = int( cdum(1) )
531 *
532  IF( n.GE.mnthr1 ) THEN
533  IF( wntqn ) THEN
534 *
535 * Path 1t (N >> M, JOBZ='N')
536 *
537  maxwrk = m + lwork_zgelqf_mn
538  maxwrk = max( maxwrk, 2*m + lwork_zgebrd_mm )
539  minwrk = 3*m
540  ELSE IF( wntqo ) THEN
541 *
542 * Path 2t (N >> M, JOBZ='O')
543 *
544  wrkbl = m + lwork_zgelqf_mn
545  wrkbl = max( wrkbl, m + lwork_zunglq_mn )
546  wrkbl = max( wrkbl, 2*m + lwork_zgebrd_mm )
547  wrkbl = max( wrkbl, 2*m + lwork_zunmbr_qln_mm )
548  wrkbl = max( wrkbl, 2*m + lwork_zunmbr_prc_mm )
549  maxwrk = m*n + m*m + wrkbl
550  minwrk = 2*m*m + 3*m
551  ELSE IF( wntqs ) THEN
552 *
553 * Path 3t (N >> M, JOBZ='S')
554 *
555  wrkbl = m + lwork_zgelqf_mn
556  wrkbl = max( wrkbl, m + lwork_zunglq_mn )
557  wrkbl = max( wrkbl, 2*m + lwork_zgebrd_mm )
558  wrkbl = max( wrkbl, 2*m + lwork_zunmbr_qln_mm )
559  wrkbl = max( wrkbl, 2*m + lwork_zunmbr_prc_mm )
560  maxwrk = m*m + wrkbl
561  minwrk = m*m + 3*m
562  ELSE IF( wntqa ) THEN
563 *
564 * Path 4t (N >> M, JOBZ='A')
565 *
566  wrkbl = m + lwork_zgelqf_mn
567  wrkbl = max( wrkbl, m + lwork_zunglq_nn )
568  wrkbl = max( wrkbl, 2*m + lwork_zgebrd_mm )
569  wrkbl = max( wrkbl, 2*m + lwork_zunmbr_qln_mm )
570  wrkbl = max( wrkbl, 2*m + lwork_zunmbr_prc_mm )
571  maxwrk = m*m + wrkbl
572  minwrk = m*m + max( 3*m, m + n )
573  END IF
574  ELSE IF( n.GE.mnthr2 ) THEN
575 *
576 * Path 5t (N >> M, but not as much as MNTHR1)
577 *
578  maxwrk = 2*m + lwork_zgebrd_mn
579  minwrk = 2*m + n
580  IF( wntqo ) THEN
581 * Path 5to (N >> M, JOBZ='O')
582  maxwrk = max( maxwrk, 2*m + lwork_zungbr_q_mm )
583  maxwrk = max( maxwrk, 2*m + lwork_zungbr_p_mn )
584  maxwrk = maxwrk + m*n
585  minwrk = minwrk + m*m
586  ELSE IF( wntqs ) THEN
587 * Path 5ts (N >> M, JOBZ='S')
588  maxwrk = max( maxwrk, 2*m + lwork_zungbr_q_mm )
589  maxwrk = max( maxwrk, 2*m + lwork_zungbr_p_mn )
590  ELSE IF( wntqa ) THEN
591 * Path 5ta (N >> M, JOBZ='A')
592  maxwrk = max( maxwrk, 2*m + lwork_zungbr_q_mm )
593  maxwrk = max( maxwrk, 2*m + lwork_zungbr_p_nn )
594  END IF
595  ELSE
596 *
597 * Path 6t (N > M, but not much larger)
598 *
599  maxwrk = 2*m + lwork_zgebrd_mn
600  minwrk = 2*m + n
601  IF( wntqo ) THEN
602 * Path 6to (N > M, JOBZ='O')
603  maxwrk = max( maxwrk, 2*m + lwork_zunmbr_qln_mm )
604  maxwrk = max( maxwrk, 2*m + lwork_zunmbr_prc_mn )
605  maxwrk = maxwrk + m*n
606  minwrk = minwrk + m*m
607  ELSE IF( wntqs ) THEN
608 * Path 6ts (N > M, JOBZ='S')
609  maxwrk = max( maxwrk, 2*m + lwork_zunmbr_qln_mm )
610  maxwrk = max( maxwrk, 2*m + lwork_zunmbr_prc_mn )
611  ELSE IF( wntqa ) THEN
612 * Path 6ta (N > M, JOBZ='A')
613  maxwrk = max( maxwrk, 2*m + lwork_zunmbr_qln_mm )
614  maxwrk = max( maxwrk, 2*m + lwork_zunmbr_prc_nn )
615  END IF
616  END IF
617  END IF
618  maxwrk = max( maxwrk, minwrk )
619  END IF
620  IF( info.EQ.0 ) THEN
621  work( 1 ) = droundup_lwork( maxwrk )
622  IF( lwork.LT.minwrk .AND. .NOT. lquery ) THEN
623  info = -12
624  END IF
625  END IF
626 *
627  IF( info.NE.0 ) THEN
628  CALL xerbla( 'ZGESDD', -info )
629  RETURN
630  ELSE IF( lquery ) THEN
631  RETURN
632  END IF
633 *
634 * Quick return if possible
635 *
636  IF( m.EQ.0 .OR. n.EQ.0 ) THEN
637  RETURN
638  END IF
639 *
640 * Get machine constants
641 *
642  eps = dlamch( 'P' )
643  smlnum = sqrt( dlamch( 'S' ) ) / eps
644  bignum = one / smlnum
645 *
646 * Scale A if max element outside range [SMLNUM,BIGNUM]
647 *
648  anrm = zlange( 'M', m, n, a, lda, dum )
649  IF( disnan( anrm ) ) THEN
650  info = -4
651  RETURN
652  END IF
653  iscl = 0
654  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
655  iscl = 1
656  CALL zlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, ierr )
657  ELSE IF( anrm.GT.bignum ) THEN
658  iscl = 1
659  CALL zlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, ierr )
660  END IF
661 *
662  IF( m.GE.n ) THEN
663 *
664 * A has at least as many rows as columns. If A has sufficiently
665 * more rows than columns, first reduce using the QR
666 * decomposition (if sufficient workspace available)
667 *
668  IF( m.GE.mnthr1 ) THEN
669 *
670  IF( wntqn ) THEN
671 *
672 * Path 1 (M >> N, JOBZ='N')
673 * No singular vectors to be computed
674 *
675  itau = 1
676  nwork = itau + n
677 *
678 * Compute A=Q*R
679 * CWorkspace: need N [tau] + N [work]
680 * CWorkspace: prefer N [tau] + N*NB [work]
681 * RWorkspace: need 0
682 *
683  CALL zgeqrf( m, n, a, lda, work( itau ), work( nwork ),
684  $ lwork-nwork+1, ierr )
685 *
686 * Zero out below R
687 *
688  CALL zlaset( 'L', n-1, n-1, czero, czero, a( 2, 1 ),
689  $ lda )
690  ie = 1
691  itauq = 1
692  itaup = itauq + n
693  nwork = itaup + n
694 *
695 * Bidiagonalize R in A
696 * CWorkspace: need 2*N [tauq, taup] + N [work]
697 * CWorkspace: prefer 2*N [tauq, taup] + 2*N*NB [work]
698 * RWorkspace: need N [e]
699 *
700  CALL zgebrd( n, n, a, lda, s, rwork( ie ), work( itauq ),
701  $ work( itaup ), work( nwork ), lwork-nwork+1,
702  $ ierr )
703  nrwork = ie + n
704 *
705 * Perform bidiagonal SVD, compute singular values only
706 * CWorkspace: need 0
707 * RWorkspace: need N [e] + BDSPAC
708 *
709  CALL dbdsdc( 'U', 'N', n, s, rwork( ie ), dum,1,dum,1,
710  $ dum, idum, rwork( nrwork ), iwork, info )
711 *
712  ELSE IF( wntqo ) THEN
713 *
714 * Path 2 (M >> N, JOBZ='O')
715 * N left singular vectors to be overwritten on A and
716 * N right singular vectors to be computed in VT
717 *
718  iu = 1
719 *
720 * WORK(IU) is N by N
721 *
722  ldwrku = n
723  ir = iu + ldwrku*n
724  IF( lwork .GE. m*n + n*n + 3*n ) THEN
725 *
726 * WORK(IR) is M by N
727 *
728  ldwrkr = m
729  ELSE
730  ldwrkr = ( lwork - n*n - 3*n ) / n
731  END IF
732  itau = ir + ldwrkr*n
733  nwork = itau + n
734 *
735 * Compute A=Q*R
736 * CWorkspace: need N*N [U] + N*N [R] + N [tau] + N [work]
737 * CWorkspace: prefer N*N [U] + N*N [R] + N [tau] + N*NB [work]
738 * RWorkspace: need 0
739 *
740  CALL zgeqrf( m, n, a, lda, work( itau ), work( nwork ),
741  $ lwork-nwork+1, ierr )
742 *
743 * Copy R to WORK( IR ), zeroing out below it
744 *
745  CALL zlacpy( 'U', n, n, a, lda, work( ir ), ldwrkr )
746  CALL zlaset( 'L', n-1, n-1, czero, czero, work( ir+1 ),
747  $ ldwrkr )
748 *
749 * Generate Q in A
750 * CWorkspace: need N*N [U] + N*N [R] + N [tau] + N [work]
751 * CWorkspace: prefer N*N [U] + N*N [R] + N [tau] + N*NB [work]
752 * RWorkspace: need 0
753 *
754  CALL zungqr( m, n, n, a, lda, work( itau ),
755  $ work( nwork ), lwork-nwork+1, ierr )
756  ie = 1
757  itauq = itau
758  itaup = itauq + n
759  nwork = itaup + n
760 *
761 * Bidiagonalize R in WORK(IR)
762 * CWorkspace: need N*N [U] + N*N [R] + 2*N [tauq, taup] + N [work]
763 * CWorkspace: prefer N*N [U] + N*N [R] + 2*N [tauq, taup] + 2*N*NB [work]
764 * RWorkspace: need N [e]
765 *
766  CALL zgebrd( n, n, work( ir ), ldwrkr, s, rwork( ie ),
767  $ work( itauq ), work( itaup ), work( nwork ),
768  $ lwork-nwork+1, ierr )
769 *
770 * Perform bidiagonal SVD, computing left singular vectors
771 * of R in WORK(IRU) and computing right singular vectors
772 * of R in WORK(IRVT)
773 * CWorkspace: need 0
774 * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
775 *
776  iru = ie + n
777  irvt = iru + n*n
778  nrwork = irvt + n*n
779  CALL dbdsdc( 'U', 'I', n, s, rwork( ie ), rwork( iru ),
780  $ n, rwork( irvt ), n, dum, idum,
781  $ rwork( nrwork ), iwork, info )
782 *
783 * Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
784 * Overwrite WORK(IU) by the left singular vectors of R
785 * CWorkspace: need N*N [U] + N*N [R] + 2*N [tauq, taup] + N [work]
786 * CWorkspace: prefer N*N [U] + N*N [R] + 2*N [tauq, taup] + N*NB [work]
787 * RWorkspace: need 0
788 *
789  CALL zlacp2( 'F', n, n, rwork( iru ), n, work( iu ),
790  $ ldwrku )
791  CALL zunmbr( 'Q', 'L', 'N', n, n, n, work( ir ), ldwrkr,
792  $ work( itauq ), work( iu ), ldwrku,
793  $ work( nwork ), lwork-nwork+1, ierr )
794 *
795 * Copy real matrix RWORK(IRVT) to complex matrix VT
796 * Overwrite VT by the right singular vectors of R
797 * CWorkspace: need N*N [U] + N*N [R] + 2*N [tauq, taup] + N [work]
798 * CWorkspace: prefer N*N [U] + N*N [R] + 2*N [tauq, taup] + N*NB [work]
799 * RWorkspace: need 0
800 *
801  CALL zlacp2( 'F', n, n, rwork( irvt ), n, vt, ldvt )
802  CALL zunmbr( 'P', 'R', 'C', n, n, n, work( ir ), ldwrkr,
803  $ work( itaup ), vt, ldvt, work( nwork ),
804  $ lwork-nwork+1, ierr )
805 *
806 * Multiply Q in A by left singular vectors of R in
807 * WORK(IU), storing result in WORK(IR) and copying to A
808 * CWorkspace: need N*N [U] + N*N [R]
809 * CWorkspace: prefer N*N [U] + M*N [R]
810 * RWorkspace: need 0
811 *
812  DO 10 i = 1, m, ldwrkr
813  chunk = min( m-i+1, ldwrkr )
814  CALL zgemm( 'N', 'N', chunk, n, n, cone, a( i, 1 ),
815  $ lda, work( iu ), ldwrku, czero,
816  $ work( ir ), ldwrkr )
817  CALL zlacpy( 'F', chunk, n, work( ir ), ldwrkr,
818  $ a( i, 1 ), lda )
819  10 CONTINUE
820 *
821  ELSE IF( wntqs ) THEN
822 *
823 * Path 3 (M >> N, JOBZ='S')
824 * N left singular vectors to be computed in U and
825 * N right singular vectors to be computed in VT
826 *
827  ir = 1
828 *
829 * WORK(IR) is N by N
830 *
831  ldwrkr = n
832  itau = ir + ldwrkr*n
833  nwork = itau + n
834 *
835 * Compute A=Q*R
836 * CWorkspace: need N*N [R] + N [tau] + N [work]
837 * CWorkspace: prefer N*N [R] + N [tau] + N*NB [work]
838 * RWorkspace: need 0
839 *
840  CALL zgeqrf( m, n, a, lda, work( itau ), work( nwork ),
841  $ lwork-nwork+1, ierr )
842 *
843 * Copy R to WORK(IR), zeroing out below it
844 *
845  CALL zlacpy( 'U', n, n, a, lda, work( ir ), ldwrkr )
846  CALL zlaset( 'L', n-1, n-1, czero, czero, work( ir+1 ),
847  $ ldwrkr )
848 *
849 * Generate Q in A
850 * CWorkspace: need N*N [R] + N [tau] + N [work]
851 * CWorkspace: prefer N*N [R] + N [tau] + N*NB [work]
852 * RWorkspace: need 0
853 *
854  CALL zungqr( m, n, n, a, lda, work( itau ),
855  $ work( nwork ), lwork-nwork+1, ierr )
856  ie = 1
857  itauq = itau
858  itaup = itauq + n
859  nwork = itaup + n
860 *
861 * Bidiagonalize R in WORK(IR)
862 * CWorkspace: need N*N [R] + 2*N [tauq, taup] + N [work]
863 * CWorkspace: prefer N*N [R] + 2*N [tauq, taup] + 2*N*NB [work]
864 * RWorkspace: need N [e]
865 *
866  CALL zgebrd( n, n, work( ir ), ldwrkr, s, rwork( ie ),
867  $ work( itauq ), work( itaup ), work( nwork ),
868  $ lwork-nwork+1, ierr )
869 *
870 * Perform bidiagonal SVD, computing left singular vectors
871 * of bidiagonal matrix in RWORK(IRU) and computing right
872 * singular vectors of bidiagonal matrix in RWORK(IRVT)
873 * CWorkspace: need 0
874 * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
875 *
876  iru = ie + n
877  irvt = iru + n*n
878  nrwork = irvt + n*n
879  CALL dbdsdc( 'U', 'I', n, s, rwork( ie ), rwork( iru ),
880  $ n, rwork( irvt ), n, dum, idum,
881  $ rwork( nrwork ), iwork, info )
882 *
883 * Copy real matrix RWORK(IRU) to complex matrix U
884 * Overwrite U by left singular vectors of R
885 * CWorkspace: need N*N [R] + 2*N [tauq, taup] + N [work]
886 * CWorkspace: prefer N*N [R] + 2*N [tauq, taup] + N*NB [work]
887 * RWorkspace: need 0
888 *
889  CALL zlacp2( 'F', n, n, rwork( iru ), n, u, ldu )
890  CALL zunmbr( 'Q', 'L', 'N', n, n, n, work( ir ), ldwrkr,
891  $ work( itauq ), u, ldu, work( nwork ),
892  $ lwork-nwork+1, ierr )
893 *
894 * Copy real matrix RWORK(IRVT) to complex matrix VT
895 * Overwrite VT by right singular vectors of R
896 * CWorkspace: need N*N [R] + 2*N [tauq, taup] + N [work]
897 * CWorkspace: prefer N*N [R] + 2*N [tauq, taup] + N*NB [work]
898 * RWorkspace: need 0
899 *
900  CALL zlacp2( 'F', n, n, rwork( irvt ), n, vt, ldvt )
901  CALL zunmbr( 'P', 'R', 'C', n, n, n, work( ir ), ldwrkr,
902  $ work( itaup ), vt, ldvt, work( nwork ),
903  $ lwork-nwork+1, ierr )
904 *
905 * Multiply Q in A by left singular vectors of R in
906 * WORK(IR), storing result in U
907 * CWorkspace: need N*N [R]
908 * RWorkspace: need 0
909 *
910  CALL zlacpy( 'F', n, n, u, ldu, work( ir ), ldwrkr )
911  CALL zgemm( 'N', 'N', m, n, n, cone, a, lda, work( ir ),
912  $ ldwrkr, czero, u, ldu )
913 *
914  ELSE IF( wntqa ) THEN
915 *
916 * Path 4 (M >> N, JOBZ='A')
917 * M left singular vectors to be computed in U and
918 * N right singular vectors to be computed in VT
919 *
920  iu = 1
921 *
922 * WORK(IU) is N by N
923 *
924  ldwrku = n
925  itau = iu + ldwrku*n
926  nwork = itau + n
927 *
928 * Compute A=Q*R, copying result to U
929 * CWorkspace: need N*N [U] + N [tau] + N [work]
930 * CWorkspace: prefer N*N [U] + N [tau] + N*NB [work]
931 * RWorkspace: need 0
932 *
933  CALL zgeqrf( m, n, a, lda, work( itau ), work( nwork ),
934  $ lwork-nwork+1, ierr )
935  CALL zlacpy( 'L', m, n, a, lda, u, ldu )
936 *
937 * Generate Q in U
938 * CWorkspace: need N*N [U] + N [tau] + M [work]
939 * CWorkspace: prefer N*N [U] + N [tau] + M*NB [work]
940 * RWorkspace: need 0
941 *
942  CALL zungqr( m, m, n, u, ldu, work( itau ),
943  $ work( nwork ), lwork-nwork+1, ierr )
944 *
945 * Produce R in A, zeroing out below it
946 *
947  CALL zlaset( 'L', n-1, n-1, czero, czero, a( 2, 1 ),
948  $ lda )
949  ie = 1
950  itauq = itau
951  itaup = itauq + n
952  nwork = itaup + n
953 *
954 * Bidiagonalize R in A
955 * CWorkspace: need N*N [U] + 2*N [tauq, taup] + N [work]
956 * CWorkspace: prefer N*N [U] + 2*N [tauq, taup] + 2*N*NB [work]
957 * RWorkspace: need N [e]
958 *
959  CALL zgebrd( n, n, a, lda, s, rwork( ie ), work( itauq ),
960  $ work( itaup ), work( nwork ), lwork-nwork+1,
961  $ ierr )
962  iru = ie + n
963  irvt = iru + n*n
964  nrwork = irvt + n*n
965 *
966 * Perform bidiagonal SVD, computing left singular vectors
967 * of bidiagonal matrix in RWORK(IRU) and computing right
968 * singular vectors of bidiagonal matrix in RWORK(IRVT)
969 * CWorkspace: need 0
970 * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
971 *
972  CALL dbdsdc( 'U', 'I', n, s, rwork( ie ), rwork( iru ),
973  $ n, rwork( irvt ), n, dum, idum,
974  $ rwork( nrwork ), iwork, info )
975 *
976 * Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
977 * Overwrite WORK(IU) by left singular vectors of R
978 * CWorkspace: need N*N [U] + 2*N [tauq, taup] + N [work]
979 * CWorkspace: prefer N*N [U] + 2*N [tauq, taup] + N*NB [work]
980 * RWorkspace: need 0
981 *
982  CALL zlacp2( 'F', n, n, rwork( iru ), n, work( iu ),
983  $ ldwrku )
984  CALL zunmbr( 'Q', 'L', 'N', n, n, n, a, lda,
985  $ work( itauq ), work( iu ), ldwrku,
986  $ work( nwork ), lwork-nwork+1, ierr )
987 *
988 * Copy real matrix RWORK(IRVT) to complex matrix VT
989 * Overwrite VT by right singular vectors of R
990 * CWorkspace: need N*N [U] + 2*N [tauq, taup] + N [work]
991 * CWorkspace: prefer N*N [U] + 2*N [tauq, taup] + N*NB [work]
992 * RWorkspace: need 0
993 *
994  CALL zlacp2( 'F', n, n, rwork( irvt ), n, vt, ldvt )
995  CALL zunmbr( 'P', 'R', 'C', n, n, n, a, lda,
996  $ work( itaup ), vt, ldvt, work( nwork ),
997  $ lwork-nwork+1, ierr )
998 *
999 * Multiply Q in U by left singular vectors of R in
1000 * WORK(IU), storing result in A
1001 * CWorkspace: need N*N [U]
1002 * RWorkspace: need 0
1003 *
1004  CALL zgemm( 'N', 'N', m, n, n, cone, u, ldu, work( iu ),
1005  $ ldwrku, czero, a, lda )
1006 *
1007 * Copy left singular vectors of A from A to U
1008 *
1009  CALL zlacpy( 'F', m, n, a, lda, u, ldu )
1010 *
1011  END IF
1012 *
1013  ELSE IF( m.GE.mnthr2 ) THEN
1014 *
1015 * MNTHR2 <= M < MNTHR1
1016 *
1017 * Path 5 (M >> N, but not as much as MNTHR1)
1018 * Reduce to bidiagonal form without QR decomposition, use
1019 * ZUNGBR and matrix multiplication to compute singular vectors
1020 *
1021  ie = 1
1022  nrwork = ie + n
1023  itauq = 1
1024  itaup = itauq + n
1025  nwork = itaup + n
1026 *
1027 * Bidiagonalize A
1028 * CWorkspace: need 2*N [tauq, taup] + M [work]
1029 * CWorkspace: prefer 2*N [tauq, taup] + (M+N)*NB [work]
1030 * RWorkspace: need N [e]
1031 *
1032  CALL zgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
1033  $ work( itaup ), work( nwork ), lwork-nwork+1,
1034  $ ierr )
1035  IF( wntqn ) THEN
1036 *
1037 * Path 5n (M >> N, JOBZ='N')
1038 * Compute singular values only
1039 * CWorkspace: need 0
1040 * RWorkspace: need N [e] + BDSPAC
1041 *
1042  CALL dbdsdc( 'U', 'N', n, s, rwork( ie ), dum, 1,dum,1,
1043  $ dum, idum, rwork( nrwork ), iwork, info )
1044  ELSE IF( wntqo ) THEN
1045  iu = nwork
1046  iru = nrwork
1047  irvt = iru + n*n
1048  nrwork = irvt + n*n
1049 *
1050 * Path 5o (M >> N, JOBZ='O')
1051 * Copy A to VT, generate P**H
1052 * CWorkspace: need 2*N [tauq, taup] + N [work]
1053 * CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1054 * RWorkspace: need 0
1055 *
1056  CALL zlacpy( 'U', n, n, a, lda, vt, ldvt )
1057  CALL zungbr( 'P', n, n, n, vt, ldvt, work( itaup ),
1058  $ work( nwork ), lwork-nwork+1, ierr )
1059 *
1060 * Generate Q in A
1061 * CWorkspace: need 2*N [tauq, taup] + N [work]
1062 * CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1063 * RWorkspace: need 0
1064 *
1065  CALL zungbr( 'Q', m, n, n, a, lda, work( itauq ),
1066  $ work( nwork ), lwork-nwork+1, ierr )
1067 *
1068  IF( lwork .GE. m*n + 3*n ) THEN
1069 *
1070 * WORK( IU ) is M by N
1071 *
1072  ldwrku = m
1073  ELSE
1074 *
1075 * WORK(IU) is LDWRKU by N
1076 *
1077  ldwrku = ( lwork - 3*n ) / n
1078  END IF
1079  nwork = iu + ldwrku*n
1080 *
1081 * Perform bidiagonal SVD, computing left singular vectors
1082 * of bidiagonal matrix in RWORK(IRU) and computing right
1083 * singular vectors of bidiagonal matrix in RWORK(IRVT)
1084 * CWorkspace: need 0
1085 * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
1086 *
1087  CALL dbdsdc( 'U', 'I', n, s, rwork( ie ), rwork( iru ),
1088  $ n, rwork( irvt ), n, dum, idum,
1089  $ rwork( nrwork ), iwork, info )
1090 *
1091 * Multiply real matrix RWORK(IRVT) by P**H in VT,
1092 * storing the result in WORK(IU), copying to VT
1093 * CWorkspace: need 2*N [tauq, taup] + N*N [U]
1094 * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + 2*N*N [rwork]
1095 *
1096  CALL zlarcm( n, n, rwork( irvt ), n, vt, ldvt,
1097  $ work( iu ), ldwrku, rwork( nrwork ) )
1098  CALL zlacpy( 'F', n, n, work( iu ), ldwrku, vt, ldvt )
1099 *
1100 * Multiply Q in A by real matrix RWORK(IRU), storing the
1101 * result in WORK(IU), copying to A
1102 * CWorkspace: need 2*N [tauq, taup] + N*N [U]
1103 * CWorkspace: prefer 2*N [tauq, taup] + M*N [U]
1104 * RWorkspace: need N [e] + N*N [RU] + 2*N*N [rwork]
1105 * RWorkspace: prefer N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here
1106 *
1107  nrwork = irvt
1108  DO 20 i = 1, m, ldwrku
1109  chunk = min( m-i+1, ldwrku )
1110  CALL zlacrm( chunk, n, a( i, 1 ), lda, rwork( iru ),
1111  $ n, work( iu ), ldwrku, rwork( nrwork ) )
1112  CALL zlacpy( 'F', chunk, n, work( iu ), ldwrku,
1113  $ a( i, 1 ), lda )
1114  20 CONTINUE
1115 *
1116  ELSE IF( wntqs ) THEN
1117 *
1118 * Path 5s (M >> N, JOBZ='S')
1119 * Copy A to VT, generate P**H
1120 * CWorkspace: need 2*N [tauq, taup] + N [work]
1121 * CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1122 * RWorkspace: need 0
1123 *
1124  CALL zlacpy( 'U', n, n, a, lda, vt, ldvt )
1125  CALL zungbr( 'P', n, n, n, vt, ldvt, work( itaup ),
1126  $ work( nwork ), lwork-nwork+1, ierr )
1127 *
1128 * Copy A to U, generate Q
1129 * CWorkspace: need 2*N [tauq, taup] + N [work]
1130 * CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1131 * RWorkspace: need 0
1132 *
1133  CALL zlacpy( 'L', m, n, a, lda, u, ldu )
1134  CALL zungbr( 'Q', m, n, n, u, ldu, work( itauq ),
1135  $ work( nwork ), lwork-nwork+1, ierr )
1136 *
1137 * Perform bidiagonal SVD, computing left singular vectors
1138 * of bidiagonal matrix in RWORK(IRU) and computing right
1139 * singular vectors of bidiagonal matrix in RWORK(IRVT)
1140 * CWorkspace: need 0
1141 * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
1142 *
1143  iru = nrwork
1144  irvt = iru + n*n
1145  nrwork = irvt + n*n
1146  CALL dbdsdc( 'U', 'I', n, s, rwork( ie ), rwork( iru ),
1147  $ n, rwork( irvt ), n, dum, idum,
1148  $ rwork( nrwork ), iwork, info )
1149 *
1150 * Multiply real matrix RWORK(IRVT) by P**H in VT,
1151 * storing the result in A, copying to VT
1152 * CWorkspace: need 0
1153 * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + 2*N*N [rwork]
1154 *
1155  CALL zlarcm( n, n, rwork( irvt ), n, vt, ldvt, a, lda,
1156  $ rwork( nrwork ) )
1157  CALL zlacpy( 'F', n, n, a, lda, vt, ldvt )
1158 *
1159 * Multiply Q in U by real matrix RWORK(IRU), storing the
1160 * result in A, copying to U
1161 * CWorkspace: need 0
1162 * RWorkspace: need N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here
1163 *
1164  nrwork = irvt
1165  CALL zlacrm( m, n, u, ldu, rwork( iru ), n, a, lda,
1166  $ rwork( nrwork ) )
1167  CALL zlacpy( 'F', m, n, a, lda, u, ldu )
1168  ELSE
1169 *
1170 * Path 5a (M >> N, JOBZ='A')
1171 * Copy A to VT, generate P**H
1172 * CWorkspace: need 2*N [tauq, taup] + N [work]
1173 * CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1174 * RWorkspace: need 0
1175 *
1176  CALL zlacpy( 'U', n, n, a, lda, vt, ldvt )
1177  CALL zungbr( 'P', n, n, n, vt, ldvt, work( itaup ),
1178  $ work( nwork ), lwork-nwork+1, ierr )
1179 *
1180 * Copy A to U, generate Q
1181 * CWorkspace: need 2*N [tauq, taup] + M [work]
1182 * CWorkspace: prefer 2*N [tauq, taup] + M*NB [work]
1183 * RWorkspace: need 0
1184 *
1185  CALL zlacpy( 'L', m, n, a, lda, u, ldu )
1186  CALL zungbr( 'Q', m, m, n, u, ldu, work( itauq ),
1187  $ work( nwork ), lwork-nwork+1, ierr )
1188 *
1189 * Perform bidiagonal SVD, computing left singular vectors
1190 * of bidiagonal matrix in RWORK(IRU) and computing right
1191 * singular vectors of bidiagonal matrix in RWORK(IRVT)
1192 * CWorkspace: need 0
1193 * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
1194 *
1195  iru = nrwork
1196  irvt = iru + n*n
1197  nrwork = irvt + n*n
1198  CALL dbdsdc( 'U', 'I', n, s, rwork( ie ), rwork( iru ),
1199  $ n, rwork( irvt ), n, dum, idum,
1200  $ rwork( nrwork ), iwork, info )
1201 *
1202 * Multiply real matrix RWORK(IRVT) by P**H in VT,
1203 * storing the result in A, copying to VT
1204 * CWorkspace: need 0
1205 * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + 2*N*N [rwork]
1206 *
1207  CALL zlarcm( n, n, rwork( irvt ), n, vt, ldvt, a, lda,
1208  $ rwork( nrwork ) )
1209  CALL zlacpy( 'F', n, n, a, lda, vt, ldvt )
1210 *
1211 * Multiply Q in U by real matrix RWORK(IRU), storing the
1212 * result in A, copying to U
1213 * CWorkspace: need 0
1214 * RWorkspace: need N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here
1215 *
1216  nrwork = irvt
1217  CALL zlacrm( m, n, u, ldu, rwork( iru ), n, a, lda,
1218  $ rwork( nrwork ) )
1219  CALL zlacpy( 'F', m, n, a, lda, u, ldu )
1220  END IF
1221 *
1222  ELSE
1223 *
1224 * M .LT. MNTHR2
1225 *
1226 * Path 6 (M >= N, but not much larger)
1227 * Reduce to bidiagonal form without QR decomposition
1228 * Use ZUNMBR to compute singular vectors
1229 *
1230  ie = 1
1231  nrwork = ie + n
1232  itauq = 1
1233  itaup = itauq + n
1234  nwork = itaup + n
1235 *
1236 * Bidiagonalize A
1237 * CWorkspace: need 2*N [tauq, taup] + M [work]
1238 * CWorkspace: prefer 2*N [tauq, taup] + (M+N)*NB [work]
1239 * RWorkspace: need N [e]
1240 *
1241  CALL zgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
1242  $ work( itaup ), work( nwork ), lwork-nwork+1,
1243  $ ierr )
1244  IF( wntqn ) THEN
1245 *
1246 * Path 6n (M >= N, JOBZ='N')
1247 * Compute singular values only
1248 * CWorkspace: need 0
1249 * RWorkspace: need N [e] + BDSPAC
1250 *
1251  CALL dbdsdc( 'U', 'N', n, s, rwork( ie ), dum,1,dum,1,
1252  $ dum, idum, rwork( nrwork ), iwork, info )
1253  ELSE IF( wntqo ) THEN
1254  iu = nwork
1255  iru = nrwork
1256  irvt = iru + n*n
1257  nrwork = irvt + n*n
1258  IF( lwork .GE. m*n + 3*n ) THEN
1259 *
1260 * WORK( IU ) is M by N
1261 *
1262  ldwrku = m
1263  ELSE
1264 *
1265 * WORK( IU ) is LDWRKU by N
1266 *
1267  ldwrku = ( lwork - 3*n ) / n
1268  END IF
1269  nwork = iu + ldwrku*n
1270 *
1271 * Path 6o (M >= N, JOBZ='O')
1272 * Perform bidiagonal SVD, computing left singular vectors
1273 * of bidiagonal matrix in RWORK(IRU) and computing right
1274 * singular vectors of bidiagonal matrix in RWORK(IRVT)
1275 * CWorkspace: need 0
1276 * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
1277 *
1278  CALL dbdsdc( 'U', 'I', n, s, rwork( ie ), rwork( iru ),
1279  $ n, rwork( irvt ), n, dum, idum,
1280  $ rwork( nrwork ), iwork, info )
1281 *
1282 * Copy real matrix RWORK(IRVT) to complex matrix VT
1283 * Overwrite VT by right singular vectors of A
1284 * CWorkspace: need 2*N [tauq, taup] + N*N [U] + N [work]
1285 * CWorkspace: prefer 2*N [tauq, taup] + N*N [U] + N*NB [work]
1286 * RWorkspace: need N [e] + N*N [RU] + N*N [RVT]
1287 *
1288  CALL zlacp2( 'F', n, n, rwork( irvt ), n, vt, ldvt )
1289  CALL zunmbr( 'P', 'R', 'C', n, n, n, a, lda,
1290  $ work( itaup ), vt, ldvt, work( nwork ),
1291  $ lwork-nwork+1, ierr )
1292 *
1293  IF( lwork .GE. m*n + 3*n ) THEN
1294 *
1295 * Path 6o-fast
1296 * Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
1297 * Overwrite WORK(IU) by left singular vectors of A, copying
1298 * to A
1299 * CWorkspace: need 2*N [tauq, taup] + M*N [U] + N [work]
1300 * CWorkspace: prefer 2*N [tauq, taup] + M*N [U] + N*NB [work]
1301 * RWorkspace: need N [e] + N*N [RU]
1302 *
1303  CALL zlaset( 'F', m, n, czero, czero, work( iu ),
1304  $ ldwrku )
1305  CALL zlacp2( 'F', n, n, rwork( iru ), n, work( iu ),
1306  $ ldwrku )
1307  CALL zunmbr( 'Q', 'L', 'N', m, n, n, a, lda,
1308  $ work( itauq ), work( iu ), ldwrku,
1309  $ work( nwork ), lwork-nwork+1, ierr )
1310  CALL zlacpy( 'F', m, n, work( iu ), ldwrku, a, lda )
1311  ELSE
1312 *
1313 * Path 6o-slow
1314 * Generate Q in A
1315 * CWorkspace: need 2*N [tauq, taup] + N*N [U] + N [work]
1316 * CWorkspace: prefer 2*N [tauq, taup] + N*N [U] + N*NB [work]
1317 * RWorkspace: need 0
1318 *
1319  CALL zungbr( 'Q', m, n, n, a, lda, work( itauq ),
1320  $ work( nwork ), lwork-nwork+1, ierr )
1321 *
1322 * Multiply Q in A by real matrix RWORK(IRU), storing the
1323 * result in WORK(IU), copying to A
1324 * CWorkspace: need 2*N [tauq, taup] + N*N [U]
1325 * CWorkspace: prefer 2*N [tauq, taup] + M*N [U]
1326 * RWorkspace: need N [e] + N*N [RU] + 2*N*N [rwork]
1327 * RWorkspace: prefer N [e] + N*N [RU] + 2*M*N [rwork] < N + 5*N*N since M < 2*N here
1328 *
1329  nrwork = irvt
1330  DO 30 i = 1, m, ldwrku
1331  chunk = min( m-i+1, ldwrku )
1332  CALL zlacrm( chunk, n, a( i, 1 ), lda,
1333  $ rwork( iru ), n, work( iu ), ldwrku,
1334  $ rwork( nrwork ) )
1335  CALL zlacpy( 'F', chunk, n, work( iu ), ldwrku,
1336  $ a( i, 1 ), lda )
1337  30 CONTINUE
1338  END IF
1339 *
1340  ELSE IF( wntqs ) THEN
1341 *
1342 * Path 6s (M >= N, JOBZ='S')
1343 * Perform bidiagonal SVD, computing left singular vectors
1344 * of bidiagonal matrix in RWORK(IRU) and computing right
1345 * singular vectors of bidiagonal matrix in RWORK(IRVT)
1346 * CWorkspace: need 0
1347 * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
1348 *
1349  iru = nrwork
1350  irvt = iru + n*n
1351  nrwork = irvt + n*n
1352  CALL dbdsdc( 'U', 'I', n, s, rwork( ie ), rwork( iru ),
1353  $ n, rwork( irvt ), n, dum, idum,
1354  $ rwork( nrwork ), iwork, info )
1355 *
1356 * Copy real matrix RWORK(IRU) to complex matrix U
1357 * Overwrite U by left singular vectors of A
1358 * CWorkspace: need 2*N [tauq, taup] + N [work]
1359 * CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1360 * RWorkspace: need N [e] + N*N [RU] + N*N [RVT]
1361 *
1362  CALL zlaset( 'F', m, n, czero, czero, u, ldu )
1363  CALL zlacp2( 'F', n, n, rwork( iru ), n, u, ldu )
1364  CALL zunmbr( 'Q', 'L', 'N', m, n, n, a, lda,
1365  $ work( itauq ), u, ldu, work( nwork ),
1366  $ lwork-nwork+1, ierr )
1367 *
1368 * Copy real matrix RWORK(IRVT) to complex matrix VT
1369 * Overwrite VT by right singular vectors of A
1370 * CWorkspace: need 2*N [tauq, taup] + N [work]
1371 * CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1372 * RWorkspace: need N [e] + N*N [RU] + N*N [RVT]
1373 *
1374  CALL zlacp2( 'F', n, n, rwork( irvt ), n, vt, ldvt )
1375  CALL zunmbr( 'P', 'R', 'C', n, n, n, a, lda,
1376  $ work( itaup ), vt, ldvt, work( nwork ),
1377  $ lwork-nwork+1, ierr )
1378  ELSE
1379 *
1380 * Path 6a (M >= N, JOBZ='A')
1381 * Perform bidiagonal SVD, computing left singular vectors
1382 * of bidiagonal matrix in RWORK(IRU) and computing right
1383 * singular vectors of bidiagonal matrix in RWORK(IRVT)
1384 * CWorkspace: need 0
1385 * RWorkspace: need N [e] + N*N [RU] + N*N [RVT] + BDSPAC
1386 *
1387  iru = nrwork
1388  irvt = iru + n*n
1389  nrwork = irvt + n*n
1390  CALL dbdsdc( 'U', 'I', n, s, rwork( ie ), rwork( iru ),
1391  $ n, rwork( irvt ), n, dum, idum,
1392  $ rwork( nrwork ), iwork, info )
1393 *
1394 * Set the right corner of U to identity matrix
1395 *
1396  CALL zlaset( 'F', m, m, czero, czero, u, ldu )
1397  IF( m.GT.n ) THEN
1398  CALL zlaset( 'F', m-n, m-n, czero, cone,
1399  $ u( n+1, n+1 ), ldu )
1400  END IF
1401 *
1402 * Copy real matrix RWORK(IRU) to complex matrix U
1403 * Overwrite U by left singular vectors of A
1404 * CWorkspace: need 2*N [tauq, taup] + M [work]
1405 * CWorkspace: prefer 2*N [tauq, taup] + M*NB [work]
1406 * RWorkspace: need N [e] + N*N [RU] + N*N [RVT]
1407 *
1408  CALL zlacp2( 'F', n, n, rwork( iru ), n, u, ldu )
1409  CALL zunmbr( 'Q', 'L', 'N', m, m, n, a, lda,
1410  $ work( itauq ), u, ldu, work( nwork ),
1411  $ lwork-nwork+1, ierr )
1412 *
1413 * Copy real matrix RWORK(IRVT) to complex matrix VT
1414 * Overwrite VT by right singular vectors of A
1415 * CWorkspace: need 2*N [tauq, taup] + N [work]
1416 * CWorkspace: prefer 2*N [tauq, taup] + N*NB [work]
1417 * RWorkspace: need N [e] + N*N [RU] + N*N [RVT]
1418 *
1419  CALL zlacp2( 'F', n, n, rwork( irvt ), n, vt, ldvt )
1420  CALL zunmbr( 'P', 'R', 'C', n, n, n, a, lda,
1421  $ work( itaup ), vt, ldvt, work( nwork ),
1422  $ lwork-nwork+1, ierr )
1423  END IF
1424 *
1425  END IF
1426 *
1427  ELSE
1428 *
1429 * A has more columns than rows. If A has sufficiently more
1430 * columns than rows, first reduce using the LQ decomposition (if
1431 * sufficient workspace available)
1432 *
1433  IF( n.GE.mnthr1 ) THEN
1434 *
1435  IF( wntqn ) THEN
1436 *
1437 * Path 1t (N >> M, JOBZ='N')
1438 * No singular vectors to be computed
1439 *
1440  itau = 1
1441  nwork = itau + m
1442 *
1443 * Compute A=L*Q
1444 * CWorkspace: need M [tau] + M [work]
1445 * CWorkspace: prefer M [tau] + M*NB [work]
1446 * RWorkspace: need 0
1447 *
1448  CALL zgelqf( m, n, a, lda, work( itau ), work( nwork ),
1449  $ lwork-nwork+1, ierr )
1450 *
1451 * Zero out above L
1452 *
1453  CALL zlaset( 'U', m-1, m-1, czero, czero, a( 1, 2 ),
1454  $ lda )
1455  ie = 1
1456  itauq = 1
1457  itaup = itauq + m
1458  nwork = itaup + m
1459 *
1460 * Bidiagonalize L in A
1461 * CWorkspace: need 2*M [tauq, taup] + M [work]
1462 * CWorkspace: prefer 2*M [tauq, taup] + 2*M*NB [work]
1463 * RWorkspace: need M [e]
1464 *
1465  CALL zgebrd( m, m, a, lda, s, rwork( ie ), work( itauq ),
1466  $ work( itaup ), work( nwork ), lwork-nwork+1,
1467  $ ierr )
1468  nrwork = ie + m
1469 *
1470 * Perform bidiagonal SVD, compute singular values only
1471 * CWorkspace: need 0
1472 * RWorkspace: need M [e] + BDSPAC
1473 *
1474  CALL dbdsdc( 'U', 'N', m, s, rwork( ie ), dum,1,dum,1,
1475  $ dum, idum, rwork( nrwork ), iwork, info )
1476 *
1477  ELSE IF( wntqo ) THEN
1478 *
1479 * Path 2t (N >> M, JOBZ='O')
1480 * M right singular vectors to be overwritten on A and
1481 * M left singular vectors to be computed in U
1482 *
1483  ivt = 1
1484  ldwkvt = m
1485 *
1486 * WORK(IVT) is M by M
1487 *
1488  il = ivt + ldwkvt*m
1489  IF( lwork .GE. m*n + m*m + 3*m ) THEN
1490 *
1491 * WORK(IL) M by N
1492 *
1493  ldwrkl = m
1494  chunk = n
1495  ELSE
1496 *
1497 * WORK(IL) is M by CHUNK
1498 *
1499  ldwrkl = m
1500  chunk = ( lwork - m*m - 3*m ) / m
1501  END IF
1502  itau = il + ldwrkl*chunk
1503  nwork = itau + m
1504 *
1505 * Compute A=L*Q
1506 * CWorkspace: need M*M [VT] + M*M [L] + M [tau] + M [work]
1507 * CWorkspace: prefer M*M [VT] + M*M [L] + M [tau] + M*NB [work]
1508 * RWorkspace: need 0
1509 *
1510  CALL zgelqf( m, n, a, lda, work( itau ), work( nwork ),
1511  $ lwork-nwork+1, ierr )
1512 *
1513 * Copy L to WORK(IL), zeroing about above it
1514 *
1515  CALL zlacpy( 'L', m, m, a, lda, work( il ), ldwrkl )
1516  CALL zlaset( 'U', m-1, m-1, czero, czero,
1517  $ work( il+ldwrkl ), ldwrkl )
1518 *
1519 * Generate Q in A
1520 * CWorkspace: need M*M [VT] + M*M [L] + M [tau] + M [work]
1521 * CWorkspace: prefer M*M [VT] + M*M [L] + M [tau] + M*NB [work]
1522 * RWorkspace: need 0
1523 *
1524  CALL zunglq( m, n, m, a, lda, work( itau ),
1525  $ work( nwork ), lwork-nwork+1, ierr )
1526  ie = 1
1527  itauq = itau
1528  itaup = itauq + m
1529  nwork = itaup + m
1530 *
1531 * Bidiagonalize L in WORK(IL)
1532 * CWorkspace: need M*M [VT] + M*M [L] + 2*M [tauq, taup] + M [work]
1533 * CWorkspace: prefer M*M [VT] + M*M [L] + 2*M [tauq, taup] + 2*M*NB [work]
1534 * RWorkspace: need M [e]
1535 *
1536  CALL zgebrd( m, m, work( il ), ldwrkl, s, rwork( ie ),
1537  $ work( itauq ), work( itaup ), work( nwork ),
1538  $ lwork-nwork+1, ierr )
1539 *
1540 * Perform bidiagonal SVD, computing left singular vectors
1541 * of bidiagonal matrix in RWORK(IRU) and computing right
1542 * singular vectors of bidiagonal matrix in RWORK(IRVT)
1543 * CWorkspace: need 0
1544 * RWorkspace: need M [e] + M*M [RU] + M*M [RVT] + BDSPAC
1545 *
1546  iru = ie + m
1547  irvt = iru + m*m
1548  nrwork = irvt + m*m
1549  CALL dbdsdc( 'U', 'I', m, s, rwork( ie ), rwork( iru ),
1550  $ m, rwork( irvt ), m, dum, idum,
1551  $ rwork( nrwork ), iwork, info )
1552 *
1553 * Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
1554 * Overwrite WORK(IU) by the left singular vectors of L
1555 * CWorkspace: need M*M [VT] + M*M [L] + 2*M [tauq, taup] + M [work]
1556 * CWorkspace: prefer M*M [VT] + M*M [L] + 2*M [tauq, taup] + M*NB [work]
1557 * RWorkspace: need 0
1558 *
1559  CALL zlacp2( 'F', m, m, rwork( iru ), m, u, ldu )
1560  CALL zunmbr( 'Q', 'L', 'N', m, m, m, work( il ), ldwrkl,
1561  $ work( itauq ), u, ldu, work( nwork ),
1562  $ lwork-nwork+1, ierr )
1563 *
1564 * Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
1565 * Overwrite WORK(IVT) by the right singular vectors of L
1566 * CWorkspace: need M*M [VT] + M*M [L] + 2*M [tauq, taup] + M [work]
1567 * CWorkspace: prefer M*M [VT] + M*M [L] + 2*M [tauq, taup] + M*NB [work]
1568 * RWorkspace: need 0
1569 *
1570  CALL zlacp2( 'F', m, m, rwork( irvt ), m, work( ivt ),
1571  $ ldwkvt )
1572  CALL zunmbr( 'P', 'R', 'C', m, m, m, work( il ), ldwrkl,
1573  $ work( itaup ), work( ivt ), ldwkvt,
1574  $ work( nwork ), lwork-nwork+1, ierr )
1575 *
1576 * Multiply right singular vectors of L in WORK(IL) by Q
1577 * in A, storing result in WORK(IL) and copying to A
1578 * CWorkspace: need M*M [VT] + M*M [L]
1579 * CWorkspace: prefer M*M [VT] + M*N [L]
1580 * RWorkspace: need 0
1581 *
1582  DO 40 i = 1, n, chunk
1583  blk = min( n-i+1, chunk )
1584  CALL zgemm( 'N', 'N', m, blk, m, cone, work( ivt ), m,
1585  $ a( 1, i ), lda, czero, work( il ),
1586  $ ldwrkl )
1587  CALL zlacpy( 'F', m, blk, work( il ), ldwrkl,
1588  $ a( 1, i ), lda )
1589  40 CONTINUE
1590 *
1591  ELSE IF( wntqs ) THEN
1592 *
1593 * Path 3t (N >> M, JOBZ='S')
1594 * M right singular vectors to be computed in VT and
1595 * M left singular vectors to be computed in U
1596 *
1597  il = 1
1598 *
1599 * WORK(IL) is M by M
1600 *
1601  ldwrkl = m
1602  itau = il + ldwrkl*m
1603  nwork = itau + m
1604 *
1605 * Compute A=L*Q
1606 * CWorkspace: need M*M [L] + M [tau] + M [work]
1607 * CWorkspace: prefer M*M [L] + M [tau] + M*NB [work]
1608 * RWorkspace: need 0
1609 *
1610  CALL zgelqf( m, n, a, lda, work( itau ), work( nwork ),
1611  $ lwork-nwork+1, ierr )
1612 *
1613 * Copy L to WORK(IL), zeroing out above it
1614 *
1615  CALL zlacpy( 'L', m, m, a, lda, work( il ), ldwrkl )
1616  CALL zlaset( 'U', m-1, m-1, czero, czero,
1617  $ work( il+ldwrkl ), ldwrkl )
1618 *
1619 * Generate Q in A
1620 * CWorkspace: need M*M [L] + M [tau] + M [work]
1621 * CWorkspace: prefer M*M [L] + M [tau] + M*NB [work]
1622 * RWorkspace: need 0
1623 *
1624  CALL zunglq( m, n, m, a, lda, work( itau ),
1625  $ work( nwork ), lwork-nwork+1, ierr )
1626  ie = 1
1627  itauq = itau
1628  itaup = itauq + m
1629  nwork = itaup + m
1630 *
1631 * Bidiagonalize L in WORK(IL)
1632 * CWorkspace: need M*M [L] + 2*M [tauq, taup] + M [work]
1633 * CWorkspace: prefer M*M [L] + 2*M [tauq, taup] + 2*M*NB [work]
1634 * RWorkspace: need M [e]
1635 *
1636  CALL zgebrd( m, m, work( il ), ldwrkl, s, rwork( ie ),
1637  $ work( itauq ), work( itaup ), work( nwork ),
1638  $ lwork-nwork+1, ierr )
1639 *
1640 * Perform bidiagonal SVD, computing left singular vectors
1641 * of bidiagonal matrix in RWORK(IRU) and computing right
1642 * singular vectors of bidiagonal matrix in RWORK(IRVT)
1643 * CWorkspace: need 0
1644 * RWorkspace: need M [e] + M*M [RU] + M*M [RVT] + BDSPAC
1645 *
1646  iru = ie + m
1647  irvt = iru + m*m
1648  nrwork = irvt + m*m
1649  CALL dbdsdc( 'U', 'I', m, s, rwork( ie ), rwork( iru ),
1650  $ m, rwork( irvt ), m, dum, idum,
1651  $ rwork( nrwork ), iwork, info )
1652 *
1653 * Copy real matrix RWORK(IRU) to complex matrix U
1654 * Overwrite U by left singular vectors of L
1655 * CWorkspace: need M*M [L] + 2*M [tauq, taup] + M [work]
1656 * CWorkspace: prefer M*M [L] + 2*M [tauq, taup] + M*NB [work]
1657 * RWorkspace: need 0
1658 *
1659  CALL zlacp2( 'F', m, m, rwork( iru ), m, u, ldu )
1660  CALL zunmbr( 'Q', 'L', 'N', m, m, m, work( il ), ldwrkl,
1661  $ work( itauq ), u, ldu, work( nwork ),
1662  $ lwork-nwork+1, ierr )
1663 *
1664 * Copy real matrix RWORK(IRVT) to complex matrix VT
1665 * Overwrite VT by left singular vectors of L
1666 * CWorkspace: need M*M [L] + 2*M [tauq, taup] + M [work]
1667 * CWorkspace: prefer M*M [L] + 2*M [tauq, taup] + M*NB [work]
1668 * RWorkspace: need 0
1669 *
1670  CALL zlacp2( 'F', m, m, rwork( irvt ), m, vt, ldvt )
1671  CALL zunmbr( 'P', 'R', 'C', m, m, m, work( il ), ldwrkl,
1672  $ work( itaup ), vt, ldvt, work( nwork ),
1673  $ lwork-nwork+1, ierr )
1674 *
1675 * Copy VT to WORK(IL), multiply right singular vectors of L
1676 * in WORK(IL) by Q in A, storing result in VT
1677 * CWorkspace: need M*M [L]
1678 * RWorkspace: need 0
1679 *
1680  CALL zlacpy( 'F', m, m, vt, ldvt, work( il ), ldwrkl )
1681  CALL zgemm( 'N', 'N', m, n, m, cone, work( il ), ldwrkl,
1682  $ a, lda, czero, vt, ldvt )
1683 *
1684  ELSE IF( wntqa ) THEN
1685 *
1686 * Path 4t (N >> M, JOBZ='A')
1687 * N right singular vectors to be computed in VT and
1688 * M left singular vectors to be computed in U
1689 *
1690  ivt = 1
1691 *
1692 * WORK(IVT) is M by M
1693 *
1694  ldwkvt = m
1695  itau = ivt + ldwkvt*m
1696  nwork = itau + m
1697 *
1698 * Compute A=L*Q, copying result to VT
1699 * CWorkspace: need M*M [VT] + M [tau] + M [work]
1700 * CWorkspace: prefer M*M [VT] + M [tau] + M*NB [work]
1701 * RWorkspace: need 0
1702 *
1703  CALL zgelqf( m, n, a, lda, work( itau ), work( nwork ),
1704  $ lwork-nwork+1, ierr )
1705  CALL zlacpy( 'U', m, n, a, lda, vt, ldvt )
1706 *
1707 * Generate Q in VT
1708 * CWorkspace: need M*M [VT] + M [tau] + N [work]
1709 * CWorkspace: prefer M*M [VT] + M [tau] + N*NB [work]
1710 * RWorkspace: need 0
1711 *
1712  CALL zunglq( n, n, m, vt, ldvt, work( itau ),
1713  $ work( nwork ), lwork-nwork+1, ierr )
1714 *
1715 * Produce L in A, zeroing out above it
1716 *
1717  CALL zlaset( 'U', m-1, m-1, czero, czero, a( 1, 2 ),
1718  $ lda )
1719  ie = 1
1720  itauq = itau
1721  itaup = itauq + m
1722  nwork = itaup + m
1723 *
1724 * Bidiagonalize L in A
1725 * CWorkspace: need M*M [VT] + 2*M [tauq, taup] + M [work]
1726 * CWorkspace: prefer M*M [VT] + 2*M [tauq, taup] + 2*M*NB [work]
1727 * RWorkspace: need M [e]
1728 *
1729  CALL zgebrd( m, m, a, lda, s, rwork( ie ), work( itauq ),
1730  $ work( itaup ), work( nwork ), lwork-nwork+1,
1731  $ ierr )
1732 *
1733 * Perform bidiagonal SVD, computing left singular vectors
1734 * of bidiagonal matrix in RWORK(IRU) and computing right
1735 * singular vectors of bidiagonal matrix in RWORK(IRVT)
1736 * CWorkspace: need 0
1737 * RWorkspace: need M [e] + M*M [RU] + M*M [RVT] + BDSPAC
1738 *
1739  iru = ie + m
1740  irvt = iru + m*m
1741  nrwork = irvt + m*m
1742  CALL dbdsdc( 'U', 'I', m, s, rwork( ie ), rwork( iru ),
1743  $ m, rwork( irvt ), m, dum, idum,
1744  $ rwork( nrwork ), iwork, info )
1745 *
1746 * Copy real matrix RWORK(IRU) to complex matrix U
1747 * Overwrite U by left singular vectors of L
1748 * CWorkspace: need M*M [VT] + 2*M [tauq, taup] + M [work]
1749 * CWorkspace: prefer M*M [VT] + 2*M [tauq, taup] + M*NB [work]
1750 * RWorkspace: need 0
1751 *
1752  CALL zlacp2( 'F', m, m, rwork( iru ), m, u, ldu )
1753  CALL zunmbr( 'Q', 'L', 'N', m, m, m, a, lda,
1754  $ work( itauq ), u, ldu, work( nwork ),
1755  $ lwork-nwork+1, ierr )
1756 *
1757 * Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
1758 * Overwrite WORK(IVT) by right singular vectors of L
1759 * CWorkspace: need M*M [VT] + 2*M [tauq, taup] + M [work]
1760 * CWorkspace: prefer M*M [VT] + 2*M [tauq, taup] + M*NB [work]
1761 * RWorkspace: need 0
1762 *
1763  CALL zlacp2( 'F', m, m, rwork( irvt ), m, work( ivt ),
1764  $ ldwkvt )
1765  CALL zunmbr( 'P', 'R', 'C', m, m, m, a, lda,
1766  $ work( itaup ), work( ivt ), ldwkvt,
1767  $ work( nwork ), lwork-nwork+1, ierr )
1768 *
1769 * Multiply right singular vectors of L in WORK(IVT) by
1770 * Q in VT, storing result in A
1771 * CWorkspace: need M*M [VT]
1772 * RWorkspace: need 0
1773 *
1774  CALL zgemm( 'N', 'N', m, n, m, cone, work( ivt ), ldwkvt,
1775  $ vt, ldvt, czero, a, lda )
1776 *
1777 * Copy right singular vectors of A from A to VT
1778 *
1779  CALL zlacpy( 'F', m, n, a, lda, vt, ldvt )
1780 *
1781  END IF
1782 *
1783  ELSE IF( n.GE.mnthr2 ) THEN
1784 *
1785 * MNTHR2 <= N < MNTHR1
1786 *
1787 * Path 5t (N >> M, but not as much as MNTHR1)
1788 * Reduce to bidiagonal form without QR decomposition, use
1789 * ZUNGBR and matrix multiplication to compute singular vectors
1790 *
1791  ie = 1
1792  nrwork = ie + m
1793  itauq = 1
1794  itaup = itauq + m
1795  nwork = itaup + m
1796 *
1797 * Bidiagonalize A
1798 * CWorkspace: need 2*M [tauq, taup] + N [work]
1799 * CWorkspace: prefer 2*M [tauq, taup] + (M+N)*NB [work]
1800 * RWorkspace: need M [e]
1801 *
1802  CALL zgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
1803  $ work( itaup ), work( nwork ), lwork-nwork+1,
1804  $ ierr )
1805 *
1806  IF( wntqn ) THEN
1807 *
1808 * Path 5tn (N >> M, JOBZ='N')
1809 * Compute singular values only
1810 * CWorkspace: need 0
1811 * RWorkspace: need M [e] + BDSPAC
1812 *
1813  CALL dbdsdc( 'L', 'N', m, s, rwork( ie ), dum,1,dum,1,
1814  $ dum, idum, rwork( nrwork ), iwork, info )
1815  ELSE IF( wntqo ) THEN
1816  irvt = nrwork
1817  iru = irvt + m*m
1818  nrwork = iru + m*m
1819  ivt = nwork
1820 *
1821 * Path 5to (N >> M, JOBZ='O')
1822 * Copy A to U, generate Q
1823 * CWorkspace: need 2*M [tauq, taup] + M [work]
1824 * CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
1825 * RWorkspace: need 0
1826 *
1827  CALL zlacpy( 'L', m, m, a, lda, u, ldu )
1828  CALL zungbr( 'Q', m, m, n, u, ldu, work( itauq ),
1829  $ work( nwork ), lwork-nwork+1, ierr )
1830 *
1831 * Generate P**H in A
1832 * CWorkspace: need 2*M [tauq, taup] + M [work]
1833 * CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
1834 * RWorkspace: need 0
1835 *
1836  CALL zungbr( 'P', m, n, m, a, lda, work( itaup ),
1837  $ work( nwork ), lwork-nwork+1, ierr )
1838 *
1839  ldwkvt = m
1840  IF( lwork .GE. m*n + 3*m ) THEN
1841 *
1842 * WORK( IVT ) is M by N
1843 *
1844  nwork = ivt + ldwkvt*n
1845  chunk = n
1846  ELSE
1847 *
1848 * WORK( IVT ) is M by CHUNK
1849 *
1850  chunk = ( lwork - 3*m ) / m
1851  nwork = ivt + ldwkvt*chunk
1852  END IF
1853 *
1854 * Perform bidiagonal SVD, computing left singular vectors
1855 * of bidiagonal matrix in RWORK(IRU) and computing right
1856 * singular vectors of bidiagonal matrix in RWORK(IRVT)
1857 * CWorkspace: need 0
1858 * RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC
1859 *
1860  CALL dbdsdc( 'L', 'I', m, s, rwork( ie ), rwork( iru ),
1861  $ m, rwork( irvt ), m, dum, idum,
1862  $ rwork( nrwork ), iwork, info )
1863 *
1864 * Multiply Q in U by real matrix RWORK(IRVT)
1865 * storing the result in WORK(IVT), copying to U
1866 * CWorkspace: need 2*M [tauq, taup] + M*M [VT]
1867 * RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + 2*M*M [rwork]
1868 *
1869  CALL zlacrm( m, m, u, ldu, rwork( iru ), m, work( ivt ),
1870  $ ldwkvt, rwork( nrwork ) )
1871  CALL zlacpy( 'F', m, m, work( ivt ), ldwkvt, u, ldu )
1872 *
1873 * Multiply RWORK(IRVT) by P**H in A, storing the
1874 * result in WORK(IVT), copying to A
1875 * CWorkspace: need 2*M [tauq, taup] + M*M [VT]
1876 * CWorkspace: prefer 2*M [tauq, taup] + M*N [VT]
1877 * RWorkspace: need M [e] + M*M [RVT] + 2*M*M [rwork]
1878 * RWorkspace: prefer M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here
1879 *
1880  nrwork = iru
1881  DO 50 i = 1, n, chunk
1882  blk = min( n-i+1, chunk )
1883  CALL zlarcm( m, blk, rwork( irvt ), m, a( 1, i ), lda,
1884  $ work( ivt ), ldwkvt, rwork( nrwork ) )
1885  CALL zlacpy( 'F', m, blk, work( ivt ), ldwkvt,
1886  $ a( 1, i ), lda )
1887  50 CONTINUE
1888  ELSE IF( wntqs ) THEN
1889 *
1890 * Path 5ts (N >> M, JOBZ='S')
1891 * Copy A to U, generate Q
1892 * CWorkspace: need 2*M [tauq, taup] + M [work]
1893 * CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
1894 * RWorkspace: need 0
1895 *
1896  CALL zlacpy( 'L', m, m, a, lda, u, ldu )
1897  CALL zungbr( 'Q', m, m, n, u, ldu, work( itauq ),
1898  $ work( nwork ), lwork-nwork+1, ierr )
1899 *
1900 * Copy A to VT, generate P**H
1901 * CWorkspace: need 2*M [tauq, taup] + M [work]
1902 * CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
1903 * RWorkspace: need 0
1904 *
1905  CALL zlacpy( 'U', m, n, a, lda, vt, ldvt )
1906  CALL zungbr( 'P', m, n, m, vt, ldvt, work( itaup ),
1907  $ work( nwork ), lwork-nwork+1, ierr )
1908 *
1909 * Perform bidiagonal SVD, computing left singular vectors
1910 * of bidiagonal matrix in RWORK(IRU) and computing right
1911 * singular vectors of bidiagonal matrix in RWORK(IRVT)
1912 * CWorkspace: need 0
1913 * RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC
1914 *
1915  irvt = nrwork
1916  iru = irvt + m*m
1917  nrwork = iru + m*m
1918  CALL dbdsdc( 'L', 'I', m, s, rwork( ie ), rwork( iru ),
1919  $ m, rwork( irvt ), m, dum, idum,
1920  $ rwork( nrwork ), iwork, info )
1921 *
1922 * Multiply Q in U by real matrix RWORK(IRU), storing the
1923 * result in A, copying to U
1924 * CWorkspace: need 0
1925 * RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + 2*M*M [rwork]
1926 *
1927  CALL zlacrm( m, m, u, ldu, rwork( iru ), m, a, lda,
1928  $ rwork( nrwork ) )
1929  CALL zlacpy( 'F', m, m, a, lda, u, ldu )
1930 *
1931 * Multiply real matrix RWORK(IRVT) by P**H in VT,
1932 * storing the result in A, copying to VT
1933 * CWorkspace: need 0
1934 * RWorkspace: need M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here
1935 *
1936  nrwork = iru
1937  CALL zlarcm( m, n, rwork( irvt ), m, vt, ldvt, a, lda,
1938  $ rwork( nrwork ) )
1939  CALL zlacpy( 'F', m, n, a, lda, vt, ldvt )
1940  ELSE
1941 *
1942 * Path 5ta (N >> M, JOBZ='A')
1943 * Copy A to U, generate Q
1944 * CWorkspace: need 2*M [tauq, taup] + M [work]
1945 * CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
1946 * RWorkspace: need 0
1947 *
1948  CALL zlacpy( 'L', m, m, a, lda, u, ldu )
1949  CALL zungbr( 'Q', m, m, n, u, ldu, work( itauq ),
1950  $ work( nwork ), lwork-nwork+1, ierr )
1951 *
1952 * Copy A to VT, generate P**H
1953 * CWorkspace: need 2*M [tauq, taup] + N [work]
1954 * CWorkspace: prefer 2*M [tauq, taup] + N*NB [work]
1955 * RWorkspace: need 0
1956 *
1957  CALL zlacpy( 'U', m, n, a, lda, vt, ldvt )
1958  CALL zungbr( 'P', n, n, m, vt, ldvt, work( itaup ),
1959  $ work( nwork ), lwork-nwork+1, ierr )
1960 *
1961 * Perform bidiagonal SVD, computing left singular vectors
1962 * of bidiagonal matrix in RWORK(IRU) and computing right
1963 * singular vectors of bidiagonal matrix in RWORK(IRVT)
1964 * CWorkspace: need 0
1965 * RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC
1966 *
1967  irvt = nrwork
1968  iru = irvt + m*m
1969  nrwork = iru + m*m
1970  CALL dbdsdc( 'L', 'I', m, s, rwork( ie ), rwork( iru ),
1971  $ m, rwork( irvt ), m, dum, idum,
1972  $ rwork( nrwork ), iwork, info )
1973 *
1974 * Multiply Q in U by real matrix RWORK(IRU), storing the
1975 * result in A, copying to U
1976 * CWorkspace: need 0
1977 * RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + 2*M*M [rwork]
1978 *
1979  CALL zlacrm( m, m, u, ldu, rwork( iru ), m, a, lda,
1980  $ rwork( nrwork ) )
1981  CALL zlacpy( 'F', m, m, a, lda, u, ldu )
1982 *
1983 * Multiply real matrix RWORK(IRVT) by P**H in VT,
1984 * storing the result in A, copying to VT
1985 * CWorkspace: need 0
1986 * RWorkspace: need M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here
1987 *
1988  nrwork = iru
1989  CALL zlarcm( m, n, rwork( irvt ), m, vt, ldvt, a, lda,
1990  $ rwork( nrwork ) )
1991  CALL zlacpy( 'F', m, n, a, lda, vt, ldvt )
1992  END IF
1993 *
1994  ELSE
1995 *
1996 * N .LT. MNTHR2
1997 *
1998 * Path 6t (N > M, but not much larger)
1999 * Reduce to bidiagonal form without LQ decomposition
2000 * Use ZUNMBR to compute singular vectors
2001 *
2002  ie = 1
2003  nrwork = ie + m
2004  itauq = 1
2005  itaup = itauq + m
2006  nwork = itaup + m
2007 *
2008 * Bidiagonalize A
2009 * CWorkspace: need 2*M [tauq, taup] + N [work]
2010 * CWorkspace: prefer 2*M [tauq, taup] + (M+N)*NB [work]
2011 * RWorkspace: need M [e]
2012 *
2013  CALL zgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
2014  $ work( itaup ), work( nwork ), lwork-nwork+1,
2015  $ ierr )
2016  IF( wntqn ) THEN
2017 *
2018 * Path 6tn (N > M, JOBZ='N')
2019 * Compute singular values only
2020 * CWorkspace: need 0
2021 * RWorkspace: need M [e] + BDSPAC
2022 *
2023  CALL dbdsdc( 'L', 'N', m, s, rwork( ie ), dum,1,dum,1,
2024  $ dum, idum, rwork( nrwork ), iwork, info )
2025  ELSE IF( wntqo ) THEN
2026 * Path 6to (N > M, JOBZ='O')
2027  ldwkvt = m
2028  ivt = nwork
2029  IF( lwork .GE. m*n + 3*m ) THEN
2030 *
2031 * WORK( IVT ) is M by N
2032 *
2033  CALL zlaset( 'F', m, n, czero, czero, work( ivt ),
2034  $ ldwkvt )
2035  nwork = ivt + ldwkvt*n
2036  ELSE
2037 *
2038 * WORK( IVT ) is M by CHUNK
2039 *
2040  chunk = ( lwork - 3*m ) / m
2041  nwork = ivt + ldwkvt*chunk
2042  END IF
2043 *
2044 * Perform bidiagonal SVD, computing left singular vectors
2045 * of bidiagonal matrix in RWORK(IRU) and computing right
2046 * singular vectors of bidiagonal matrix in RWORK(IRVT)
2047 * CWorkspace: need 0
2048 * RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC
2049 *
2050  irvt = nrwork
2051  iru = irvt + m*m
2052  nrwork = iru + m*m
2053  CALL dbdsdc( 'L', 'I', m, s, rwork( ie ), rwork( iru ),
2054  $ m, rwork( irvt ), m, dum, idum,
2055  $ rwork( nrwork ), iwork, info )
2056 *
2057 * Copy real matrix RWORK(IRU) to complex matrix U
2058 * Overwrite U by left singular vectors of A
2059 * CWorkspace: need 2*M [tauq, taup] + M*M [VT] + M [work]
2060 * CWorkspace: prefer 2*M [tauq, taup] + M*M [VT] + M*NB [work]
2061 * RWorkspace: need M [e] + M*M [RVT] + M*M [RU]
2062 *
2063  CALL zlacp2( 'F', m, m, rwork( iru ), m, u, ldu )
2064  CALL zunmbr( 'Q', 'L', 'N', m, m, n, a, lda,
2065  $ work( itauq ), u, ldu, work( nwork ),
2066  $ lwork-nwork+1, ierr )
2067 *
2068  IF( lwork .GE. m*n + 3*m ) THEN
2069 *
2070 * Path 6to-fast
2071 * Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
2072 * Overwrite WORK(IVT) by right singular vectors of A,
2073 * copying to A
2074 * CWorkspace: need 2*M [tauq, taup] + M*N [VT] + M [work]
2075 * CWorkspace: prefer 2*M [tauq, taup] + M*N [VT] + M*NB [work]
2076 * RWorkspace: need M [e] + M*M [RVT]
2077 *
2078  CALL zlacp2( 'F', m, m, rwork( irvt ), m, work( ivt ),
2079  $ ldwkvt )
2080  CALL zunmbr( 'P', 'R', 'C', m, n, m, a, lda,
2081  $ work( itaup ), work( ivt ), ldwkvt,
2082  $ work( nwork ), lwork-nwork+1, ierr )
2083  CALL zlacpy( 'F', m, n, work( ivt ), ldwkvt, a, lda )
2084  ELSE
2085 *
2086 * Path 6to-slow
2087 * Generate P**H in A
2088 * CWorkspace: need 2*M [tauq, taup] + M*M [VT] + M [work]
2089 * CWorkspace: prefer 2*M [tauq, taup] + M*M [VT] + M*NB [work]
2090 * RWorkspace: need 0
2091 *
2092  CALL zungbr( 'P', m, n, m, a, lda, work( itaup ),
2093  $ work( nwork ), lwork-nwork+1, ierr )
2094 *
2095 * Multiply Q in A by real matrix RWORK(IRU), storing the
2096 * result in WORK(IU), copying to A
2097 * CWorkspace: need 2*M [tauq, taup] + M*M [VT]
2098 * CWorkspace: prefer 2*M [tauq, taup] + M*N [VT]
2099 * RWorkspace: need M [e] + M*M [RVT] + 2*M*M [rwork]
2100 * RWorkspace: prefer M [e] + M*M [RVT] + 2*M*N [rwork] < M + 5*M*M since N < 2*M here
2101 *
2102  nrwork = iru
2103  DO 60 i = 1, n, chunk
2104  blk = min( n-i+1, chunk )
2105  CALL zlarcm( m, blk, rwork( irvt ), m, a( 1, i ),
2106  $ lda, work( ivt ), ldwkvt,
2107  $ rwork( nrwork ) )
2108  CALL zlacpy( 'F', m, blk, work( ivt ), ldwkvt,
2109  $ a( 1, i ), lda )
2110  60 CONTINUE
2111  END IF
2112  ELSE IF( wntqs ) THEN
2113 *
2114 * Path 6ts (N > M, JOBZ='S')
2115 * Perform bidiagonal SVD, computing left singular vectors
2116 * of bidiagonal matrix in RWORK(IRU) and computing right
2117 * singular vectors of bidiagonal matrix in RWORK(IRVT)
2118 * CWorkspace: need 0
2119 * RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC
2120 *
2121  irvt = nrwork
2122  iru = irvt + m*m
2123  nrwork = iru + m*m
2124  CALL dbdsdc( 'L', 'I', m, s, rwork( ie ), rwork( iru ),
2125  $ m, rwork( irvt ), m, dum, idum,
2126  $ rwork( nrwork ), iwork, info )
2127 *
2128 * Copy real matrix RWORK(IRU) to complex matrix U
2129 * Overwrite U by left singular vectors of A
2130 * CWorkspace: need 2*M [tauq, taup] + M [work]
2131 * CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
2132 * RWorkspace: need M [e] + M*M [RVT] + M*M [RU]
2133 *
2134  CALL zlacp2( 'F', m, m, rwork( iru ), m, u, ldu )
2135  CALL zunmbr( 'Q', 'L', 'N', m, m, n, a, lda,
2136  $ work( itauq ), u, ldu, work( nwork ),
2137  $ lwork-nwork+1, ierr )
2138 *
2139 * Copy real matrix RWORK(IRVT) to complex matrix VT
2140 * Overwrite VT by right singular vectors of A
2141 * CWorkspace: need 2*M [tauq, taup] + M [work]
2142 * CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
2143 * RWorkspace: need M [e] + M*M [RVT]
2144 *
2145  CALL zlaset( 'F', m, n, czero, czero, vt, ldvt )
2146  CALL zlacp2( 'F', m, m, rwork( irvt ), m, vt, ldvt )
2147  CALL zunmbr( 'P', 'R', 'C', m, n, m, a, lda,
2148  $ work( itaup ), vt, ldvt, work( nwork ),
2149  $ lwork-nwork+1, ierr )
2150  ELSE
2151 *
2152 * Path 6ta (N > M, JOBZ='A')
2153 * Perform bidiagonal SVD, computing left singular vectors
2154 * of bidiagonal matrix in RWORK(IRU) and computing right
2155 * singular vectors of bidiagonal matrix in RWORK(IRVT)
2156 * CWorkspace: need 0
2157 * RWorkspace: need M [e] + M*M [RVT] + M*M [RU] + BDSPAC
2158 *
2159  irvt = nrwork
2160  iru = irvt + m*m
2161  nrwork = iru + m*m
2162 *
2163  CALL dbdsdc( 'L', 'I', m, s, rwork( ie ), rwork( iru ),
2164  $ m, rwork( irvt ), m, dum, idum,
2165  $ rwork( nrwork ), iwork, info )
2166 *
2167 * Copy real matrix RWORK(IRU) to complex matrix U
2168 * Overwrite U by left singular vectors of A
2169 * CWorkspace: need 2*M [tauq, taup] + M [work]
2170 * CWorkspace: prefer 2*M [tauq, taup] + M*NB [work]
2171 * RWorkspace: need M [e] + M*M [RVT] + M*M [RU]
2172 *
2173  CALL zlacp2( 'F', m, m, rwork( iru ), m, u, ldu )
2174  CALL zunmbr( 'Q', 'L', 'N', m, m, n, a, lda,
2175  $ work( itauq ), u, ldu, work( nwork ),
2176  $ lwork-nwork+1, ierr )
2177 *
2178 * Set all of VT to identity matrix
2179 *
2180  CALL zlaset( 'F', n, n, czero, cone, vt, ldvt )
2181 *
2182 * Copy real matrix RWORK(IRVT) to complex matrix VT
2183 * Overwrite VT by right singular vectors of A
2184 * CWorkspace: need 2*M [tauq, taup] + N [work]
2185 * CWorkspace: prefer 2*M [tauq, taup] + N*NB [work]
2186 * RWorkspace: need M [e] + M*M [RVT]
2187 *
2188  CALL zlacp2( 'F', m, m, rwork( irvt ), m, vt, ldvt )
2189  CALL zunmbr( 'P', 'R', 'C', n, n, m, a, lda,
2190  $ work( itaup ), vt, ldvt, work( nwork ),
2191  $ lwork-nwork+1, ierr )
2192  END IF
2193 *
2194  END IF
2195 *
2196  END IF
2197 *
2198 * Undo scaling if necessary
2199 *
2200  IF( iscl.EQ.1 ) THEN
2201  IF( anrm.GT.bignum )
2202  $ CALL dlascl( 'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
2203  $ ierr )
2204  IF( info.NE.0 .AND. anrm.GT.bignum )
2205  $ CALL dlascl( 'G', 0, 0, bignum, anrm, minmn-1, 1,
2206  $ rwork( ie ), minmn, ierr )
2207  IF( anrm.LT.smlnum )
2208  $ CALL dlascl( 'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
2209  $ ierr )
2210  IF( info.NE.0 .AND. anrm.LT.smlnum )
2211  $ CALL dlascl( 'G', 0, 0, smlnum, anrm, minmn-1, 1,
2212  $ rwork( ie ), minmn, ierr )
2213  END IF
2214 *
2215 * Return optimal workspace in WORK(1)
2216 *
2217  work( 1 ) = droundup_lwork( maxwrk )
2218 *
2219  RETURN
2220 *
2221 * End of ZGESDD
2222 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
double precision function droundup_lwork(LWORK)
DROUNDUP_LWORK
logical function disnan(DIN)
DISNAN tests input for NaN.
Definition: disnan.f:59
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: dlascl.f:143
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine dbdsdc(UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ, WORK, IWORK, INFO)
DBDSDC
Definition: dbdsdc.f:205
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:187
subroutine zungbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGBR
Definition: zungbr.f:157
double precision function zlange(NORM, M, N, A, LDA, WORK)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: zlange.f:115
subroutine zgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGELQF
Definition: zgelqf.f:143
subroutine zgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
ZGEBRD
Definition: zgebrd.f:205
subroutine zlarcm(M, N, A, LDA, B, LDB, C, LDC, RWORK)
ZLARCM copies all or part of a real two-dimensional array to a complex array.
Definition: zlarcm.f:114
subroutine zlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: zlascl.f:143
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:103
subroutine zlacp2(UPLO, M, N, A, LDA, B, LDB)
ZLACP2 copies all or part of a real two-dimensional array to a complex array.
Definition: zlacp2.f:104
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: zlaset.f:106
subroutine zlacrm(M, N, A, LDA, B, LDB, C, LDC, RWORK)
ZLACRM multiplies a complex matrix by a square real matrix.
Definition: zlacrm.f:114
subroutine zunglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGLQ
Definition: zunglq.f:127
subroutine zungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGQR
Definition: zungqr.f:128
subroutine zunmbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMBR
Definition: zunmbr.f:196
subroutine zgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGEQRF VARIANT: left-looking Level 3 BLAS of the algorithm.
Definition: zgeqrf.f:151
Here is the call graph for this function:
Here is the caller graph for this function: