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