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