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