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