LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
zgesvd.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine zgesvd (JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, RWORK, INFO)
  ZGESVD computes the singular value decomposition (SVD) for GE matrices More...
 

Function/Subroutine Documentation

subroutine zgesvd ( character  JOBU,
character  JOBVT,
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  INFO 
)

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

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

Purpose:
 ZGESVD 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*16 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 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)
          (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*16 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*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 >=  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 DOUBLE PRECISION 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 ZBDSQR 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 zgesvd.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: