LAPACK  3.9.0
LAPACK: Linear Algebra PACKage
cgesvdx.f
Go to the documentation of this file.
1 *> \brief <b> CGESVDX computes the singular value decomposition (SVD) for GE matrices</b>
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CGESVDX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgesvdx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgesvdx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgesvdx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CGESVDX( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
22 * $ IL, IU, NS, S, U, LDU, VT, LDVT, WORK,
23 * $ LWORK, RWORK, IWORK, INFO )
24 *
25 *
26 * .. Scalar Arguments ..
27 * CHARACTER JOBU, JOBVT, RANGE
28 * INTEGER IL, INFO, IU, LDA, LDU, LDVT, LWORK, M, N, NS
29 * REAL VL, VU
30 * ..
31 * .. Array Arguments ..
32 * INTEGER IWORK( * )
33 * REAL S( * ), RWORK( * )
34 * COMPLEX A( LDA, * ), U( LDU, * ), VT( LDVT, * ),
35 * $ WORK( * )
36 * ..
37 *
38 *
39 *> \par Purpose:
40 * =============
41 *>
42 *> \verbatim
43 *>
44 *> CGESVDX computes the singular value decomposition (SVD) of a complex
45 *> M-by-N matrix A, optionally computing the left and/or right singular
46 *> vectors. The SVD is written
47 *>
48 *> A = U * SIGMA * transpose(V)
49 *>
50 *> where SIGMA is an M-by-N matrix which is zero except for its
51 *> min(m,n) diagonal elements, U is an M-by-M unitary matrix, and
52 *> V is an N-by-N unitary matrix. The diagonal elements of SIGMA
53 *> are the singular values of A; they are real and non-negative, and
54 *> are returned in descending order. The first min(m,n) columns of
55 *> U and V are the left and right singular vectors of A.
56 *>
57 *> CGESVDX uses an eigenvalue problem for obtaining the SVD, which
58 *> allows for the computation of a subset of singular values and
59 *> vectors. See SBDSVDX for details.
60 *>
61 *> Note that the routine returns V**T, not V.
62 *> \endverbatim
63 *
64 * Arguments:
65 * ==========
66 *
67 *> \param[in] JOBU
68 *> \verbatim
69 *> JOBU is CHARACTER*1
70 *> Specifies options for computing all or part of the matrix U:
71 *> = 'V': the first min(m,n) columns of U (the left singular
72 *> vectors) or as specified by RANGE are returned in
73 *> the array U;
74 *> = 'N': no columns of U (no left singular vectors) are
75 *> computed.
76 *> \endverbatim
77 *>
78 *> \param[in] JOBVT
79 *> \verbatim
80 *> JOBVT is CHARACTER*1
81 *> Specifies options for computing all or part of the matrix
82 *> V**T:
83 *> = 'V': the first min(m,n) rows of V**T (the right singular
84 *> vectors) or as specified by RANGE are returned in
85 *> the array VT;
86 *> = 'N': no rows of V**T (no right singular vectors) are
87 *> computed.
88 *> \endverbatim
89 *>
90 *> \param[in] RANGE
91 *> \verbatim
92 *> RANGE is CHARACTER*1
93 *> = 'A': all singular values will be found.
94 *> = 'V': all singular values in the half-open interval (VL,VU]
95 *> will be found.
96 *> = 'I': the IL-th through IU-th singular values will be found.
97 *> \endverbatim
98 *>
99 *> \param[in] M
100 *> \verbatim
101 *> M is INTEGER
102 *> The number of rows of the input matrix A. M >= 0.
103 *> \endverbatim
104 *>
105 *> \param[in] N
106 *> \verbatim
107 *> N is INTEGER
108 *> The number of columns of the input matrix A. N >= 0.
109 *> \endverbatim
110 *>
111 *> \param[in,out] A
112 *> \verbatim
113 *> A is COMPLEX array, dimension (LDA,N)
114 *> On entry, the M-by-N matrix A.
115 *> On exit, the contents of A are destroyed.
116 *> \endverbatim
117 *>
118 *> \param[in] LDA
119 *> \verbatim
120 *> LDA is INTEGER
121 *> The leading dimension of the array A. LDA >= max(1,M).
122 *> \endverbatim
123 *>
124 *> \param[in] VL
125 *> \verbatim
126 *> VL is REAL
127 *> If RANGE='V', the lower bound of the interval to
128 *> be searched for singular values. VU > VL.
129 *> Not referenced if RANGE = 'A' or 'I'.
130 *> \endverbatim
131 *>
132 *> \param[in] VU
133 *> \verbatim
134 *> VU is REAL
135 *> If RANGE='V', the upper bound of the interval to
136 *> be searched for singular values. VU > VL.
137 *> Not referenced if RANGE = 'A' or 'I'.
138 *> \endverbatim
139 *>
140 *> \param[in] IL
141 *> \verbatim
142 *> IL is INTEGER
143 *> If RANGE='I', the index of the
144 *> smallest singular value to be returned.
145 *> 1 <= IL <= IU <= min(M,N), if min(M,N) > 0.
146 *> Not referenced if RANGE = 'A' or 'V'.
147 *> \endverbatim
148 *>
149 *> \param[in] IU
150 *> \verbatim
151 *> IU is INTEGER
152 *> If RANGE='I', the index of the
153 *> largest singular value to be returned.
154 *> 1 <= IL <= IU <= min(M,N), if min(M,N) > 0.
155 *> Not referenced if RANGE = 'A' or 'V'.
156 *> \endverbatim
157 *>
158 *> \param[out] NS
159 *> \verbatim
160 *> NS is INTEGER
161 *> The total number of singular values found,
162 *> 0 <= NS <= min(M,N).
163 *> If RANGE = 'A', NS = min(M,N); if RANGE = 'I', NS = IU-IL+1.
164 *> \endverbatim
165 *>
166 *> \param[out] S
167 *> \verbatim
168 *> S is REAL array, dimension (min(M,N))
169 *> The singular values of A, sorted so that S(i) >= S(i+1).
170 *> \endverbatim
171 *>
172 *> \param[out] U
173 *> \verbatim
174 *> U is COMPLEX array, dimension (LDU,UCOL)
175 *> If JOBU = 'V', U contains columns of U (the left singular
176 *> vectors, stored columnwise) as specified by RANGE; if
177 *> JOBU = 'N', U is not referenced.
178 *> Note: The user must ensure that UCOL >= NS; if RANGE = 'V',
179 *> the exact value of NS is not known in advance and an upper
180 *> bound must be used.
181 *> \endverbatim
182 *>
183 *> \param[in] LDU
184 *> \verbatim
185 *> LDU is INTEGER
186 *> The leading dimension of the array U. LDU >= 1; if
187 *> JOBU = 'V', LDU >= M.
188 *> \endverbatim
189 *>
190 *> \param[out] VT
191 *> \verbatim
192 *> VT is COMPLEX array, dimension (LDVT,N)
193 *> If JOBVT = 'V', VT contains the rows of V**T (the right singular
194 *> vectors, stored rowwise) as specified by RANGE; if JOBVT = 'N',
195 *> VT is not referenced.
196 *> Note: The user must ensure that LDVT >= NS; if RANGE = 'V',
197 *> the exact value of NS is not known in advance and an upper
198 *> bound must be used.
199 *> \endverbatim
200 *>
201 *> \param[in] LDVT
202 *> \verbatim
203 *> LDVT is INTEGER
204 *> The leading dimension of the array VT. LDVT >= 1; if
205 *> JOBVT = 'V', LDVT >= NS (see above).
206 *> \endverbatim
207 *>
208 *> \param[out] WORK
209 *> \verbatim
210 *> WORK is COMPLEX array, dimension (MAX(1,LWORK))
211 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK;
212 *> \endverbatim
213 *>
214 *> \param[in] LWORK
215 *> \verbatim
216 *> LWORK is INTEGER
217 *> The dimension of the array WORK.
218 *> LWORK >= MAX(1,MIN(M,N)*(MIN(M,N)+4)) for the paths (see
219 *> comments inside the code):
220 *> - PATH 1 (M much larger than N)
221 *> - PATH 1t (N much larger than M)
222 *> LWORK >= MAX(1,MIN(M,N)*2+MAX(M,N)) for the other paths.
223 *> For good performance, LWORK should generally be larger.
224 *>
225 *> If LWORK = -1, then a workspace query is assumed; the routine
226 *> only calculates the optimal size of the WORK array, returns
227 *> this value as the first entry of the WORK array, and no error
228 *> message related to LWORK is issued by XERBLA.
229 *> \endverbatim
230 *>
231 *> \param[out] RWORK
232 *> \verbatim
233 *> RWORK is REAL array, dimension (MAX(1,LRWORK))
234 *> LRWORK >= MIN(M,N)*(MIN(M,N)*2+15*MIN(M,N)).
235 *> \endverbatim
236 *>
237 *> \param[out] IWORK
238 *> \verbatim
239 *> IWORK is INTEGER array, dimension (12*MIN(M,N))
240 *> If INFO = 0, the first NS elements of IWORK are zero. If INFO > 0,
241 *> then IWORK contains the indices of the eigenvectors that failed
242 *> to converge in SBDSVDX/SSTEVX.
243 *> \endverbatim
244 *>
245 *> \param[out] INFO
246 *> \verbatim
247 *> INFO is INTEGER
248 *> = 0: successful exit
249 *> < 0: if INFO = -i, the i-th argument had an illegal value
250 *> > 0: if INFO = i, then i eigenvectors failed to converge
251 *> in SBDSVDX/SSTEVX.
252 *> if INFO = N*2 + 1, an internal error occurred in
253 *> SBDSVDX
254 *> \endverbatim
255 *
256 * Authors:
257 * ========
258 *
259 *> \author Univ. of Tennessee
260 *> \author Univ. of California Berkeley
261 *> \author Univ. of Colorado Denver
262 *> \author NAG Ltd.
263 *
264 *> \date June 2016
265 *
266 *> \ingroup complexGEsing
267 *
268 * =====================================================================
269  SUBROUTINE cgesvdx( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
270  $ IL, IU, NS, S, U, LDU, VT, LDVT, WORK,
271  $ LWORK, RWORK, IWORK, INFO )
272 *
273 * -- LAPACK driver routine (version 3.8.0) --
274 * -- LAPACK is a software package provided by Univ. of Tennessee, --
275 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
276 * June 2016
277 *
278 * .. Scalar Arguments ..
279  CHARACTER JOBU, JOBVT, RANGE
280  INTEGER IL, INFO, IU, LDA, LDU, LDVT, LWORK, M, N, NS
281  REAL VL, VU
282 * ..
283 * .. Array Arguments ..
284  INTEGER IWORK( * )
285  REAL S( * ), RWORK( * )
286  COMPLEX A( LDA, * ), U( LDU, * ), VT( LDVT, * ),
287  $ work( * )
288 * ..
289 *
290 * =====================================================================
291 *
292 * .. Parameters ..
293  COMPLEX CZERO, CONE
294  PARAMETER ( CZERO = ( 0.0e0, 0.0e0 ),
295  $ cone = ( 1.0e0, 0.0e0 ) )
296  REAL ZERO, ONE
297  PARAMETER ( ZERO = 0.0e0, one = 1.0e0 )
298 * ..
299 * .. Local Scalars ..
300  CHARACTER JOBZ, RNGTGK
301  LOGICAL ALLS, INDS, LQUERY, VALS, WANTU, WANTVT
302  INTEGER I, ID, IE, IERR, ILQF, ILTGK, IQRF, ISCL,
303  $ itau, itaup, itauq, itemp, itempr, itgkz,
304  $ iutgk, j, k, maxwrk, minmn, minwrk, mnthr
305  REAL ABSTOL, ANRM, BIGNUM, EPS, SMLNUM
306 * ..
307 * .. Local Arrays ..
308  REAL DUM( 1 )
309 * ..
310 * .. External Subroutines ..
311  EXTERNAL cgebrd, cgelqf, cgeqrf, clascl, claset,
312  $ cunmbr, cunmqr, cunmlq, clacpy,
313  $ sbdsvdx, slascl, xerbla
314 * ..
315 * .. External Functions ..
316  LOGICAL LSAME
317  INTEGER ILAENV
318  REAL SLAMCH, CLANGE
319  EXTERNAL lsame, ilaenv, slamch, clange
320 * ..
321 * .. Intrinsic Functions ..
322  INTRINSIC max, min, sqrt
323 * ..
324 * .. Executable Statements ..
325 *
326 * Test the input arguments.
327 *
328  ns = 0
329  info = 0
330  abstol = 2*slamch('S')
331  lquery = ( lwork.EQ.-1 )
332  minmn = min( m, n )
333 
334  wantu = lsame( jobu, 'V' )
335  wantvt = lsame( jobvt, 'V' )
336  IF( wantu .OR. wantvt ) THEN
337  jobz = 'V'
338  ELSE
339  jobz = 'N'
340  END IF
341  alls = lsame( range, 'A' )
342  vals = lsame( range, 'V' )
343  inds = lsame( range, 'I' )
344 *
345  info = 0
346  IF( .NOT.lsame( jobu, 'V' ) .AND.
347  $ .NOT.lsame( jobu, 'N' ) ) THEN
348  info = -1
349  ELSE IF( .NOT.lsame( jobvt, 'V' ) .AND.
350  $ .NOT.lsame( jobvt, 'N' ) ) THEN
351  info = -2
352  ELSE IF( .NOT.( alls .OR. vals .OR. inds ) ) THEN
353  info = -3
354  ELSE IF( m.LT.0 ) THEN
355  info = -4
356  ELSE IF( n.LT.0 ) THEN
357  info = -5
358  ELSE IF( m.GT.lda ) THEN
359  info = -7
360  ELSE IF( minmn.GT.0 ) THEN
361  IF( vals ) THEN
362  IF( vl.LT.zero ) THEN
363  info = -8
364  ELSE IF( vu.LE.vl ) THEN
365  info = -9
366  END IF
367  ELSE IF( inds ) THEN
368  IF( il.LT.1 .OR. il.GT.max( 1, minmn ) ) THEN
369  info = -10
370  ELSE IF( iu.LT.min( minmn, il ) .OR. iu.GT.minmn ) THEN
371  info = -11
372  END IF
373  END IF
374  IF( info.EQ.0 ) THEN
375  IF( wantu .AND. ldu.LT.m ) THEN
376  info = -15
377  ELSE IF( wantvt ) THEN
378  IF( inds ) THEN
379  IF( ldvt.LT.iu-il+1 ) THEN
380  info = -17
381  END IF
382  ELSE IF( ldvt.LT.minmn ) THEN
383  info = -17
384  END IF
385  END IF
386  END IF
387  END IF
388 *
389 * Compute workspace
390 * (Note: Comments in the code beginning "Workspace:" describe the
391 * minimal amount of workspace needed at that point in the code,
392 * as well as the preferred amount for good performance.
393 * NB refers to the optimal block size for the immediately
394 * following subroutine, as returned by ILAENV.)
395 *
396  IF( info.EQ.0 ) THEN
397  minwrk = 1
398  maxwrk = 1
399  IF( minmn.GT.0 ) THEN
400  IF( m.GE.n ) THEN
401  mnthr = ilaenv( 6, 'CGESVD', jobu // jobvt, m, n, 0, 0 )
402  IF( m.GE.mnthr ) THEN
403 *
404 * Path 1 (M much larger than N)
405 *
406  minwrk = n*(n+5)
407  maxwrk = n + n*ilaenv(1,'CGEQRF',' ',m,n,-1,-1)
408  maxwrk = max(maxwrk,
409  $ n*n+2*n+2*n*ilaenv(1,'CGEBRD',' ',n,n,-1,-1))
410  IF (wantu .OR. wantvt) THEN
411  maxwrk = max(maxwrk,
412  $ n*n+2*n+n*ilaenv(1,'CUNMQR','LN',n,n,n,-1))
413  END IF
414  ELSE
415 *
416 * Path 2 (M at least N, but not much larger)
417 *
418  minwrk = 3*n + m
419  maxwrk = 2*n + (m+n)*ilaenv(1,'CGEBRD',' ',m,n,-1,-1)
420  IF (wantu .OR. wantvt) THEN
421  maxwrk = max(maxwrk,
422  $ 2*n+n*ilaenv(1,'CUNMQR','LN',n,n,n,-1))
423  END IF
424  END IF
425  ELSE
426  mnthr = ilaenv( 6, 'CGESVD', jobu // jobvt, m, n, 0, 0 )
427  IF( n.GE.mnthr ) THEN
428 *
429 * Path 1t (N much larger than M)
430 *
431  minwrk = m*(m+5)
432  maxwrk = m + m*ilaenv(1,'CGELQF',' ',m,n,-1,-1)
433  maxwrk = max(maxwrk,
434  $ m*m+2*m+2*m*ilaenv(1,'CGEBRD',' ',m,m,-1,-1))
435  IF (wantu .OR. wantvt) THEN
436  maxwrk = max(maxwrk,
437  $ m*m+2*m+m*ilaenv(1,'CUNMQR','LN',m,m,m,-1))
438  END IF
439  ELSE
440 *
441 * Path 2t (N greater than M, but not much larger)
442 *
443 *
444  minwrk = 3*m + n
445  maxwrk = 2*m + (m+n)*ilaenv(1,'CGEBRD',' ',m,n,-1,-1)
446  IF (wantu .OR. wantvt) THEN
447  maxwrk = max(maxwrk,
448  $ 2*m+m*ilaenv(1,'CUNMQR','LN',m,m,m,-1))
449  END IF
450  END IF
451  END IF
452  END IF
453  maxwrk = max( maxwrk, minwrk )
454  work( 1 ) = cmplx( real( maxwrk ), zero )
455 *
456  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
457  info = -19
458  END IF
459  END IF
460 *
461  IF( info.NE.0 ) THEN
462  CALL xerbla( 'CGESVDX', -info )
463  RETURN
464  ELSE IF( lquery ) THEN
465  RETURN
466  END IF
467 *
468 * Quick return if possible
469 *
470  IF( m.EQ.0 .OR. n.EQ.0 ) THEN
471  RETURN
472  END IF
473 *
474 * Set singular values indices accord to RANGE='A'.
475 *
476  IF( alls ) THEN
477  rngtgk = 'I'
478  iltgk = 1
479  iutgk = min( m, n )
480  ELSE IF( inds ) THEN
481  rngtgk = 'I'
482  iltgk = il
483  iutgk = iu
484  ELSE
485  rngtgk = 'V'
486  iltgk = 0
487  iutgk = 0
488  END IF
489 *
490 * Get machine constants
491 *
492  eps = slamch( 'P' )
493  smlnum = sqrt( slamch( 'S' ) ) / eps
494  bignum = one / smlnum
495 *
496 * Scale A if max element outside range [SMLNUM,BIGNUM]
497 *
498  anrm = clange( 'M', m, n, a, lda, dum )
499  iscl = 0
500  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
501  iscl = 1
502  CALL clascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
503  ELSE IF( anrm.GT.bignum ) THEN
504  iscl = 1
505  CALL clascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
506  END IF
507 *
508  IF( m.GE.n ) THEN
509 *
510 * A has at least as many rows as columns. If A has sufficiently
511 * more rows than columns, first reduce A using the QR
512 * decomposition.
513 *
514  IF( m.GE.mnthr ) THEN
515 *
516 * Path 1 (M much larger than N):
517 * A = Q * R = Q * ( QB * B * PB**T )
518 * = Q * ( QB * ( UB * S * VB**T ) * PB**T )
519 * U = Q * QB * UB; V**T = VB**T * PB**T
520 *
521 * Compute A=Q*R
522 * (Workspace: need 2*N, prefer N+N*NB)
523 *
524  itau = 1
525  itemp = itau + n
526  CALL cgeqrf( m, n, a, lda, work( itau ), work( itemp ),
527  $ lwork-itemp+1, info )
528 *
529 * Copy R into WORK and bidiagonalize it:
530 * (Workspace: need N*N+3*N, prefer N*N+N+2*N*NB)
531 *
532  iqrf = itemp
533  itauq = itemp + n*n
534  itaup = itauq + n
535  itemp = itaup + n
536  id = 1
537  ie = id + n
538  itgkz = ie + n
539  CALL clacpy( 'U', n, n, a, lda, work( iqrf ), n )
540  CALL claset( 'L', n-1, n-1, czero, czero,
541  $ work( iqrf+1 ), n )
542  CALL cgebrd( n, n, work( iqrf ), n, rwork( id ),
543  $ rwork( ie ), work( itauq ), work( itaup ),
544  $ work( itemp ), lwork-itemp+1, info )
545  itempr = itgkz + n*(n*2+1)
546 *
547 * Solve eigenvalue problem TGK*Z=Z*S.
548 * (Workspace: need 2*N*N+14*N)
549 *
550  CALL sbdsvdx( 'U', jobz, rngtgk, n, rwork( id ),
551  $ rwork( ie ), vl, vu, iltgk, iutgk, ns, s,
552  $ rwork( itgkz ), n*2, rwork( itempr ),
553  $ iwork, info)
554 *
555 * If needed, compute left singular vectors.
556 *
557  IF( wantu ) THEN
558  k = itgkz
559  DO i = 1, ns
560  DO j = 1, n
561  u( j, i ) = cmplx( rwork( k ), zero )
562  k = k + 1
563  END DO
564  k = k + n
565  END DO
566  CALL claset( 'A', m-n, ns, czero, czero, u( n+1,1 ), ldu)
567 *
568 * Call CUNMBR to compute QB*UB.
569 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
570 *
571  CALL cunmbr( 'Q', 'L', 'N', n, ns, n, work( iqrf ), n,
572  $ work( itauq ), u, ldu, work( itemp ),
573  $ lwork-itemp+1, info )
574 *
575 * Call CUNMQR to compute Q*(QB*UB).
576 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
577 *
578  CALL cunmqr( 'L', 'N', m, ns, n, a, lda,
579  $ work( itau ), u, ldu, work( itemp ),
580  $ lwork-itemp+1, info )
581  END IF
582 *
583 * If needed, compute right singular vectors.
584 *
585  IF( wantvt) THEN
586  k = itgkz + n
587  DO i = 1, ns
588  DO j = 1, n
589  vt( i, j ) = cmplx( rwork( k ), zero )
590  k = k + 1
591  END DO
592  k = k + n
593  END DO
594 *
595 * Call CUNMBR to compute VB**T * PB**T
596 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
597 *
598  CALL cunmbr( 'P', 'R', 'C', ns, n, n, work( iqrf ), n,
599  $ work( itaup ), vt, ldvt, work( itemp ),
600  $ lwork-itemp+1, info )
601  END IF
602  ELSE
603 *
604 * Path 2 (M at least N, but not much larger)
605 * Reduce A to bidiagonal form without QR decomposition
606 * A = QB * B * PB**T = QB * ( UB * S * VB**T ) * PB**T
607 * U = QB * UB; V**T = VB**T * PB**T
608 *
609 * Bidiagonalize A
610 * (Workspace: need 2*N+M, prefer 2*N+(M+N)*NB)
611 *
612  itauq = 1
613  itaup = itauq + n
614  itemp = itaup + n
615  id = 1
616  ie = id + n
617  itgkz = ie + n
618  CALL cgebrd( m, n, a, lda, rwork( id ), rwork( ie ),
619  $ work( itauq ), work( itaup ), work( itemp ),
620  $ lwork-itemp+1, info )
621  itempr = itgkz + n*(n*2+1)
622 *
623 * Solve eigenvalue problem TGK*Z=Z*S.
624 * (Workspace: need 2*N*N+14*N)
625 *
626  CALL sbdsvdx( 'U', jobz, rngtgk, n, rwork( id ),
627  $ rwork( ie ), vl, vu, iltgk, iutgk, ns, s,
628  $ rwork( itgkz ), n*2, rwork( itempr ),
629  $ iwork, info)
630 *
631 * If needed, compute left singular vectors.
632 *
633  IF( wantu ) THEN
634  k = itgkz
635  DO i = 1, ns
636  DO j = 1, n
637  u( j, i ) = cmplx( rwork( k ), zero )
638  k = k + 1
639  END DO
640  k = k + n
641  END DO
642  CALL claset( 'A', m-n, ns, czero, czero, u( n+1,1 ), ldu)
643 *
644 * Call CUNMBR to compute QB*UB.
645 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
646 *
647  CALL cunmbr( 'Q', 'L', 'N', m, ns, n, a, lda,
648  $ work( itauq ), u, ldu, work( itemp ),
649  $ lwork-itemp+1, ierr )
650  END IF
651 *
652 * If needed, compute right singular vectors.
653 *
654  IF( wantvt) THEN
655  k = itgkz + n
656  DO i = 1, ns
657  DO j = 1, n
658  vt( i, j ) = cmplx( rwork( k ), zero )
659  k = k + 1
660  END DO
661  k = k + n
662  END DO
663 *
664 * Call CUNMBR to compute VB**T * PB**T
665 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
666 *
667  CALL cunmbr( 'P', 'R', 'C', ns, n, n, a, lda,
668  $ work( itaup ), vt, ldvt, work( itemp ),
669  $ lwork-itemp+1, ierr )
670  END IF
671  END IF
672  ELSE
673 *
674 * A has more columns than rows. If A has sufficiently more
675 * columns than rows, first reduce A using the LQ decomposition.
676 *
677  IF( n.GE.mnthr ) THEN
678 *
679 * Path 1t (N much larger than M):
680 * A = L * Q = ( QB * B * PB**T ) * Q
681 * = ( QB * ( UB * S * VB**T ) * PB**T ) * Q
682 * U = QB * UB ; V**T = VB**T * PB**T * Q
683 *
684 * Compute A=L*Q
685 * (Workspace: need 2*M, prefer M+M*NB)
686 *
687  itau = 1
688  itemp = itau + m
689  CALL cgelqf( m, n, a, lda, work( itau ), work( itemp ),
690  $ lwork-itemp+1, info )
691 
692 * Copy L into WORK and bidiagonalize it:
693 * (Workspace in WORK( ITEMP ): need M*M+3*M, prefer M*M+M+2*M*NB)
694 *
695  ilqf = itemp
696  itauq = ilqf + m*m
697  itaup = itauq + m
698  itemp = itaup + m
699  id = 1
700  ie = id + m
701  itgkz = ie + m
702  CALL clacpy( 'L', m, m, a, lda, work( ilqf ), m )
703  CALL claset( 'U', m-1, m-1, czero, czero,
704  $ work( ilqf+m ), m )
705  CALL cgebrd( m, m, work( ilqf ), m, rwork( id ),
706  $ rwork( ie ), work( itauq ), work( itaup ),
707  $ work( itemp ), lwork-itemp+1, info )
708  itempr = itgkz + m*(m*2+1)
709 *
710 * Solve eigenvalue problem TGK*Z=Z*S.
711 * (Workspace: need 2*M*M+14*M)
712 *
713  CALL sbdsvdx( 'U', jobz, rngtgk, m, rwork( id ),
714  $ rwork( ie ), vl, vu, iltgk, iutgk, ns, s,
715  $ rwork( itgkz ), m*2, rwork( itempr ),
716  $ iwork, info)
717 *
718 * If needed, compute left singular vectors.
719 *
720  IF( wantu ) THEN
721  k = itgkz
722  DO i = 1, ns
723  DO j = 1, m
724  u( j, i ) = cmplx( rwork( k ), zero )
725  k = k + 1
726  END DO
727  k = k + m
728  END DO
729 *
730 * Call CUNMBR to compute QB*UB.
731 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
732 *
733  CALL cunmbr( 'Q', 'L', 'N', m, ns, m, work( ilqf ), m,
734  $ work( itauq ), u, ldu, work( itemp ),
735  $ lwork-itemp+1, info )
736  END IF
737 *
738 * If needed, compute right singular vectors.
739 *
740  IF( wantvt) THEN
741  k = itgkz + m
742  DO i = 1, ns
743  DO j = 1, m
744  vt( i, j ) = cmplx( rwork( k ), zero )
745  k = k + 1
746  END DO
747  k = k + m
748  END DO
749  CALL claset( 'A', ns, n-m, czero, czero,
750  $ vt( 1,m+1 ), ldvt )
751 *
752 * Call CUNMBR to compute (VB**T)*(PB**T)
753 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
754 *
755  CALL cunmbr( 'P', 'R', 'C', ns, m, m, work( ilqf ), m,
756  $ work( itaup ), vt, ldvt, work( itemp ),
757  $ lwork-itemp+1, info )
758 *
759 * Call CUNMLQ to compute ((VB**T)*(PB**T))*Q.
760 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
761 *
762  CALL cunmlq( 'R', 'N', ns, n, m, a, lda,
763  $ work( itau ), vt, ldvt, work( itemp ),
764  $ lwork-itemp+1, info )
765  END IF
766  ELSE
767 *
768 * Path 2t (N greater than M, but not much larger)
769 * Reduce to bidiagonal form without LQ decomposition
770 * A = QB * B * PB**T = QB * ( UB * S * VB**T ) * PB**T
771 * U = QB * UB; V**T = VB**T * PB**T
772 *
773 * Bidiagonalize A
774 * (Workspace: need 2*M+N, prefer 2*M+(M+N)*NB)
775 *
776  itauq = 1
777  itaup = itauq + m
778  itemp = itaup + m
779  id = 1
780  ie = id + m
781  itgkz = ie + m
782  CALL cgebrd( m, n, a, lda, rwork( id ), rwork( ie ),
783  $ work( itauq ), work( itaup ), work( itemp ),
784  $ lwork-itemp+1, info )
785  itempr = itgkz + m*(m*2+1)
786 *
787 * Solve eigenvalue problem TGK*Z=Z*S.
788 * (Workspace: need 2*M*M+14*M)
789 *
790  CALL sbdsvdx( 'L', jobz, rngtgk, m, rwork( id ),
791  $ rwork( ie ), vl, vu, iltgk, iutgk, ns, s,
792  $ rwork( itgkz ), m*2, rwork( itempr ),
793  $ iwork, info)
794 *
795 * If needed, compute left singular vectors.
796 *
797  IF( wantu ) THEN
798  k = itgkz
799  DO i = 1, ns
800  DO j = 1, m
801  u( j, i ) = cmplx( rwork( k ), zero )
802  k = k + 1
803  END DO
804  k = k + m
805  END DO
806 *
807 * Call CUNMBR to compute QB*UB.
808 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
809 *
810  CALL cunmbr( 'Q', 'L', 'N', m, ns, n, a, lda,
811  $ work( itauq ), u, ldu, work( itemp ),
812  $ lwork-itemp+1, info )
813  END IF
814 *
815 * If needed, compute right singular vectors.
816 *
817  IF( wantvt) THEN
818  k = itgkz + m
819  DO i = 1, ns
820  DO j = 1, m
821  vt( i, j ) = cmplx( rwork( k ), zero )
822  k = k + 1
823  END DO
824  k = k + m
825  END DO
826  CALL claset( 'A', ns, n-m, czero, czero,
827  $ vt( 1,m+1 ), ldvt )
828 *
829 * Call CUNMBR to compute VB**T * PB**T
830 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
831 *
832  CALL cunmbr( 'P', 'R', 'C', ns, n, m, a, lda,
833  $ work( itaup ), vt, ldvt, work( itemp ),
834  $ lwork-itemp+1, info )
835  END IF
836  END IF
837  END IF
838 *
839 * Undo scaling if necessary
840 *
841  IF( iscl.EQ.1 ) THEN
842  IF( anrm.GT.bignum )
843  $ CALL slascl( 'G', 0, 0, bignum, anrm, minmn, 1,
844  $ s, minmn, info )
845  IF( anrm.LT.smlnum )
846  $ CALL slascl( 'G', 0, 0, smlnum, anrm, minmn, 1,
847  $ s, minmn, info )
848  END IF
849 *
850 * Return optimal workspace in WORK(1)
851 *
852  work( 1 ) = cmplx( real( maxwrk ), zero )
853 *
854  RETURN
855 *
856 * End of CGESVDX
857 *
858  END
cunmbr
subroutine cunmbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMBR
Definition: cunmbr.f:199
cunmqr
subroutine cunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMQR
Definition: cunmqr.f:170
cunmlq
subroutine cunmlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMLQ
Definition: cunmlq.f:170
clacpy
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:105
slascl
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: slascl.f:145
cgesvdx
subroutine cgesvdx(JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU, IL, IU, NS, S, U, LDU, VT, LDVT, WORK, LWORK, RWORK, IWORK, INFO)
CGESVDX computes the singular value decomposition (SVD) for GE matrices
Definition: cgesvdx.f:272
cgebrd
subroutine cgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
CGEBRD
Definition: cgebrd.f:208
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
cgelqf
subroutine cgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGELQF
Definition: cgelqf.f:145
claset
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: claset.f:108
cgeqrf
subroutine cgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGEQRF
Definition: cgeqrf.f:147
sbdsvdx
subroutine sbdsvdx(UPLO, JOBZ, RANGE, N, D, E, VL, VU, IL, IU, NS, S, Z, LDZ, WORK, IWORK, INFO)
SBDSVDX
Definition: sbdsvdx.f:228
clascl
subroutine clascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: clascl.f:145