Home
Browse Files
   1       SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK,
   2      $                   LWORK, IWORK, INFO )
   3 *
   4 *  -- LAPACK driver routine (version 3.2.1)                                  --
   5 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
   6 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
   7 *     March 2009
   8 *
   9 *     .. Scalar Arguments ..
  10       CHARACTER          JOBZ
  11       INTEGER            INFO, LDA, LDU, LDVT, LWORK, M, N
  12 *     ..
  13 *     .. Array Arguments ..
  14       INTEGER            IWORK( * )
  15       DOUBLE PRECISION   A( LDA, * ), S( * ), U( LDU, * ),
  16      $                   VT( LDVT, * ), WORK( * )
  17 *     ..
  18 *
  19 *  Purpose
  20 *  =======
  21 *
  22 *  DGESDD computes the singular value decomposition (SVD) of a real
  23 *  M-by-N matrix A, optionally computing the left and right singular
  24 *  vectors.  If singular vectors are desired, it uses a
  25 *  divide-and-conquer algorithm.
  26 *
  27 *  The SVD is written
  28 *
  29 *       A = U * SIGMA * transpose(V)
  30 *
  31 *  where SIGMA is an M-by-N matrix which is zero except for its
  32 *  min(m,n) diagonal elements, U is an M-by-M orthogonal matrix, and
  33 *  V is an N-by-N orthogonal matrix.  The diagonal elements of SIGMA
  34 *  are the singular values of A; they are real and non-negative, and
  35 *  are returned in descending order.  The first min(m,n) columns of
  36 *  U and V are the left and right singular vectors of A.
  37 *
  38 *  Note that the routine returns VT = V**T, not V.
  39 *
  40 *  The divide and conquer algorithm makes very mild assumptions about
  41 *  floating point arithmetic. It will work on machines with a guard
  42 *  digit in add/subtract, or on those binary machines without guard
  43 *  digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
  44 *  Cray-2. It could conceivably fail on hexadecimal or decimal machines
  45 *  without guard digits, but we know of none.
  46 *
  47 *  Arguments
  48 *  =========
  49 *
  50 *  JOBZ    (input) CHARACTER*1
  51 *          Specifies options for computing all or part of the matrix U:
  52 *          = 'A':  all M columns of U and all N rows of V**T are
  53 *                  returned in the arrays U and VT;
  54 *          = 'S':  the first min(M,N) columns of U and the first
  55 *                  min(M,N) rows of V**T are returned in the arrays U
  56 *                  and VT;
  57 *          = 'O':  If M >= N, the first N columns of U are overwritten
  58 *                  on the array A and all rows of V**T are returned in
  59 *                  the array VT;
  60 *                  otherwise, all columns of U are returned in the
  61 *                  array U and the first M rows of V**T are overwritten
  62 *                  in the array A;
  63 *          = 'N':  no columns of U or rows of V**T are computed.
  64 *
  65 *  M       (input) INTEGER
  66 *          The number of rows of the input matrix A.  M >= 0.
  67 *
  68 *  N       (input) INTEGER
  69 *          The number of columns of the input matrix A.  N >= 0.
  70 *
  71 *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
  72 *          On entry, the M-by-N matrix A.
  73 *          On exit,
  74 *          if JOBZ = 'O',  A is overwritten with the first N columns
  75 *                          of U (the left singular vectors, stored
  76 *                          columnwise) if M >= N;
  77 *                          A is overwritten with the first M rows
  78 *                          of V**T (the right singular vectors, stored
  79 *                          rowwise) otherwise.
  80 *          if JOBZ .ne. 'O', the contents of A are destroyed.
  81 *
  82 *  LDA     (input) INTEGER
  83 *          The leading dimension of the array A.  LDA >= max(1,M).
  84 *
  85 *  S       (output) DOUBLE PRECISION array, dimension (min(M,N))
  86 *          The singular values of A, sorted so that S(i) >= S(i+1).
  87 *
  88 *  U       (output) DOUBLE PRECISION array, dimension (LDU,UCOL)
  89 *          UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N;
  90 *          UCOL = min(M,N) if JOBZ = 'S'.
  91 *          If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M
  92 *          orthogonal matrix U;
  93 *          if JOBZ = 'S', U contains the first min(M,N) columns of U
  94 *          (the left singular vectors, stored columnwise);
  95 *          if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced.
  96 *
  97 *  LDU     (input) INTEGER
  98 *          The leading dimension of the array U.  LDU >= 1; if
  99 *          JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M.
 100 *
 101 *  VT      (output) DOUBLE PRECISION array, dimension (LDVT,N)
 102 *          If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the
 103 *          N-by-N orthogonal matrix V**T;
 104 *          if JOBZ = 'S', VT contains the first min(M,N) rows of
 105 *          V**T (the right singular vectors, stored rowwise);
 106 *          if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced.
 107 *
 108 *  LDVT    (input) INTEGER
 109 *          The leading dimension of the array VT.  LDVT >= 1; if
 110 *          JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N;
 111 *          if JOBZ = 'S', LDVT >= min(M,N).
 112 *
 113 *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 114 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
 115 *
 116 *  LWORK   (input) INTEGER
 117 *          The dimension of the array WORK. LWORK >= 1.
 118 *          If JOBZ = 'N',
 119 *            LWORK >= 3*min(M,N) + max(max(M,N),7*min(M,N)).
 120 *          If JOBZ = 'O',
 121 *            LWORK >= 3*min(M,N) + 
 122 *                     max(max(M,N),5*min(M,N)*min(M,N)+4*min(M,N)).
 123 *          If JOBZ = 'S' or 'A'
 124 *            LWORK >= 3*min(M,N) +
 125 *                     max(max(M,N),4*min(M,N)*min(M,N)+4*min(M,N)).
 126 *          For good performance, LWORK should generally be larger.
 127 *          If LWORK = -1 but other input arguments are legal, WORK(1)
 128 *          returns the optimal LWORK.
 129 *
 130 *  IWORK   (workspace) INTEGER array, dimension (8*min(M,N))
 131 *
 132 *  INFO    (output) INTEGER
 133 *          = 0:  successful exit.
 134 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
 135 *          > 0:  DBDSDC did not converge, updating process failed.
 136 *
 137 *  Further Details
 138 *  ===============
 139 *
 140 *  Based on contributions by
 141 *     Ming Gu and Huan Ren, Computer Science Division, University of
 142 *     California at Berkeley, USA
 143 *
 144 *  =====================================================================
 145 *
 146 *     .. Parameters ..
 147       DOUBLE PRECISION   ZERO, ONE
 148       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
 149 *     ..
 150 *     .. Local Scalars ..
 151       LOGICAL            LQUERY, WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS
 152       INTEGER            BDSPAC, BLK, CHUNK, I, IE, IERR, IL,
 153      $                   IR, ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT,
 154      $                   LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK,
 155      $                   MNTHR, NWORK, WRKBL
 156       DOUBLE PRECISION   ANRM, BIGNUM, EPS, SMLNUM
 157 *     ..
 158 *     .. Local Arrays ..
 159       INTEGER            IDUM( 1 )
 160       DOUBLE PRECISION   DUM( 1 )
 161 *     ..
 162 *     .. External Subroutines ..
 163       EXTERNAL           DBDSDC, DGEBRD, DGELQF, DGEMM, DGEQRF, DLACPY,
 164      $                   DLASCL, DLASET, DORGBR, DORGLQ, DORGQR, DORMBR,
 165      $                   XERBLA
 166 *     ..
 167 *     .. External Functions ..
 168       LOGICAL            LSAME
 169       INTEGER            ILAENV
 170       DOUBLE PRECISION   DLAMCH, DLANGE
 171       EXTERNAL           DLAMCH, DLANGE, ILAENV, LSAME
 172 *     ..
 173 *     .. Intrinsic Functions ..
 174       INTRINSIC          INT, MAX, MIN, SQRT
 175 *     ..
 176 *     .. Executable Statements ..
 177 *
 178 *     Test the input arguments
 179 *
 180       INFO = 0
 181       MINMN = MIN( M, N )
 182       WNTQA = LSAME( JOBZ, 'A' )
 183       WNTQS = LSAME( JOBZ, 'S' )
 184       WNTQAS = WNTQA .OR. WNTQS
 185       WNTQO = LSAME( JOBZ, 'O' )
 186       WNTQN = LSAME( JOBZ, 'N' )
 187       LQUERY = ( LWORK.EQ.-1 )
 188 *
 189       IF( .NOT.( WNTQA .OR. WNTQS .OR. WNTQO .OR. WNTQN ) ) THEN
 190          INFO = -1
 191       ELSE IF( M.LT.0 ) THEN
 192          INFO = -2
 193       ELSE IF( N.LT.0 ) THEN
 194          INFO = -3
 195       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
 196          INFO = -5
 197       ELSE IF( LDU.LT.1 .OR. ( WNTQAS .AND. LDU.LT.M ) .OR.
 198      $         ( WNTQO .AND. M.LT.N .AND. LDU.LT.M ) ) THEN
 199          INFO = -8
 200       ELSE IF( LDVT.LT.1 .OR. ( WNTQA .AND. LDVT.LT.N ) .OR.
 201      $         ( WNTQS .AND. LDVT.LT.MINMN ) .OR.
 202      $         ( WNTQO .AND. M.GE.N .AND. LDVT.LT.N ) ) THEN
 203          INFO = -10
 204       END IF
 205 *
 206 *     Compute workspace
 207 *      (Note: Comments in the code beginning "Workspace:" describe the
 208 *       minimal amount of workspace needed at that point in the code,
 209 *       as well as the preferred amount for good performance.
 210 *       NB refers to the optimal block size for the immediately
 211 *       following subroutine, as returned by ILAENV.)
 212 *
 213       IF( INFO.EQ.0 ) THEN
 214          MINWRK = 1
 215          MAXWRK = 1
 216          IF( M.GE.N .AND. MINMN.GT.0 ) THEN
 217 *
 218 *           Compute space needed for DBDSDC
 219 *
 220             MNTHR = INT( MINMN*11.0D0 / 6.0D0 )
 221             IF( WNTQN ) THEN
 222                BDSPAC = 7*N
 223             ELSE
 224                BDSPAC = 3*N*N + 4*N
 225             END IF
 226             IF( M.GE.MNTHR ) THEN
 227                IF( WNTQN ) THEN
 228 *
 229 *                 Path 1 (M much larger than N, JOBZ='N')
 230 *
 231                   WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1,
 232      $                    -1 )
 233                   WRKBL = MAX( WRKBL, 3*N+2*N*
 234      $                    ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
 235                   MAXWRK = MAX( WRKBL, BDSPAC+N )
 236                   MINWRK = BDSPAC + N
 237                ELSE IF( WNTQO ) THEN
 238 *
 239 *                 Path 2 (M much larger than N, JOBZ='O')
 240 *
 241                   WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
 242                   WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M,
 243      $                    N, N, -1 ) )
 244                   WRKBL = MAX( WRKBL, 3*N+2*N*
 245      $                    ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
 246                   WRKBL = MAX( WRKBL, 3*N+N*
 247      $                    ILAENV( 1, 'DORMBR', 'QLN', N, N, N, -1 ) )
 248                   WRKBL = MAX( WRKBL, 3*N+N*
 249      $                    ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
 250                   WRKBL = MAX( WRKBL, BDSPAC+3*N )
 251                   MAXWRK = WRKBL + 2*N*N
 252                   MINWRK = BDSPAC + 2*N*N + 3*N
 253                ELSE IF( WNTQS ) THEN
 254 *
 255 *                 Path 3 (M much larger than N, JOBZ='S')
 256 *
 257                   WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
 258                   WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'DORGQR', ' ', M,
 259      $                    N, N, -1 ) )
 260                   WRKBL = MAX( WRKBL, 3*N+2*N*
 261      $                    ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
 262                   WRKBL = MAX( WRKBL, 3*N+N*
 263      $                    ILAENV( 1, 'DORMBR', 'QLN', N, N, N, -1 ) )
 264                   WRKBL = MAX( WRKBL, 3*N+N*
 265      $                    ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
 266                   WRKBL = MAX( WRKBL, BDSPAC+3*N )
 267                   MAXWRK = WRKBL + N*N
 268                   MINWRK = BDSPAC + N*N + 3*N
 269                ELSE IF( WNTQA ) THEN
 270 *
 271 *                 Path 4 (M much larger than N, JOBZ='A')
 272 *
 273                   WRKBL = N + N*ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
 274                   WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'DORGQR', ' ', M,
 275      $                    M, N, -1 ) )
 276                   WRKBL = MAX( WRKBL, 3*N+2*N*
 277      $                    ILAENV( 1, 'DGEBRD', ' ', N, N, -1, -1 ) )
 278                   WRKBL = MAX( WRKBL, 3*N+N*
 279      $                    ILAENV( 1, 'DORMBR', 'QLN', N, N, N, -1 ) )
 280                   WRKBL = MAX( WRKBL, 3*N+N*
 281      $                    ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
 282                   WRKBL = MAX( WRKBL, BDSPAC+3*N )
 283                   MAXWRK = WRKBL + N*N
 284                   MINWRK = BDSPAC + N*N + 3*N
 285                END IF
 286             ELSE
 287 *
 288 *              Path 5 (M at least N, but not much larger)
 289 *
 290                WRKBL = 3*N + ( M+N )*ILAENV( 1, 'DGEBRD', ' ', M, N, -1,
 291      $                 -1 )
 292                IF( WNTQN ) THEN
 293                   MAXWRK = MAX( WRKBL, BDSPAC+3*N )
 294                   MINWRK = 3*N + MAX( M, BDSPAC )
 295                ELSE IF( WNTQO ) THEN
 296                   WRKBL = MAX( WRKBL, 3*N+N*
 297      $                    ILAENV( 1, 'DORMBR', 'QLN', M, N, N, -1 ) )
 298                   WRKBL = MAX( WRKBL, 3*N+N*
 299      $                    ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
 300                   WRKBL = MAX( WRKBL, BDSPAC+3*N )
 301                   MAXWRK = WRKBL + M*N
 302                   MINWRK = 3*N + MAX( M, N*N+BDSPAC )
 303                ELSE IF( WNTQS ) THEN
 304                   WRKBL = MAX( WRKBL, 3*N+N*
 305      $                    ILAENV( 1, 'DORMBR', 'QLN', M, N, N, -1 ) )
 306                   WRKBL = MAX( WRKBL, 3*N+N*
 307      $                    ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
 308                   MAXWRK = MAX( WRKBL, BDSPAC+3*N )
 309                   MINWRK = 3*N + MAX( M, BDSPAC )
 310                ELSE IF( WNTQA ) THEN
 311                   WRKBL = MAX( WRKBL, 3*N+M*
 312      $                    ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) )
 313                   WRKBL = MAX( WRKBL, 3*N+N*
 314      $                    ILAENV( 1, 'DORMBR', 'PRT', N, N, N, -1 ) )
 315                   MAXWRK = MAX( MAXWRK, BDSPAC+3*N )
 316                   MINWRK = 3*N + MAX( M, BDSPAC )
 317                END IF
 318             END IF
 319          ELSE IF( MINMN.GT.0 ) THEN
 320 *
 321 *           Compute space needed for DBDSDC
 322 *
 323             MNTHR = INT( MINMN*11.0D0 / 6.0D0 )
 324             IF( WNTQN ) THEN
 325                BDSPAC = 7*M
 326             ELSE
 327                BDSPAC = 3*M*M + 4*M
 328             END IF
 329             IF( N.GE.MNTHR ) THEN
 330                IF( WNTQN ) THEN
 331 *
 332 *                 Path 1t (N much larger than M, JOBZ='N')
 333 *
 334                   WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1,
 335      $                    -1 )
 336                   WRKBL = MAX( WRKBL, 3*M+2*M*
 337      $                    ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
 338                   MAXWRK = MAX( WRKBL, BDSPAC+M )
 339                   MINWRK = BDSPAC + M
 340                ELSE IF( WNTQO ) THEN
 341 *
 342 *                 Path 2t (N much larger than M, JOBZ='O')
 343 *
 344                   WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
 345                   WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M,
 346      $                    N, M, -1 ) )
 347                   WRKBL = MAX( WRKBL, 3*M+2*M*
 348      $                    ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
 349                   WRKBL = MAX( WRKBL, 3*M+M*
 350      $                    ILAENV( 1, 'DORMBR', 'QLN', M, M, M, -1 ) )
 351                   WRKBL = MAX( WRKBL, 3*M+M*
 352      $                    ILAENV( 1, 'DORMBR', 'PRT', M, M, M, -1 ) )
 353                   WRKBL = MAX( WRKBL, BDSPAC+3*M )
 354                   MAXWRK = WRKBL + 2*M*M
 355                   MINWRK = BDSPAC + 2*M*M + 3*M
 356                ELSE IF( WNTQS ) THEN
 357 *
 358 *                 Path 3t (N much larger than M, JOBZ='S')
 359 *
 360                   WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
 361                   WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'DORGLQ', ' ', M,
 362      $                    N, M, -1 ) )
 363                   WRKBL = MAX( WRKBL, 3*M+2*M*
 364      $                    ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
 365                   WRKBL = MAX( WRKBL, 3*M+M*
 366      $                    ILAENV( 1, 'DORMBR', 'QLN', M, M, M, -1 ) )
 367                   WRKBL = MAX( WRKBL, 3*M+M*
 368      $                    ILAENV( 1, 'DORMBR', 'PRT', M, M, M, -1 ) )
 369                   WRKBL = MAX( WRKBL, BDSPAC+3*M )
 370                   MAXWRK = WRKBL + M*M
 371                   MINWRK = BDSPAC + M*M + 3*M
 372                ELSE IF( WNTQA ) THEN
 373 *
 374 *                 Path 4t (N much larger than M, JOBZ='A')
 375 *
 376                   WRKBL = M + M*ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
 377                   WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'DORGLQ', ' ', N,
 378      $                    N, M, -1 ) )
 379                   WRKBL = MAX( WRKBL, 3*M+2*M*
 380      $                    ILAENV( 1, 'DGEBRD', ' ', M, M, -1, -1 ) )
 381                   WRKBL = MAX( WRKBL, 3*M+M*
 382      $                    ILAENV( 1, 'DORMBR', 'QLN', M, M, M, -1 ) )
 383                   WRKBL = MAX( WRKBL, 3*M+M*
 384      $                    ILAENV( 1, 'DORMBR', 'PRT', M, M, M, -1 ) )
 385                   WRKBL = MAX( WRKBL, BDSPAC+3*M )
 386                   MAXWRK = WRKBL + M*M
 387                   MINWRK = BDSPAC + M*M + 3*M
 388                END IF
 389             ELSE
 390 *
 391 *              Path 5t (N greater than M, but not much larger)
 392 *
 393                WRKBL = 3*M + ( M+N )*ILAENV( 1, 'DGEBRD', ' ', M, N, -1,
 394      $                 -1 )
 395                IF( WNTQN ) THEN
 396                   MAXWRK = MAX( WRKBL, BDSPAC+3*M )
 397                   MINWRK = 3*M + MAX( N, BDSPAC )
 398                ELSE IF( WNTQO ) THEN
 399                   WRKBL = MAX( WRKBL, 3*M+M*
 400      $                    ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) )
 401                   WRKBL = MAX( WRKBL, 3*M+M*
 402      $                    ILAENV( 1, 'DORMBR', 'PRT', M, N, M, -1 ) )
 403                   WRKBL = MAX( WRKBL, BDSPAC+3*M )
 404                   MAXWRK = WRKBL + M*N
 405                   MINWRK = 3*M + MAX( N, M*M+BDSPAC )
 406                ELSE IF( WNTQS ) THEN
 407                   WRKBL = MAX( WRKBL, 3*M+M*
 408      $                    ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) )
 409                   WRKBL = MAX( WRKBL, 3*M+M*
 410      $                    ILAENV( 1, 'DORMBR', 'PRT', M, N, M, -1 ) )
 411                   MAXWRK = MAX( WRKBL, BDSPAC+3*M )
 412                   MINWRK = 3*M + MAX( N, BDSPAC )
 413                ELSE IF( WNTQA ) THEN
 414                   WRKBL = MAX( WRKBL, 3*M+M*
 415      $                    ILAENV( 1, 'DORMBR', 'QLN', M, M, N, -1 ) )
 416                   WRKBL = MAX( WRKBL, 3*M+M*
 417      $                    ILAENV( 1, 'DORMBR', 'PRT', N, N, M, -1 ) )
 418                   MAXWRK = MAX( WRKBL, BDSPAC+3*M )
 419                   MINWRK = 3*M + MAX( N, BDSPAC )
 420                END IF
 421             END IF
 422          END IF
 423          MAXWRK = MAX( MAXWRK, MINWRK )
 424          WORK( 1 ) = MAXWRK
 425 *
 426          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
 427             INFO = -12
 428          END IF
 429       END IF
 430 *
 431       IF( INFO.NE.0 ) THEN
 432          CALL XERBLA( 'DGESDD', -INFO )
 433          RETURN
 434       ELSE IF( LQUERY ) THEN
 435          RETURN
 436       END IF
 437 *
 438 *     Quick return if possible
 439 *
 440       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
 441          RETURN
 442       END IF
 443 *
 444 *     Get machine constants
 445 *
 446       EPS = DLAMCH( 'P' )
 447       SMLNUM = SQRT( DLAMCH( 'S' ) ) / EPS
 448       BIGNUM = ONE / SMLNUM
 449 *
 450 *     Scale A if max element outside range [SMLNUM,BIGNUM]
 451 *
 452       ANRM = DLANGE( 'M', M, N, A, LDA, DUM )
 453       ISCL = 0
 454       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
 455          ISCL = 1
 456          CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, IERR )
 457       ELSE IF( ANRM.GT.BIGNUM ) THEN
 458          ISCL = 1
 459          CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, IERR )
 460       END IF
 461 *
 462       IF( M.GE.N ) THEN
 463 *
 464 *        A has at least as many rows as columns. If A has sufficiently
 465 *        more rows than columns, first reduce using the QR
 466 *        decomposition (if sufficient workspace available)
 467 *
 468          IF( M.GE.MNTHR ) THEN
 469 *
 470             IF( WNTQN ) THEN
 471 *
 472 *              Path 1 (M much larger than N, JOBZ='N')
 473 *              No singular vectors to be computed
 474 *
 475                ITAU = 1
 476                NWORK = ITAU + N
 477 *
 478 *              Compute A=Q*R
 479 *              (Workspace: need 2*N, prefer N+N*NB)
 480 *
 481                CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
 482      $                      LWORK-NWORK+1, IERR )
 483 *
 484 *              Zero out below R
 485 *
 486                CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
 487                IE = 1
 488                ITAUQ = IE + N
 489                ITAUP = ITAUQ + N
 490                NWORK = ITAUP + N
 491 *
 492 *              Bidiagonalize R in A
 493 *              (Workspace: need 4*N, prefer 3*N+2*N*NB)
 494 *
 495                CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
 496      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
 497      $                      IERR )
 498                NWORK = IE + N
 499 *
 500 *              Perform bidiagonal SVD, computing singular values only
 501 *              (Workspace: need N+BDSPAC)
 502 *
 503                CALL DBDSDC( 'U', 'N', N, S, WORK( IE ), DUM, 1, DUM, 1,
 504      $                      DUM, IDUM, WORK( NWORK ), IWORK, INFO )
 505 *
 506             ELSE IF( WNTQO ) THEN
 507 *
 508 *              Path 2 (M much larger than N, JOBZ = 'O')
 509 *              N left singular vectors to be overwritten on A and
 510 *              N right singular vectors to be computed in VT
 511 *
 512                IR = 1
 513 *
 514 *              WORK(IR) is LDWRKR by N
 515 *
 516                IF( LWORK.GE.LDA*N+N*N+3*N+BDSPAC ) THEN
 517                   LDWRKR = LDA
 518                ELSE
 519                   LDWRKR = ( LWORK-N*N-3*N-BDSPAC ) / N
 520                END IF
 521                ITAU = IR + LDWRKR*N
 522                NWORK = ITAU + N
 523 *
 524 *              Compute A=Q*R
 525 *              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
 526 *
 527                CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
 528      $                      LWORK-NWORK+1, IERR )
 529 *
 530 *              Copy R to WORK(IR), zeroing out below it
 531 *
 532                CALL DLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
 533                CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, WORK( IR+1 ),
 534      $                      LDWRKR )
 535 *
 536 *              Generate Q in A
 537 *              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
 538 *
 539                CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
 540      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
 541                IE = ITAU
 542                ITAUQ = IE + N
 543                ITAUP = ITAUQ + N
 544                NWORK = ITAUP + N
 545 *
 546 *              Bidiagonalize R in VT, copying result to WORK(IR)
 547 *              (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
 548 *
 549                CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ),
 550      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
 551      $                      LWORK-NWORK+1, IERR )
 552 *
 553 *              WORK(IU) is N by N
 554 *
 555                IU = NWORK
 556                NWORK = IU + N*N
 557 *
 558 *              Perform bidiagonal SVD, computing left singular vectors
 559 *              of bidiagonal matrix in WORK(IU) and computing right
 560 *              singular vectors of bidiagonal matrix in VT
 561 *              (Workspace: need N+N*N+BDSPAC)
 562 *
 563                CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), N,
 564      $                      VT, LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
 565      $                      INFO )
 566 *
 567 *              Overwrite WORK(IU) by left singular vectors of R
 568 *              and VT by right singular vectors of R
 569 *              (Workspace: need 2*N*N+3*N, prefer 2*N*N+2*N+N*NB)
 570 *
 571                CALL DORMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
 572      $                      WORK( ITAUQ ), WORK( IU ), N, WORK( NWORK ),
 573      $                      LWORK-NWORK+1, IERR )
 574                CALL DORMBR( 'P', 'R', 'T', N, N, N, WORK( IR ), LDWRKR,
 575      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
 576      $                      LWORK-NWORK+1, IERR )
 577 *
 578 *              Multiply Q in A by left singular vectors of R in
 579 *              WORK(IU), storing result in WORK(IR) and copying to A
 580 *              (Workspace: need 2*N*N, prefer N*N+M*N)
 581 *
 582                DO 10 I = 1, M, LDWRKR
 583                   CHUNK = MIN( M-I+1, LDWRKR )
 584                   CALL DGEMM( 'N', 'N', CHUNK, N, N, ONE, A( I, 1 ),
 585      $                        LDA, WORK( IU ), N, ZERO, WORK( IR ),
 586      $                        LDWRKR )
 587                   CALL DLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR,
 588      $                         A( I, 1 ), LDA )
 589    10          CONTINUE
 590 *
 591             ELSE IF( WNTQS ) THEN
 592 *
 593 *              Path 3 (M much larger than N, JOBZ='S')
 594 *              N left singular vectors to be computed in U and
 595 *              N right singular vectors to be computed in VT
 596 *
 597                IR = 1
 598 *
 599 *              WORK(IR) is N by N
 600 *
 601                LDWRKR = N
 602                ITAU = IR + LDWRKR*N
 603                NWORK = ITAU + N
 604 *
 605 *              Compute A=Q*R
 606 *              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
 607 *
 608                CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
 609      $                      LWORK-NWORK+1, IERR )
 610 *
 611 *              Copy R to WORK(IR), zeroing out below it
 612 *
 613                CALL DLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
 614                CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, WORK( IR+1 ),
 615      $                      LDWRKR )
 616 *
 617 *              Generate Q in A
 618 *              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
 619 *
 620                CALL DORGQR( M, N, N, A, LDA, WORK( ITAU ),
 621      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
 622                IE = ITAU
 623                ITAUQ = IE + N
 624                ITAUP = ITAUQ + N
 625                NWORK = ITAUP + N
 626 *
 627 *              Bidiagonalize R in WORK(IR)
 628 *              (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
 629 *
 630                CALL DGEBRD( N, N, WORK( IR ), LDWRKR, S, WORK( IE ),
 631      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
 632      $                      LWORK-NWORK+1, IERR )
 633 *
 634 *              Perform bidiagonal SVD, computing left singular vectors
 635 *              of bidiagoal matrix in U and computing right singular
 636 *              vectors of bidiagonal matrix in VT
 637 *              (Workspace: need N+BDSPAC)
 638 *
 639                CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
 640      $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
 641      $                      INFO )
 642 *
 643 *              Overwrite U by left singular vectors of R and VT
 644 *              by right singular vectors of R
 645 *              (Workspace: need N*N+3*N, prefer N*N+2*N+N*NB)
 646 *
 647                CALL DORMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
 648      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
 649      $                      LWORK-NWORK+1, IERR )
 650 *
 651                CALL DORMBR( 'P', 'R', 'T', N, N, N, WORK( IR ), LDWRKR,
 652      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
 653      $                      LWORK-NWORK+1, IERR )
 654 *
 655 *              Multiply Q in A by left singular vectors of R in
 656 *              WORK(IR), storing result in U
 657 *              (Workspace: need N*N)
 658 *
 659                CALL DLACPY( 'F', N, N, U, LDU, WORK( IR ), LDWRKR )
 660                CALL DGEMM( 'N', 'N', M, N, N, ONE, A, LDA, WORK( IR ),
 661      $                     LDWRKR, ZERO, U, LDU )
 662 *
 663             ELSE IF( WNTQA ) THEN
 664 *
 665 *              Path 4 (M much larger than N, JOBZ='A')
 666 *              M left singular vectors to be computed in U and
 667 *              N right singular vectors to be computed in VT
 668 *
 669                IU = 1
 670 *
 671 *              WORK(IU) is N by N
 672 *
 673                LDWRKU = N
 674                ITAU = IU + LDWRKU*N
 675                NWORK = ITAU + N
 676 *
 677 *              Compute A=Q*R, copying result to U
 678 *              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
 679 *
 680                CALL DGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
 681      $                      LWORK-NWORK+1, IERR )
 682                CALL DLACPY( 'L', M, N, A, LDA, U, LDU )
 683 *
 684 *              Generate Q in U
 685 *              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
 686                CALL DORGQR( M, M, N, U, LDU, WORK( ITAU ),
 687      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
 688 *
 689 *              Produce R in A, zeroing out other entries
 690 *
 691                CALL DLASET( 'L', N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
 692                IE = ITAU
 693                ITAUQ = IE + N
 694                ITAUP = ITAUQ + N
 695                NWORK = ITAUP + N
 696 *
 697 *              Bidiagonalize R in A
 698 *              (Workspace: need N*N+4*N, prefer N*N+3*N+2*N*NB)
 699 *
 700                CALL DGEBRD( N, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
 701      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
 702      $                      IERR )
 703 *
 704 *              Perform bidiagonal SVD, computing left singular vectors
 705 *              of bidiagonal matrix in WORK(IU) and computing right
 706 *              singular vectors of bidiagonal matrix in VT
 707 *              (Workspace: need N+N*N+BDSPAC)
 708 *
 709                CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ), N,
 710      $                      VT, LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
 711      $                      INFO )
 712 *
 713 *              Overwrite WORK(IU) by left singular vectors of R and VT
 714 *              by right singular vectors of R
 715 *              (Workspace: need N*N+3*N, prefer N*N+2*N+N*NB)
 716 *
 717                CALL DORMBR( 'Q', 'L', 'N', N, N, N, A, LDA,
 718      $                      WORK( ITAUQ ), WORK( IU ), LDWRKU,
 719      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
 720                CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
 721      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
 722      $                      LWORK-NWORK+1, IERR )
 723 *
 724 *              Multiply Q in U by left singular vectors of R in
 725 *              WORK(IU), storing result in A
 726 *              (Workspace: need N*N)
 727 *
 728                CALL DGEMM( 'N', 'N', M, N, N, ONE, U, LDU, WORK( IU ),
 729      $                     LDWRKU, ZERO, A, LDA )
 730 *
 731 *              Copy left singular vectors of A from A to U
 732 *
 733                CALL DLACPY( 'F', M, N, A, LDA, U, LDU )
 734 *
 735             END IF
 736 *
 737          ELSE
 738 *
 739 *           M .LT. MNTHR
 740 *
 741 *           Path 5 (M at least N, but not much larger)
 742 *           Reduce to bidiagonal form without QR decomposition
 743 *
 744             IE = 1
 745             ITAUQ = IE + N
 746             ITAUP = ITAUQ + N
 747             NWORK = ITAUP + N
 748 *
 749 *           Bidiagonalize A
 750 *           (Workspace: need 3*N+M, prefer 3*N+(M+N)*NB)
 751 *
 752             CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
 753      $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
 754      $                   IERR )
 755             IF( WNTQN ) THEN
 756 *
 757 *              Perform bidiagonal SVD, only computing singular values
 758 *              (Workspace: need N+BDSPAC)
 759 *
 760                CALL DBDSDC( 'U', 'N', N, S, WORK( IE ), DUM, 1, DUM, 1,
 761      $                      DUM, IDUM, WORK( NWORK ), IWORK, INFO )
 762             ELSE IF( WNTQO ) THEN
 763                IU = NWORK
 764                IF( LWORK.GE.M*N+3*N+BDSPAC ) THEN
 765 *
 766 *                 WORK( IU ) is M by N
 767 *
 768                   LDWRKU = M
 769                   NWORK = IU + LDWRKU*N
 770                   CALL DLASET( 'F', M, N, ZERO, ZERO, WORK( IU ),
 771      $                         LDWRKU )
 772                ELSE
 773 *
 774 *                 WORK( IU ) is N by N
 775 *
 776                   LDWRKU = N
 777                   NWORK = IU + LDWRKU*N
 778 *
 779 *                 WORK(IR) is LDWRKR by N
 780 *
 781                   IR = NWORK
 782                   LDWRKR = ( LWORK-N*N-3*N ) / N
 783                END IF
 784                NWORK = IU + LDWRKU*N
 785 *
 786 *              Perform bidiagonal SVD, computing left singular vectors
 787 *              of bidiagonal matrix in WORK(IU) and computing right
 788 *              singular vectors of bidiagonal matrix in VT
 789 *              (Workspace: need N+N*N+BDSPAC)
 790 *
 791                CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), WORK( IU ),
 792      $                      LDWRKU, VT, LDVT, DUM, IDUM, WORK( NWORK ),
 793      $                      IWORK, INFO )
 794 *
 795 *              Overwrite VT by right singular vectors of A
 796 *              (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
 797 *
 798                CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
 799      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
 800      $                      LWORK-NWORK+1, IERR )
 801 *
 802                IF( LWORK.GE.M*N+3*N+BDSPAC ) THEN
 803 *
 804 *                 Overwrite WORK(IU) by left singular vectors of A
 805 *                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
 806 *
 807                   CALL DORMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
 808      $                         WORK( ITAUQ ), WORK( IU ), LDWRKU,
 809      $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
 810 *
 811 *                 Copy left singular vectors of A from WORK(IU) to A
 812 *
 813                   CALL DLACPY( 'F', M, N, WORK( IU ), LDWRKU, A, LDA )
 814                ELSE
 815 *
 816 *                 Generate Q in A
 817 *                 (Workspace: need N*N+2*N, prefer N*N+N+N*NB)
 818 *
 819                   CALL DORGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
 820      $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
 821 *
 822 *                 Multiply Q in A by left singular vectors of
 823 *                 bidiagonal matrix in WORK(IU), storing result in
 824 *                 WORK(IR) and copying to A
 825 *                 (Workspace: need 2*N*N, prefer N*N+M*N)
 826 *
 827                   DO 20 I = 1, M, LDWRKR
 828                      CHUNK = MIN( M-I+1, LDWRKR )
 829                      CALL DGEMM( 'N', 'N', CHUNK, N, N, ONE, A( I, 1 ),
 830      $                           LDA, WORK( IU ), LDWRKU, ZERO,
 831      $                           WORK( IR ), LDWRKR )
 832                      CALL DLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR,
 833      $                            A( I, 1 ), LDA )
 834    20             CONTINUE
 835                END IF
 836 *
 837             ELSE IF( WNTQS ) THEN
 838 *
 839 *              Perform bidiagonal SVD, computing left singular vectors
 840 *              of bidiagonal matrix in U and computing right singular
 841 *              vectors of bidiagonal matrix in VT
 842 *              (Workspace: need N+BDSPAC)
 843 *
 844                CALL DLASET( 'F', M, N, ZERO, ZERO, U, LDU )
 845                CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
 846      $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
 847      $                      INFO )
 848 *
 849 *              Overwrite U by left singular vectors of A and VT
 850 *              by right singular vectors of A
 851 *              (Workspace: need 3*N, prefer 2*N+N*NB)
 852 *
 853                CALL DORMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
 854      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
 855      $                      LWORK-NWORK+1, IERR )
 856                CALL DORMBR( 'P', 'R', 'T', N, N, N, A, LDA,
 857      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
 858      $                      LWORK-NWORK+1, IERR )
 859             ELSE IF( WNTQA ) THEN
 860 *
 861 *              Perform bidiagonal SVD, computing left singular vectors
 862 *              of bidiagonal matrix in U and computing right singular
 863 *              vectors of bidiagonal matrix in VT
 864 *              (Workspace: need N+BDSPAC)
 865 *
 866                CALL DLASET( 'F', M, M, ZERO, ZERO, U, LDU )
 867                CALL DBDSDC( 'U', 'I', N, S, WORK( IE ), U, LDU, VT,
 868      $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
 869      $                      INFO )
 870 *
 871 *              Set the right corner of U to identity matrix
 872 *
 873                IF( M.GT.N ) THEN
 874                   CALL DLASET( 'F', M-N, M-N, ZERO, ONE, U( N+1, N+1 ),
 875      $                         LDU )
 876                END IF
 877 *
 878 *              Overwrite U by left singular vectors of A and VT
 879 *              by right singular vectors of A
 880 *              (Workspace: need N*N+2*N+M, prefer N*N+2*N+M*NB)
 881 *
 882                CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
 883      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
 884      $                      LWORK-NWORK+1, IERR )
 885                CALL DORMBR( 'P', 'R', 'T', N, N, M, A, LDA,
 886      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
 887      $                      LWORK-NWORK+1, IERR )
 888             END IF
 889 *
 890          END IF
 891 *
 892       ELSE
 893 *
 894 *        A has more columns than rows. If A has sufficiently more
 895 *        columns than rows, first reduce using the LQ decomposition (if
 896 *        sufficient workspace available)
 897 *
 898          IF( N.GE.MNTHR ) THEN
 899 *
 900             IF( WNTQN ) THEN
 901 *
 902 *              Path 1t (N much larger than M, JOBZ='N')
 903 *              No singular vectors to be computed
 904 *
 905                ITAU = 1
 906                NWORK = ITAU + M
 907 *
 908 *              Compute A=L*Q
 909 *              (Workspace: need 2*M, prefer M+M*NB)
 910 *
 911                CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
 912      $                      LWORK-NWORK+1, IERR )
 913 *
 914 *              Zero out above L
 915 *
 916                CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ), LDA )
 917                IE = 1
 918                ITAUQ = IE + M
 919                ITAUP = ITAUQ + M
 920                NWORK = ITAUP + M
 921 *
 922 *              Bidiagonalize L in A
 923 *              (Workspace: need 4*M, prefer 3*M+2*M*NB)
 924 *
 925                CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
 926      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
 927      $                      IERR )
 928                NWORK = IE + M
 929 *
 930 *              Perform bidiagonal SVD, computing singular values only
 931 *              (Workspace: need M+BDSPAC)
 932 *
 933                CALL DBDSDC( 'U', 'N', M, S, WORK( IE ), DUM, 1, DUM, 1,
 934      $                      DUM, IDUM, WORK( NWORK ), IWORK, INFO )
 935 *
 936             ELSE IF( WNTQO ) THEN
 937 *
 938 *              Path 2t (N much larger than M, JOBZ='O')
 939 *              M right singular vectors to be overwritten on A and
 940 *              M left singular vectors to be computed in U
 941 *
 942                IVT = 1
 943 *
 944 *              IVT is M by M
 945 *
 946                IL = IVT + M*M
 947                IF( LWORK.GE.M*N+M*M+3*M+BDSPAC ) THEN
 948 *
 949 *                 WORK(IL) is M by N
 950 *
 951                   LDWRKL = M
 952                   CHUNK = N
 953                ELSE
 954                   LDWRKL = M
 955                   CHUNK = ( LWORK-M*M ) / M
 956                END IF
 957                ITAU = IL + LDWRKL*M
 958                NWORK = ITAU + M
 959 *
 960 *              Compute A=L*Q
 961 *              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
 962 *
 963                CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
 964      $                      LWORK-NWORK+1, IERR )
 965 *
 966 *              Copy L to WORK(IL), zeroing about above it
 967 *
 968                CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
 969                CALL DLASET( 'U', M-1, M-1, ZERO, ZERO,
 970      $                      WORK( IL+LDWRKL ), LDWRKL )
 971 *
 972 *              Generate Q in A
 973 *              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
 974 *
 975                CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
 976      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
 977                IE = ITAU
 978                ITAUQ = IE + M
 979                ITAUP = ITAUQ + M
 980                NWORK = ITAUP + M
 981 *
 982 *              Bidiagonalize L in WORK(IL)
 983 *              (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
 984 *
 985                CALL DGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ),
 986      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
 987      $                      LWORK-NWORK+1, IERR )
 988 *
 989 *              Perform bidiagonal SVD, computing left singular vectors
 990 *              of bidiagonal matrix in U, and computing right singular
 991 *              vectors of bidiagonal matrix in WORK(IVT)
 992 *              (Workspace: need M+M*M+BDSPAC)
 993 *
 994                CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU,
 995      $                      WORK( IVT ), M, DUM, IDUM, WORK( NWORK ),
 996      $                      IWORK, INFO )
 997 *
 998 *              Overwrite U by left singular vectors of L and WORK(IVT)
 999 *              by right singular vectors of L
1000 *              (Workspace: need 2*M*M+3*M, prefer 2*M*M+2*M+M*NB)
1001 *
1002                CALL DORMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
1003      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1004      $                      LWORK-NWORK+1, IERR )
1005                CALL DORMBR( 'P', 'R', 'T', M, M, M, WORK( IL ), LDWRKL,
1006      $                      WORK( ITAUP ), WORK( IVT ), M,
1007      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1008 *
1009 *              Multiply right singular vectors of L in WORK(IVT) by Q
1010 *              in A, storing result in WORK(IL) and copying to A
1011 *              (Workspace: need 2*M*M, prefer M*M+M*N)
1012 *
1013                DO 30 I = 1, N, CHUNK
1014                   BLK = MIN( N-I+1, CHUNK )
1015                   CALL DGEMM( 'N', 'N', M, BLK, M, ONE, WORK( IVT ), M,
1016      $                        A( 1, I ), LDA, ZERO, WORK( IL ), LDWRKL )
1017                   CALL DLACPY( 'F', M, BLK, WORK( IL ), LDWRKL,
1018      $                         A( 1, I ), LDA )
1019    30          CONTINUE
1020 *
1021             ELSE IF( WNTQS ) THEN
1022 *
1023 *              Path 3t (N much larger than M, JOBZ='S')
1024 *              M right singular vectors to be computed in VT and
1025 *              M left singular vectors to be computed in U
1026 *
1027                IL = 1
1028 *
1029 *              WORK(IL) is M by M
1030 *
1031                LDWRKL = M
1032                ITAU = IL + LDWRKL*M
1033                NWORK = ITAU + M
1034 *
1035 *              Compute A=L*Q
1036 *              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1037 *
1038                CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1039      $                      LWORK-NWORK+1, IERR )
1040 *
1041 *              Copy L to WORK(IL), zeroing out above it
1042 *
1043                CALL DLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
1044                CALL DLASET( 'U', M-1, M-1, ZERO, ZERO,
1045      $                      WORK( IL+LDWRKL ), LDWRKL )
1046 *
1047 *              Generate Q in A
1048 *              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1049 *
1050                CALL DORGLQ( M, N, M, A, LDA, WORK( ITAU ),
1051      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1052                IE = ITAU
1053                ITAUQ = IE + M
1054                ITAUP = ITAUQ + M
1055                NWORK = ITAUP + M
1056 *
1057 *              Bidiagonalize L in WORK(IU), copying result to U
1058 *              (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
1059 *
1060                CALL DGEBRD( M, M, WORK( IL ), LDWRKL, S, WORK( IE ),
1061      $                      WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
1062      $                      LWORK-NWORK+1, IERR )
1063 *
1064 *              Perform bidiagonal SVD, computing left singular vectors
1065 *              of bidiagonal matrix in U and computing right singular
1066 *              vectors of bidiagonal matrix in VT
1067 *              (Workspace: need M+BDSPAC)
1068 *
1069                CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU, VT,
1070      $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
1071      $                      INFO )
1072 *
1073 *              Overwrite U by left singular vectors of L and VT
1074 *              by right singular vectors of L
1075 *              (Workspace: need M*M+3*M, prefer M*M+2*M+M*NB)
1076 *
1077                CALL DORMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
1078      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1079      $                      LWORK-NWORK+1, IERR )
1080                CALL DORMBR( 'P', 'R', 'T', M, M, M, WORK( IL ), LDWRKL,
1081      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1082      $                      LWORK-NWORK+1, IERR )
1083 *
1084 *              Multiply right singular vectors of L in WORK(IL) by
1085 *              Q in A, storing result in VT
1086 *              (Workspace: need M*M)
1087 *
1088                CALL DLACPY( 'F', M, M, VT, LDVT, WORK( IL ), LDWRKL )
1089                CALL DGEMM( 'N', 'N', M, N, M, ONE, WORK( IL ), LDWRKL,
1090      $                     A, LDA, ZERO, VT, LDVT )
1091 *
1092             ELSE IF( WNTQA ) THEN
1093 *
1094 *              Path 4t (N much larger than M, JOBZ='A')
1095 *              N right singular vectors to be computed in VT and
1096 *              M left singular vectors to be computed in U
1097 *
1098                IVT = 1
1099 *
1100 *              WORK(IVT) is M by M
1101 *
1102                LDWKVT = M
1103                ITAU = IVT + LDWKVT*M
1104                NWORK = ITAU + M
1105 *
1106 *              Compute A=L*Q, copying result to VT
1107 *              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1108 *
1109                CALL DGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1110      $                      LWORK-NWORK+1, IERR )
1111                CALL DLACPY( 'U', M, N, A, LDA, VT, LDVT )
1112 *
1113 *              Generate Q in VT
1114 *              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1115 *
1116                CALL DORGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
1117      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1118 *
1119 *              Produce L in A, zeroing out other entries
1120 *
1121                CALL DLASET( 'U', M-1, M-1, ZERO, ZERO, A( 1, 2 ), LDA )
1122                IE = ITAU
1123                ITAUQ = IE + M
1124                ITAUP = ITAUQ + M
1125                NWORK = ITAUP + M
1126 *
1127 *              Bidiagonalize L in A
1128 *              (Workspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
1129 *
1130                CALL DGEBRD( M, M, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
1131      $                      WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1132      $                      IERR )
1133 *
1134 *              Perform bidiagonal SVD, computing left singular vectors
1135 *              of bidiagonal matrix in U and computing right singular
1136 *              vectors of bidiagonal matrix in WORK(IVT)
1137 *              (Workspace: need M+M*M+BDSPAC)
1138 *
1139                CALL DBDSDC( 'U', 'I', M, S, WORK( IE ), U, LDU,
1140      $                      WORK( IVT ), LDWKVT, DUM, IDUM,
1141      $                      WORK( NWORK ), IWORK, INFO )
1142 *
1143 *              Overwrite U by left singular vectors of L and WORK(IVT)
1144 *              by right singular vectors of L
1145 *              (Workspace: need M*M+3*M, prefer M*M+2*M+M*NB)
1146 *
1147                CALL DORMBR( 'Q', 'L', 'N', M, M, M, A, LDA,
1148      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1149      $                      LWORK-NWORK+1, IERR )
1150                CALL DORMBR( 'P', 'R', 'T', M, M, M, A, LDA,
1151      $                      WORK( ITAUP ), WORK( IVT ), LDWKVT,
1152      $                      WORK( NWORK ), LWORK-NWORK+1, IERR )
1153 *
1154 *              Multiply right singular vectors of L in WORK(IVT) by
1155 *              Q in VT, storing result in A
1156 *              (Workspace: need M*M)
1157 *
1158                CALL DGEMM( 'N', 'N', M, N, M, ONE, WORK( IVT ), LDWKVT,
1159      $                     VT, LDVT, ZERO, A, LDA )
1160 *
1161 *              Copy right singular vectors of A from A to VT
1162 *
1163                CALL DLACPY( 'F', M, N, A, LDA, VT, LDVT )
1164 *
1165             END IF
1166 *
1167          ELSE
1168 *
1169 *           N .LT. MNTHR
1170 *
1171 *           Path 5t (N greater than M, but not much larger)
1172 *           Reduce to bidiagonal form without LQ decomposition
1173 *
1174             IE = 1
1175             ITAUQ = IE + M
1176             ITAUP = ITAUQ + M
1177             NWORK = ITAUP + M
1178 *
1179 *           Bidiagonalize A
1180 *           (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB)
1181 *
1182             CALL DGEBRD( M, N, A, LDA, S, WORK( IE ), WORK( ITAUQ ),
1183      $                   WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1184      $                   IERR )
1185             IF( WNTQN ) THEN
1186 *
1187 *              Perform bidiagonal SVD, only computing singular values
1188 *              (Workspace: need M+BDSPAC)
1189 *
1190                CALL DBDSDC( 'L', 'N', M, S, WORK( IE ), DUM, 1, DUM, 1,
1191      $                      DUM, IDUM, WORK( NWORK ), IWORK, INFO )
1192             ELSE IF( WNTQO ) THEN
1193                LDWKVT = M
1194                IVT = NWORK
1195                IF( LWORK.GE.M*N+3*M+BDSPAC ) THEN
1196 *
1197 *                 WORK( IVT ) is M by N
1198 *
1199                   CALL DLASET( 'F', M, N, ZERO, ZERO, WORK( IVT ),
1200      $                         LDWKVT )
1201                   NWORK = IVT + LDWKVT*N
1202                ELSE
1203 *
1204 *                 WORK( IVT ) is M by M
1205 *
1206                   NWORK = IVT + LDWKVT*M
1207                   IL = NWORK
1208 *
1209 *                 WORK(IL) is M by CHUNK
1210 *
1211                   CHUNK = ( LWORK-M*M-3*M ) / M
1212                END IF
1213 *
1214 *              Perform bidiagonal SVD, computing left singular vectors
1215 *              of bidiagonal matrix in U and computing right singular
1216 *              vectors of bidiagonal matrix in WORK(IVT)
1217 *              (Workspace: need M*M+BDSPAC)
1218 *
1219                CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU,
1220      $                      WORK( IVT ), LDWKVT, DUM, IDUM,
1221      $                      WORK( NWORK ), IWORK, INFO )
1222 *
1223 *              Overwrite U by left singular vectors of A
1224 *              (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1225 *
1226                CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1227      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1228      $                      LWORK-NWORK+1, IERR )
1229 *
1230                IF( LWORK.GE.M*N+3*M+BDSPAC ) THEN
1231 *
1232 *                 Overwrite WORK(IVT) by left singular vectors of A
1233 *                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1234 *
1235                   CALL DORMBR( 'P', 'R', 'T', M, N, M, A, LDA,
1236      $                         WORK( ITAUP ), WORK( IVT ), LDWKVT,
1237      $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
1238 *
1239 *                 Copy right singular vectors of A from WORK(IVT) to A
1240 *
1241                   CALL DLACPY( 'F', M, N, WORK( IVT ), LDWKVT, A, LDA )
1242                ELSE
1243 *
1244 *                 Generate P**T in A
1245 *                 (Workspace: need M*M+2*M, prefer M*M+M+M*NB)
1246 *
1247                   CALL DORGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
1248      $                         WORK( NWORK ), LWORK-NWORK+1, IERR )
1249 *
1250 *                 Multiply Q in A by right singular vectors of
1251 *                 bidiagonal matrix in WORK(IVT), storing result in
1252 *                 WORK(IL) and copying to A
1253 *                 (Workspace: need 2*M*M, prefer M*M+M*N)
1254 *
1255                   DO 40 I = 1, N, CHUNK
1256                      BLK = MIN( N-I+1, CHUNK )
1257                      CALL DGEMM( 'N', 'N', M, BLK, M, ONE, WORK( IVT ),
1258      $                           LDWKVT, A( 1, I ), LDA, ZERO,
1259      $                           WORK( IL ), M )
1260                      CALL DLACPY( 'F', M, BLK, WORK( IL ), M, A( 1, I ),
1261      $                            LDA )
1262    40             CONTINUE
1263                END IF
1264             ELSE IF( WNTQS ) THEN
1265 *
1266 *              Perform bidiagonal SVD, computing left singular vectors
1267 *              of bidiagonal matrix in U and computing right singular
1268 *              vectors of bidiagonal matrix in VT
1269 *              (Workspace: need M+BDSPAC)
1270 *
1271                CALL DLASET( 'F', M, N, ZERO, ZERO, VT, LDVT )
1272                CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, VT,
1273      $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
1274      $                      INFO )
1275 *
1276 *              Overwrite U by left singular vectors of A and VT
1277 *              by right singular vectors of A
1278 *              (Workspace: need 3*M, prefer 2*M+M*NB)
1279 *
1280                CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1281      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1282      $                      LWORK-NWORK+1, IERR )
1283                CALL DORMBR( 'P', 'R', 'T', M, N, M, A, LDA,
1284      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1285      $                      LWORK-NWORK+1, IERR )
1286             ELSE IF( WNTQA ) THEN
1287 *
1288 *              Perform bidiagonal SVD, computing left singular vectors
1289 *              of bidiagonal matrix in U and computing right singular
1290 *              vectors of bidiagonal matrix in VT
1291 *              (Workspace: need M+BDSPAC)
1292 *
1293                CALL DLASET( 'F', N, N, ZERO, ZERO, VT, LDVT )
1294                CALL DBDSDC( 'L', 'I', M, S, WORK( IE ), U, LDU, VT,
1295      $                      LDVT, DUM, IDUM, WORK( NWORK ), IWORK,
1296      $                      INFO )
1297 *
1298 *              Set the right corner of VT to identity matrix
1299 *
1300                IF( N.GT.M ) THEN
1301                   CALL DLASET( 'F', N-M, N-M, ZERO, ONE, VT( M+1, M+1 ),
1302      $                         LDVT )
1303                END IF
1304 *
1305 *              Overwrite U by left singular vectors of A and VT
1306 *              by right singular vectors of A
1307 *              (Workspace: need 2*M+N, prefer 2*M+N*NB)
1308 *
1309                CALL DORMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1310      $                      WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1311      $                      LWORK-NWORK+1, IERR )
1312                CALL DORMBR( 'P', 'R', 'T', N, N, M, A, LDA,
1313      $                      WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1314      $                      LWORK-NWORK+1, IERR )
1315             END IF
1316 *
1317          END IF
1318 *
1319       END IF
1320 *
1321 *     Undo scaling if necessary
1322 *
1323       IF( ISCL.EQ.1 ) THEN
1324          IF( ANRM.GT.BIGNUM )
1325      $      CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
1326      $                   IERR )
1327          IF( ANRM.LT.SMLNUM )
1328      $      CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
1329      $                   IERR )
1330       END IF
1331 *
1332 *     Return optimal workspace in WORK(1)
1333 *
1334       WORK( 1 ) = MAXWRK
1335 *
1336       RETURN
1337 *
1338 *     End of DGESDD
1339 *
1340       END
Home