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