LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
sgesvdx.f
Go to the documentation of this file.
1*> \brief <b> SGESVDX 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 SGESVDX + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgesvdx.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgesvdx.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgesvdx.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE SGESVDX( 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* REAL VL, VU
28* ..
29* .. Array Arguments ..
30* INTEGER IWORK( * )
31* REAL A( LDA, * ), S( * ), U( LDU, * ),
32* $ VT( LDVT, * ), WORK( * )
33* ..
34*
35*
36*> \par Purpose:
37* =============
38*>
39*> \verbatim
40*>
41*> SGESVDX 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*> SGESVDX uses an eigenvalue problem for obtaining the SVD, which
55*> allows for the computation of a subset of singular values and
56*> vectors. See SBDSVDX 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 REAL 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 REAL
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 REAL
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 REAL 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 REAL 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 REAL 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 REAL 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 SBDSVDX/SSTEVX.
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 SBDSVDX/SSTEVX.
243*> if INFO = N*2 + 1, an internal error occurred in
244*> SBDSVDX
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 sgesvdx( 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 REAL VL, VU
270* ..
271* .. Array Arguments ..
272 INTEGER IWORK( * )
273 REAL A( LDA, * ), S( * ), U( LDU, * ),
274 $ vt( ldvt, * ), work( * )
275* ..
276*
277* =====================================================================
278*
279* .. Parameters ..
280 REAL ZERO, ONE
281 PARAMETER ( ZERO = 0.0e0, one = 1.0e0 )
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 REAL ABSTOL, ANRM, BIGNUM, EPS, SMLNUM
290* ..
291* .. Local Arrays ..
292 REAL DUM( 1 )
293* ..
294* .. External Subroutines ..
295 EXTERNAL sbdsvdx, sgebrd, sgelqf, sgeqrf,
296 $ slacpy,
297 $ slascl, slaset, sormbr, sormlq, sormqr,
298 $ scopy, xerbla
299* ..
300* .. External Functions ..
301 LOGICAL LSAME
302 INTEGER ILAENV
303 REAL SLAMCH, SLANGE, SROUNDUP_LWORK
304 EXTERNAL lsame, ilaenv, slamch, slange,
305 $ sroundup_lwork
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*slamch('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, 'SGESVD', jobu // jobvt, m, n, 0,
388 $ 0 )
389 IF( m.GE.mnthr ) THEN
390*
391* Path 1 (M much larger than N)
392*
393 maxwrk = n +
394 $ n*ilaenv( 1, 'SGEQRF', ' ', m, n, -1, -1 )
395 maxwrk = max( maxwrk, n*(n+5) + 2*n*
396 $ ilaenv( 1, 'SGEBRD', ' ', n, n, -1, -1 ) )
397 IF (wantu) THEN
398 maxwrk = max(maxwrk,n*(n*3+6)+n*
399 $ ilaenv( 1, 'SORMQR', ' ', n, n, -1, -1 ) )
400 END IF
401 IF (wantvt) THEN
402 maxwrk = max(maxwrk,n*(n*3+6)+n*
403 $ ilaenv( 1, 'SORMLQ', ' ', n, n, -1, -1 ) )
404 END IF
405 minwrk = n*(n*3+20)
406 ELSE
407*
408* Path 2 (M at least N, but not much larger)
409*
410 maxwrk = 4*n + ( m+n )*
411 $ ilaenv( 1, 'SGEBRD', ' ', m, n, -1, -1 )
412 IF (wantu) THEN
413 maxwrk = max(maxwrk,n*(n*2+5)+n*
414 $ ilaenv( 1, 'SORMQR', ' ', n, n, -1, -1 ) )
415 END IF
416 IF (wantvt) THEN
417 maxwrk = max(maxwrk,n*(n*2+5)+n*
418 $ ilaenv( 1, 'SORMLQ', ' ', n, n, -1, -1 ) )
419 END IF
420 minwrk = max(n*(n*2+19),4*n+m)
421 END IF
422 ELSE
423 mnthr = ilaenv( 6, 'SGESVD', jobu // jobvt, m, n, 0,
424 $ 0 )
425 IF( n.GE.mnthr ) THEN
426*
427* Path 1t (N much larger than M)
428*
429 maxwrk = m +
430 $ m*ilaenv( 1, 'SGELQF', ' ', m, n, -1, -1 )
431 maxwrk = max( maxwrk, m*(m+5) + 2*m*
432 $ ilaenv( 1, 'SGEBRD', ' ', m, m, -1, -1 ) )
433 IF (wantu) THEN
434 maxwrk = max(maxwrk,m*(m*3+6)+m*
435 $ ilaenv( 1, 'SORMQR', ' ', m, m, -1, -1 ) )
436 END IF
437 IF (wantvt) THEN
438 maxwrk = max(maxwrk,m*(m*3+6)+m*
439 $ ilaenv( 1, 'SORMLQ', ' ', m, m, -1, -1 ) )
440 END IF
441 minwrk = m*(m*3+20)
442 ELSE
443*
444* Path 2t (N at least M, but not much larger)
445*
446 maxwrk = 4*m + ( m+n )*
447 $ ilaenv( 1, 'SGEBRD', ' ', m, n, -1, -1 )
448 IF (wantu) THEN
449 maxwrk = max(maxwrk,m*(m*2+5)+m*
450 $ ilaenv( 1, 'SORMQR', ' ', m, m, -1, -1 ) )
451 END IF
452 IF (wantvt) THEN
453 maxwrk = max(maxwrk,m*(m*2+5)+m*
454 $ ilaenv( 1, 'SORMLQ', ' ', m, m, -1, -1 ) )
455 END IF
456 minwrk = max(m*(m*2+19),4*m+n)
457 END IF
458 END IF
459 END IF
460 maxwrk = max( maxwrk, minwrk )
461 work( 1 ) = sroundup_lwork( maxwrk )
462*
463 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
464 info = -19
465 END IF
466 END IF
467*
468 IF( info.NE.0 ) THEN
469 CALL xerbla( 'SGESVDX', -info )
470 RETURN
471 ELSE IF( lquery ) THEN
472 RETURN
473 END IF
474*
475* Quick return if possible
476*
477 IF( m.EQ.0 .OR. n.EQ.0 ) THEN
478 RETURN
479 END IF
480*
481* Set singular values indices accord to RANGE.
482*
483 IF( alls ) THEN
484 rngtgk = 'I'
485 iltgk = 1
486 iutgk = min( m, n )
487 ELSE IF( inds ) THEN
488 rngtgk = 'I'
489 iltgk = il
490 iutgk = iu
491 ELSE
492 rngtgk = 'V'
493 iltgk = 0
494 iutgk = 0
495 END IF
496*
497* Get machine constants
498*
499 eps = slamch( 'P' )
500 smlnum = sqrt( slamch( 'S' ) ) / eps
501 bignum = one / smlnum
502*
503* Scale A if max element outside range [SMLNUM,BIGNUM]
504*
505 anrm = slange( 'M', m, n, a, lda, dum )
506 iscl = 0
507 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
508 iscl = 1
509 CALL slascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
510 ELSE IF( anrm.GT.bignum ) THEN
511 iscl = 1
512 CALL slascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
513 END IF
514*
515 IF( m.GE.n ) THEN
516*
517* A has at least as many rows as columns. If A has sufficiently
518* more rows than columns, first reduce A using the QR
519* decomposition.
520*
521 IF( m.GE.mnthr ) THEN
522*
523* Path 1 (M much larger than N):
524* A = Q * R = Q * ( QB * B * PB**T )
525* = Q * ( QB * ( UB * S * VB**T ) * PB**T )
526* U = Q * QB * UB; V**T = VB**T * PB**T
527*
528* Compute A=Q*R
529* (Workspace: need 2*N, prefer N+N*NB)
530*
531 itau = 1
532 itemp = itau + n
533 CALL sgeqrf( m, n, a, lda, work( itau ), work( itemp ),
534 $ lwork-itemp+1, info )
535*
536* Copy R into WORK and bidiagonalize it:
537* (Workspace: need N*N+5*N, prefer N*N+4*N+2*N*NB)
538*
539 iqrf = itemp
540 id = iqrf + n*n
541 ie = id + n
542 itauq = ie + n
543 itaup = itauq + n
544 itemp = itaup + n
545 CALL slacpy( 'U', n, n, a, lda, work( iqrf ), n )
546 CALL slaset( 'L', n-1, n-1, zero, zero, work( iqrf+1 ),
547 $ n )
548 CALL sgebrd( n, n, work( iqrf ), n, work( id ),
549 $ work( ie ),
550 $ work( itauq ), work( itaup ), work( itemp ),
551 $ lwork-itemp+1, info )
552*
553* Solve eigenvalue problem TGK*Z=Z*S.
554* (Workspace: need 14*N + 2*N*(N+1))
555*
556 itgkz = itemp
557 itemp = itgkz + n*(n*2+1)
558 CALL sbdsvdx( 'U', jobz, rngtgk, n, work( id ),
559 $ work( ie ),
560 $ vl, vu, iltgk, iutgk, ns, s, work( itgkz ),
561 $ n*2, work( itemp ), iwork, info)
562*
563* If needed, compute left singular vectors.
564*
565 IF( wantu ) THEN
566 j = itgkz
567 DO i = 1, ns
568 CALL scopy( n, work( j ), 1, u( 1,i ), 1 )
569 j = j + n*2
570 END DO
571 CALL slaset( 'A', m-n, ns, zero, zero, u( n+1,1 ),
572 $ ldu )
573*
574* Call SORMBR to compute QB*UB.
575* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
576*
577 CALL sormbr( 'Q', 'L', 'N', n, ns, n, work( iqrf ), n,
578 $ work( itauq ), u, ldu, work( itemp ),
579 $ lwork-itemp+1, info )
580*
581* Call SORMQR to compute Q*(QB*UB).
582* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
583*
584 CALL sormqr( 'L', 'N', m, ns, n, a, lda,
585 $ work( itau ), u, ldu, work( itemp ),
586 $ lwork-itemp+1, info )
587 END IF
588*
589* If needed, compute right singular vectors.
590*
591 IF( wantvt) THEN
592 j = itgkz + n
593 DO i = 1, ns
594 CALL scopy( n, work( j ), 1, vt( i,1 ), ldvt )
595 j = j + n*2
596 END DO
597*
598* Call SORMBR to compute VB**T * PB**T
599* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
600*
601 CALL sormbr( 'P', 'R', 'T', ns, n, n, work( iqrf ), n,
602 $ work( itaup ), vt, ldvt, work( itemp ),
603 $ lwork-itemp+1, info )
604 END IF
605 ELSE
606*
607* Path 2 (M at least N, but not much larger)
608* Reduce A to bidiagonal form without QR decomposition
609* A = QB * B * PB**T = QB * ( UB * S * VB**T ) * PB**T
610* U = QB * UB; V**T = VB**T * PB**T
611*
612* Bidiagonalize A
613* (Workspace: need 4*N+M, prefer 4*N+(M+N)*NB)
614*
615 id = 1
616 ie = id + n
617 itauq = ie + n
618 itaup = itauq + n
619 itemp = itaup + n
620 CALL sgebrd( m, n, a, lda, work( id ), work( ie ),
621 $ work( itauq ), work( itaup ), work( itemp ),
622 $ lwork-itemp+1, info )
623*
624* Solve eigenvalue problem TGK*Z=Z*S.
625* (Workspace: need 14*N + 2*N*(N+1))
626*
627 itgkz = itemp
628 itemp = itgkz + n*(n*2+1)
629 CALL sbdsvdx( 'U', jobz, rngtgk, n, work( id ),
630 $ work( ie ),
631 $ vl, vu, iltgk, iutgk, ns, s, work( itgkz ),
632 $ n*2, work( itemp ), iwork, info)
633*
634* If needed, compute left singular vectors.
635*
636 IF( wantu ) THEN
637 j = itgkz
638 DO i = 1, ns
639 CALL scopy( n, work( j ), 1, u( 1,i ), 1 )
640 j = j + n*2
641 END DO
642 CALL slaset( 'A', m-n, ns, zero, zero, u( n+1,1 ),
643 $ ldu )
644*
645* Call SORMBR to compute QB*UB.
646* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
647*
648 CALL sormbr( 'Q', 'L', 'N', m, ns, n, a, lda,
649 $ work( itauq ), u, ldu, work( itemp ),
650 $ lwork-itemp+1, ierr )
651 END IF
652*
653* If needed, compute right singular vectors.
654*
655 IF( wantvt) THEN
656 j = itgkz + n
657 DO i = 1, ns
658 CALL scopy( n, work( j ), 1, vt( i,1 ), ldvt )
659 j = j + n*2
660 END DO
661*
662* Call SORMBR to compute VB**T * PB**T
663* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
664*
665 CALL sormbr( 'P', 'R', 'T', ns, n, n, a, lda,
666 $ work( itaup ), vt, ldvt, work( itemp ),
667 $ lwork-itemp+1, ierr )
668 END IF
669 END IF
670 ELSE
671*
672* A has more columns than rows. If A has sufficiently more
673* columns than rows, first reduce A using the LQ decomposition.
674*
675 IF( n.GE.mnthr ) THEN
676*
677* Path 1t (N much larger than M):
678* A = L * Q = ( QB * B * PB**T ) * Q
679* = ( QB * ( UB * S * VB**T ) * PB**T ) * Q
680* U = QB * UB ; V**T = VB**T * PB**T * Q
681*
682* Compute A=L*Q
683* (Workspace: need 2*M, prefer M+M*NB)
684*
685 itau = 1
686 itemp = itau + m
687 CALL sgelqf( m, n, a, lda, work( itau ), work( itemp ),
688 $ lwork-itemp+1, info )
689
690* Copy L into WORK and bidiagonalize it:
691* (Workspace in WORK( ITEMP ): need M*M+5*N, prefer M*M+4*M+2*M*NB)
692*
693 ilqf = itemp
694 id = ilqf + m*m
695 ie = id + m
696 itauq = ie + m
697 itaup = itauq + m
698 itemp = itaup + m
699 CALL slacpy( 'L', m, m, a, lda, work( ilqf ), m )
700 CALL slaset( 'U', m-1, m-1, zero, zero, work( ilqf+m ),
701 $ m )
702 CALL sgebrd( m, m, work( ilqf ), m, work( id ),
703 $ work( ie ),
704 $ work( itauq ), work( itaup ), work( itemp ),
705 $ lwork-itemp+1, info )
706*
707* Solve eigenvalue problem TGK*Z=Z*S.
708* (Workspace: need 2*M*M+14*M)
709*
710 itgkz = itemp
711 itemp = itgkz + m*(m*2+1)
712 CALL sbdsvdx( 'U', jobz, rngtgk, m, work( id ),
713 $ work( ie ),
714 $ vl, vu, iltgk, iutgk, ns, s, work( itgkz ),
715 $ m*2, work( itemp ), iwork, info)
716*
717* If needed, compute left singular vectors.
718*
719 IF( wantu ) THEN
720 j = itgkz
721 DO i = 1, ns
722 CALL scopy( m, work( j ), 1, u( 1,i ), 1 )
723 j = j + m*2
724 END DO
725*
726* Call SORMBR to compute QB*UB.
727* (Workspace in WORK( ITEMP ): need M, prefer M*NB)
728*
729 CALL sormbr( '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 j = itgkz + m
738 DO i = 1, ns
739 CALL scopy( m, work( j ), 1, vt( i,1 ), ldvt )
740 j = j + m*2
741 END DO
742 CALL slaset( 'A', ns, n-m, zero, zero, vt( 1,m+1 ),
743 $ ldvt)
744*
745* Call SORMBR to compute (VB**T)*(PB**T)
746* (Workspace in WORK( ITEMP ): need M, prefer M*NB)
747*
748 CALL sormbr( 'P', 'R', 'T', ns, m, m, work( ilqf ), m,
749 $ work( itaup ), vt, ldvt, work( itemp ),
750 $ lwork-itemp+1, info )
751*
752* Call SORMLQ to compute ((VB**T)*(PB**T))*Q.
753* (Workspace in WORK( ITEMP ): need M, prefer M*NB)
754*
755 CALL sormlq( 'R', 'N', ns, n, m, a, lda,
756 $ work( itau ), vt, ldvt, work( itemp ),
757 $ lwork-itemp+1, info )
758 END IF
759 ELSE
760*
761* Path 2t (N greater than M, but not much larger)
762* Reduce to bidiagonal form without LQ decomposition
763* A = QB * B * PB**T = QB * ( UB * S * VB**T ) * PB**T
764* U = QB * UB; V**T = VB**T * PB**T
765*
766* Bidiagonalize A
767* (Workspace: need 4*M+N, prefer 4*M+(M+N)*NB)
768*
769 id = 1
770 ie = id + m
771 itauq = ie + m
772 itaup = itauq + m
773 itemp = itaup + m
774 CALL sgebrd( m, n, a, lda, work( id ), work( ie ),
775 $ work( itauq ), work( itaup ), work( itemp ),
776 $ lwork-itemp+1, info )
777*
778* Solve eigenvalue problem TGK*Z=Z*S.
779* (Workspace: need 2*M*M+14*M)
780*
781 itgkz = itemp
782 itemp = itgkz + m*(m*2+1)
783 CALL sbdsvdx( 'L', jobz, rngtgk, m, work( id ),
784 $ work( ie ),
785 $ vl, vu, iltgk, iutgk, ns, s, work( itgkz ),
786 $ m*2, work( itemp ), iwork, info)
787*
788* If needed, compute left singular vectors.
789*
790 IF( wantu ) THEN
791 j = itgkz
792 DO i = 1, ns
793 CALL scopy( m, work( j ), 1, u( 1,i ), 1 )
794 j = j + m*2
795 END DO
796*
797* Call SORMBR to compute QB*UB.
798* (Workspace in WORK( ITEMP ): need M, prefer M*NB)
799*
800 CALL sormbr( 'Q', 'L', 'N', m, ns, n, a, lda,
801 $ work( itauq ), u, ldu, work( itemp ),
802 $ lwork-itemp+1, info )
803 END IF
804*
805* If needed, compute right singular vectors.
806*
807 IF( wantvt) THEN
808 j = itgkz + m
809 DO i = 1, ns
810 CALL scopy( m, work( j ), 1, vt( i,1 ), ldvt )
811 j = j + m*2
812 END DO
813 CALL slaset( 'A', ns, n-m, zero, zero, vt( 1,m+1 ),
814 $ ldvt)
815*
816* Call SORMBR to compute VB**T * PB**T
817* (Workspace in WORK( ITEMP ): need M, prefer M*NB)
818*
819 CALL sormbr( 'P', 'R', 'T', ns, n, m, a, lda,
820 $ work( itaup ), vt, ldvt, work( itemp ),
821 $ lwork-itemp+1, info )
822 END IF
823 END IF
824 END IF
825*
826* Undo scaling if necessary
827*
828 IF( iscl.EQ.1 ) THEN
829 IF( anrm.GT.bignum )
830 $ CALL slascl( 'G', 0, 0, bignum, anrm, minmn, 1,
831 $ s, minmn, info )
832 IF( anrm.LT.smlnum )
833 $ CALL slascl( 'G', 0, 0, smlnum, anrm, minmn, 1,
834 $ s, minmn, info )
835 END IF
836*
837* Return optimal workspace in WORK(1)
838*
839 work( 1 ) = sroundup_lwork( maxwrk )
840*
841 RETURN
842*
843* End of SGESVDX
844*
845 END
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
sbdsvdx
subroutine sbdsvdx(uplo, jobz, range, n, d, e, vl, vu, il, iu, ns, s, z, ldz, work, iwork, info)
SBDSVDX
Definition sbdsvdx.f:224
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
sgebrd
subroutine sgebrd(m, n, a, lda, d, e, tauq, taup, work, lwork, info)
SGEBRD
Definition sgebrd.f:204
sgelqf
subroutine sgelqf(m, n, a, lda, tau, work, lwork, info)
SGELQF
Definition sgelqf.f:142
sgeqrf
subroutine sgeqrf(m, n, a, lda, tau, work, lwork, info)
SGEQRF
Definition sgeqrf.f:144
sgesvdx
subroutine sgesvdx(jobu, jobvt, range, m, n, a, lda, vl, vu, il, iu, ns, s, u, ldu, vt, ldvt, work, lwork, iwork, info)
SGESVDX computes the singular value decomposition (SVD) for GE matrices
Definition sgesvdx.f:261
slacpy
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:101
slascl
subroutine slascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition slascl.f:142
slaset
subroutine slaset(uplo, m, n, alpha, beta, a, lda)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition slaset.f:108
sormbr
subroutine sormbr(vect, side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
SORMBR
Definition sormbr.f:194
sormlq
subroutine sormlq(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
SORMLQ
Definition sormlq.f:166
sormqr
subroutine sormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
SORMQR
Definition sormqr.f:166