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