Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
ElmerCSC
GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/mathlibs/src/lapack/cgesdd.f
5191 views
1
SUBROUTINE CGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK,
2
$ LWORK, RWORK, IWORK, 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 JOBZ
11
INTEGER INFO, LDA, LDU, LDVT, LWORK, M, N
12
* ..
13
* .. Array Arguments ..
14
INTEGER IWORK( * )
15
REAL RWORK( * ), S( * )
16
COMPLEX A( LDA, * ), U( LDU, * ), VT( LDVT, * ),
17
$ WORK( * )
18
* ..
19
*
20
* Purpose
21
* =======
22
*
23
* CGESDD computes the singular value decomposition (SVD) of a complex
24
* M-by-N matrix A, optionally computing the left and/or right singular
25
* vectors, by using divide-and-conquer method. The SVD is written
26
*
27
* A = U * SIGMA * conjugate-transpose(V)
28
*
29
* where SIGMA is an M-by-N matrix which is zero except for its
30
* min(m,n) diagonal elements, U is an M-by-M unitary matrix, and
31
* V is an N-by-N unitary matrix. The diagonal elements of SIGMA
32
* are the singular values of A; they are real and non-negative, and
33
* are returned in descending order. The first min(m,n) columns of
34
* U and V are the left and right singular vectors of A.
35
*
36
* Note that the routine returns VT = V**H, not V.
37
*
38
* The divide and conquer algorithm makes very mild assumptions about
39
* floating point arithmetic. It will work on machines with a guard
40
* digit in add/subtract, or on those binary machines without guard
41
* digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or
42
* Cray-2. It could conceivably fail on hexadecimal or decimal machines
43
* without guard digits, but we know of none.
44
*
45
* Arguments
46
* =========
47
*
48
* JOBZ (input) CHARACTER*1
49
* Specifies options for computing all or part of the matrix U:
50
* = 'A': all M columns of U and all N rows of V**H are
51
* returned in the arrays U and VT;
52
* = 'S': the first min(M,N) columns of U and the first
53
* min(M,N) rows of V**H are returned in the arrays U
54
* and VT;
55
* = 'O': If M >= N, the first N columns of U are overwritten
56
* on the array A and all rows of V**H are returned in
57
* the array VT;
58
* otherwise, all columns of U are returned in the
59
* array U and the first M rows of V**H are overwritten
60
* in the array VT;
61
* = 'N': no columns of U or rows of V**H are computed.
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 JOBZ = 'O', A is overwritten with the first N columns
73
* of U (the left singular vectors, stored
74
* columnwise) if M >= N;
75
* A is overwritten with the first M rows
76
* of V**H (the right singular vectors, stored
77
* rowwise) otherwise.
78
* if JOBZ .ne. 'O', the contents of A are destroyed.
79
*
80
* LDA (input) INTEGER
81
* The leading dimension of the array A. LDA >= max(1,M).
82
*
83
* S (output) REAL array, dimension (min(M,N))
84
* The singular values of A, sorted so that S(i) >= S(i+1).
85
*
86
* U (output) COMPLEX array, dimension (LDU,UCOL)
87
* UCOL = M if JOBZ = 'A' or JOBZ = 'O' and M < N;
88
* UCOL = min(M,N) if JOBZ = 'S'.
89
* If JOBZ = 'A' or JOBZ = 'O' and M < N, U contains the M-by-M
90
* unitary matrix U;
91
* if JOBZ = 'S', U contains the first min(M,N) columns of U
92
* (the left singular vectors, stored columnwise);
93
* if JOBZ = 'O' and M >= N, or JOBZ = 'N', U is not referenced.
94
*
95
* LDU (input) INTEGER
96
* The leading dimension of the array U. LDU >= 1; if
97
* JOBZ = 'S' or 'A' or JOBZ = 'O' and M < N, LDU >= M.
98
*
99
* VT (output) COMPLEX array, dimension (LDVT,N)
100
* If JOBZ = 'A' or JOBZ = 'O' and M >= N, VT contains the
101
* N-by-N unitary matrix V**H;
102
* if JOBZ = 'S', VT contains the first min(M,N) rows of
103
* V**H (the right singular vectors, stored rowwise);
104
* if JOBZ = 'O' and M < N, or JOBZ = 'N', VT is not referenced.
105
*
106
* LDVT (input) INTEGER
107
* The leading dimension of the array VT. LDVT >= 1; if
108
* JOBZ = 'A' or JOBZ = 'O' and M >= N, LDVT >= N;
109
* if JOBZ = 'S', LDVT >= min(M,N).
110
*
111
* WORK (workspace/output) COMPLEX array, dimension (LWORK)
112
* On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
113
*
114
* LWORK (input) INTEGER
115
* The dimension of the array WORK. LWORK >= 1.
116
* if JOBZ = 'N', LWORK >= 2*min(M,N)+max(M,N).
117
* if JOBZ = 'O',
118
* LWORK >= 2*min(M,N)*min(M,N)+2*min(M,N)+max(M,N).
119
* if JOBZ = 'S' or 'A',
120
* LWORK >= min(M,N)*min(M,N)+2*min(M,N)+max(M,N).
121
* For good performance, LWORK should generally be larger.
122
* If LWORK < 0 but other input arguments are legal, WORK(1)
123
* returns the optimal LWORK.
124
*
125
* RWORK (workspace) REAL array, dimension (LRWORK)
126
* If JOBZ = 'N', LRWORK >= 7*min(M,N).
127
* Otherwise, LRWORK >= 5*min(M,N)*min(M,N) + 5*min(M,N)
128
*
129
* IWORK (workspace) INTEGER array, dimension (8*min(M,N))
130
*
131
* INFO (output) INTEGER
132
* = 0: successful exit.
133
* < 0: if INFO = -i, the i-th argument had an illegal value.
134
* > 0: The updating process of SBDSDC did not converge.
135
*
136
* Further Details
137
* ===============
138
*
139
* Based on contributions by
140
* Ming Gu and Huan Ren, Computer Science Division, University of
141
* California at Berkeley, USA
142
*
143
* =====================================================================
144
*
145
* .. Parameters ..
146
COMPLEX CZERO, CONE
147
PARAMETER ( CZERO = ( 0.0E0, 0.0E0 ),
148
$ CONE = ( 1.0E0, 0.0E0 ) )
149
REAL ZERO, ONE
150
PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 )
151
* ..
152
* .. Local Scalars ..
153
LOGICAL LQUERY, WNTQA, WNTQAS, WNTQN, WNTQO, WNTQS
154
INTEGER BLK, CHUNK, I, IE, IERR, IL, IR, IRU, IRVT,
155
$ ISCL, ITAU, ITAUP, ITAUQ, IU, IVT, LDWKVT,
156
$ LDWRKL, LDWRKR, LDWRKU, MAXWRK, MINMN, MINWRK,
157
$ MNTHR1, MNTHR2, NRWORK, NWORK, WRKBL
158
REAL ANRM, BIGNUM, EPS, SMLNUM
159
* ..
160
* .. Local Arrays ..
161
INTEGER IDUM( 1 )
162
REAL DUM( 1 )
163
* ..
164
* .. External Subroutines ..
165
EXTERNAL CGEBRD, CGELQF, CGEMM, CGEQRF, CLACP2, CLACPY,
166
$ CLACRM, CLARCM, CLASCL, CLASET, CUNGBR, CUNGLQ,
167
$ CUNGQR, CUNMBR, SBDSDC, SLASCL, XERBLA
168
* ..
169
* .. External Functions ..
170
LOGICAL LSAME
171
INTEGER ILAENV
172
REAL CLANGE, SLAMCH
173
EXTERNAL CLANGE, ILAENV, LSAME, SLAMCH
174
* ..
175
* .. Intrinsic Functions ..
176
INTRINSIC INT, MAX, MIN, SQRT
177
* ..
178
* .. Executable Statements ..
179
*
180
* Test the input arguments
181
*
182
INFO = 0
183
MINMN = MIN( M, N )
184
MNTHR1 = INT( MINMN*17.0E0 / 9.0E0 )
185
MNTHR2 = INT( MINMN*5.0E0 / 3.0E0 )
186
WNTQA = LSAME( JOBZ, 'A' )
187
WNTQS = LSAME( JOBZ, 'S' )
188
WNTQAS = WNTQA .OR. WNTQS
189
WNTQO = LSAME( JOBZ, 'O' )
190
WNTQN = LSAME( JOBZ, 'N' )
191
MINWRK = 1
192
MAXWRK = 1
193
LQUERY = ( LWORK.EQ.-1 )
194
*
195
IF( .NOT.( WNTQA .OR. WNTQS .OR. WNTQO .OR. WNTQN ) ) THEN
196
INFO = -1
197
ELSE IF( M.LT.0 ) THEN
198
INFO = -2
199
ELSE IF( N.LT.0 ) THEN
200
INFO = -3
201
ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
202
INFO = -5
203
ELSE IF( LDU.LT.1 .OR. ( WNTQAS .AND. LDU.LT.M ) .OR.
204
$ ( WNTQO .AND. M.LT.N .AND. LDU.LT.M ) ) THEN
205
INFO = -8
206
ELSE IF( LDVT.LT.1 .OR. ( WNTQA .AND. LDVT.LT.N ) .OR.
207
$ ( WNTQS .AND. LDVT.LT.MINMN ) .OR.
208
$ ( WNTQO .AND. M.GE.N .AND. LDVT.LT.N ) ) THEN
209
INFO = -10
210
END IF
211
*
212
* Compute workspace
213
* (Note: Comments in the code beginning "Workspace:" describe the
214
* minimal amount of workspace needed at that point in the code,
215
* as well as the preferred amount for good performance.
216
* CWorkspace refers to complex workspace, and RWorkspace to
217
* real workspace. NB refers to the optimal block size for the
218
* immediately following subroutine, as returned by ILAENV.)
219
*
220
IF( INFO.EQ.0 .AND. M.GT.0 .AND. N.GT.0 ) THEN
221
IF( M.GE.N ) THEN
222
*
223
* There is no complex work space needed for bidiagonal SVD
224
* The real work space needed for bidiagonal SVD is BDSPAC,
225
* BDSPAC = 3*N*N + 4*N
226
*
227
IF( M.GE.MNTHR1 ) THEN
228
IF( WNTQN ) THEN
229
*
230
* Path 1 (M much larger than N, JOBZ='N')
231
*
232
WRKBL = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1,
233
$ -1 )
234
WRKBL = MAX( WRKBL, 2*N+2*N*
235
$ ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) )
236
MAXWRK = WRKBL
237
MINWRK = 3*N
238
ELSE IF( WNTQO ) THEN
239
*
240
* Path 2 (M much larger than N, JOBZ='O')
241
*
242
WRKBL = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 )
243
WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'CUNGQR', ' ', M,
244
$ N, N, -1 ) )
245
WRKBL = MAX( WRKBL, 2*N+2*N*
246
$ ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) )
247
WRKBL = MAX( WRKBL, 2*N+N*
248
$ ILAENV( 1, 'CUNMBR', 'QLN', N, N, N, -1 ) )
249
WRKBL = MAX( WRKBL, 2*N+N*
250
$ ILAENV( 1, 'CUNMBR', 'PRC', N, N, N, -1 ) )
251
MAXWRK = M*N + N*N + WRKBL
252
MINWRK = 2*N*N + 3*N
253
ELSE IF( WNTQS ) THEN
254
*
255
* Path 3 (M much larger than N, JOBZ='S')
256
*
257
WRKBL = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 )
258
WRKBL = MAX( WRKBL, N+N*ILAENV( 1, 'CUNGQR', ' ', M,
259
$ N, N, -1 ) )
260
WRKBL = MAX( WRKBL, 2*N+2*N*
261
$ ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) )
262
WRKBL = MAX( WRKBL, 2*N+N*
263
$ ILAENV( 1, 'CUNMBR', 'QLN', N, N, N, -1 ) )
264
WRKBL = MAX( WRKBL, 2*N+N*
265
$ ILAENV( 1, 'CUNMBR', 'PRC', N, N, N, -1 ) )
266
MAXWRK = N*N + WRKBL
267
MINWRK = N*N + 3*N
268
ELSE IF( WNTQA ) THEN
269
*
270
* Path 4 (M much larger than N, JOBZ='A')
271
*
272
WRKBL = N + N*ILAENV( 1, 'CGEQRF', ' ', M, N, -1, -1 )
273
WRKBL = MAX( WRKBL, N+M*ILAENV( 1, 'CUNGQR', ' ', M,
274
$ M, N, -1 ) )
275
WRKBL = MAX( WRKBL, 2*N+2*N*
276
$ ILAENV( 1, 'CGEBRD', ' ', N, N, -1, -1 ) )
277
WRKBL = MAX( WRKBL, 2*N+N*
278
$ ILAENV( 1, 'CUNMBR', 'QLN', N, N, N, -1 ) )
279
WRKBL = MAX( WRKBL, 2*N+N*
280
$ ILAENV( 1, 'CUNMBR', 'PRC', N, N, N, -1 ) )
281
MAXWRK = N*N + WRKBL
282
MINWRK = N*N + 2*N + M
283
END IF
284
ELSE IF( M.GE.MNTHR2 ) THEN
285
*
286
* Path 5 (M much larger than N, but not as much as MNTHR1)
287
*
288
MAXWRK = 2*N + ( M+N )*ILAENV( 1, 'CGEBRD', ' ', M, N,
289
$ -1, -1 )
290
MINWRK = 2*N + M
291
IF( WNTQO ) THEN
292
MAXWRK = MAX( MAXWRK, 2*N+N*
293
$ ILAENV( 1, 'CUNGBR', 'P', N, N, N, -1 ) )
294
MAXWRK = MAX( MAXWRK, 2*N+N*
295
$ ILAENV( 1, 'CUNGBR', 'Q', M, N, N, -1 ) )
296
MAXWRK = MAXWRK + M*N
297
MINWRK = MINWRK + N*N
298
ELSE IF( WNTQS ) THEN
299
MAXWRK = MAX( MAXWRK, 2*N+N*
300
$ ILAENV( 1, 'CUNGBR', 'P', N, N, N, -1 ) )
301
MAXWRK = MAX( MAXWRK, 2*N+N*
302
$ ILAENV( 1, 'CUNGBR', 'Q', M, N, N, -1 ) )
303
ELSE IF( WNTQA ) THEN
304
MAXWRK = MAX( MAXWRK, 2*N+N*
305
$ ILAENV( 1, 'CUNGBR', 'P', N, N, N, -1 ) )
306
MAXWRK = MAX( MAXWRK, 2*N+M*
307
$ ILAENV( 1, 'CUNGBR', 'Q', M, M, N, -1 ) )
308
END IF
309
ELSE
310
*
311
* Path 6 (M at least N, but not much larger)
312
*
313
MAXWRK = 2*N + ( M+N )*ILAENV( 1, 'CGEBRD', ' ', M, N,
314
$ -1, -1 )
315
MINWRK = 2*N + M
316
IF( WNTQO ) THEN
317
MAXWRK = MAX( MAXWRK, 2*N+N*
318
$ ILAENV( 1, 'CUNMBR', 'PRC', N, N, N, -1 ) )
319
MAXWRK = MAX( MAXWRK, 2*N+N*
320
$ ILAENV( 1, 'CUNMBR', 'QLN', M, N, N, -1 ) )
321
MAXWRK = MAXWRK + M*N
322
MINWRK = MINWRK + N*N
323
ELSE IF( WNTQS ) THEN
324
MAXWRK = MAX( MAXWRK, 2*N+N*
325
$ ILAENV( 1, 'CUNMBR', 'PRC', N, N, N, -1 ) )
326
MAXWRK = MAX( MAXWRK, 2*N+N*
327
$ ILAENV( 1, 'CUNMBR', 'QLN', M, N, N, -1 ) )
328
ELSE IF( WNTQA ) THEN
329
MAXWRK = MAX( MAXWRK, 2*N+N*
330
$ ILAENV( 1, 'CUNGBR', 'PRC', N, N, N, -1 ) )
331
MAXWRK = MAX( MAXWRK, 2*N+M*
332
$ ILAENV( 1, 'CUNGBR', 'QLN', M, M, N, -1 ) )
333
END IF
334
END IF
335
ELSE
336
*
337
* There is no complex work space needed for bidiagonal SVD
338
* The real work space needed for bidiagonal SVD is BDSPAC,
339
* BDSPAC = 3*M*M + 4*M
340
*
341
IF( N.GE.MNTHR1 ) THEN
342
IF( WNTQN ) THEN
343
*
344
* Path 1t (N much larger than M, JOBZ='N')
345
*
346
MAXWRK = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1,
347
$ -1 )
348
MAXWRK = MAX( MAXWRK, 2*M+2*M*
349
$ ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) )
350
MINWRK = 3*M
351
ELSE IF( WNTQO ) THEN
352
*
353
* Path 2t (N much larger than M, JOBZ='O')
354
*
355
WRKBL = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 )
356
WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'CUNGLQ', ' ', M,
357
$ N, M, -1 ) )
358
WRKBL = MAX( WRKBL, 2*M+2*M*
359
$ ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) )
360
WRKBL = MAX( WRKBL, 2*M+M*
361
$ ILAENV( 1, 'CUNMBR', 'PRC', M, M, M, -1 ) )
362
WRKBL = MAX( WRKBL, 2*M+M*
363
$ ILAENV( 1, 'CUNMBR', 'QLN', M, M, M, -1 ) )
364
MAXWRK = M*N + M*M + WRKBL
365
MINWRK = 2*M*M + 3*M
366
ELSE IF( WNTQS ) THEN
367
*
368
* Path 3t (N much larger than M, JOBZ='S')
369
*
370
WRKBL = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 )
371
WRKBL = MAX( WRKBL, M+M*ILAENV( 1, 'CUNGLQ', ' ', M,
372
$ N, M, -1 ) )
373
WRKBL = MAX( WRKBL, 2*M+2*M*
374
$ ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) )
375
WRKBL = MAX( WRKBL, 2*M+M*
376
$ ILAENV( 1, 'CUNMBR', 'PRC', M, M, M, -1 ) )
377
WRKBL = MAX( WRKBL, 2*M+M*
378
$ ILAENV( 1, 'CUNMBR', 'QLN', M, M, M, -1 ) )
379
MAXWRK = M*M + WRKBL
380
MINWRK = M*M + 3*M
381
ELSE IF( WNTQA ) THEN
382
*
383
* Path 4t (N much larger than M, JOBZ='A')
384
*
385
WRKBL = M + M*ILAENV( 1, 'CGELQF', ' ', M, N, -1, -1 )
386
WRKBL = MAX( WRKBL, M+N*ILAENV( 1, 'CUNGLQ', ' ', N,
387
$ N, M, -1 ) )
388
WRKBL = MAX( WRKBL, 2*M+2*M*
389
$ ILAENV( 1, 'CGEBRD', ' ', M, M, -1, -1 ) )
390
WRKBL = MAX( WRKBL, 2*M+M*
391
$ ILAENV( 1, 'CUNMBR', 'PRC', M, M, M, -1 ) )
392
WRKBL = MAX( WRKBL, 2*M+M*
393
$ ILAENV( 1, 'CUNMBR', 'QLN', M, M, M, -1 ) )
394
MAXWRK = M*M + WRKBL
395
MINWRK = M*M + 2*M + N
396
END IF
397
ELSE IF( N.GE.MNTHR2 ) THEN
398
*
399
* Path 5t (N much larger than M, but not as much as MNTHR1)
400
*
401
MAXWRK = 2*M + ( M+N )*ILAENV( 1, 'CGEBRD', ' ', M, N,
402
$ -1, -1 )
403
MINWRK = 2*M + N
404
IF( WNTQO ) THEN
405
MAXWRK = MAX( MAXWRK, 2*M+M*
406
$ ILAENV( 1, 'CUNGBR', 'P', M, N, M, -1 ) )
407
MAXWRK = MAX( MAXWRK, 2*M+M*
408
$ ILAENV( 1, 'CUNGBR', 'Q', M, M, N, -1 ) )
409
MAXWRK = MAXWRK + M*N
410
MINWRK = MINWRK + M*M
411
ELSE IF( WNTQS ) THEN
412
MAXWRK = MAX( MAXWRK, 2*M+M*
413
$ ILAENV( 1, 'CUNGBR', 'P', M, N, M, -1 ) )
414
MAXWRK = MAX( MAXWRK, 2*M+M*
415
$ ILAENV( 1, 'CUNGBR', 'Q', M, M, N, -1 ) )
416
ELSE IF( WNTQA ) THEN
417
MAXWRK = MAX( MAXWRK, 2*M+N*
418
$ ILAENV( 1, 'CUNGBR', 'P', N, N, M, -1 ) )
419
MAXWRK = MAX( MAXWRK, 2*M+M*
420
$ ILAENV( 1, 'CUNGBR', 'Q', M, M, N, -1 ) )
421
END IF
422
ELSE
423
*
424
* Path 6t (N greater than M, but not much larger)
425
*
426
MAXWRK = 2*M + ( M+N )*ILAENV( 1, 'CGEBRD', ' ', M, N,
427
$ -1, -1 )
428
MINWRK = 2*M + N
429
IF( WNTQO ) THEN
430
MAXWRK = MAX( MAXWRK, 2*M+M*
431
$ ILAENV( 1, 'CUNMBR', 'PRC', M, N, M, -1 ) )
432
MAXWRK = MAX( MAXWRK, 2*M+M*
433
$ ILAENV( 1, 'CUNMBR', 'QLN', M, M, N, -1 ) )
434
MAXWRK = MAXWRK + M*N
435
MINWRK = MINWRK + M*M
436
ELSE IF( WNTQS ) THEN
437
MAXWRK = MAX( MAXWRK, 2*M+M*
438
$ ILAENV( 1, 'CUNGBR', 'PRC', M, N, M, -1 ) )
439
MAXWRK = MAX( MAXWRK, 2*M+M*
440
$ ILAENV( 1, 'CUNGBR', 'QLN', M, M, N, -1 ) )
441
ELSE IF( WNTQA ) THEN
442
MAXWRK = MAX( MAXWRK, 2*M+N*
443
$ ILAENV( 1, 'CUNGBR', 'PRC', N, N, M, -1 ) )
444
MAXWRK = MAX( MAXWRK, 2*M+M*
445
$ ILAENV( 1, 'CUNGBR', 'QLN', M, M, N, -1 ) )
446
END IF
447
END IF
448
END IF
449
MAXWRK = MAX( MAXWRK, MINWRK )
450
WORK( 1 ) = MAXWRK
451
END IF
452
*
453
IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
454
INFO = -13
455
END IF
456
IF( INFO.NE.0 ) THEN
457
CALL XERBLA( 'CGESDD', -INFO )
458
RETURN
459
ELSE IF( LQUERY ) THEN
460
RETURN
461
END IF
462
*
463
* Quick return if possible
464
*
465
IF( M.EQ.0 .OR. N.EQ.0 ) THEN
466
IF( LWORK.GE.1 )
467
$ WORK( 1 ) = ONE
468
RETURN
469
END IF
470
*
471
* Get machine constants
472
*
473
EPS = SLAMCH( 'P' )
474
SMLNUM = SQRT( SLAMCH( 'S' ) ) / EPS
475
BIGNUM = ONE / SMLNUM
476
*
477
* Scale A if max element outside range [SMLNUM,BIGNUM]
478
*
479
ANRM = CLANGE( 'M', M, N, A, LDA, DUM )
480
ISCL = 0
481
IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
482
ISCL = 1
483
CALL CLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, IERR )
484
ELSE IF( ANRM.GT.BIGNUM ) THEN
485
ISCL = 1
486
CALL CLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, IERR )
487
END IF
488
*
489
IF( M.GE.N ) THEN
490
*
491
* A has at least as many rows as columns. If A has sufficiently
492
* more rows than columns, first reduce using the QR
493
* decomposition (if sufficient workspace available)
494
*
495
IF( M.GE.MNTHR1 ) THEN
496
*
497
IF( WNTQN ) THEN
498
*
499
* Path 1 (M much larger than N, JOBZ='N')
500
* No singular vectors to be computed
501
*
502
ITAU = 1
503
NWORK = ITAU + N
504
*
505
* Compute A=Q*R
506
* (CWorkspace: need 2*N, prefer N+N*NB)
507
* (RWorkspace: need 0)
508
*
509
CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
510
$ LWORK-NWORK+1, IERR )
511
*
512
* Zero out below R
513
*
514
CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ),
515
$ LDA )
516
IE = 1
517
ITAUQ = 1
518
ITAUP = ITAUQ + N
519
NWORK = ITAUP + N
520
*
521
* Bidiagonalize R in A
522
* (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
523
* (RWorkspace: need N)
524
*
525
CALL CGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
526
$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
527
$ IERR )
528
NRWORK = IE + N
529
*
530
* Perform bidiagonal SVD, compute singular values only
531
* (CWorkspace: 0)
532
* (RWorkspace: need BDSPAC)
533
*
534
CALL SBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1,
535
$ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
536
*
537
ELSE IF( WNTQO ) THEN
538
*
539
* Path 2 (M much larger than N, JOBZ='O')
540
* N left singular vectors to be overwritten on A and
541
* N right singular vectors to be computed in VT
542
*
543
IU = 1
544
*
545
* WORK(IU) is N by N
546
*
547
LDWRKU = N
548
IR = IU + LDWRKU*N
549
IF( LWORK.GE.M*N+N*N+3*N ) THEN
550
*
551
* WORK(IR) is M by N
552
*
553
LDWRKR = M
554
ELSE
555
LDWRKR = ( LWORK-N*N-3*N ) / N
556
END IF
557
ITAU = IR + LDWRKR*N
558
NWORK = ITAU + N
559
*
560
* Compute A=Q*R
561
* (CWorkspace: need N*N+2*N, prefer M*N+N+N*NB)
562
* (RWorkspace: 0)
563
*
564
CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
565
$ LWORK-NWORK+1, IERR )
566
*
567
* Copy R to WORK( IR ), zeroing out below it
568
*
569
CALL CLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
570
CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, WORK( IR+1 ),
571
$ LDWRKR )
572
*
573
* Generate Q in A
574
* (CWorkspace: need 2*N, prefer N+N*NB)
575
* (RWorkspace: 0)
576
*
577
CALL CUNGQR( M, N, N, A, LDA, WORK( ITAU ),
578
$ WORK( NWORK ), LWORK-NWORK+1, IERR )
579
IE = 1
580
ITAUQ = ITAU
581
ITAUP = ITAUQ + N
582
NWORK = ITAUP + N
583
*
584
* Bidiagonalize R in WORK(IR)
585
* (CWorkspace: need N*N+3*N, prefer M*N+2*N+2*N*NB)
586
* (RWorkspace: need N)
587
*
588
CALL CGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ),
589
$ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
590
$ LWORK-NWORK+1, IERR )
591
*
592
* Perform bidiagonal SVD, computing left singular vectors
593
* of R in WORK(IRU) and computing right singular vectors
594
* of R in WORK(IRVT)
595
* (CWorkspace: need 0)
596
* (RWorkspace: need BDSPAC)
597
*
598
IRU = IE + N
599
IRVT = IRU + N*N
600
NRWORK = IRVT + N*N
601
CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
602
$ N, RWORK( IRVT ), N, DUM, IDUM,
603
$ RWORK( NRWORK ), IWORK, INFO )
604
*
605
* Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
606
* Overwrite WORK(IU) by the left singular vectors of R
607
* (CWorkspace: need 2*N*N+3*N, prefer M*N+N*N+2*N+N*NB)
608
* (RWorkspace: 0)
609
*
610
CALL CLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ),
611
$ LDWRKU )
612
CALL CUNMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
613
$ WORK( ITAUQ ), WORK( IU ), LDWRKU,
614
$ WORK( NWORK ), LWORK-NWORK+1, IERR )
615
*
616
* Copy real matrix RWORK(IRVT) to complex matrix VT
617
* Overwrite VT by the right singular vectors of R
618
* (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB)
619
* (RWorkspace: 0)
620
*
621
CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
622
CALL CUNMBR( 'P', 'R', 'C', N, N, N, WORK( IR ), LDWRKR,
623
$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
624
$ LWORK-NWORK+1, IERR )
625
*
626
* Multiply Q in A by left singular vectors of R in
627
* WORK(IU), storing result in WORK(IR) and copying to A
628
* (CWorkspace: need 2*N*N, prefer N*N+M*N)
629
* (RWorkspace: 0)
630
*
631
DO 10 I = 1, M, LDWRKR
632
CHUNK = MIN( M-I+1, LDWRKR )
633
CALL CGEMM( 'N', 'N', CHUNK, N, N, CONE, A( I, 1 ),
634
$ LDA, WORK( IU ), LDWRKU, CZERO,
635
$ WORK( IR ), LDWRKR )
636
CALL CLACPY( 'F', CHUNK, N, WORK( IR ), LDWRKR,
637
$ A( I, 1 ), LDA )
638
10 CONTINUE
639
*
640
ELSE IF( WNTQS ) THEN
641
*
642
* Path 3 (M much larger than N, JOBZ='S')
643
* N left singular vectors to be computed in U and
644
* N right singular vectors to be computed in VT
645
*
646
IR = 1
647
*
648
* WORK(IR) is N by N
649
*
650
LDWRKR = N
651
ITAU = IR + LDWRKR*N
652
NWORK = ITAU + N
653
*
654
* Compute A=Q*R
655
* (CWorkspace: need N*N+2*N, prefer N*N+N+N*NB)
656
* (RWorkspace: 0)
657
*
658
CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
659
$ LWORK-NWORK+1, IERR )
660
*
661
* Copy R to WORK(IR), zeroing out below it
662
*
663
CALL CLACPY( 'U', N, N, A, LDA, WORK( IR ), LDWRKR )
664
CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, WORK( IR+1 ),
665
$ LDWRKR )
666
*
667
* Generate Q in A
668
* (CWorkspace: need 2*N, prefer N+N*NB)
669
* (RWorkspace: 0)
670
*
671
CALL CUNGQR( M, N, N, A, LDA, WORK( ITAU ),
672
$ WORK( NWORK ), LWORK-NWORK+1, IERR )
673
IE = 1
674
ITAUQ = ITAU
675
ITAUP = ITAUQ + N
676
NWORK = ITAUP + N
677
*
678
* Bidiagonalize R in WORK(IR)
679
* (CWorkspace: need N*N+3*N, prefer N*N+2*N+2*N*NB)
680
* (RWorkspace: need N)
681
*
682
CALL CGEBRD( N, N, WORK( IR ), LDWRKR, S, RWORK( IE ),
683
$ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
684
$ LWORK-NWORK+1, IERR )
685
*
686
* Perform bidiagonal SVD, computing left singular vectors
687
* of bidiagonal matrix in RWORK(IRU) and computing right
688
* singular vectors of bidiagonal matrix in RWORK(IRVT)
689
* (CWorkspace: need 0)
690
* (RWorkspace: need BDSPAC)
691
*
692
IRU = IE + N
693
IRVT = IRU + N*N
694
NRWORK = IRVT + N*N
695
CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
696
$ N, RWORK( IRVT ), N, DUM, IDUM,
697
$ RWORK( NRWORK ), IWORK, INFO )
698
*
699
* Copy real matrix RWORK(IRU) to complex matrix U
700
* Overwrite U by left singular vectors of R
701
* (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
702
* (RWorkspace: 0)
703
*
704
CALL CLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )
705
CALL CUNMBR( 'Q', 'L', 'N', N, N, N, WORK( IR ), LDWRKR,
706
$ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
707
$ LWORK-NWORK+1, IERR )
708
*
709
* Copy real matrix RWORK(IRVT) to complex matrix VT
710
* Overwrite VT by right singular vectors of R
711
* (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
712
* (RWorkspace: 0)
713
*
714
CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
715
CALL CUNMBR( 'P', 'R', 'C', N, N, N, WORK( IR ), LDWRKR,
716
$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
717
$ LWORK-NWORK+1, IERR )
718
*
719
* Multiply Q in A by left singular vectors of R in
720
* WORK(IR), storing result in U
721
* (CWorkspace: need N*N)
722
* (RWorkspace: 0)
723
*
724
CALL CLACPY( 'F', N, N, U, LDU, WORK( IR ), LDWRKR )
725
CALL CGEMM( 'N', 'N', M, N, N, CONE, A, LDA, WORK( IR ),
726
$ LDWRKR, CZERO, U, LDU )
727
*
728
ELSE IF( WNTQA ) THEN
729
*
730
* Path 4 (M much larger than N, JOBZ='A')
731
* M left singular vectors to be computed in U and
732
* N right singular vectors to be computed in VT
733
*
734
IU = 1
735
*
736
* WORK(IU) is N by N
737
*
738
LDWRKU = N
739
ITAU = IU + LDWRKU*N
740
NWORK = ITAU + N
741
*
742
* Compute A=Q*R, copying result to U
743
* (CWorkspace: need 2*N, prefer N+N*NB)
744
* (RWorkspace: 0)
745
*
746
CALL CGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
747
$ LWORK-NWORK+1, IERR )
748
CALL CLACPY( 'L', M, N, A, LDA, U, LDU )
749
*
750
* Generate Q in U
751
* (CWorkspace: need N+M, prefer N+M*NB)
752
* (RWorkspace: 0)
753
*
754
CALL CUNGQR( M, M, N, U, LDU, WORK( ITAU ),
755
$ WORK( NWORK ), LWORK-NWORK+1, IERR )
756
*
757
* Produce R in A, zeroing out below it
758
*
759
CALL CLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ),
760
$ LDA )
761
IE = 1
762
ITAUQ = ITAU
763
ITAUP = ITAUQ + N
764
NWORK = ITAUP + N
765
*
766
* Bidiagonalize R in A
767
* (CWorkspace: need 3*N, prefer 2*N+2*N*NB)
768
* (RWorkspace: need N)
769
*
770
CALL CGEBRD( N, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
771
$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
772
$ IERR )
773
IRU = IE + N
774
IRVT = IRU + N*N
775
NRWORK = IRVT + N*N
776
*
777
* Perform bidiagonal SVD, computing left singular vectors
778
* of bidiagonal matrix in RWORK(IRU) and computing right
779
* singular vectors of bidiagonal matrix in RWORK(IRVT)
780
* (CWorkspace: need 0)
781
* (RWorkspace: need BDSPAC)
782
*
783
CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
784
$ N, RWORK( IRVT ), N, DUM, IDUM,
785
$ RWORK( NRWORK ), IWORK, INFO )
786
*
787
* Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
788
* Overwrite WORK(IU) by left singular vectors of R
789
* (CWorkspace: need N*N+3*N, prefer N*N+2*N+N*NB)
790
* (RWorkspace: 0)
791
*
792
CALL CLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ),
793
$ LDWRKU )
794
CALL CUNMBR( 'Q', 'L', 'N', N, N, N, A, LDA,
795
$ WORK( ITAUQ ), WORK( IU ), LDWRKU,
796
$ WORK( NWORK ), LWORK-NWORK+1, IERR )
797
*
798
* Copy real matrix RWORK(IRVT) to complex matrix VT
799
* Overwrite VT by right singular vectors of R
800
* (CWorkspace: need 3*N, prefer 2*N+N*NB)
801
* (RWorkspace: 0)
802
*
803
CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
804
CALL CUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,
805
$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
806
$ LWORK-NWORK+1, IERR )
807
*
808
* Multiply Q in U by left singular vectors of R in
809
* WORK(IU), storing result in A
810
* (CWorkspace: need N*N)
811
* (RWorkspace: 0)
812
*
813
CALL CGEMM( 'N', 'N', M, N, N, CONE, U, LDU, WORK( IU ),
814
$ LDWRKU, CZERO, A, LDA )
815
*
816
* Copy left singular vectors of A from A to U
817
*
818
CALL CLACPY( 'F', M, N, A, LDA, U, LDU )
819
*
820
END IF
821
*
822
ELSE IF( M.GE.MNTHR2 ) THEN
823
*
824
* MNTHR2 <= M < MNTHR1
825
*
826
* Path 5 (M much larger than N, but not as much as MNTHR1)
827
* Reduce to bidiagonal form without QR decomposition, use
828
* CUNGBR and matrix multiplication to compute singular vectors
829
*
830
IE = 1
831
NRWORK = IE + N
832
ITAUQ = 1
833
ITAUP = ITAUQ + N
834
NWORK = ITAUP + N
835
*
836
* Bidiagonalize A
837
* (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB)
838
* (RWorkspace: need N)
839
*
840
CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
841
$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
842
$ IERR )
843
IF( WNTQN ) THEN
844
*
845
* Compute singular values only
846
* (Cworkspace: 0)
847
* (Rworkspace: need BDSPAC)
848
*
849
CALL SBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1,
850
$ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
851
ELSE IF( WNTQO ) THEN
852
IU = NWORK
853
IRU = NRWORK
854
IRVT = IRU + N*N
855
NRWORK = IRVT + N*N
856
*
857
* Copy A to VT, generate P**H
858
* (Cworkspace: need 2*N, prefer N+N*NB)
859
* (Rworkspace: 0)
860
*
861
CALL CLACPY( 'U', N, N, A, LDA, VT, LDVT )
862
CALL CUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
863
$ WORK( NWORK ), LWORK-NWORK+1, IERR )
864
*
865
* Generate Q in A
866
* (CWorkspace: need 2*N, prefer N+N*NB)
867
* (RWorkspace: 0)
868
*
869
CALL CUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
870
$ WORK( NWORK ), LWORK-NWORK+1, IERR )
871
*
872
IF( LWORK.GE.M*N+3*N ) THEN
873
*
874
* WORK( IU ) is M by N
875
*
876
LDWRKU = M
877
ELSE
878
*
879
* WORK(IU) is LDWRKU by N
880
*
881
LDWRKU = ( LWORK-3*N ) / N
882
END IF
883
NWORK = IU + LDWRKU*N
884
*
885
* Perform bidiagonal SVD, computing left singular vectors
886
* of bidiagonal matrix in RWORK(IRU) and computing right
887
* singular vectors of bidiagonal matrix in RWORK(IRVT)
888
* (CWorkspace: need 0)
889
* (RWorkspace: need BDSPAC)
890
*
891
CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
892
$ N, RWORK( IRVT ), N, DUM, IDUM,
893
$ RWORK( NRWORK ), IWORK, INFO )
894
*
895
* Multiply real matrix RWORK(IRVT) by P**H in VT,
896
* storing the result in WORK(IU), copying to VT
897
* (Cworkspace: need 0)
898
* (Rworkspace: need 3*N*N)
899
*
900
CALL CLARCM( N, N, RWORK( IRVT ), N, VT, LDVT,
901
$ WORK( IU ), LDWRKU, RWORK( NRWORK ) )
902
CALL CLACPY( 'F', N, N, WORK( IU ), LDWRKU, VT, LDVT )
903
*
904
* Multiply Q in A by real matrix RWORK(IRU), storing the
905
* result in WORK(IU), copying to A
906
* (CWorkspace: need N*N, prefer M*N)
907
* (Rworkspace: need 3*N*N, prefer N*N+2*M*N)
908
*
909
NRWORK = IRVT
910
DO 20 I = 1, M, LDWRKU
911
CHUNK = MIN( M-I+1, LDWRKU )
912
CALL CLACRM( CHUNK, N, A( I, 1 ), LDA, RWORK( IRU ),
913
$ N, WORK( IU ), LDWRKU, RWORK( NRWORK ) )
914
CALL CLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU,
915
$ A( I, 1 ), LDA )
916
20 CONTINUE
917
*
918
ELSE IF( WNTQS ) THEN
919
*
920
* Copy A to VT, generate P**H
921
* (Cworkspace: need 2*N, prefer N+N*NB)
922
* (Rworkspace: 0)
923
*
924
CALL CLACPY( 'U', N, N, A, LDA, VT, LDVT )
925
CALL CUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
926
$ WORK( NWORK ), LWORK-NWORK+1, IERR )
927
*
928
* Copy A to U, generate Q
929
* (Cworkspace: need 2*N, prefer N+N*NB)
930
* (Rworkspace: 0)
931
*
932
CALL CLACPY( 'L', M, N, A, LDA, U, LDU )
933
CALL CUNGBR( 'Q', M, N, N, U, LDU, WORK( ITAUQ ),
934
$ WORK( NWORK ), LWORK-NWORK+1, IERR )
935
*
936
* Perform bidiagonal SVD, computing left singular vectors
937
* of bidiagonal matrix in RWORK(IRU) and computing right
938
* singular vectors of bidiagonal matrix in RWORK(IRVT)
939
* (CWorkspace: need 0)
940
* (RWorkspace: need BDSPAC)
941
*
942
IRU = NRWORK
943
IRVT = IRU + N*N
944
NRWORK = IRVT + N*N
945
CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
946
$ N, RWORK( IRVT ), N, DUM, IDUM,
947
$ RWORK( NRWORK ), IWORK, INFO )
948
*
949
* Multiply real matrix RWORK(IRVT) by P**H in VT,
950
* storing the result in A, copying to VT
951
* (Cworkspace: need 0)
952
* (Rworkspace: need 3*N*N)
953
*
954
CALL CLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, A, LDA,
955
$ RWORK( NRWORK ) )
956
CALL CLACPY( 'F', N, N, A, LDA, VT, LDVT )
957
*
958
* Multiply Q in U by real matrix RWORK(IRU), storing the
959
* result in A, copying to U
960
* (CWorkspace: need 0)
961
* (Rworkspace: need N*N+2*M*N)
962
*
963
NRWORK = IRVT
964
CALL CLACRM( M, N, U, LDU, RWORK( IRU ), N, A, LDA,
965
$ RWORK( NRWORK ) )
966
CALL CLACPY( 'F', M, N, A, LDA, U, LDU )
967
ELSE
968
*
969
* Copy A to VT, generate P**H
970
* (Cworkspace: need 2*N, prefer N+N*NB)
971
* (Rworkspace: 0)
972
*
973
CALL CLACPY( 'U', N, N, A, LDA, VT, LDVT )
974
CALL CUNGBR( 'P', N, N, N, VT, LDVT, WORK( ITAUP ),
975
$ WORK( NWORK ), LWORK-NWORK+1, IERR )
976
*
977
* Copy A to U, generate Q
978
* (Cworkspace: need 2*N, prefer N+N*NB)
979
* (Rworkspace: 0)
980
*
981
CALL CLACPY( 'L', M, N, A, LDA, U, LDU )
982
CALL CUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
983
$ WORK( NWORK ), LWORK-NWORK+1, IERR )
984
*
985
* Perform bidiagonal SVD, computing left singular vectors
986
* of bidiagonal matrix in RWORK(IRU) and computing right
987
* singular vectors of bidiagonal matrix in RWORK(IRVT)
988
* (CWorkspace: need 0)
989
* (RWorkspace: need BDSPAC)
990
*
991
IRU = NRWORK
992
IRVT = IRU + N*N
993
NRWORK = IRVT + N*N
994
CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
995
$ N, RWORK( IRVT ), N, DUM, IDUM,
996
$ RWORK( NRWORK ), IWORK, INFO )
997
*
998
* Multiply real matrix RWORK(IRVT) by P**H in VT,
999
* storing the result in A, copying to VT
1000
* (Cworkspace: need 0)
1001
* (Rworkspace: need 3*N*N)
1002
*
1003
CALL CLARCM( N, N, RWORK( IRVT ), N, VT, LDVT, A, LDA,
1004
$ RWORK( NRWORK ) )
1005
CALL CLACPY( 'F', N, N, A, LDA, VT, LDVT )
1006
*
1007
* Multiply Q in U by real matrix RWORK(IRU), storing the
1008
* result in A, copying to U
1009
* (CWorkspace: 0)
1010
* (Rworkspace: need 3*N*N)
1011
*
1012
NRWORK = IRVT
1013
CALL CLACRM( M, N, U, LDU, RWORK( IRU ), N, A, LDA,
1014
$ RWORK( NRWORK ) )
1015
CALL CLACPY( 'F', M, N, A, LDA, U, LDU )
1016
END IF
1017
*
1018
ELSE
1019
*
1020
* M .LT. MNTHR2
1021
*
1022
* Path 6 (M at least N, but not much larger)
1023
* Reduce to bidiagonal form without QR decomposition
1024
* Use CUNMBR to compute singular vectors
1025
*
1026
IE = 1
1027
NRWORK = IE + N
1028
ITAUQ = 1
1029
ITAUP = ITAUQ + N
1030
NWORK = ITAUP + N
1031
*
1032
* Bidiagonalize A
1033
* (CWorkspace: need 2*N+M, prefer 2*N+(M+N)*NB)
1034
* (RWorkspace: need N)
1035
*
1036
CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
1037
$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1038
$ IERR )
1039
IF( WNTQN ) THEN
1040
*
1041
* Compute singular values only
1042
* (Cworkspace: 0)
1043
* (Rworkspace: need BDSPAC)
1044
*
1045
CALL SBDSDC( 'U', 'N', N, S, RWORK( IE ), DUM, 1, DUM, 1,
1046
$ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
1047
ELSE IF( WNTQO ) THEN
1048
IU = NWORK
1049
IRU = NRWORK
1050
IRVT = IRU + N*N
1051
NRWORK = IRVT + N*N
1052
IF( LWORK.GE.M*N+3*N ) THEN
1053
*
1054
* WORK( IU ) is M by N
1055
*
1056
LDWRKU = M
1057
ELSE
1058
*
1059
* WORK( IU ) is LDWRKU by N
1060
*
1061
LDWRKU = ( LWORK-3*N ) / N
1062
END IF
1063
NWORK = IU + LDWRKU*N
1064
*
1065
* Perform bidiagonal SVD, computing left singular vectors
1066
* of bidiagonal matrix in RWORK(IRU) and computing right
1067
* singular vectors of bidiagonal matrix in RWORK(IRVT)
1068
* (CWorkspace: need 0)
1069
* (RWorkspace: need BDSPAC)
1070
*
1071
CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
1072
$ N, RWORK( IRVT ), N, DUM, IDUM,
1073
$ RWORK( NRWORK ), IWORK, INFO )
1074
*
1075
* Copy real matrix RWORK(IRVT) to complex matrix VT
1076
* Overwrite VT by right singular vectors of A
1077
* (Cworkspace: need 2*N, prefer N+N*NB)
1078
* (Rworkspace: need 0)
1079
*
1080
CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
1081
CALL CUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,
1082
$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1083
$ LWORK-NWORK+1, IERR )
1084
*
1085
IF( LWORK.GE.M*N+3*N ) THEN
1086
*
1087
* Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
1088
* Overwrite WORK(IU) by left singular vectors of A, copying
1089
* to A
1090
* (Cworkspace: need M*N+2*N, prefer M*N+N+N*NB)
1091
* (Rworkspace: need 0)
1092
*
1093
CALL CLASET( 'F', M, N, CZERO, CZERO, WORK( IU ),
1094
$ LDWRKU )
1095
CALL CLACP2( 'F', N, N, RWORK( IRU ), N, WORK( IU ),
1096
$ LDWRKU )
1097
CALL CUNMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
1098
$ WORK( ITAUQ ), WORK( IU ), LDWRKU,
1099
$ WORK( NWORK ), LWORK-NWORK+1, IERR )
1100
CALL CLACPY( 'F', M, N, WORK( IU ), LDWRKU, A, LDA )
1101
ELSE
1102
*
1103
* Generate Q in A
1104
* (Cworkspace: need 2*N, prefer N+N*NB)
1105
* (Rworkspace: need 0)
1106
*
1107
CALL CUNGBR( 'Q', M, N, N, A, LDA, WORK( ITAUQ ),
1108
$ WORK( NWORK ), LWORK-NWORK+1, IERR )
1109
*
1110
* Multiply Q in A by real matrix RWORK(IRU), storing the
1111
* result in WORK(IU), copying to A
1112
* (CWorkspace: need N*N, prefer M*N)
1113
* (Rworkspace: need 3*N*N, prefer N*N+2*M*N)
1114
*
1115
NRWORK = IRVT
1116
DO 30 I = 1, M, LDWRKU
1117
CHUNK = MIN( M-I+1, LDWRKU )
1118
CALL CLACRM( CHUNK, N, A( I, 1 ), LDA,
1119
$ RWORK( IRU ), N, WORK( IU ), LDWRKU,
1120
$ RWORK( NRWORK ) )
1121
CALL CLACPY( 'F', CHUNK, N, WORK( IU ), LDWRKU,
1122
$ A( I, 1 ), LDA )
1123
30 CONTINUE
1124
END IF
1125
*
1126
ELSE IF( WNTQS ) THEN
1127
*
1128
* Perform bidiagonal SVD, computing left singular vectors
1129
* of bidiagonal matrix in RWORK(IRU) and computing right
1130
* singular vectors of bidiagonal matrix in RWORK(IRVT)
1131
* (CWorkspace: need 0)
1132
* (RWorkspace: need BDSPAC)
1133
*
1134
IRU = NRWORK
1135
IRVT = IRU + N*N
1136
NRWORK = IRVT + N*N
1137
CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
1138
$ N, RWORK( IRVT ), N, DUM, IDUM,
1139
$ RWORK( NRWORK ), IWORK, INFO )
1140
*
1141
* Copy real matrix RWORK(IRU) to complex matrix U
1142
* Overwrite U by left singular vectors of A
1143
* (CWorkspace: need 3*N, prefer 2*N+N*NB)
1144
* (RWorkspace: 0)
1145
*
1146
CALL CLASET( 'F', M, N, CZERO, CZERO, U, LDU )
1147
CALL CLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )
1148
CALL CUNMBR( 'Q', 'L', 'N', M, N, N, A, LDA,
1149
$ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1150
$ LWORK-NWORK+1, IERR )
1151
*
1152
* Copy real matrix RWORK(IRVT) to complex matrix VT
1153
* Overwrite VT by right singular vectors of A
1154
* (CWorkspace: need 3*N, prefer 2*N+N*NB)
1155
* (RWorkspace: 0)
1156
*
1157
CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
1158
CALL CUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,
1159
$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1160
$ LWORK-NWORK+1, IERR )
1161
ELSE
1162
*
1163
* Perform bidiagonal SVD, computing left singular vectors
1164
* of bidiagonal matrix in RWORK(IRU) and computing right
1165
* singular vectors of bidiagonal matrix in RWORK(IRVT)
1166
* (CWorkspace: need 0)
1167
* (RWorkspace: need BDSPAC)
1168
*
1169
IRU = NRWORK
1170
IRVT = IRU + N*N
1171
NRWORK = IRVT + N*N
1172
CALL SBDSDC( 'U', 'I', N, S, RWORK( IE ), RWORK( IRU ),
1173
$ N, RWORK( IRVT ), N, DUM, IDUM,
1174
$ RWORK( NRWORK ), IWORK, INFO )
1175
*
1176
* Set the right corner of U to identity matrix
1177
*
1178
CALL CLASET( 'F', M, M, CZERO, CZERO, U, LDU )
1179
CALL CLASET( 'F', M-N, M-N, CZERO, CONE, U( N+1, N+1 ),
1180
$ LDU )
1181
*
1182
* Copy real matrix RWORK(IRU) to complex matrix U
1183
* Overwrite U by left singular vectors of A
1184
* (CWorkspace: need 2*N+M, prefer 2*N+M*NB)
1185
* (RWorkspace: 0)
1186
*
1187
CALL CLACP2( 'F', N, N, RWORK( IRU ), N, U, LDU )
1188
CALL CUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1189
$ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1190
$ LWORK-NWORK+1, IERR )
1191
*
1192
* Copy real matrix RWORK(IRVT) to complex matrix VT
1193
* Overwrite VT by right singular vectors of A
1194
* (CWorkspace: need 3*N, prefer 2*N+N*NB)
1195
* (RWorkspace: 0)
1196
*
1197
CALL CLACP2( 'F', N, N, RWORK( IRVT ), N, VT, LDVT )
1198
CALL CUNMBR( 'P', 'R', 'C', N, N, N, A, LDA,
1199
$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1200
$ LWORK-NWORK+1, IERR )
1201
END IF
1202
*
1203
END IF
1204
*
1205
ELSE
1206
*
1207
* A has more columns than rows. If A has sufficiently more
1208
* columns than rows, first reduce using the LQ decomposition
1209
* (if sufficient workspace available)
1210
*
1211
IF( N.GE.MNTHR1 ) THEN
1212
*
1213
IF( WNTQN ) THEN
1214
*
1215
* Path 1t (N much larger than M, JOBZ='N')
1216
* No singular vectors to be computed
1217
*
1218
ITAU = 1
1219
NWORK = ITAU + M
1220
*
1221
* Compute A=L*Q
1222
* (CWorkspace: need 2*M, prefer M+M*NB)
1223
* (RWorkspace: 0)
1224
*
1225
CALL CGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1226
$ LWORK-NWORK+1, IERR )
1227
*
1228
* Zero out above L
1229
*
1230
CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, A( 1, 2 ),
1231
$ LDA )
1232
IE = 1
1233
ITAUQ = 1
1234
ITAUP = ITAUQ + M
1235
NWORK = ITAUP + M
1236
*
1237
* Bidiagonalize L in A
1238
* (CWorkspace: need 3*M, prefer 2*M+2*M*NB)
1239
* (RWorkspace: need M)
1240
*
1241
CALL CGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
1242
$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1243
$ IERR )
1244
NRWORK = IE + M
1245
*
1246
* Perform bidiagonal SVD, compute singular values only
1247
* (CWorkspace: 0)
1248
* (RWorkspace: need BDSPAC)
1249
*
1250
CALL SBDSDC( 'U', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1,
1251
$ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
1252
*
1253
ELSE IF( WNTQO ) THEN
1254
*
1255
* Path 2t (N much larger than M, JOBZ='O')
1256
* M right singular vectors to be overwritten on A and
1257
* M left singular vectors to be computed in U
1258
*
1259
IVT = 1
1260
LDWKVT = M
1261
*
1262
* WORK(IVT) is M by M
1263
*
1264
IL = IVT + LDWKVT*M
1265
IF( LWORK.GE.M*N+M*M+3*M ) THEN
1266
*
1267
* WORK(IL) M by N
1268
*
1269
LDWRKL = M
1270
CHUNK = N
1271
ELSE
1272
*
1273
* WORK(IL) is M by CHUNK
1274
*
1275
LDWRKL = M
1276
CHUNK = ( LWORK-M*M-3*M ) / M
1277
END IF
1278
ITAU = IL + LDWRKL*CHUNK
1279
NWORK = ITAU + M
1280
*
1281
* Compute A=L*Q
1282
* (CWorkspace: need 2*M, prefer M+M*NB)
1283
* (RWorkspace: 0)
1284
*
1285
CALL CGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1286
$ LWORK-NWORK+1, IERR )
1287
*
1288
* Copy L to WORK(IL), zeroing about above it
1289
*
1290
CALL CLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
1291
CALL CLASET( 'U', M-1, M-1, CZERO, CZERO,
1292
$ WORK( IL+LDWRKL ), LDWRKL )
1293
*
1294
* Generate Q in A
1295
* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
1296
* (RWorkspace: 0)
1297
*
1298
CALL CUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
1299
$ WORK( NWORK ), LWORK-NWORK+1, IERR )
1300
IE = 1
1301
ITAUQ = ITAU
1302
ITAUP = ITAUQ + M
1303
NWORK = ITAUP + M
1304
*
1305
* Bidiagonalize L in WORK(IL)
1306
* (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
1307
* (RWorkspace: need M)
1308
*
1309
CALL CGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ),
1310
$ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
1311
$ LWORK-NWORK+1, IERR )
1312
*
1313
* Perform bidiagonal SVD, computing left singular vectors
1314
* of bidiagonal matrix in RWORK(IRU) and computing right
1315
* singular vectors of bidiagonal matrix in RWORK(IRVT)
1316
* (CWorkspace: need 0)
1317
* (RWorkspace: need BDSPAC)
1318
*
1319
IRU = IE + M
1320
IRVT = IRU + M*M
1321
NRWORK = IRVT + M*M
1322
CALL SBDSDC( 'U', 'I', M, S, RWORK( IE ), RWORK( IRU ),
1323
$ M, RWORK( IRVT ), M, DUM, IDUM,
1324
$ RWORK( NRWORK ), IWORK, INFO )
1325
*
1326
* Copy real matrix RWORK(IRU) to complex matrix WORK(IU)
1327
* Overwrite WORK(IU) by the left singular vectors of L
1328
* (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB)
1329
* (RWorkspace: 0)
1330
*
1331
CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
1332
CALL CUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
1333
$ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1334
$ LWORK-NWORK+1, IERR )
1335
*
1336
* Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
1337
* Overwrite WORK(IVT) by the right singular vectors of L
1338
* (CWorkspace: need N*N+3*N, prefer M*N+2*N+N*NB)
1339
* (RWorkspace: 0)
1340
*
1341
CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ),
1342
$ LDWKVT )
1343
CALL CUNMBR( 'P', 'R', 'C', M, M, M, WORK( IL ), LDWRKL,
1344
$ WORK( ITAUP ), WORK( IVT ), LDWKVT,
1345
$ WORK( NWORK ), LWORK-NWORK+1, IERR )
1346
*
1347
* Multiply right singular vectors of L in WORK(IL) by Q
1348
* in A, storing result in WORK(IL) and copying to A
1349
* (CWorkspace: need 2*M*M, prefer M*M+M*N))
1350
* (RWorkspace: 0)
1351
*
1352
DO 40 I = 1, N, CHUNK
1353
BLK = MIN( N-I+1, CHUNK )
1354
CALL CGEMM( 'N', 'N', M, BLK, M, CONE, WORK( IVT ), M,
1355
$ A( 1, I ), LDA, CZERO, WORK( IL ),
1356
$ LDWRKL )
1357
CALL CLACPY( 'F', M, BLK, WORK( IL ), LDWRKL,
1358
$ A( 1, I ), LDA )
1359
40 CONTINUE
1360
*
1361
ELSE IF( WNTQS ) THEN
1362
*
1363
* Path 3t (N much larger than M, JOBZ='S')
1364
* M right singular vectors to be computed in VT and
1365
* M left singular vectors to be computed in U
1366
*
1367
IL = 1
1368
*
1369
* WORK(IL) is M by M
1370
*
1371
LDWRKL = M
1372
ITAU = IL + LDWRKL*M
1373
NWORK = ITAU + M
1374
*
1375
* Compute A=L*Q
1376
* (CWorkspace: need 2*M, prefer M+M*NB)
1377
* (RWorkspace: 0)
1378
*
1379
CALL CGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1380
$ LWORK-NWORK+1, IERR )
1381
*
1382
* Copy L to WORK(IL), zeroing out above it
1383
*
1384
CALL CLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWRKL )
1385
CALL CLASET( 'U', M-1, M-1, CZERO, CZERO,
1386
$ WORK( IL+LDWRKL ), LDWRKL )
1387
*
1388
* Generate Q in A
1389
* (CWorkspace: need M*M+2*M, prefer M*M+M+M*NB)
1390
* (RWorkspace: 0)
1391
*
1392
CALL CUNGLQ( M, N, M, A, LDA, WORK( ITAU ),
1393
$ WORK( NWORK ), LWORK-NWORK+1, IERR )
1394
IE = 1
1395
ITAUQ = ITAU
1396
ITAUP = ITAUQ + M
1397
NWORK = ITAUP + M
1398
*
1399
* Bidiagonalize L in WORK(IL)
1400
* (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
1401
* (RWorkspace: need M)
1402
*
1403
CALL CGEBRD( M, M, WORK( IL ), LDWRKL, S, RWORK( IE ),
1404
$ WORK( ITAUQ ), WORK( ITAUP ), WORK( NWORK ),
1405
$ LWORK-NWORK+1, IERR )
1406
*
1407
* Perform bidiagonal SVD, computing left singular vectors
1408
* of bidiagonal matrix in RWORK(IRU) and computing right
1409
* singular vectors of bidiagonal matrix in RWORK(IRVT)
1410
* (CWorkspace: need 0)
1411
* (RWorkspace: need BDSPAC)
1412
*
1413
IRU = IE + M
1414
IRVT = IRU + M*M
1415
NRWORK = IRVT + M*M
1416
CALL SBDSDC( 'U', 'I', M, S, RWORK( IE ), RWORK( IRU ),
1417
$ M, RWORK( IRVT ), M, DUM, IDUM,
1418
$ RWORK( NRWORK ), IWORK, INFO )
1419
*
1420
* Copy real matrix RWORK(IRU) to complex matrix U
1421
* Overwrite U by left singular vectors of L
1422
* (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
1423
* (RWorkspace: 0)
1424
*
1425
CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
1426
CALL CUNMBR( 'Q', 'L', 'N', M, M, M, WORK( IL ), LDWRKL,
1427
$ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1428
$ LWORK-NWORK+1, IERR )
1429
*
1430
* Copy real matrix RWORK(IRVT) to complex matrix VT
1431
* Overwrite VT by left singular vectors of L
1432
* (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
1433
* (RWorkspace: 0)
1434
*
1435
CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT )
1436
CALL CUNMBR( 'P', 'R', 'C', M, M, M, WORK( IL ), LDWRKL,
1437
$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1438
$ LWORK-NWORK+1, IERR )
1439
*
1440
* Copy VT to WORK(IL), multiply right singular vectors of L
1441
* in WORK(IL) by Q in A, storing result in VT
1442
* (CWorkspace: need M*M)
1443
* (RWorkspace: 0)
1444
*
1445
CALL CLACPY( 'F', M, M, VT, LDVT, WORK( IL ), LDWRKL )
1446
CALL CGEMM( 'N', 'N', M, N, M, CONE, WORK( IL ), LDWRKL,
1447
$ A, LDA, CZERO, VT, LDVT )
1448
*
1449
ELSE IF( WNTQA ) THEN
1450
*
1451
* Path 9t (N much larger than M, JOBZ='A')
1452
* N right singular vectors to be computed in VT and
1453
* M left singular vectors to be computed in U
1454
*
1455
IVT = 1
1456
*
1457
* WORK(IVT) is M by M
1458
*
1459
LDWKVT = M
1460
ITAU = IVT + LDWKVT*M
1461
NWORK = ITAU + M
1462
*
1463
* Compute A=L*Q, copying result to VT
1464
* (CWorkspace: need 2*M, prefer M+M*NB)
1465
* (RWorkspace: 0)
1466
*
1467
CALL CGELQF( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
1468
$ LWORK-NWORK+1, IERR )
1469
CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT )
1470
*
1471
* Generate Q in VT
1472
* (CWorkspace: need M+N, prefer M+N*NB)
1473
* (RWorkspace: 0)
1474
*
1475
CALL CUNGLQ( N, N, M, VT, LDVT, WORK( ITAU ),
1476
$ WORK( NWORK ), LWORK-NWORK+1, IERR )
1477
*
1478
* Produce L in A, zeroing out above it
1479
*
1480
CALL CLASET( 'U', M-1, M-1, CZERO, CZERO, A( 1, 2 ),
1481
$ LDA )
1482
IE = 1
1483
ITAUQ = ITAU
1484
ITAUP = ITAUQ + M
1485
NWORK = ITAUP + M
1486
*
1487
* Bidiagonalize L in A
1488
* (CWorkspace: need M*M+3*M, prefer M*M+2*M+2*M*NB)
1489
* (RWorkspace: need M)
1490
*
1491
CALL CGEBRD( M, M, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
1492
$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1493
$ IERR )
1494
*
1495
* Perform bidiagonal SVD, computing left singular vectors
1496
* of bidiagonal matrix in RWORK(IRU) and computing right
1497
* singular vectors of bidiagonal matrix in RWORK(IRVT)
1498
* (CWorkspace: need 0)
1499
* (RWorkspace: need BDSPAC)
1500
*
1501
IRU = IE + M
1502
IRVT = IRU + M*M
1503
NRWORK = IRVT + M*M
1504
CALL SBDSDC( 'U', 'I', M, S, RWORK( IE ), RWORK( IRU ),
1505
$ M, RWORK( IRVT ), M, DUM, IDUM,
1506
$ RWORK( NRWORK ), IWORK, INFO )
1507
*
1508
* Copy real matrix RWORK(IRU) to complex matrix U
1509
* Overwrite U by left singular vectors of L
1510
* (CWorkspace: need 3*M, prefer 2*M+M*NB)
1511
* (RWorkspace: 0)
1512
*
1513
CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
1514
CALL CUNMBR( 'Q', 'L', 'N', M, M, M, A, LDA,
1515
$ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1516
$ LWORK-NWORK+1, IERR )
1517
*
1518
* Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
1519
* Overwrite WORK(IVT) by right singular vectors of L
1520
* (CWorkspace: need M*M+3*M, prefer M*M+2*M+M*NB)
1521
* (RWorkspace: 0)
1522
*
1523
CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ),
1524
$ LDWKVT )
1525
CALL CUNMBR( 'P', 'R', 'C', M, M, M, A, LDA,
1526
$ WORK( ITAUP ), WORK( IVT ), LDWKVT,
1527
$ WORK( NWORK ), LWORK-NWORK+1, IERR )
1528
*
1529
* Multiply right singular vectors of L in WORK(IVT) by
1530
* Q in VT, storing result in A
1531
* (CWorkspace: need M*M)
1532
* (RWorkspace: 0)
1533
*
1534
CALL CGEMM( 'N', 'N', M, N, M, CONE, WORK( IVT ), LDWKVT,
1535
$ VT, LDVT, CZERO, A, LDA )
1536
*
1537
* Copy right singular vectors of A from A to VT
1538
*
1539
CALL CLACPY( 'F', M, N, A, LDA, VT, LDVT )
1540
*
1541
END IF
1542
*
1543
ELSE IF( N.GE.MNTHR2 ) THEN
1544
*
1545
* MNTHR2 <= N < MNTHR1
1546
*
1547
* Path 5t (N much larger than M, but not as much as MNTHR1)
1548
* Reduce to bidiagonal form without QR decomposition, use
1549
* CUNGBR and matrix multiplication to compute singular vectors
1550
*
1551
*
1552
IE = 1
1553
NRWORK = IE + M
1554
ITAUQ = 1
1555
ITAUP = ITAUQ + M
1556
NWORK = ITAUP + M
1557
*
1558
* Bidiagonalize A
1559
* (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)
1560
* (RWorkspace: M)
1561
*
1562
CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
1563
$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1564
$ IERR )
1565
*
1566
IF( WNTQN ) THEN
1567
*
1568
* Compute singular values only
1569
* (Cworkspace: 0)
1570
* (Rworkspace: need BDSPAC)
1571
*
1572
CALL SBDSDC( 'L', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1,
1573
$ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
1574
ELSE IF( WNTQO ) THEN
1575
IRVT = NRWORK
1576
IRU = IRVT + M*M
1577
NRWORK = IRU + M*M
1578
IVT = NWORK
1579
*
1580
* Copy A to U, generate Q
1581
* (Cworkspace: need 2*M, prefer M+M*NB)
1582
* (Rworkspace: 0)
1583
*
1584
CALL CLACPY( 'L', M, M, A, LDA, U, LDU )
1585
CALL CUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
1586
$ WORK( NWORK ), LWORK-NWORK+1, IERR )
1587
*
1588
* Generate P**H in A
1589
* (Cworkspace: need 2*M, prefer M+M*NB)
1590
* (Rworkspace: 0)
1591
*
1592
CALL CUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
1593
$ WORK( NWORK ), LWORK-NWORK+1, IERR )
1594
*
1595
LDWKVT = M
1596
IF( LWORK.GE.M*N+3*M ) THEN
1597
*
1598
* WORK( IVT ) is M by N
1599
*
1600
NWORK = IVT + LDWKVT*N
1601
CHUNK = N
1602
ELSE
1603
*
1604
* WORK( IVT ) is M by CHUNK
1605
*
1606
CHUNK = ( LWORK-3*M ) / M
1607
NWORK = IVT + LDWKVT*CHUNK
1608
END IF
1609
*
1610
* Perform bidiagonal SVD, computing left singular vectors
1611
* of bidiagonal matrix in RWORK(IRU) and computing right
1612
* singular vectors of bidiagonal matrix in RWORK(IRVT)
1613
* (CWorkspace: need 0)
1614
* (RWorkspace: need BDSPAC)
1615
*
1616
CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
1617
$ M, RWORK( IRVT ), M, DUM, IDUM,
1618
$ RWORK( NRWORK ), IWORK, INFO )
1619
*
1620
* Multiply Q in U by real matrix RWORK(IRVT)
1621
* storing the result in WORK(IVT), copying to U
1622
* (Cworkspace: need 0)
1623
* (Rworkspace: need 2*M*M)
1624
*
1625
CALL CLACRM( M, M, U, LDU, RWORK( IRU ), M, WORK( IVT ),
1626
$ LDWKVT, RWORK( NRWORK ) )
1627
CALL CLACPY( 'F', M, M, WORK( IVT ), LDWKVT, U, LDU )
1628
*
1629
* Multiply RWORK(IRVT) by P**H in A, storing the
1630
* result in WORK(IVT), copying to A
1631
* (CWorkspace: need M*M, prefer M*N)
1632
* (Rworkspace: need 2*M*M, prefer 2*M*N)
1633
*
1634
NRWORK = IRU
1635
DO 50 I = 1, N, CHUNK
1636
BLK = MIN( N-I+1, CHUNK )
1637
CALL CLARCM( M, BLK, RWORK( IRVT ), M, A( 1, I ), LDA,
1638
$ WORK( IVT ), LDWKVT, RWORK( NRWORK ) )
1639
CALL CLACPY( 'F', M, BLK, WORK( IVT ), LDWKVT,
1640
$ A( 1, I ), LDA )
1641
50 CONTINUE
1642
ELSE IF( WNTQS ) THEN
1643
*
1644
* Copy A to U, generate Q
1645
* (Cworkspace: need 2*M, prefer M+M*NB)
1646
* (Rworkspace: 0)
1647
*
1648
CALL CLACPY( 'L', M, M, A, LDA, U, LDU )
1649
CALL CUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
1650
$ WORK( NWORK ), LWORK-NWORK+1, IERR )
1651
*
1652
* Copy A to VT, generate P**H
1653
* (Cworkspace: need 2*M, prefer M+M*NB)
1654
* (Rworkspace: 0)
1655
*
1656
CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT )
1657
CALL CUNGBR( 'P', M, N, M, VT, LDVT, WORK( ITAUP ),
1658
$ WORK( NWORK ), LWORK-NWORK+1, IERR )
1659
*
1660
* Perform bidiagonal SVD, computing left singular vectors
1661
* of bidiagonal matrix in RWORK(IRU) and computing right
1662
* singular vectors of bidiagonal matrix in RWORK(IRVT)
1663
* (CWorkspace: need 0)
1664
* (RWorkspace: need BDSPAC)
1665
*
1666
IRVT = NRWORK
1667
IRU = IRVT + M*M
1668
NRWORK = IRU + M*M
1669
CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
1670
$ M, RWORK( IRVT ), M, DUM, IDUM,
1671
$ RWORK( NRWORK ), IWORK, INFO )
1672
*
1673
* Multiply Q in U by real matrix RWORK(IRU), storing the
1674
* result in A, copying to U
1675
* (CWorkspace: need 0)
1676
* (Rworkspace: need 3*M*M)
1677
*
1678
CALL CLACRM( M, M, U, LDU, RWORK( IRU ), M, A, LDA,
1679
$ RWORK( NRWORK ) )
1680
CALL CLACPY( 'F', M, M, A, LDA, U, LDU )
1681
*
1682
* Multiply real matrix RWORK(IRVT) by P**H in VT,
1683
* storing the result in A, copying to VT
1684
* (Cworkspace: need 0)
1685
* (Rworkspace: need M*M+2*M*N)
1686
*
1687
NRWORK = IRU
1688
CALL CLARCM( M, N, RWORK( IRVT ), M, VT, LDVT, A, LDA,
1689
$ RWORK( NRWORK ) )
1690
CALL CLACPY( 'F', M, N, A, LDA, VT, LDVT )
1691
ELSE
1692
*
1693
* Copy A to U, generate Q
1694
* (Cworkspace: need 2*M, prefer M+M*NB)
1695
* (Rworkspace: 0)
1696
*
1697
CALL CLACPY( 'L', M, M, A, LDA, U, LDU )
1698
CALL CUNGBR( 'Q', M, M, N, U, LDU, WORK( ITAUQ ),
1699
$ WORK( NWORK ), LWORK-NWORK+1, IERR )
1700
*
1701
* Copy A to VT, generate P**H
1702
* (Cworkspace: need 2*M, prefer M+M*NB)
1703
* (Rworkspace: 0)
1704
*
1705
CALL CLACPY( 'U', M, N, A, LDA, VT, LDVT )
1706
CALL CUNGBR( 'P', N, N, M, VT, LDVT, WORK( ITAUP ),
1707
$ WORK( NWORK ), LWORK-NWORK+1, IERR )
1708
*
1709
* Perform bidiagonal SVD, computing left singular vectors
1710
* of bidiagonal matrix in RWORK(IRU) and computing right
1711
* singular vectors of bidiagonal matrix in RWORK(IRVT)
1712
* (CWorkspace: need 0)
1713
* (RWorkspace: need BDSPAC)
1714
*
1715
IRVT = NRWORK
1716
IRU = IRVT + M*M
1717
NRWORK = IRU + M*M
1718
CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
1719
$ M, RWORK( IRVT ), M, DUM, IDUM,
1720
$ RWORK( NRWORK ), IWORK, INFO )
1721
*
1722
* Multiply Q in U by real matrix RWORK(IRU), storing the
1723
* result in A, copying to U
1724
* (CWorkspace: need 0)
1725
* (Rworkspace: need 3*M*M)
1726
*
1727
CALL CLACRM( M, M, U, LDU, RWORK( IRU ), M, A, LDA,
1728
$ RWORK( NRWORK ) )
1729
CALL CLACPY( 'F', M, M, A, LDA, U, LDU )
1730
*
1731
* Multiply real matrix RWORK(IRVT) by P**H in VT,
1732
* storing the result in A, copying to VT
1733
* (Cworkspace: need 0)
1734
* (Rworkspace: need M*M+2*M*N)
1735
*
1736
CALL CLARCM( M, N, RWORK( IRVT ), M, VT, LDVT, A, LDA,
1737
$ RWORK( NRWORK ) )
1738
CALL CLACPY( 'F', M, N, A, LDA, VT, LDVT )
1739
END IF
1740
*
1741
ELSE
1742
*
1743
* N .LT. MNTHR2
1744
*
1745
* Path 6t (N greater than M, but not much larger)
1746
* Reduce to bidiagonal form without LQ decomposition
1747
* Use CUNMBR to compute singular vectors
1748
*
1749
IE = 1
1750
NRWORK = IE + M
1751
ITAUQ = 1
1752
ITAUP = ITAUQ + M
1753
NWORK = ITAUP + M
1754
*
1755
* Bidiagonalize A
1756
* (CWorkspace: need 2*M+N, prefer 2*M+(M+N)*NB)
1757
* (RWorkspace: M)
1758
*
1759
CALL CGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
1760
$ WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
1761
$ IERR )
1762
IF( WNTQN ) THEN
1763
*
1764
* Compute singular values only
1765
* (Cworkspace: 0)
1766
* (Rworkspace: need BDSPAC)
1767
*
1768
CALL SBDSDC( 'L', 'N', M, S, RWORK( IE ), DUM, 1, DUM, 1,
1769
$ DUM, IDUM, RWORK( NRWORK ), IWORK, INFO )
1770
ELSE IF( WNTQO ) THEN
1771
LDWKVT = M
1772
IVT = NWORK
1773
IF( LWORK.GE.M*N+3*M ) THEN
1774
*
1775
* WORK( IVT ) is M by N
1776
*
1777
CALL CLASET( 'F', M, N, CZERO, CZERO, WORK( IVT ),
1778
$ LDWKVT )
1779
NWORK = IVT + LDWKVT*N
1780
ELSE
1781
*
1782
* WORK( IVT ) is M by CHUNK
1783
*
1784
CHUNK = ( LWORK-3*M ) / M
1785
NWORK = IVT + LDWKVT*CHUNK
1786
END IF
1787
*
1788
* Perform bidiagonal SVD, computing left singular vectors
1789
* of bidiagonal matrix in RWORK(IRU) and computing right
1790
* singular vectors of bidiagonal matrix in RWORK(IRVT)
1791
* (CWorkspace: need 0)
1792
* (RWorkspace: need BDSPAC)
1793
*
1794
IRVT = NRWORK
1795
IRU = IRVT + M*M
1796
NRWORK = IRU + M*M
1797
CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
1798
$ M, RWORK( IRVT ), M, DUM, IDUM,
1799
$ RWORK( NRWORK ), IWORK, INFO )
1800
*
1801
* Copy real matrix RWORK(IRU) to complex matrix U
1802
* Overwrite U by left singular vectors of A
1803
* (Cworkspace: need 2*M, prefer M+M*NB)
1804
* (Rworkspace: need 0)
1805
*
1806
CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
1807
CALL CUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1808
$ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1809
$ LWORK-NWORK+1, IERR )
1810
*
1811
IF( LWORK.GE.M*N+3*M ) THEN
1812
*
1813
* Copy real matrix RWORK(IRVT) to complex matrix WORK(IVT)
1814
* Overwrite WORK(IVT) by right singular vectors of A,
1815
* copying to A
1816
* (Cworkspace: need M*N+2*M, prefer M*N+M+M*NB)
1817
* (Rworkspace: need 0)
1818
*
1819
CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, WORK( IVT ),
1820
$ LDWKVT )
1821
CALL CUNMBR( 'P', 'R', 'C', M, N, M, A, LDA,
1822
$ WORK( ITAUP ), WORK( IVT ), LDWKVT,
1823
$ WORK( NWORK ), LWORK-NWORK+1, IERR )
1824
CALL CLACPY( 'F', M, N, WORK( IVT ), LDWKVT, A, LDA )
1825
ELSE
1826
*
1827
* Generate P**H in A
1828
* (Cworkspace: need 2*M, prefer M+M*NB)
1829
* (Rworkspace: need 0)
1830
*
1831
CALL CUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
1832
$ WORK( NWORK ), LWORK-NWORK+1, IERR )
1833
*
1834
* Multiply Q in A by real matrix RWORK(IRU), storing the
1835
* result in WORK(IU), copying to A
1836
* (CWorkspace: need M*M, prefer M*N)
1837
* (Rworkspace: need 3*M*M, prefer M*M+2*M*N)
1838
*
1839
NRWORK = IRU
1840
DO 60 I = 1, N, CHUNK
1841
BLK = MIN( N-I+1, CHUNK )
1842
CALL CLARCM( M, BLK, RWORK( IRVT ), M, A( 1, I ),
1843
$ LDA, WORK( IVT ), LDWKVT,
1844
$ RWORK( NRWORK ) )
1845
CALL CLACPY( 'F', M, BLK, WORK( IVT ), LDWKVT,
1846
$ A( 1, I ), LDA )
1847
60 CONTINUE
1848
END IF
1849
ELSE IF( WNTQS ) THEN
1850
*
1851
* Perform bidiagonal SVD, computing left singular vectors
1852
* of bidiagonal matrix in RWORK(IRU) and computing right
1853
* singular vectors of bidiagonal matrix in RWORK(IRVT)
1854
* (CWorkspace: need 0)
1855
* (RWorkspace: need BDSPAC)
1856
*
1857
IRVT = NRWORK
1858
IRU = IRVT + M*M
1859
NRWORK = IRU + M*M
1860
CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
1861
$ M, RWORK( IRVT ), M, DUM, IDUM,
1862
$ RWORK( NRWORK ), IWORK, INFO )
1863
*
1864
* Copy real matrix RWORK(IRU) to complex matrix U
1865
* Overwrite U by left singular vectors of A
1866
* (CWorkspace: need 3*M, prefer 2*M+M*NB)
1867
* (RWorkspace: M*M)
1868
*
1869
CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
1870
CALL CUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1871
$ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1872
$ LWORK-NWORK+1, IERR )
1873
*
1874
* Copy real matrix RWORK(IRVT) to complex matrix VT
1875
* Overwrite VT by right singular vectors of A
1876
* (CWorkspace: need 3*M, prefer 2*M+M*NB)
1877
* (RWorkspace: M*M)
1878
*
1879
CALL CLASET( 'F', M, N, CZERO, CZERO, VT, LDVT )
1880
CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT )
1881
CALL CUNMBR( 'P', 'R', 'C', M, N, M, A, LDA,
1882
$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1883
$ LWORK-NWORK+1, IERR )
1884
ELSE
1885
*
1886
* Perform bidiagonal SVD, computing left singular vectors
1887
* of bidiagonal matrix in RWORK(IRU) and computing right
1888
* singular vectors of bidiagonal matrix in RWORK(IRVT)
1889
* (CWorkspace: need 0)
1890
* (RWorkspace: need BDSPAC)
1891
*
1892
IRVT = NRWORK
1893
IRU = IRVT + M*M
1894
NRWORK = IRU + M*M
1895
*
1896
CALL SBDSDC( 'L', 'I', M, S, RWORK( IE ), RWORK( IRU ),
1897
$ M, RWORK( IRVT ), M, DUM, IDUM,
1898
$ RWORK( NRWORK ), IWORK, INFO )
1899
*
1900
* Copy real matrix RWORK(IRU) to complex matrix U
1901
* Overwrite U by left singular vectors of A
1902
* (CWorkspace: need 3*M, prefer 2*M+M*NB)
1903
* (RWorkspace: M*M)
1904
*
1905
CALL CLACP2( 'F', M, M, RWORK( IRU ), M, U, LDU )
1906
CALL CUNMBR( 'Q', 'L', 'N', M, M, N, A, LDA,
1907
$ WORK( ITAUQ ), U, LDU, WORK( NWORK ),
1908
$ LWORK-NWORK+1, IERR )
1909
*
1910
* Set the right corner of VT to identity matrix
1911
*
1912
CALL CLASET( 'F', N-M, N-M, CZERO, CONE, VT( M+1, M+1 ),
1913
$ LDVT )
1914
*
1915
* Copy real matrix RWORK(IRVT) to complex matrix VT
1916
* Overwrite VT by right singular vectors of A
1917
* (CWorkspace: need 2*M+N, prefer 2*M+N*NB)
1918
* (RWorkspace: M*M)
1919
*
1920
CALL CLASET( 'F', N, N, CZERO, CZERO, VT, LDVT )
1921
CALL CLACP2( 'F', M, M, RWORK( IRVT ), M, VT, LDVT )
1922
CALL CUNMBR( 'P', 'R', 'C', N, N, M, A, LDA,
1923
$ WORK( ITAUP ), VT, LDVT, WORK( NWORK ),
1924
$ LWORK-NWORK+1, IERR )
1925
END IF
1926
*
1927
END IF
1928
*
1929
END IF
1930
*
1931
* Undo scaling if necessary
1932
*
1933
IF( ISCL.EQ.1 ) THEN
1934
IF( ANRM.GT.BIGNUM )
1935
$ CALL SLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
1936
$ IERR )
1937
IF( ANRM.LT.SMLNUM )
1938
$ CALL SLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
1939
$ IERR )
1940
END IF
1941
*
1942
* Return optimal workspace in WORK(1)
1943
*
1944
WORK( 1 ) = MAXWRK
1945
*
1946
RETURN
1947
*
1948
* End of CGESDD
1949
*
1950
END
1951
1952