CoCalc -- cgesvd.f

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

3620
Product

Resources

Company
