LAPACK  3.11.0
LAPACK: Linear Algebra PACKage
cgesvdx.f
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 *> \ingroup gesvdx
265 *
266 * =====================================================================
267  SUBROUTINE cgesvdx( JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU,
268  $ IL, IU, NS, S, U, LDU, VT, LDVT, WORK,
269  $ LWORK, RWORK, IWORK, INFO )
270 *
271 * -- LAPACK driver routine --
272 * -- LAPACK is a software package provided by Univ. of Tennessee, --
273 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
274 *
275 * .. Scalar Arguments ..
276  CHARACTER JOBU, JOBVT, RANGE
277  INTEGER IL, INFO, IU, LDA, LDU, LDVT, LWORK, M, N, NS
278  REAL VL, VU
279 * ..
280 * .. Array Arguments ..
281  INTEGER IWORK( * )
282  REAL S( * ), RWORK( * )
283  COMPLEX A( lda, * ), U( ldu, * ), VT( ldvt, * ),
284  $ work( * )
285 * ..
286 *
287 * =====================================================================
288 *
289 * .. Parameters ..
290  COMPLEX CZERO, CONE
291  parameter( czero = ( 0.0e0, 0.0e0 ),
292  $ cone = ( 1.0e0, 0.0e0 ) )
293  REAL ZERO, ONE
294  parameter( zero = 0.0e0, one = 1.0e0 )
295 * ..
296 * .. Local Scalars ..
297  CHARACTER JOBZ, RNGTGK
298  LOGICAL ALLS, INDS, LQUERY, VALS, WANTU, WANTVT
299  INTEGER I, ID, IE, IERR, ILQF, ILTGK, IQRF, ISCL,
300  $ itau, itaup, itauq, itemp, itempr, itgkz,
301  $ iutgk, j, k, maxwrk, minmn, minwrk, mnthr
302  REAL ABSTOL, ANRM, BIGNUM, EPS, SMLNUM
303 * ..
304 * .. Local Arrays ..
305  REAL DUM( 1 )
306 * ..
307 * .. External Subroutines ..
308  EXTERNAL cgebrd, cgelqf, cgeqrf, clascl, claset,
309  $ cunmbr, cunmqr, cunmlq, clacpy,
310  $ sbdsvdx, slascl, xerbla
311 * ..
312 * .. External Functions ..
313  LOGICAL LSAME
314  INTEGER ILAENV
315  REAL SLAMCH, CLANGE
316  EXTERNAL lsame, ilaenv, slamch, clange
317 * ..
318 * .. Intrinsic Functions ..
319  INTRINSIC max, min, sqrt
320 * ..
321 * .. Executable Statements ..
322 *
323 * Test the input arguments.
324 *
325  ns = 0
326  info = 0
327  abstol = 2*slamch('S')
328  lquery = ( lwork.EQ.-1 )
329  minmn = min( m, n )
330 
331  wantu = lsame( jobu, 'V' )
332  wantvt = lsame( jobvt, 'V' )
333  IF( wantu .OR. wantvt ) THEN
334  jobz = 'V'
335  ELSE
336  jobz = 'N'
337  END IF
338  alls = lsame( range, 'A' )
339  vals = lsame( range, 'V' )
340  inds = lsame( range, 'I' )
341 *
342  info = 0
343  IF( .NOT.lsame( jobu, 'V' ) .AND.
344  $ .NOT.lsame( jobu, 'N' ) ) THEN
345  info = -1
346  ELSE IF( .NOT.lsame( jobvt, 'V' ) .AND.
347  $ .NOT.lsame( jobvt, 'N' ) ) THEN
348  info = -2
349  ELSE IF( .NOT.( alls .OR. vals .OR. inds ) ) THEN
350  info = -3
351  ELSE IF( m.LT.0 ) THEN
352  info = -4
353  ELSE IF( n.LT.0 ) THEN
354  info = -5
355  ELSE IF( m.GT.lda ) THEN
356  info = -7
357  ELSE IF( minmn.GT.0 ) THEN
358  IF( vals ) THEN
359  IF( vl.LT.zero ) THEN
360  info = -8
361  ELSE IF( vu.LE.vl ) THEN
362  info = -9
363  END IF
364  ELSE IF( inds ) THEN
365  IF( il.LT.1 .OR. il.GT.max( 1, minmn ) ) THEN
366  info = -10
367  ELSE IF( iu.LT.min( minmn, il ) .OR. iu.GT.minmn ) THEN
368  info = -11
369  END IF
370  END IF
371  IF( info.EQ.0 ) THEN
372  IF( wantu .AND. ldu.LT.m ) THEN
373  info = -15
374  ELSE IF( wantvt ) THEN
375  IF( inds ) THEN
376  IF( ldvt.LT.iu-il+1 ) THEN
377  info = -17
378  END IF
379  ELSE IF( ldvt.LT.minmn ) THEN
380  info = -17
381  END IF
382  END IF
383  END IF
384  END IF
385 *
386 * Compute workspace
387 * (Note: Comments in the code beginning "Workspace:" describe the
388 * minimal amount of workspace needed at that point in the code,
389 * as well as the preferred amount for good performance.
390 * NB refers to the optimal block size for the immediately
391 * following subroutine, as returned by ILAENV.)
392 *
393  IF( info.EQ.0 ) THEN
394  minwrk = 1
395  maxwrk = 1
396  IF( minmn.GT.0 ) THEN
397  IF( m.GE.n ) THEN
398  mnthr = ilaenv( 6, 'CGESVD', jobu // jobvt, m, n, 0, 0 )
399  IF( m.GE.mnthr ) THEN
400 *
401 * Path 1 (M much larger than N)
402 *
403  minwrk = n*(n+5)
404  maxwrk = n + n*ilaenv(1,'CGEQRF',' ',m,n,-1,-1)
405  maxwrk = max(maxwrk,
406  $ n*n+2*n+2*n*ilaenv(1,'CGEBRD',' ',n,n,-1,-1))
407  IF (wantu .OR. wantvt) THEN
408  maxwrk = max(maxwrk,
409  $ n*n+2*n+n*ilaenv(1,'CUNMQR','LN',n,n,n,-1))
410  END IF
411  ELSE
412 *
413 * Path 2 (M at least N, but not much larger)
414 *
415  minwrk = 3*n + m
416  maxwrk = 2*n + (m+n)*ilaenv(1,'CGEBRD',' ',m,n,-1,-1)
417  IF (wantu .OR. wantvt) THEN
418  maxwrk = max(maxwrk,
419  $ 2*n+n*ilaenv(1,'CUNMQR','LN',n,n,n,-1))
420  END IF
421  END IF
422  ELSE
423  mnthr = ilaenv( 6, 'CGESVD', jobu // jobvt, m, n, 0, 0 )
424  IF( n.GE.mnthr ) THEN
425 *
426 * Path 1t (N much larger than M)
427 *
428  minwrk = m*(m+5)
429  maxwrk = m + m*ilaenv(1,'CGELQF',' ',m,n,-1,-1)
430  maxwrk = max(maxwrk,
431  $ m*m+2*m+2*m*ilaenv(1,'CGEBRD',' ',m,m,-1,-1))
432  IF (wantu .OR. wantvt) THEN
433  maxwrk = max(maxwrk,
434  $ m*m+2*m+m*ilaenv(1,'CUNMQR','LN',m,m,m,-1))
435  END IF
436  ELSE
437 *
438 * Path 2t (N greater than M, but not much larger)
439 *
440 *
441  minwrk = 3*m + n
442  maxwrk = 2*m + (m+n)*ilaenv(1,'CGEBRD',' ',m,n,-1,-1)
443  IF (wantu .OR. wantvt) THEN
444  maxwrk = max(maxwrk,
445  $ 2*m+m*ilaenv(1,'CUNMQR','LN',m,m,m,-1))
446  END IF
447  END IF
448  END IF
449  END IF
450  maxwrk = max( maxwrk, minwrk )
451  work( 1 ) = cmplx( REAL( MAXWRK ), ZERO )
452 *
453  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
454  info = -19
455  END IF
456  END IF
457 *
458  IF( info.NE.0 ) THEN
459  CALL xerbla( 'CGESVDX', -info )
460  RETURN
461  ELSE IF( lquery ) THEN
462  RETURN
463  END IF
464 *
465 * Quick return if possible
466 *
467  IF( m.EQ.0 .OR. n.EQ.0 ) THEN
468  RETURN
469  END IF
470 *
471 * Set singular values indices accord to RANGE='A'.
472 *
473  IF( alls ) THEN
474  rngtgk = 'I'
475  iltgk = 1
476  iutgk = min( m, n )
477  ELSE IF( inds ) THEN
478  rngtgk = 'I'
479  iltgk = il
480  iutgk = iu
481  ELSE
482  rngtgk = 'V'
483  iltgk = 0
484  iutgk = 0
485  END IF
486 *
487 * Get machine constants
488 *
489  eps = slamch( 'P' )
490  smlnum = sqrt( slamch( 'S' ) ) / eps
491  bignum = one / smlnum
492 *
493 * Scale A if max element outside range [SMLNUM,BIGNUM]
494 *
495  anrm = clange( 'M', m, n, a, lda, dum )
496  iscl = 0
497  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
498  iscl = 1
499  CALL clascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
500  ELSE IF( anrm.GT.bignum ) THEN
501  iscl = 1
502  CALL clascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
503  END IF
504 *
505  IF( m.GE.n ) THEN
506 *
507 * A has at least as many rows as columns. If A has sufficiently
508 * more rows than columns, first reduce A using the QR
509 * decomposition.
510 *
511  IF( m.GE.mnthr ) THEN
512 *
513 * Path 1 (M much larger than N):
514 * A = Q * R = Q * ( QB * B * PB**T )
515 * = Q * ( QB * ( UB * S * VB**T ) * PB**T )
516 * U = Q * QB * UB; V**T = VB**T * PB**T
517 *
518 * Compute A=Q*R
519 * (Workspace: need 2*N, prefer N+N*NB)
520 *
521  itau = 1
522  itemp = itau + n
523  CALL cgeqrf( m, n, a, lda, work( itau ), work( itemp ),
524  $ lwork-itemp+1, info )
525 *
526 * Copy R into WORK and bidiagonalize it:
527 * (Workspace: need N*N+3*N, prefer N*N+N+2*N*NB)
528 *
529  iqrf = itemp
530  itauq = itemp + n*n
531  itaup = itauq + n
532  itemp = itaup + n
533  id = 1
534  ie = id + n
535  itgkz = ie + n
536  CALL clacpy( 'U', n, n, a, lda, work( iqrf ), n )
537  CALL claset( 'L', n-1, n-1, czero, czero,
538  $ work( iqrf+1 ), n )
539  CALL cgebrd( n, n, work( iqrf ), n, rwork( id ),
540  $ rwork( ie ), work( itauq ), work( itaup ),
541  $ work( itemp ), lwork-itemp+1, info )
542  itempr = itgkz + n*(n*2+1)
543 *
544 * Solve eigenvalue problem TGK*Z=Z*S.
545 * (Workspace: need 2*N*N+14*N)
546 *
547  CALL sbdsvdx( 'U', jobz, rngtgk, n, rwork( id ),
548  $ rwork( ie ), vl, vu, iltgk, iutgk, ns, s,
549  $ rwork( itgkz ), n*2, rwork( itempr ),
550  $ iwork, info)
551 *
552 * If needed, compute left singular vectors.
553 *
554  IF( wantu ) THEN
555  k = itgkz
556  DO i = 1, ns
557  DO j = 1, n
558  u( j, i ) = cmplx( rwork( k ), zero )
559  k = k + 1
560  END DO
561  k = k + n
562  END DO
563  CALL claset( 'A', m-n, ns, czero, czero, u( n+1,1 ), ldu)
564 *
565 * Call CUNMBR to compute QB*UB.
566 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
567 *
568  CALL cunmbr( 'Q', 'L', 'N', n, ns, n, work( iqrf ), n,
569  $ work( itauq ), u, ldu, work( itemp ),
570  $ lwork-itemp+1, info )
571 *
572 * Call CUNMQR to compute Q*(QB*UB).
573 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
574 *
575  CALL cunmqr( 'L', 'N', m, ns, n, a, lda,
576  $ work( itau ), u, ldu, work( itemp ),
577  $ lwork-itemp+1, info )
578  END IF
579 *
580 * If needed, compute right singular vectors.
581 *
582  IF( wantvt) THEN
583  k = itgkz + n
584  DO i = 1, ns
585  DO j = 1, n
586  vt( i, j ) = cmplx( rwork( k ), zero )
587  k = k + 1
588  END DO
589  k = k + n
590  END DO
591 *
592 * Call CUNMBR to compute VB**T * PB**T
593 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
594 *
595  CALL cunmbr( 'P', 'R', 'C', ns, n, n, work( iqrf ), n,
596  $ work( itaup ), vt, ldvt, work( itemp ),
597  $ lwork-itemp+1, info )
598  END IF
599  ELSE
600 *
601 * Path 2 (M at least N, but not much larger)
602 * Reduce A to bidiagonal form without QR decomposition
603 * A = QB * B * PB**T = QB * ( UB * S * VB**T ) * PB**T
604 * U = QB * UB; V**T = VB**T * PB**T
605 *
606 * Bidiagonalize A
607 * (Workspace: need 2*N+M, prefer 2*N+(M+N)*NB)
608 *
609  itauq = 1
610  itaup = itauq + n
611  itemp = itaup + n
612  id = 1
613  ie = id + n
614  itgkz = ie + n
615  CALL cgebrd( m, n, a, lda, rwork( id ), rwork( ie ),
616  $ work( itauq ), work( itaup ), work( itemp ),
617  $ lwork-itemp+1, info )
618  itempr = itgkz + n*(n*2+1)
619 *
620 * Solve eigenvalue problem TGK*Z=Z*S.
621 * (Workspace: need 2*N*N+14*N)
622 *
623  CALL sbdsvdx( 'U', jobz, rngtgk, n, rwork( id ),
624  $ rwork( ie ), vl, vu, iltgk, iutgk, ns, s,
625  $ rwork( itgkz ), n*2, rwork( itempr ),
626  $ iwork, info)
627 *
628 * If needed, compute left singular vectors.
629 *
630  IF( wantu ) THEN
631  k = itgkz
632  DO i = 1, ns
633  DO j = 1, n
634  u( j, i ) = cmplx( rwork( k ), zero )
635  k = k + 1
636  END DO
637  k = k + n
638  END DO
639  CALL claset( 'A', m-n, ns, czero, czero, u( n+1,1 ), ldu)
640 *
641 * Call CUNMBR to compute QB*UB.
642 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
643 *
644  CALL cunmbr( 'Q', 'L', 'N', m, ns, n, a, lda,
645  $ work( itauq ), u, ldu, work( itemp ),
646  $ lwork-itemp+1, ierr )
647  END IF
648 *
649 * If needed, compute right singular vectors.
650 *
651  IF( wantvt) THEN
652  k = itgkz + n
653  DO i = 1, ns
654  DO j = 1, n
655  vt( i, j ) = cmplx( rwork( k ), zero )
656  k = k + 1
657  END DO
658  k = k + n
659  END DO
660 *
661 * Call CUNMBR to compute VB**T * PB**T
662 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
663 *
664  CALL cunmbr( 'P', 'R', 'C', ns, n, n, a, lda,
665  $ work( itaup ), vt, ldvt, work( itemp ),
666  $ lwork-itemp+1, ierr )
667  END IF
668  END IF
669  ELSE
670 *
671 * A has more columns than rows. If A has sufficiently more
672 * columns than rows, first reduce A using the LQ decomposition.
673 *
674  IF( n.GE.mnthr ) THEN
675 *
676 * Path 1t (N much larger than M):
677 * A = L * Q = ( QB * B * PB**T ) * Q
678 * = ( QB * ( UB * S * VB**T ) * PB**T ) * Q
679 * U = QB * UB ; V**T = VB**T * PB**T * Q
680 *
681 * Compute A=L*Q
682 * (Workspace: need 2*M, prefer M+M*NB)
683 *
684  itau = 1
685  itemp = itau + m
686  CALL cgelqf( m, n, a, lda, work( itau ), work( itemp ),
687  $ lwork-itemp+1, info )
688 
689 * Copy L into WORK and bidiagonalize it:
690 * (Workspace in WORK( ITEMP ): need M*M+3*M, prefer M*M+M+2*M*NB)
691 *
692  ilqf = itemp
693  itauq = ilqf + m*m
694  itaup = itauq + m
695  itemp = itaup + m
696  id = 1
697  ie = id + m
698  itgkz = ie + m
699  CALL clacpy( 'L', m, m, a, lda, work( ilqf ), m )
700  CALL claset( 'U', m-1, m-1, czero, czero,
701  $ work( ilqf+m ), m )
702  CALL cgebrd( m, m, work( ilqf ), m, rwork( id ),
703  $ rwork( ie ), work( itauq ), work( itaup ),
704  $ work( itemp ), lwork-itemp+1, info )
705  itempr = itgkz + m*(m*2+1)
706 *
707 * Solve eigenvalue problem TGK*Z=Z*S.
708 * (Workspace: need 2*M*M+14*M)
709 *
710  CALL sbdsvdx( 'U', jobz, rngtgk, m, rwork( id ),
711  $ rwork( ie ), vl, vu, iltgk, iutgk, ns, s,
712  $ rwork( itgkz ), m*2, rwork( itempr ),
713  $ iwork, info)
714 *
715 * If needed, compute left singular vectors.
716 *
717  IF( wantu ) THEN
718  k = itgkz
719  DO i = 1, ns
720  DO j = 1, m
721  u( j, i ) = cmplx( rwork( k ), zero )
722  k = k + 1
723  END DO
724  k = k + m
725  END DO
726 *
727 * Call CUNMBR to compute QB*UB.
728 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
729 *
730  CALL cunmbr( 'Q', 'L', 'N', m, ns, m, work( ilqf ), m,
731  $ work( itauq ), u, ldu, work( itemp ),
732  $ lwork-itemp+1, info )
733  END IF
734 *
735 * If needed, compute right singular vectors.
736 *
737  IF( wantvt) THEN
738  k = itgkz + m
739  DO i = 1, ns
740  DO j = 1, m
741  vt( i, j ) = cmplx( rwork( k ), zero )
742  k = k + 1
743  END DO
744  k = k + m
745  END DO
746  CALL claset( 'A', ns, n-m, czero, czero,
747  $ vt( 1,m+1 ), ldvt )
748 *
749 * Call CUNMBR to compute (VB**T)*(PB**T)
750 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
751 *
752  CALL cunmbr( 'P', 'R', 'C', ns, m, m, work( ilqf ), m,
753  $ work( itaup ), vt, ldvt, work( itemp ),
754  $ lwork-itemp+1, info )
755 *
756 * Call CUNMLQ to compute ((VB**T)*(PB**T))*Q.
757 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
758 *
759  CALL cunmlq( 'R', 'N', ns, n, m, a, lda,
760  $ work( itau ), vt, ldvt, work( itemp ),
761  $ lwork-itemp+1, info )
762  END IF
763  ELSE
764 *
765 * Path 2t (N greater than M, but not much larger)
766 * Reduce to bidiagonal form without LQ decomposition
767 * A = QB * B * PB**T = QB * ( UB * S * VB**T ) * PB**T
768 * U = QB * UB; V**T = VB**T * PB**T
769 *
770 * Bidiagonalize A
771 * (Workspace: need 2*M+N, prefer 2*M+(M+N)*NB)
772 *
773  itauq = 1
774  itaup = itauq + m
775  itemp = itaup + m
776  id = 1
777  ie = id + m
778  itgkz = ie + m
779  CALL cgebrd( m, n, a, lda, rwork( id ), rwork( ie ),
780  $ work( itauq ), work( itaup ), work( itemp ),
781  $ lwork-itemp+1, info )
782  itempr = itgkz + m*(m*2+1)
783 *
784 * Solve eigenvalue problem TGK*Z=Z*S.
785 * (Workspace: need 2*M*M+14*M)
786 *
787  CALL sbdsvdx( 'L', jobz, rngtgk, m, rwork( id ),
788  $ rwork( ie ), vl, vu, iltgk, iutgk, ns, s,
789  $ rwork( itgkz ), m*2, rwork( itempr ),
790  $ iwork, info)
791 *
792 * If needed, compute left singular vectors.
793 *
794  IF( wantu ) THEN
795  k = itgkz
796  DO i = 1, ns
797  DO j = 1, m
798  u( j, i ) = cmplx( rwork( k ), zero )
799  k = k + 1
800  END DO
801  k = k + m
802  END DO
803 *
804 * Call CUNMBR to compute QB*UB.
805 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
806 *
807  CALL cunmbr( 'Q', 'L', 'N', m, ns, n, a, lda,
808  $ work( itauq ), u, ldu, work( itemp ),
809  $ lwork-itemp+1, info )
810  END IF
811 *
812 * If needed, compute right singular vectors.
813 *
814  IF( wantvt) THEN
815  k = itgkz + m
816  DO i = 1, ns
817  DO j = 1, m
818  vt( i, j ) = cmplx( rwork( k ), zero )
819  k = k + 1
820  END DO
821  k = k + m
822  END DO
823  CALL claset( 'A', ns, n-m, czero, czero,
824  $ vt( 1,m+1 ), ldvt )
825 *
826 * Call CUNMBR to compute VB**T * PB**T
827 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
828 *
829  CALL cunmbr( 'P', 'R', 'C', ns, n, m, a, lda,
830  $ work( itaup ), vt, ldvt, work( itemp ),
831  $ lwork-itemp+1, info )
832  END IF
833  END IF
834  END IF
835 *
836 * Undo scaling if necessary
837 *
838  IF( iscl.EQ.1 ) THEN
839  IF( anrm.GT.bignum )
840  $ CALL slascl( 'G', 0, 0, bignum, anrm, minmn, 1,
841  $ s, minmn, info )
842  IF( anrm.LT.smlnum )
843  $ CALL slascl( 'G', 0, 0, smlnum, anrm, minmn, 1,
844  $ s, minmn, info )
845  END IF
846 *
847 * Return optimal workspace in WORK(1)
848 *
849  work( 1 ) = cmplx( REAL( MAXWRK ), ZERO )
850 *
851  RETURN
852 *
853 * End of CGESVDX
854 *
855  END
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:103
sbdsvdx
subroutine sbdsvdx(UPLO, JOBZ, RANGE, N, D, E, VL, VU, IL, IU, NS, S, Z, LDZ, WORK, IWORK, INFO)
SBDSVDX
Definition: sbdsvdx.f:226
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
cunmbr
subroutine cunmbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMBR
Definition: cunmbr.f:197
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:270
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:106
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:143
cgelqf
subroutine cgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGELQF
Definition: cgelqf.f:143
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:143
cunmqr
subroutine cunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMQR
Definition: cunmqr.f:168
cgebrd
subroutine cgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
CGEBRD
Definition: cgebrd.f:206
cgeqrf
subroutine cgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGEQRF
Definition: cgeqrf.f:146
cunmlq
subroutine cunmlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMLQ
Definition: cunmlq.f:168