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