LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine dgesvd ( character  JOBU,
character  JOBVT,
integer  M,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  S,
double precision, dimension( ldu, * )  U,
integer  LDU,
double precision, dimension( ldvt, * )  VT,
integer  LDVT,
double precision, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

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

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

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

      A = U * SIGMA * 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 orthogonal matrix, and
 V is an N-by-N orthogonal 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**T, 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**T:
          = 'A':  all N rows of V**T are returned in the array VT;
          = 'S':  the first min(m,n) rows of V**T (the right singular
                  vectors) are returned in the array VT;
          = 'O':  the first min(m,n) rows of V**T (the right singular
                  vectors) are overwritten on the array A;
          = 'N':  no rows of V**T (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 DOUBLE PRECISION 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**T (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 DOUBLE PRECISION 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 orthogonal 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 DOUBLE PRECISION array, dimension (LDVT,N)
          If JOBVT = 'A', VT contains the N-by-N orthogonal matrix
          V**T;
          if JOBVT = 'S', VT contains the first min(m,n) rows of
          V**T (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 DOUBLE PRECISION array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
          if INFO > 0, WORK(2:MIN(M,N)) 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.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.
          LWORK >= MAX(1,5*MIN(M,N)) for the paths (see comments inside code):
             - PATH 1  (M much larger than N, JOBU='N') 
             - PATH 1t (N much larger than M, JOBVT='N')
          LWORK >= MAX(1,3*MIN(M,N) + MAX(M,N),5*MIN(M,N)) for the other paths
          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]INFO
          INFO is INTEGER
          = 0:  successful exit.
          < 0:  if INFO = -i, the i-th argument had an illegal value.
          > 0:  if DBDSQR did not converge, INFO specifies how many
                superdiagonals of an intermediate bidiagonal form B
                did not converge to zero. See the description of WORK
                above for details.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
April 2012

Definition at line 213 of file dgesvd.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: