LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zgesvdx.f
Go to the documentation of this file.
1*> \brief <b> ZGESVDX 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 ZGESVDX + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgesvdx.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgesvdx.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgesvdx.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE ZGESVDX( 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* DOUBLE PRECISION VL, VU
28* ..
29* .. Array Arguments ..
30* INTEGER IWORK( * )
31* DOUBLE PRECISION S( * ), RWORK( * )
32* COMPLEX*16 A( LDA, * ), U( LDU, * ), VT( LDVT, * ),
33* $ WORK( * )
34* ..
35*
36*
37*> \par Purpose:
38* =============
39*>
40*> \verbatim
41*>
42*> ZGESVDX 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*> ZGESVDX uses an eigenvalue problem for obtaining the SVD, which
56*> allows for the computation of a subset of singular values and
57*> vectors. See DBDSVDX 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*16 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 DOUBLE PRECISION
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 DOUBLE PRECISION
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 DOUBLE PRECISION 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*16 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*16 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*16 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 DOUBLE PRECISION 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 DBDSVDX/DSTEVX.
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 DBDSVDX/DSTEVX.
250*> if INFO = N*2 + 1, an internal error occurred in
251*> DBDSVDX
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 zgesvdx( 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 DOUBLE PRECISION VL, VU
277* ..
278* .. Array Arguments ..
279 INTEGER IWORK( * )
280 DOUBLE PRECISION S( * ), RWORK( * )
281 COMPLEX*16 A( LDA, * ), U( LDU, * ), VT( LDVT, * ),
282 $ work( * )
283* ..
284*
285* =====================================================================
286*
287* .. Parameters ..
288 COMPLEX*16 CZERO, CONE
289 PARAMETER ( CZERO = ( 0.0d0, 0.0d0 ),
290 $ cone = ( 1.0d0, 0.0d0 ) )
291 DOUBLE PRECISION ZERO, ONE
292 PARAMETER ( ZERO = 0.0d0, one = 1.0d0 )
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 DOUBLE PRECISION ABSTOL, ANRM, BIGNUM, EPS, SMLNUM
301* ..
302* .. Local Arrays ..
303 DOUBLE PRECISION DUM( 1 )
304* ..
305* .. External Subroutines ..
306 EXTERNAL zgebrd, zgelqf, zgeqrf, zlascl, zlaset,
307 $ zlacpy,
308 $ zunmlq, zunmbr, zunmqr, dbdsvdx, dlascl, xerbla
309* ..
310* .. External Functions ..
311 LOGICAL LSAME
312 INTEGER ILAENV
313 DOUBLE PRECISION DLAMCH, ZLANGE
314 EXTERNAL lsame, ilaenv, dlamch, zlange
315* ..
316* .. Intrinsic Functions ..
317 INTRINSIC max, min, sqrt
318* ..
319* .. Executable Statements ..
320*
321* Test the input arguments.
322*
323 ns = 0
324 info = 0
325 abstol = 2*dlamch('S')
326 lquery = ( lwork.EQ.-1 )
327 minmn = min( m, n )
328
329 wantu = lsame( jobu, 'V' )
330 wantvt = lsame( jobvt, 'V' )
331 IF( wantu .OR. wantvt ) THEN
332 jobz = 'V'
333 ELSE
334 jobz = 'N'
335 END IF
336 alls = lsame( range, 'A' )
337 vals = lsame( range, 'V' )
338 inds = lsame( range, 'I' )
339*
340 info = 0
341 IF( .NOT.lsame( jobu, 'V' ) .AND.
342 $ .NOT.lsame( jobu, 'N' ) ) THEN
343 info = -1
344 ELSE IF( .NOT.lsame( jobvt, 'V' ) .AND.
345 $ .NOT.lsame( jobvt, 'N' ) ) THEN
346 info = -2
347 ELSE IF( .NOT.( alls .OR. vals .OR. inds ) ) THEN
348 info = -3
349 ELSE IF( m.LT.0 ) THEN
350 info = -4
351 ELSE IF( n.LT.0 ) THEN
352 info = -5
353 ELSE IF( m.GT.lda ) THEN
354 info = -7
355 ELSE IF( minmn.GT.0 ) THEN
356 IF( vals ) THEN
357 IF( vl.LT.zero ) THEN
358 info = -8
359 ELSE IF( vu.LE.vl ) THEN
360 info = -9
361 END IF
362 ELSE IF( inds ) THEN
363 IF( il.LT.1 .OR. il.GT.max( 1, minmn ) ) THEN
364 info = -10
365 ELSE IF( iu.LT.min( minmn, il ) .OR. iu.GT.minmn ) THEN
366 info = -11
367 END IF
368 END IF
369 IF( info.EQ.0 ) THEN
370 IF( wantu .AND. ldu.LT.m ) THEN
371 info = -15
372 ELSE IF( wantvt ) THEN
373 IF( inds ) THEN
374 IF( ldvt.LT.iu-il+1 ) THEN
375 info = -17
376 END IF
377 ELSE IF( ldvt.LT.minmn ) THEN
378 info = -17
379 END IF
380 END IF
381 END IF
382 END IF
383*
384* Compute workspace
385* (Note: Comments in the code beginning "Workspace:" describe the
386* minimal amount of workspace needed at that point in the code,
387* as well as the preferred amount for good performance.
388* NB refers to the optimal block size for the immediately
389* following subroutine, as returned by ILAENV.)
390*
391 IF( info.EQ.0 ) THEN
392 minwrk = 1
393 maxwrk = 1
394 IF( minmn.GT.0 ) THEN
395 IF( m.GE.n ) THEN
396 mnthr = ilaenv( 6, 'ZGESVD', jobu // jobvt, m, n, 0,
397 $ 0 )
398 IF( m.GE.mnthr ) THEN
399*
400* Path 1 (M much larger than N)
401*
402 minwrk = n*(n+5)
403 maxwrk = n + n*ilaenv(1,'ZGEQRF',' ',m,n,-1,-1)
404 maxwrk = max(maxwrk,
405 $ n*n+2*n+2*n*ilaenv(1,'ZGEBRD',' ',n,n,-1,
406 $ -1))
407 IF (wantu .OR. wantvt) THEN
408 maxwrk = max(maxwrk,
409 $ n*n+2*n+n*ilaenv(1,'ZUNMQR','LN',n,n,n,
410 $ -1))
411 END IF
412 ELSE
413*
414* Path 2 (M at least N, but not much larger)
415*
416 minwrk = 3*n + m
417 maxwrk = 2*n + (m+n)*ilaenv(1,'ZGEBRD',' ',m,n,-1,
418 $ -1)
419 IF (wantu .OR. wantvt) THEN
420 maxwrk = max(maxwrk,
421 $ 2*n+n*ilaenv(1,'ZUNMQR','LN',n,n,n,-1))
422 END IF
423 END IF
424 ELSE
425 mnthr = ilaenv( 6, 'ZGESVD', jobu // jobvt, m, n, 0,
426 $ 0 )
427 IF( n.GE.mnthr ) THEN
428*
429* Path 1t (N much larger than M)
430*
431 minwrk = m*(m+5)
432 maxwrk = m + m*ilaenv(1,'ZGELQF',' ',m,n,-1,-1)
433 maxwrk = max(maxwrk,
434 $ m*m+2*m+2*m*ilaenv(1,'ZGEBRD',' ',m,m,-1,
435 $ -1))
436 IF (wantu .OR. wantvt) THEN
437 maxwrk = max(maxwrk,
438 $ m*m+2*m+m*ilaenv(1,'ZUNMQR','LN',m,m,m,
439 $ -1))
440 END IF
441 ELSE
442*
443* Path 2t (N greater than M, but not much larger)
444*
445*
446 minwrk = 3*m + n
447 maxwrk = 2*m + (m+n)*ilaenv(1,'ZGEBRD',' ',m,n,-1,
448 $ -1)
449 IF (wantu .OR. wantvt) THEN
450 maxwrk = max(maxwrk,
451 $ 2*m+m*ilaenv(1,'ZUNMQR','LN',m,m,m,-1))
452 END IF
453 END IF
454 END IF
455 END IF
456 maxwrk = max( maxwrk, minwrk )
457 work( 1 ) = dcmplx( dble( maxwrk ), zero )
458*
459 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
460 info = -19
461 END IF
462 END IF
463*
464 IF( info.NE.0 ) THEN
465 CALL xerbla( 'ZGESVDX', -info )
466 RETURN
467 ELSE IF( lquery ) THEN
468 RETURN
469 END IF
470*
471* Quick return if possible
472*
473 IF( m.EQ.0 .OR. n.EQ.0 ) THEN
474 RETURN
475 END IF
476*
477* Set singular values indices accord to RANGE='A'.
478*
479 IF( alls ) THEN
480 rngtgk = 'I'
481 iltgk = 1
482 iutgk = min( m, n )
483 ELSE IF( inds ) THEN
484 rngtgk = 'I'
485 iltgk = il
486 iutgk = iu
487 ELSE
488 rngtgk = 'V'
489 iltgk = 0
490 iutgk = 0
491 END IF
492*
493* Get machine constants
494*
495 eps = dlamch( 'P' )
496 smlnum = sqrt( dlamch( 'S' ) ) / eps
497 bignum = one / smlnum
498*
499* Scale A if max element outside range [SMLNUM,BIGNUM]
500*
501 anrm = zlange( 'M', m, n, a, lda, dum )
502 iscl = 0
503 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
504 iscl = 1
505 CALL zlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
506 ELSE IF( anrm.GT.bignum ) THEN
507 iscl = 1
508 CALL zlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
509 END IF
510*
511 IF( m.GE.n ) THEN
512*
513* A has at least as many rows as columns. If A has sufficiently
514* more rows than columns, first reduce A using the QR
515* decomposition.
516*
517 IF( m.GE.mnthr ) THEN
518*
519* Path 1 (M much larger than N):
520* A = Q * R = Q * ( QB * B * PB**T )
521* = Q * ( QB * ( UB * S * VB**T ) * PB**T )
522* U = Q * QB * UB; V**T = VB**T * PB**T
523*
524* Compute A=Q*R
525* (Workspace: need 2*N, prefer N+N*NB)
526*
527 itau = 1
528 itemp = itau + n
529 CALL zgeqrf( m, n, a, lda, work( itau ), work( itemp ),
530 $ lwork-itemp+1, info )
531*
532* Copy R into WORK and bidiagonalize it:
533* (Workspace: need N*N+3*N, prefer N*N+N+2*N*NB)
534*
535 iqrf = itemp
536 itauq = itemp + n*n
537 itaup = itauq + n
538 itemp = itaup + n
539 id = 1
540 ie = id + n
541 itgkz = ie + n
542 CALL zlacpy( 'U', n, n, a, lda, work( iqrf ), n )
543 CALL zlaset( 'L', n-1, n-1, czero, czero,
544 $ work( iqrf+1 ), n )
545 CALL zgebrd( n, n, work( iqrf ), n, rwork( id ),
546 $ rwork( ie ), work( itauq ), work( itaup ),
547 $ work( itemp ), lwork-itemp+1, info )
548 itempr = itgkz + n*(n*2+1)
549*
550* Solve eigenvalue problem TGK*Z=Z*S.
551* (Workspace: need 2*N*N+14*N)
552*
553 CALL dbdsvdx( 'U', jobz, rngtgk, n, rwork( id ),
554 $ rwork( ie ), vl, vu, iltgk, iutgk, ns, s,
555 $ rwork( itgkz ), n*2, rwork( itempr ),
556 $ iwork, info)
557*
558* If needed, compute left singular vectors.
559*
560 IF( wantu ) THEN
561 k = itgkz
562 DO i = 1, ns
563 DO j = 1, n
564 u( j, i ) = dcmplx( rwork( k ), zero )
565 k = k + 1
566 END DO
567 k = k + n
568 END DO
569 CALL zlaset( 'A', m-n, ns, czero, czero, u( n+1,1 ),
570 $ ldu)
571*
572* Call ZUNMBR to compute QB*UB.
573* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
574*
575 CALL zunmbr( 'Q', 'L', 'N', n, ns, n, work( iqrf ), n,
576 $ work( itauq ), u, ldu, work( itemp ),
577 $ lwork-itemp+1, info )
578*
579* Call ZUNMQR to compute Q*(QB*UB).
580* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
581*
582 CALL zunmqr( 'L', 'N', m, ns, n, a, lda,
583 $ work( itau ), u, ldu, work( itemp ),
584 $ lwork-itemp+1, info )
585 END IF
586*
587* If needed, compute right singular vectors.
588*
589 IF( wantvt) THEN
590 k = itgkz + n
591 DO i = 1, ns
592 DO j = 1, n
593 vt( i, j ) = dcmplx( rwork( k ), zero )
594 k = k + 1
595 END DO
596 k = k + n
597 END DO
598*
599* Call ZUNMBR to compute VB**T * PB**T
600* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
601*
602 CALL zunmbr( 'P', 'R', 'C', ns, n, n, work( iqrf ), n,
603 $ work( itaup ), vt, ldvt, work( itemp ),
604 $ lwork-itemp+1, info )
605 END IF
606 ELSE
607*
608* Path 2 (M at least N, but not much larger)
609* Reduce A to bidiagonal form without QR decomposition
610* A = QB * B * PB**T = QB * ( UB * S * VB**T ) * PB**T
611* U = QB * UB; V**T = VB**T * PB**T
612*
613* Bidiagonalize A
614* (Workspace: need 2*N+M, prefer 2*N+(M+N)*NB)
615*
616 itauq = 1
617 itaup = itauq + n
618 itemp = itaup + n
619 id = 1
620 ie = id + n
621 itgkz = ie + n
622 CALL zgebrd( m, n, a, lda, rwork( id ), rwork( ie ),
623 $ work( itauq ), work( itaup ), work( itemp ),
624 $ lwork-itemp+1, info )
625 itempr = itgkz + n*(n*2+1)
626*
627* Solve eigenvalue problem TGK*Z=Z*S.
628* (Workspace: need 2*N*N+14*N)
629*
630 CALL dbdsvdx( 'U', jobz, rngtgk, n, rwork( id ),
631 $ rwork( ie ), vl, vu, iltgk, iutgk, ns, s,
632 $ rwork( itgkz ), n*2, rwork( itempr ),
633 $ iwork, info)
634*
635* If needed, compute left singular vectors.
636*
637 IF( wantu ) THEN
638 k = itgkz
639 DO i = 1, ns
640 DO j = 1, n
641 u( j, i ) = dcmplx( rwork( k ), zero )
642 k = k + 1
643 END DO
644 k = k + n
645 END DO
646 CALL zlaset( 'A', m-n, ns, czero, czero, u( n+1,1 ),
647 $ ldu)
648*
649* Call ZUNMBR to compute QB*UB.
650* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
651*
652 CALL zunmbr( 'Q', 'L', 'N', m, ns, n, a, lda,
653 $ work( itauq ), u, ldu, work( itemp ),
654 $ lwork-itemp+1, ierr )
655 END IF
656*
657* If needed, compute right singular vectors.
658*
659 IF( wantvt) THEN
660 k = itgkz + n
661 DO i = 1, ns
662 DO j = 1, n
663 vt( i, j ) = dcmplx( rwork( k ), zero )
664 k = k + 1
665 END DO
666 k = k + n
667 END DO
668*
669* Call ZUNMBR to compute VB**T * PB**T
670* (Workspace in WORK( ITEMP ): need N, prefer N*NB)
671*
672 CALL zunmbr( 'P', 'R', 'C', ns, n, n, a, lda,
673 $ work( itaup ), vt, ldvt, work( itemp ),
674 $ lwork-itemp+1, ierr )
675 END IF
676 END IF
677 ELSE
678*
679* A has more columns than rows. If A has sufficiently more
680* columns than rows, first reduce A using the LQ decomposition.
681*
682 IF( n.GE.mnthr ) THEN
683*
684* Path 1t (N much larger than M):
685* A = L * Q = ( QB * B * PB**T ) * Q
686* = ( QB * ( UB * S * VB**T ) * PB**T ) * Q
687* U = QB * UB ; V**T = VB**T * PB**T * Q
688*
689* Compute A=L*Q
690* (Workspace: need 2*M, prefer M+M*NB)
691*
692 itau = 1
693 itemp = itau + m
694 CALL zgelqf( m, n, a, lda, work( itau ), work( itemp ),
695 $ lwork-itemp+1, info )
696
697* Copy L into WORK and bidiagonalize it:
698* (Workspace in WORK( ITEMP ): need M*M+3*M, prefer M*M+M+2*M*NB)
699*
700 ilqf = itemp
701 itauq = ilqf + m*m
702 itaup = itauq + m
703 itemp = itaup + m
704 id = 1
705 ie = id + m
706 itgkz = ie + m
707 CALL zlacpy( 'L', m, m, a, lda, work( ilqf ), m )
708 CALL zlaset( 'U', m-1, m-1, czero, czero,
709 $ work( ilqf+m ), m )
710 CALL zgebrd( m, m, work( ilqf ), m, rwork( id ),
711 $ rwork( ie ), work( itauq ), work( itaup ),
712 $ work( itemp ), lwork-itemp+1, info )
713 itempr = itgkz + m*(m*2+1)
714*
715* Solve eigenvalue problem TGK*Z=Z*S.
716* (Workspace: need 2*M*M+14*M)
717*
718 CALL dbdsvdx( 'U', jobz, rngtgk, m, rwork( id ),
719 $ rwork( ie ), vl, vu, iltgk, iutgk, ns, s,
720 $ rwork( itgkz ), m*2, rwork( itempr ),
721 $ iwork, info)
722*
723* If needed, compute left singular vectors.
724*
725 IF( wantu ) THEN
726 k = itgkz
727 DO i = 1, ns
728 DO j = 1, m
729 u( j, i ) = dcmplx( rwork( k ), zero )
730 k = k + 1
731 END DO
732 k = k + m
733 END DO
734*
735* Call ZUNMBR to compute QB*UB.
736* (Workspace in WORK( ITEMP ): need M, prefer M*NB)
737*
738 CALL zunmbr( 'Q', 'L', 'N', m, ns, m, work( ilqf ), m,
739 $ work( itauq ), u, ldu, work( itemp ),
740 $ lwork-itemp+1, info )
741 END IF
742*
743* If needed, compute right singular vectors.
744*
745 IF( wantvt) THEN
746 k = itgkz + m
747 DO i = 1, ns
748 DO j = 1, m
749 vt( i, j ) = dcmplx( rwork( k ), zero )
750 k = k + 1
751 END DO
752 k = k + m
753 END DO
754 CALL zlaset( 'A', ns, n-m, czero, czero,
755 $ vt( 1,m+1 ), ldvt )
756*
757* Call ZUNMBR to compute (VB**T)*(PB**T)
758* (Workspace in WORK( ITEMP ): need M, prefer M*NB)
759*
760 CALL zunmbr( 'P', 'R', 'C', ns, m, m, work( ilqf ), m,
761 $ work( itaup ), vt, ldvt, work( itemp ),
762 $ lwork-itemp+1, info )
763*
764* Call ZUNMLQ to compute ((VB**T)*(PB**T))*Q.
765* (Workspace in WORK( ITEMP ): need M, prefer M*NB)
766*
767 CALL zunmlq( 'R', 'N', ns, n, m, a, lda,
768 $ work( itau ), vt, ldvt, work( itemp ),
769 $ lwork-itemp+1, info )
770 END IF
771 ELSE
772*
773* Path 2t (N greater than M, but not much larger)
774* Reduce to bidiagonal form without LQ decomposition
775* A = QB * B * PB**T = QB * ( UB * S * VB**T ) * PB**T
776* U = QB * UB; V**T = VB**T * PB**T
777*
778* Bidiagonalize A
779* (Workspace: need 2*M+N, prefer 2*M+(M+N)*NB)
780*
781 itauq = 1
782 itaup = itauq + m
783 itemp = itaup + m
784 id = 1
785 ie = id + m
786 itgkz = ie + m
787 CALL zgebrd( m, n, a, lda, rwork( id ), rwork( ie ),
788 $ work( itauq ), work( itaup ), work( itemp ),
789 $ lwork-itemp+1, info )
790 itempr = itgkz + m*(m*2+1)
791*
792* Solve eigenvalue problem TGK*Z=Z*S.
793* (Workspace: need 2*M*M+14*M)
794*
795 CALL dbdsvdx( 'L', jobz, rngtgk, m, rwork( id ),
796 $ rwork( ie ), vl, vu, iltgk, iutgk, ns, s,
797 $ rwork( itgkz ), m*2, rwork( itempr ),
798 $ iwork, info)
799*
800* If needed, compute left singular vectors.
801*
802 IF( wantu ) THEN
803 k = itgkz
804 DO i = 1, ns
805 DO j = 1, m
806 u( j, i ) = dcmplx( rwork( k ), zero )
807 k = k + 1
808 END DO
809 k = k + m
810 END DO
811*
812* Call ZUNMBR to compute QB*UB.
813* (Workspace in WORK( ITEMP ): need M, prefer M*NB)
814*
815 CALL zunmbr( 'Q', 'L', 'N', m, ns, n, a, lda,
816 $ work( itauq ), u, ldu, work( itemp ),
817 $ lwork-itemp+1, info )
818 END IF
819*
820* If needed, compute right singular vectors.
821*
822 IF( wantvt) THEN
823 k = itgkz + m
824 DO i = 1, ns
825 DO j = 1, m
826 vt( i, j ) = dcmplx( rwork( k ), zero )
827 k = k + 1
828 END DO
829 k = k + m
830 END DO
831 CALL zlaset( 'A', ns, n-m, czero, czero,
832 $ vt( 1,m+1 ), ldvt )
833*
834* Call ZUNMBR to compute VB**T * PB**T
835* (Workspace in WORK( ITEMP ): need M, prefer M*NB)
836*
837 CALL zunmbr( 'P', 'R', 'C', ns, n, m, a, lda,
838 $ work( itaup ), vt, ldvt, work( itemp ),
839 $ lwork-itemp+1, info )
840 END IF
841 END IF
842 END IF
843*
844* Undo scaling if necessary
845*
846 IF( iscl.EQ.1 ) THEN
847 IF( anrm.GT.bignum )
848 $ CALL dlascl( 'G', 0, 0, bignum, anrm, minmn, 1,
849 $ s, minmn, info )
850 IF( anrm.LT.smlnum )
851 $ CALL dlascl( 'G', 0, 0, smlnum, anrm, minmn, 1,
852 $ s, minmn, info )
853 END IF
854*
855* Return optimal workspace in WORK(1)
856*
857 work( 1 ) = dcmplx( dble( maxwrk ), zero )
858*
859 RETURN
860*
861* End of ZGESVDX
862*
863 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
zgebrd
subroutine zgebrd(m, n, a, lda, d, e, tauq, taup, work, lwork, info)
ZGEBRD
Definition zgebrd.f:204
zgelqf
subroutine zgelqf(m, n, a, lda, tau, work, lwork, info)
ZGELQF
Definition zgelqf.f:142
zgeqrf
subroutine zgeqrf(m, n, a, lda, tau, work, lwork, info)
ZGEQRF
Definition zgeqrf.f:144
zgesvdx
subroutine zgesvdx(jobu, jobvt, range, m, n, a, lda, vl, vu, il, iu, ns, s, u, ldu, vt, ldvt, work, lwork, rwork, iwork, info)
ZGESVDX computes the singular value decomposition (SVD) for GE matrices
Definition zgesvdx.f:268
zlacpy
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:101
zlascl
subroutine zlascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition zlascl.f:142
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
zlaset
subroutine zlaset(uplo, m, n, alpha, beta, a, lda)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition zlaset.f:104
zunmbr
subroutine zunmbr(vect, side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
ZUNMBR
Definition zunmbr.f:194
zunmlq
subroutine zunmlq(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
ZUNMLQ
Definition zunmlq.f:165
zunmqr
subroutine zunmqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
ZUNMQR
Definition zunmqr.f:165