LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
zgesvdx.f
Go to the documentation of this file.
1 *> \brief <b> ZGESVDX 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 ZGESVDX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgesvdx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgesvdx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgesvdx.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 * DOUBLE PRECISION VL, VU
30 * ..
31 * .. Array Arguments ..
32 * INTEGER IWORK( * )
33 * DOUBLE PRECISION S( * ), RWORK( * )
34 * COMPLEX*16 A( LDA, * ), U( LDU, * ), VT( LDVT, * ),
35 * $ WORK( * )
36 * ..
37 *
38 *
39 *> \par Purpose:
40 * =============
41 *>
42 *> \verbatim
43 *>
44 *> ZGESVDX 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 *> ZGESVDX uses an eigenvalue problem for obtaining the SVD, which
58 *> allows for the computation of a subset of singular values and
59 *> vectors. See DBDSVDX 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*16 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 DOUBLE PRECISION
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 DOUBLE PRECISION
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 DOUBLE PRECISION 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*16 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*16 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*16 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 DOUBLE PRECISION 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 DBDSVDX/DSTEVX.
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 DBDSVDX/DSTEVX.
252 *> if INFO = N*2 + 1, an internal error occurred in
253 *> DBDSVDX
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 complex16GEsing
267 *
268 * =====================================================================
269  SUBROUTINE zgesvdx( 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.6.1) --
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  DOUBLE PRECISION VL, VU
282 * ..
283 * .. Array Arguments ..
284  INTEGER IWORK( * )
285  DOUBLE PRECISION S( * ), RWORK( * )
286  COMPLEX*16 A( lda, * ), U( ldu, * ), VT( ldvt, * ),
287  $ work( * )
288 * ..
289 *
290 * =====================================================================
291 *
292 * .. Parameters ..
293  COMPLEX*16 CZERO, CONE
294  parameter ( czero = ( 0.0d0, 0.0d0 ),
295  $ cone = ( 1.0d0, 0.0d0 ) )
296  DOUBLE PRECISION ZERO, ONE
297  parameter ( zero = 0.0d0, one = 1.0d0 )
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  DOUBLE PRECISION ABSTOL, ANRM, BIGNUM, EPS, SMLNUM
306 * ..
307 * .. Local Arrays ..
308  DOUBLE PRECISION DUM( 1 )
309 * ..
310 * .. External Subroutines ..
311  EXTERNAL zgebrd, zgelqf, zgeqrf, zlascl, zlaset,
312  $ dlascl, xerbla
313 * ..
314 * .. External Functions ..
315  LOGICAL LSAME
316  INTEGER ILAENV
317  DOUBLE PRECISION DLAMCH, ZLANGE
318  EXTERNAL lsame, ilaenv, dlamch, zlange
319 * ..
320 * .. Intrinsic Functions ..
321  INTRINSIC max, min, sqrt
322 * ..
323 * .. Executable Statements ..
324 *
325 * Test the input arguments.
326 *
327  ns = 0
328  info = 0
329  abstol = 2*dlamch('S')
330  lquery = ( lwork.EQ.-1 )
331  minmn = min( m, n )
332 
333  wantu = lsame( jobu, 'V' )
334  wantvt = lsame( jobvt, 'V' )
335  IF( wantu .OR. wantvt ) THEN
336  jobz = 'V'
337  ELSE
338  jobz = 'N'
339  END IF
340  alls = lsame( range, 'A' )
341  vals = lsame( range, 'V' )
342  inds = lsame( range, 'I' )
343 *
344  info = 0
345  IF( .NOT.lsame( jobu, 'V' ) .AND.
346  $ .NOT.lsame( jobu, 'N' ) ) THEN
347  info = -1
348  ELSE IF( .NOT.lsame( jobvt, 'V' ) .AND.
349  $ .NOT.lsame( jobvt, 'N' ) ) THEN
350  info = -2
351  ELSE IF( .NOT.( alls .OR. vals .OR. inds ) ) THEN
352  info = -3
353  ELSE IF( m.LT.0 ) THEN
354  info = -4
355  ELSE IF( n.LT.0 ) THEN
356  info = -5
357  ELSE IF( m.GT.lda ) THEN
358  info = -7
359  ELSE IF( minmn.GT.0 ) THEN
360  IF( vals ) THEN
361  IF( vl.LT.zero ) THEN
362  info = -8
363  ELSE IF( vu.LE.vl ) THEN
364  info = -9
365  END IF
366  ELSE IF( inds ) THEN
367  IF( il.LT.1 .OR. il.GT.max( 1, minmn ) ) THEN
368  info = -10
369  ELSE IF( iu.LT.min( minmn, il ) .OR. iu.GT.minmn ) THEN
370  info = -11
371  END IF
372  END IF
373  IF( info.EQ.0 ) THEN
374  IF( wantu .AND. ldu.LT.m ) THEN
375  info = -15
376  ELSE IF( wantvt ) THEN
377  IF( inds ) THEN
378  IF( ldvt.LT.iu-il+1 ) THEN
379  info = -17
380  END IF
381  ELSE IF( ldvt.LT.minmn ) THEN
382  info = -17
383  END IF
384  END IF
385  END IF
386  END IF
387 *
388 * Compute workspace
389 * (Note: Comments in the code beginning "Workspace:" describe the
390 * minimal amount of workspace needed at that point in the code,
391 * as well as the preferred amount for good performance.
392 * NB refers to the optimal block size for the immediately
393 * following subroutine, as returned by ILAENV.)
394 *
395  IF( info.EQ.0 ) THEN
396  minwrk = 1
397  maxwrk = 1
398  IF( minmn.GT.0 ) THEN
399  IF( m.GE.n ) THEN
400  mnthr = ilaenv( 6, 'ZGESVD', jobu // jobvt, m, n, 0, 0 )
401  IF( m.GE.mnthr ) THEN
402 *
403 * Path 1 (M much larger than N)
404 *
405  minwrk = n*(n+5)
406  maxwrk = n + n*ilaenv(1,'ZGEQRF',' ',m,n,-1,-1)
407  maxwrk = max(maxwrk,
408  $ n*n+2*n+2*n*ilaenv(1,'ZGEBRD',' ',n,n,-1,-1))
409  IF (wantu .OR. wantvt) THEN
410  maxwrk = max(maxwrk,
411  $ n*n+2*n+n*ilaenv(1,'ZUNMQR','LN',n,n,n,-1))
412  END IF
413  ELSE
414 *
415 * Path 2 (M at least N, but not much larger)
416 *
417  minwrk = 3*n + m
418  maxwrk = 2*n + (m+n)*ilaenv(1,'ZGEBRD',' ',m,n,-1,-1)
419  IF (wantu .OR. wantvt) THEN
420  maxwrk = max(maxwrk,
421  $ 2*n+n*ilaenv(1,'ZUNMQR','LN',n,n,n,-1))
422  END IF
423  END IF
424  ELSE
425  mnthr = ilaenv( 6, 'ZGESVD', jobu // jobvt, m, n, 0, 0 )
426  IF( n.GE.mnthr ) THEN
427 *
428 * Path 1t (N much larger than M)
429 *
430  minwrk = m*(m+5)
431  maxwrk = m + m*ilaenv(1,'ZGELQF',' ',m,n,-1,-1)
432  maxwrk = max(maxwrk,
433  $ m*m+2*m+2*m*ilaenv(1,'ZGEBRD',' ',m,m,-1,-1))
434  IF (wantu .OR. wantvt) THEN
435  maxwrk = max(maxwrk,
436  $ m*m+2*m+m*ilaenv(1,'ZUNMQR','LN',m,m,m,-1))
437  END IF
438  ELSE
439 *
440 * Path 2t (N greater than M, but not much larger)
441 *
442 *
443  minwrk = 3*m + n
444  maxwrk = 2*m + (m+n)*ilaenv(1,'ZGEBRD',' ',m,n,-1,-1)
445  IF (wantu .OR. wantvt) THEN
446  maxwrk = max(maxwrk,
447  $ 2*m+m*ilaenv(1,'ZUNMQR','LN',m,m,m,-1))
448  END IF
449  END IF
450  END IF
451  END IF
452  maxwrk = max( maxwrk, minwrk )
453  work( 1 ) = dcmplx( dble( maxwrk ), zero )
454 *
455  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
456  info = -19
457  END IF
458  END IF
459 *
460  IF( info.NE.0 ) THEN
461  CALL xerbla( 'ZGESVDX', -info )
462  RETURN
463  ELSE IF( lquery ) THEN
464  RETURN
465  END IF
466 *
467 * Quick return if possible
468 *
469  IF( m.EQ.0 .OR. n.EQ.0 ) THEN
470  RETURN
471  END IF
472 *
473 * Set singular values indices accord to RANGE='A'.
474 *
475  IF( alls ) THEN
476  rngtgk = 'I'
477  iltgk = 1
478  iutgk = min( m, n )
479  ELSE IF( inds ) THEN
480  rngtgk = 'I'
481  iltgk = il
482  iutgk = iu
483  ELSE
484  rngtgk = 'V'
485  iltgk = 0
486  iutgk = 0
487  END IF
488 *
489 * Get machine constants
490 *
491  eps = dlamch( 'P' )
492  smlnum = sqrt( dlamch( 'S' ) ) / eps
493  bignum = one / smlnum
494 *
495 * Scale A if max element outside range [SMLNUM,BIGNUM]
496 *
497  anrm = zlange( 'M', m, n, a, lda, dum )
498  iscl = 0
499  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
500  iscl = 1
501  CALL zlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
502  ELSE IF( anrm.GT.bignum ) THEN
503  iscl = 1
504  CALL zlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
505  END IF
506 *
507  IF( m.GE.n ) THEN
508 *
509 * A has at least as many rows as columns. If A has sufficiently
510 * more rows than columns, first reduce A using the QR
511 * decomposition.
512 *
513  IF( m.GE.mnthr ) THEN
514 *
515 * Path 1 (M much larger than N):
516 * A = Q * R = Q * ( QB * B * PB**T )
517 * = Q * ( QB * ( UB * S * VB**T ) * PB**T )
518 * U = Q * QB * UB; V**T = VB**T * PB**T
519 *
520 * Compute A=Q*R
521 * (Workspace: need 2*N, prefer N+N*NB)
522 *
523  itau = 1
524  itemp = itau + n
525  CALL zgeqrf( m, n, a, lda, work( itau ), work( itemp ),
526  $ lwork-itemp+1, info )
527 *
528 * Copy R into WORK and bidiagonalize it:
529 * (Workspace: need N*N+3*N, prefer N*N+N+2*N*NB)
530 *
531  iqrf = itemp
532  itauq = itemp + n*n
533  itaup = itauq + n
534  itemp = itaup + n
535  id = 1
536  ie = id + n
537  itgkz = ie + n
538  CALL zlacpy( 'U', n, n, a, lda, work( iqrf ), n )
539  CALL zlaset( 'L', n-1, n-1, czero, czero,
540  $ work( iqrf+1 ), n )
541  CALL zgebrd( n, n, work( iqrf ), n, rwork( id ),
542  $ rwork( ie ), work( itauq ), work( itaup ),
543  $ work( itemp ), lwork-itemp+1, info )
544  itempr = itgkz + n*(n*2+1)
545 *
546 * Solve eigenvalue problem TGK*Z=Z*S.
547 * (Workspace: need 2*N*N+14*N)
548 *
549  CALL dbdsvdx( 'U', jobz, rngtgk, n, rwork( id ),
550  $ rwork( ie ), vl, vu, iltgk, iutgk, ns, s,
551  $ rwork( itgkz ), n*2, rwork( itempr ),
552  $ iwork, info)
553 *
554 * If needed, compute left singular vectors.
555 *
556  IF( wantu ) THEN
557  k = itgkz
558  DO i = 1, ns
559  DO j = 1, n
560  u( j, i ) = dcmplx( rwork( k ), zero )
561  k = k + 1
562  END DO
563  k = k + n
564  END DO
565  CALL zlaset( 'A', m-n, ns, czero, czero, u( n+1,1 ), ldu)
566 *
567 * Call ZUNMBR to compute QB*UB.
568 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
569 *
570  CALL zunmbr( 'Q', 'L', 'N', n, ns, n, work( iqrf ), n,
571  $ work( itauq ), u, ldu, work( itemp ),
572  $ lwork-itemp+1, info )
573 *
574 * Call ZUNMQR to compute Q*(QB*UB).
575 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
576 *
577  CALL zunmqr( 'L', 'N', m, ns, n, a, lda,
578  $ work( itau ), u, ldu, work( itemp ),
579  $ lwork-itemp+1, info )
580  END IF
581 *
582 * If needed, compute right singular vectors.
583 *
584  IF( wantvt) THEN
585  k = itgkz + n
586  DO i = 1, ns
587  DO j = 1, n
588  vt( i, j ) = dcmplx( rwork( k ), zero )
589  k = k + 1
590  END DO
591  k = k + n
592  END DO
593 *
594 * Call ZUNMBR to compute VB**T * PB**T
595 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
596 *
597  CALL zunmbr( 'P', 'R', 'C', ns, n, n, work( iqrf ), n,
598  $ work( itaup ), vt, ldvt, work( itemp ),
599  $ lwork-itemp+1, info )
600  END IF
601  ELSE
602 *
603 * Path 2 (M at least N, but not much larger)
604 * Reduce A to bidiagonal form without QR decomposition
605 * A = QB * B * PB**T = QB * ( UB * S * VB**T ) * PB**T
606 * U = QB * UB; V**T = VB**T * PB**T
607 *
608 * Bidiagonalize A
609 * (Workspace: need 2*N+M, prefer 2*N+(M+N)*NB)
610 *
611  itauq = 1
612  itaup = itauq + n
613  itemp = itaup + n
614  id = 1
615  ie = id + n
616  itgkz = ie + n
617  CALL zgebrd( m, n, a, lda, rwork( id ), rwork( ie ),
618  $ work( itauq ), work( itaup ), work( itemp ),
619  $ lwork-itemp+1, info )
620  itempr = itgkz + n*(n*2+1)
621 *
622 * Solve eigenvalue problem TGK*Z=Z*S.
623 * (Workspace: need 2*N*N+14*N)
624 *
625  CALL dbdsvdx( 'U', jobz, rngtgk, n, rwork( id ),
626  $ rwork( ie ), vl, vu, iltgk, iutgk, ns, s,
627  $ rwork( itgkz ), n*2, rwork( itempr ),
628  $ iwork, info)
629 *
630 * If needed, compute left singular vectors.
631 *
632  IF( wantu ) THEN
633  k = itgkz
634  DO i = 1, ns
635  DO j = 1, n
636  u( j, i ) = dcmplx( rwork( k ), zero )
637  k = k + 1
638  END DO
639  k = k + n
640  END DO
641  CALL zlaset( 'A', m-n, ns, czero, czero, u( n+1,1 ), ldu)
642 *
643 * Call ZUNMBR to compute QB*UB.
644 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
645 *
646  CALL zunmbr( 'Q', 'L', 'N', m, ns, n, a, lda,
647  $ work( itauq ), u, ldu, work( itemp ),
648  $ lwork-itemp+1, ierr )
649  END IF
650 *
651 * If needed, compute right singular vectors.
652 *
653  IF( wantvt) THEN
654  k = itgkz + n
655  DO i = 1, ns
656  DO j = 1, n
657  vt( i, j ) = dcmplx( rwork( k ), zero )
658  k = k + 1
659  END DO
660  k = k + n
661  END DO
662 *
663 * Call ZUNMBR to compute VB**T * PB**T
664 * (Workspace in WORK( ITEMP ): need N, prefer N*NB)
665 *
666  CALL zunmbr( 'P', 'R', 'C', ns, n, n, a, lda,
667  $ work( itaup ), vt, ldvt, work( itemp ),
668  $ lwork-itemp+1, ierr )
669  END IF
670  END IF
671  ELSE
672 *
673 * A has more columns than rows. If A has sufficiently more
674 * columns than rows, first reduce A using the LQ decomposition.
675 *
676  IF( n.GE.mnthr ) THEN
677 *
678 * Path 1t (N much larger than M):
679 * A = L * Q = ( QB * B * PB**T ) * Q
680 * = ( QB * ( UB * S * VB**T ) * PB**T ) * Q
681 * U = QB * UB ; V**T = VB**T * PB**T * Q
682 *
683 * Compute A=L*Q
684 * (Workspace: need 2*M, prefer M+M*NB)
685 *
686  itau = 1
687  itemp = itau + m
688  CALL zgelqf( m, n, a, lda, work( itau ), work( itemp ),
689  $ lwork-itemp+1, info )
690 
691 * Copy L into WORK and bidiagonalize it:
692 * (Workspace in WORK( ITEMP ): need M*M+3*M, prefer M*M+M+2*M*NB)
693 *
694  ilqf = itemp
695  itauq = ilqf + m*m
696  itaup = itauq + m
697  itemp = itaup + m
698  id = 1
699  ie = id + m
700  itgkz = ie + m
701  CALL zlacpy( 'L', m, m, a, lda, work( ilqf ), m )
702  CALL zlaset( 'U', m-1, m-1, czero, czero,
703  $ work( ilqf+m ), m )
704  CALL zgebrd( m, m, work( ilqf ), m, rwork( id ),
705  $ rwork( ie ), work( itauq ), work( itaup ),
706  $ work( itemp ), lwork-itemp+1, info )
707  itempr = itgkz + m*(m*2+1)
708 *
709 * Solve eigenvalue problem TGK*Z=Z*S.
710 * (Workspace: need 2*M*M+14*M)
711 *
712  CALL dbdsvdx( 'U', jobz, rngtgk, m, rwork( id ),
713  $ rwork( ie ), vl, vu, iltgk, iutgk, ns, s,
714  $ rwork( itgkz ), m*2, rwork( itempr ),
715  $ iwork, info)
716 *
717 * If needed, compute left singular vectors.
718 *
719  IF( wantu ) THEN
720  k = itgkz
721  DO i = 1, ns
722  DO j = 1, m
723  u( j, i ) = dcmplx( rwork( k ), zero )
724  k = k + 1
725  END DO
726  k = k + m
727  END DO
728 *
729 * Call ZUNMBR to compute QB*UB.
730 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
731 *
732  CALL zunmbr( 'Q', 'L', 'N', m, ns, m, work( ilqf ), m,
733  $ work( itauq ), u, ldu, work( itemp ),
734  $ lwork-itemp+1, info )
735  END IF
736 *
737 * If needed, compute right singular vectors.
738 *
739  IF( wantvt) THEN
740  k = itgkz + m
741  DO i = 1, ns
742  DO j = 1, m
743  vt( i, j ) = dcmplx( rwork( k ), zero )
744  k = k + 1
745  END DO
746  k = k + m
747  END DO
748  CALL zlaset( 'A', ns, n-m, czero, czero,
749  $ vt( 1,m+1 ), ldvt )
750 *
751 * Call ZUNMBR to compute (VB**T)*(PB**T)
752 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
753 *
754  CALL zunmbr( 'P', 'R', 'C', ns, m, m, work( ilqf ), m,
755  $ work( itaup ), vt, ldvt, work( itemp ),
756  $ lwork-itemp+1, info )
757 *
758 * Call ZUNMLQ to compute ((VB**T)*(PB**T))*Q.
759 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
760 *
761  CALL zunmlq( 'R', 'N', ns, n, m, a, lda,
762  $ work( itau ), vt, ldvt, work( itemp ),
763  $ lwork-itemp+1, info )
764  END IF
765  ELSE
766 *
767 * Path 2t (N greater than M, but not much larger)
768 * Reduce to bidiagonal form without LQ decomposition
769 * A = QB * B * PB**T = QB * ( UB * S * VB**T ) * PB**T
770 * U = QB * UB; V**T = VB**T * PB**T
771 *
772 * Bidiagonalize A
773 * (Workspace: need 2*M+N, prefer 2*M+(M+N)*NB)
774 *
775  itauq = 1
776  itaup = itauq + m
777  itemp = itaup + m
778  id = 1
779  ie = id + m
780  itgkz = ie + m
781  CALL zgebrd( m, n, a, lda, rwork( id ), rwork( ie ),
782  $ work( itauq ), work( itaup ), work( itemp ),
783  $ lwork-itemp+1, info )
784  itempr = itgkz + m*(m*2+1)
785 *
786 * Solve eigenvalue problem TGK*Z=Z*S.
787 * (Workspace: need 2*M*M+14*M)
788 *
789  CALL dbdsvdx( 'L', jobz, rngtgk, m, rwork( id ),
790  $ rwork( ie ), vl, vu, iltgk, iutgk, ns, s,
791  $ rwork( itgkz ), m*2, rwork( itempr ),
792  $ iwork, info)
793 *
794 * If needed, compute left singular vectors.
795 *
796  IF( wantu ) THEN
797  k = itgkz
798  DO i = 1, ns
799  DO j = 1, m
800  u( j, i ) = dcmplx( rwork( k ), zero )
801  k = k + 1
802  END DO
803  k = k + m
804  END DO
805 *
806 * Call ZUNMBR to compute QB*UB.
807 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
808 *
809  CALL zunmbr( 'Q', 'L', 'N', m, ns, n, a, lda,
810  $ work( itauq ), u, ldu, work( itemp ),
811  $ lwork-itemp+1, info )
812  END IF
813 *
814 * If needed, compute right singular vectors.
815 *
816  IF( wantvt) THEN
817  k = itgkz + m
818  DO i = 1, ns
819  DO j = 1, m
820  vt( i, j ) = dcmplx( rwork( k ), zero )
821  k = k + 1
822  END DO
823  k = k + m
824  END DO
825  CALL zlaset( 'A', ns, n-m, czero, czero,
826  $ vt( 1,m+1 ), ldvt )
827 *
828 * Call ZUNMBR to compute VB**T * PB**T
829 * (Workspace in WORK( ITEMP ): need M, prefer M*NB)
830 *
831  CALL zunmbr( 'P', 'R', 'C', ns, n, m, a, lda,
832  $ work( itaup ), vt, ldvt, work( itemp ),
833  $ lwork-itemp+1, info )
834  END IF
835  END IF
836  END IF
837 *
838 * Undo scaling if necessary
839 *
840  IF( iscl.EQ.1 ) THEN
841  IF( anrm.GT.bignum )
842  $ CALL dlascl( 'G', 0, 0, bignum, anrm, minmn, 1,
843  $ s, minmn, info )
844  IF( anrm.LT.smlnum )
845  $ CALL dlascl( 'G', 0, 0, smlnum, anrm, minmn, 1,
846  $ s, minmn, info )
847  END IF
848 *
849 * Return optimal workspace in WORK(1)
850 *
851  work( 1 ) = dcmplx( dble( maxwrk ), zero )
852 *
853  RETURN
854 *
855 * End of ZGESVDX
856 *
857  END
zlacpy
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
zgeqrf
subroutine zgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGEQRF VARIANT: left-looking Level 3 BLAS of the algorithm.
Definition: zgeqrf.f:151
zunmbr
subroutine zunmbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMBR
Definition: zunmbr.f:198
dlascl
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: dlascl.f:145
zlaset
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: zlaset.f:108
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
zgesvdx
subroutine zgesvdx(JOBU, JOBVT, RANGE, M, N, A, LDA, VL, VU, IL, IU, NS, S, U, LDU, VT, LDVT, WORK, LWORK, RWORK, IWORK, INFO)
ZGESVDX computes the singular value decomposition (SVD) for GE matrices
Definition: zgesvdx.f:272
zgebrd
subroutine zgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
ZGEBRD
Definition: zgebrd.f:207
zunmqr
subroutine zunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMQR
Definition: zunmqr.f:169
zunmlq
subroutine zunmlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMLQ
Definition: zunmlq.f:169
dbdsvdx
subroutine dbdsvdx(UPLO, JOBZ, RANGE, N, D, E, VL, VU, IL, IU, NS, S, Z, LDZ, WORK, IWORK, INFO)
DBDSVDX
Definition: dbdsvdx.f:228
zgelqf
subroutine zgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGELQF
Definition: zgelqf.f:137
zlascl
subroutine zlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: zlascl.f:145