LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dbdsvdx.f
Go to the documentation of this file.
1*> \brief \b DBDSVDX
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download DBDSVDX + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dbdsvdx.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dbdsvdx.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dbdsvdx.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE DBDSVDX( UPLO, JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
20* $ NS, S, Z, LDZ, WORK, IWORK, INFO )
21*
22* .. Scalar Arguments ..
23* CHARACTER JOBZ, RANGE, UPLO
24* INTEGER IL, INFO, IU, LDZ, N, NS
25* DOUBLE PRECISION VL, VU
26* ..
27* .. Array Arguments ..
28* INTEGER IWORK( * )
29* DOUBLE PRECISION D( * ), E( * ), S( * ), WORK( * ),
30* Z( LDZ, * )
31* ..
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> DBDSVDX computes the singular value decomposition (SVD) of a real
39*> N-by-N (upper or lower) bidiagonal matrix B, B = U * S * VT,
40*> where S is a diagonal matrix with non-negative diagonal elements
41*> (the singular values of B), and U and VT are orthogonal matrices
42*> of left and right singular vectors, respectively.
43*>
44*> Given an upper bidiagonal B with diagonal D = [ d_1 d_2 ... d_N ]
45*> and superdiagonal E = [ e_1 e_2 ... e_N-1 ], DBDSVDX computes the
46*> singular value decomposition of B through the eigenvalues and
47*> eigenvectors of the N*2-by-N*2 tridiagonal matrix
48*>
49*> | 0 d_1 |
50*> | d_1 0 e_1 |
51*> TGK = | e_1 0 d_2 |
52*> | d_2 . . |
53*> | . . . |
54*>
55*> If (s,u,v) is a singular triplet of B with ||u|| = ||v|| = 1, then
56*> (+/-s,q), ||q|| = 1, are eigenpairs of TGK, with q = P * ( u' +/-v' ) /
57*> sqrt(2) = ( v_1 u_1 v_2 u_2 ... v_n u_n ) / sqrt(2), and
58*> P = [ e_{n+1} e_{1} e_{n+2} e_{2} ... ].
59*>
60*> Given a TGK matrix, one can either a) compute -s,-v and change signs
61*> so that the singular values (and corresponding vectors) are already in
62*> descending order (as in DGESVD/DGESDD) or b) compute s,v and reorder
63*> the values (and corresponding vectors). DBDSVDX implements a) by
64*> calling DSTEVX (bisection plus inverse iteration, to be replaced
65*> with a version of the Multiple Relative Robust Representation
66*> algorithm. (See P. Willems and B. Lang, A framework for the MR^3
67*> algorithm: theory and implementation, SIAM J. Sci. Comput.,
68*> 35:740-766, 2013.)
69*> \endverbatim
70*
71* Arguments:
72* ==========
73*
74*> \param[in] UPLO
75*> \verbatim
76*> UPLO is CHARACTER*1
77*> = 'U': B is upper bidiagonal;
78*> = 'L': B is lower bidiagonal.
79*> \endverbatim
80*>
81*> \param[in] JOBZ
82*> \verbatim
83*> JOBZ is CHARACTER*1
84*> = 'N': Compute singular values only;
85*> = 'V': Compute singular values and singular vectors.
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] N
98*> \verbatim
99*> N is INTEGER
100*> The order of the bidiagonal matrix. N >= 0.
101*> \endverbatim
102*>
103*> \param[in] D
104*> \verbatim
105*> D is DOUBLE PRECISION array, dimension (N)
106*> The n diagonal elements of the bidiagonal matrix B.
107*> \endverbatim
108*>
109*> \param[in] E
110*> \verbatim
111*> E is DOUBLE PRECISION array, dimension (max(1,N-1))
112*> The (n-1) superdiagonal elements of the bidiagonal matrix
113*> B in elements 1 to N-1.
114*> \endverbatim
115*>
116*> \param[in] VL
117*> \verbatim
118*> VL is DOUBLE PRECISION
119*> If RANGE='V', the lower bound of the interval to
120*> be searched for singular values. VU > VL.
121*> Not referenced if RANGE = 'A' or 'I'.
122*> \endverbatim
123*>
124*> \param[in] VU
125*> \verbatim
126*> VU is DOUBLE PRECISION
127*> If RANGE='V', the upper 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] IL
133*> \verbatim
134*> IL is INTEGER
135*> If RANGE='I', the index of the
136*> smallest singular value to be returned.
137*> 1 <= IL <= IU <= min(M,N), if min(M,N) > 0.
138*> Not referenced if RANGE = 'A' or 'V'.
139*> \endverbatim
140*>
141*> \param[in] IU
142*> \verbatim
143*> IU is INTEGER
144*> If RANGE='I', the index of the
145*> largest singular value to be returned.
146*> 1 <= IL <= IU <= min(M,N), if min(M,N) > 0.
147*> Not referenced if RANGE = 'A' or 'V'.
148*> \endverbatim
149*>
150*> \param[out] NS
151*> \verbatim
152*> NS is INTEGER
153*> The total number of singular values found. 0 <= NS <= N.
154*> If RANGE = 'A', NS = N, and if RANGE = 'I', NS = IU-IL+1.
155*> \endverbatim
156*>
157*> \param[out] S
158*> \verbatim
159*> S is DOUBLE PRECISION array, dimension (N)
160*> The first NS elements contain the selected singular values in
161*> ascending order.
162*> \endverbatim
163*>
164*> \param[out] Z
165*> \verbatim
166*> Z is DOUBLE PRECISION array, dimension (2*N,K)
167*> If JOBZ = 'V', then if INFO = 0 the first NS columns of Z
168*> contain the singular vectors of the matrix B corresponding to
169*> the selected singular values, with U in rows 1 to N and V
170*> in rows N+1 to N*2, i.e.
171*> Z = [ U ]
172*> [ V ]
173*> If JOBZ = 'N', then Z is not referenced.
174*> Note: The user must ensure that at least K = NS+1 columns are
175*> supplied in the array Z; if RANGE = 'V', the exact value of
176*> NS is not known in advance and an upper bound must be used.
177*> \endverbatim
178*>
179*> \param[in] LDZ
180*> \verbatim
181*> LDZ is INTEGER
182*> The leading dimension of the array Z. LDZ >= 1, and if
183*> JOBZ = 'V', LDZ >= max(2,N*2).
184*> \endverbatim
185*>
186*> \param[out] WORK
187*> \verbatim
188*> WORK is DOUBLE PRECISION array, dimension (14*N)
189*> \endverbatim
190*>
191*> \param[out] IWORK
192*> \verbatim
193*> IWORK is INTEGER array, dimension (12*N)
194*> If JOBZ = 'V', then if INFO = 0, the first NS elements of
195*> IWORK are zero. If INFO > 0, then IWORK contains the indices
196*> of the eigenvectors that failed to converge in DSTEVX.
197*> \endverbatim
198*>
199*> \param[out] INFO
200*> \verbatim
201*> INFO is INTEGER
202*> = 0: successful exit
203*> < 0: if INFO = -i, the i-th argument had an illegal value
204*> > 0: if INFO = i, then i eigenvectors failed to converge
205*> in DSTEVX. The indices of the eigenvectors
206*> (as returned by DSTEVX) are stored in the
207*> array IWORK.
208*> if INFO = N*2 + 1, an internal error occurred.
209*> \endverbatim
210*
211* Authors:
212* ========
213*
214*> \author Univ. of Tennessee
215*> \author Univ. of California Berkeley
216*> \author Univ. of Colorado Denver
217*> \author NAG Ltd.
218*
219*> \ingroup bdsvdx
220*
221* =====================================================================
222 SUBROUTINE dbdsvdx( UPLO, JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
223 $ NS, S, Z, LDZ, WORK, IWORK, INFO)
224*
225* -- LAPACK driver routine --
226* -- LAPACK is a software package provided by Univ. of Tennessee, --
227* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
228*
229* .. Scalar Arguments ..
230 CHARACTER JOBZ, RANGE, UPLO
231 INTEGER IL, INFO, IU, LDZ, N, NS
232 DOUBLE PRECISION VL, VU
233* ..
234* .. Array Arguments ..
235 INTEGER IWORK( * )
236 DOUBLE PRECISION D( * ), E( * ), S( * ), WORK( * ),
237 $ z( ldz, * )
238* ..
239*
240* =====================================================================
241*
242* .. Parameters ..
243 DOUBLE PRECISION ZERO, ONE, TEN, HNDRD, MEIGTH
244 parameter( zero = 0.0d0, one = 1.0d0, ten = 10.0d0,
245 $ hndrd = 100.0d0, meigth = -0.1250d0 )
246 DOUBLE PRECISION FUDGE
247 parameter( fudge = 2.0d0 )
248* ..
249* .. Local Scalars ..
250 CHARACTER RNGVX
251 LOGICAL ALLSV, INDSV, LOWER, SPLIT, SVEQ0, VALSV, WANTZ
252 INTEGER I, ICOLZ, IDBEG, IDEND, IDTGK, IDPTR, IEPTR,
253 $ ietgk, iifail, iiwork, iltgk, irowu, irowv,
254 $ irowz, isbeg, isplt, itemp, iutgk, j, k,
255 $ ntgk, nru, nrv, nsl
256 DOUBLE PRECISION ABSTOL, EPS, EMIN, MU, NRMU, NRMV, ORTOL, SMAX,
257 $ smin, sqrt2, thresh, tol, ulp,
258 $ vltgk, vutgk, zjtji
259* ..
260* .. External Functions ..
261 LOGICAL LSAME
262 INTEGER IDAMAX
263 DOUBLE PRECISION DDOT, DLAMCH, DNRM2
264 EXTERNAL idamax, lsame, daxpy, ddot, dlamch,
265 $ dnrm2
266* ..
267* .. External Subroutines ..
268 EXTERNAL dstevx, dcopy, dlaset, dscal, dswap,
269 $ xerbla
270* ..
271* .. Intrinsic Functions ..
272 INTRINSIC abs, dble, sign, sqrt
273* ..
274* .. Executable Statements ..
275*
276* Test the input parameters.
277*
278 allsv = lsame( range, 'A' )
279 valsv = lsame( range, 'V' )
280 indsv = lsame( range, 'I' )
281 wantz = lsame( jobz, 'V' )
282 lower = lsame( uplo, 'L' )
283*
284 info = 0
285 IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lower ) THEN
286 info = -1
287 ELSE IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
288 info = -2
289 ELSE IF( .NOT.( allsv .OR. valsv .OR. indsv ) ) THEN
290 info = -3
291 ELSE IF( n.LT.0 ) THEN
292 info = -4
293 ELSE IF( n.GT.0 ) THEN
294 IF( valsv ) THEN
295 IF( vl.LT.zero ) THEN
296 info = -7
297 ELSE IF( vu.LE.vl ) THEN
298 info = -8
299 END IF
300 ELSE IF( indsv ) THEN
301 IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
302 info = -9
303 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
304 info = -10
305 END IF
306 END IF
307 END IF
308 IF( info.EQ.0 ) THEN
309 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n*2 ) ) info = -14
310 END IF
311*
312 IF( info.NE.0 ) THEN
313 CALL xerbla( 'DBDSVDX', -info )
314 RETURN
315 END IF
316*
317* Quick return if possible (N.LE.1)
318*
319 ns = 0
320 IF( n.EQ.0 ) RETURN
321*
322 IF( n.EQ.1 ) THEN
323 IF( allsv .OR. indsv ) THEN
324 ns = 1
325 s( 1 ) = abs( d( 1 ) )
326 ELSE
327 IF( vl.LT.abs( d( 1 ) ) .AND. vu.GE.abs( d( 1 ) ) ) THEN
328 ns = 1
329 s( 1 ) = abs( d( 1 ) )
330 END IF
331 END IF
332 IF( wantz ) THEN
333 z( 1, 1 ) = sign( one, d( 1 ) )
334 z( 2, 1 ) = one
335 END IF
336 RETURN
337 END IF
338*
339 abstol = 2*dlamch( 'Safe Minimum' )
340 ulp = dlamch( 'Precision' )
341 eps = dlamch( 'Epsilon' )
342 sqrt2 = sqrt( 2.0d0 )
343 ortol = sqrt( ulp )
344*
345* Criterion for splitting is taken from DBDSQR when singular
346* values are computed to relative accuracy TOL. (See J. Demmel and
347* W. Kahan, Accurate singular values of bidiagonal matrices, SIAM
348* J. Sci. and Stat. Comput., 11:873–912, 1990.)
349*
350 tol = max( ten, min( hndrd, eps**meigth ) )*eps
351*
352* Compute approximate maximum, minimum singular values.
353*
354 i = idamax( n, d, 1 )
355 smax = abs( d( i ) )
356 i = idamax( n-1, e, 1 )
357 smax = max( smax, abs( e( i ) ) )
358*
359* Compute threshold for neglecting D's and E's.
360*
361 smin = abs( d( 1 ) )
362 IF( smin.NE.zero ) THEN
363 mu = smin
364 DO i = 2, n
365 mu = abs( d( i ) )*( mu / ( mu+abs( e( i-1 ) ) ) )
366 smin = min( smin, mu )
367 IF( smin.EQ.zero ) EXIT
368 END DO
369 END IF
370 smin = smin / sqrt( dble( n ) )
371 thresh = tol*smin
372*
373* Check for zeros in D and E (splits), i.e. submatrices.
374*
375 DO i = 1, n-1
376 IF( abs( d( i ) ).LE.thresh ) d( i ) = zero
377 IF( abs( e( i ) ).LE.thresh ) e( i ) = zero
378 END DO
379 IF( abs( d( n ) ).LE.thresh ) d( n ) = zero
380*
381* Pointers for arrays used by DSTEVX.
382*
383 idtgk = 1
384 ietgk = idtgk + n*2
385 itemp = ietgk + n*2
386 iifail = 1
387 iiwork = iifail + n*2
388*
389* Set RNGVX, which corresponds to RANGE for DSTEVX in TGK mode.
390* VL,VU or IL,IU are redefined to conform to implementation a)
391* described in the leading comments.
392*
393 iltgk = 0
394 iutgk = 0
395 vltgk = zero
396 vutgk = zero
397*
398 IF( allsv ) THEN
399*
400* All singular values will be found. We aim at -s (see
401* leading comments) with RNGVX = 'I'. IL and IU are set
402* later (as ILTGK and IUTGK) according to the dimension
403* of the active submatrix.
404*
405 rngvx = 'I'
406 IF( wantz ) CALL dlaset( 'F', n*2, n+1, zero, zero, z, ldz )
407 ELSE IF( valsv ) THEN
408*
409* Find singular values in a half-open interval. We aim
410* at -s (see leading comments) and we swap VL and VU
411* (as VUTGK and VLTGK), changing their signs.
412*
413 rngvx = 'V'
414 vltgk = -vu
415 vutgk = -vl
416 work( idtgk:idtgk+2*n-1 ) = zero
417 CALL dcopy( n, d, 1, work( ietgk ), 2 )
418 CALL dcopy( n-1, e, 1, work( ietgk+1 ), 2 )
419 CALL dstevx( 'N', 'V', n*2, work( idtgk ), work( ietgk ),
420 $ vltgk, vutgk, iltgk, iltgk, abstol, ns, s,
421 $ z, ldz, work( itemp ), iwork( iiwork ),
422 $ iwork( iifail ), info )
423 IF( ns.EQ.0 ) THEN
424 RETURN
425 ELSE
426 IF( wantz ) CALL dlaset( 'F', n*2, ns, zero, zero, z,
427 $ ldz )
428 END IF
429 ELSE IF( indsv ) THEN
430*
431* Find the IL-th through the IU-th singular values. We aim
432* at -s (see leading comments) and indices are mapped into
433* values, therefore mimicking DSTEBZ, where
434*
435* GL = GL - FUDGE*TNORM*ULP*N - FUDGE*TWO*PIVMIN
436* GU = GU + FUDGE*TNORM*ULP*N + FUDGE*PIVMIN
437*
438 iltgk = il
439 iutgk = iu
440 rngvx = 'V'
441 work( idtgk:idtgk+2*n-1 ) = zero
442 CALL dcopy( n, d, 1, work( ietgk ), 2 )
443 CALL dcopy( n-1, e, 1, work( ietgk+1 ), 2 )
444 CALL dstevx( 'N', 'I', n*2, work( idtgk ), work( ietgk ),
445 $ vltgk, vltgk, iltgk, iltgk, abstol, ns, s,
446 $ z, ldz, work( itemp ), iwork( iiwork ),
447 $ iwork( iifail ), info )
448 vltgk = s( 1 ) - fudge*smax*ulp*n
449 work( idtgk:idtgk+2*n-1 ) = zero
450 CALL dcopy( n, d, 1, work( ietgk ), 2 )
451 CALL dcopy( n-1, e, 1, work( ietgk+1 ), 2 )
452 CALL dstevx( 'N', 'I', n*2, work( idtgk ), work( ietgk ),
453 $ vutgk, vutgk, iutgk, iutgk, abstol, ns, s,
454 $ z, ldz, work( itemp ), iwork( iiwork ),
455 $ iwork( iifail ), info )
456 vutgk = s( 1 ) + fudge*smax*ulp*n
457 vutgk = min( vutgk, zero )
458*
459* If VLTGK=VUTGK, DSTEVX returns an error message,
460* so if needed we change VUTGK slightly.
461*
462 IF( vltgk.EQ.vutgk ) vltgk = vltgk - tol
463*
464 IF( wantz ) CALL dlaset( 'F', n*2, iu-il+1, zero, zero, z,
465 $ ldz)
466 END IF
467*
468* Initialize variables and pointers for S, Z, and WORK.
469*
470* NRU, NRV: number of rows in U and V for the active submatrix
471* IDBEG, ISBEG: offsets for the entries of D and S
472* IROWZ, ICOLZ: offsets for the rows and columns of Z
473* IROWU, IROWV: offsets for the rows of U and V
474*
475 ns = 0
476 nru = 0
477 nrv = 0
478 idbeg = 1
479 isbeg = 1
480 irowz = 1
481 icolz = 1
482 irowu = 2
483 irowv = 1
484 split = .false.
485 sveq0 = .false.
486*
487* Form the tridiagonal TGK matrix.
488*
489 s( 1:n ) = zero
490 work( ietgk+2*n-1 ) = zero
491 work( idtgk:idtgk+2*n-1 ) = zero
492 CALL dcopy( n, d, 1, work( ietgk ), 2 )
493 CALL dcopy( n-1, e, 1, work( ietgk+1 ), 2 )
494*
495*
496* Check for splits in two levels, outer level
497* in E and inner level in D.
498*
499 DO ieptr = 2, n*2, 2
500 IF( work( ietgk+ieptr-1 ).EQ.zero ) THEN
501*
502* Split in E (this piece of B is square) or bottom
503* of the (input bidiagonal) matrix.
504*
505 isplt = idbeg
506 idend = ieptr - 1
507 DO idptr = idbeg, idend, 2
508 IF( work( ietgk+idptr-1 ).EQ.zero ) THEN
509*
510* Split in D (rectangular submatrix). Set the number
511* of rows in U and V (NRU and NRV) accordingly.
512*
513 IF( idptr.EQ.idbeg ) THEN
514*
515* D=0 at the top.
516*
517 sveq0 = .true.
518 IF( idbeg.EQ.idend) THEN
519 nru = 1
520 nrv = 1
521 END IF
522 ELSE IF( idptr.EQ.idend ) THEN
523*
524* D=0 at the bottom.
525*
526 sveq0 = .true.
527 nru = (idend-isplt)/2 + 1
528 nrv = nru
529 IF( isplt.NE.idbeg ) THEN
530 nru = nru + 1
531 END IF
532 ELSE
533 IF( isplt.EQ.idbeg ) THEN
534*
535* Split: top rectangular submatrix.
536*
537 nru = (idptr-idbeg)/2
538 nrv = nru + 1
539 ELSE
540*
541* Split: middle square submatrix.
542*
543 nru = (idptr-isplt)/2 + 1
544 nrv = nru
545 END IF
546 END IF
547 ELSE IF( idptr.EQ.idend ) THEN
548*
549* Last entry of D in the active submatrix.
550*
551 IF( isplt.EQ.idbeg ) THEN
552*
553* No split (trivial case).
554*
555 nru = (idend-idbeg)/2 + 1
556 nrv = nru
557 ELSE
558*
559* Split: bottom rectangular submatrix.
560*
561 nrv = (idend-isplt)/2 + 1
562 nru = nrv + 1
563 END IF
564 END IF
565*
566 ntgk = nru + nrv
567*
568 IF( ntgk.GT.0 ) THEN
569*
570* Compute eigenvalues/vectors of the active
571* submatrix according to RANGE:
572* if RANGE='A' (ALLSV) then RNGVX = 'I'
573* if RANGE='V' (VALSV) then RNGVX = 'V'
574* if RANGE='I' (INDSV) then RNGVX = 'V'
575*
576 iltgk = 1
577 iutgk = ntgk / 2
578 IF( allsv .OR. vutgk.EQ.zero ) THEN
579 IF( sveq0 .OR.
580 $ smin.LT.eps .OR.
581 $ mod(ntgk,2).GT.0 ) THEN
582* Special case: eigenvalue equal to zero or very
583* small, additional eigenvector is needed.
584 iutgk = iutgk + 1
585 END IF
586 END IF
587*
588* Workspace needed by DSTEVX:
589* WORK( ITEMP: ): 2*5*NTGK
590* IWORK( 1: ): 2*6*NTGK
591*
592 CALL dstevx( jobz, rngvx, ntgk,
593 $ work( idtgk+isplt-1 ),
594 $ work( ietgk+isplt-1 ), vltgk, vutgk,
595 $ iltgk, iutgk, abstol, nsl, s( isbeg ),
596 $ z( irowz,icolz ), ldz, work( itemp ),
597 $ iwork( iiwork ), iwork( iifail ),
598 $ info )
599 IF( info.NE.0 ) THEN
600* Exit with the error code from DSTEVX.
601 RETURN
602 END IF
603 emin = abs( maxval( s( isbeg:isbeg+nsl-1 ) ) )
604*
605 IF( nsl.GT.0 .AND. wantz ) THEN
606*
607* Normalize u=Z([2,4,...],:) and v=Z([1,3,...],:),
608* changing the sign of v as discussed in the leading
609* comments. The norms of u and v may be (slightly)
610* different from 1/sqrt(2) if the corresponding
611* eigenvalues are very small or too close. We check
612* those norms and, if needed, reorthogonalize the
613* vectors.
614*
615 IF( nsl.GT.1 .AND.
616 $ vutgk.EQ.zero .AND.
617 $ mod(ntgk,2).EQ.0 .AND.
618 $ emin.EQ.0 .AND. .NOT.split ) THEN
619*
620* D=0 at the top or bottom of the active submatrix:
621* one eigenvalue is equal to zero; concatenate the
622* eigenvectors corresponding to the two smallest
623* eigenvalues.
624*
625 z( irowz:irowz+ntgk-1,icolz+nsl-2 ) =
626 $ z( irowz:irowz+ntgk-1,icolz+nsl-2 ) +
627 $ z( irowz:irowz+ntgk-1,icolz+nsl-1 )
628 z( irowz:irowz+ntgk-1,icolz+nsl-1 ) =
629 $ zero
630* IF( IUTGK*2.GT.NTGK ) THEN
631* Eigenvalue equal to zero or very small.
632* NSL = NSL - 1
633* END IF
634 END IF
635*
636 DO i = 0, min( nsl-1, nru-1 )
637 nrmu = dnrm2( nru, z( irowu, icolz+i ), 2 )
638 IF( nrmu.EQ.zero ) THEN
639 info = n*2 + 1
640 RETURN
641 END IF
642 CALL dscal( nru, one/nrmu,
643 $ z( irowu,icolz+i ), 2 )
644 IF( nrmu.NE.one .AND.
645 $ abs( nrmu-ortol )*sqrt2.GT.one )
646 $ THEN
647 DO j = 0, i-1
648 zjtji = -ddot( nru, z( irowu,
649 $ icolz+j ),
650 $ 2, z( irowu, icolz+i ), 2 )
651 CALL daxpy( nru, zjtji,
652 $ z( irowu, icolz+j ), 2,
653 $ z( irowu, icolz+i ), 2 )
654 END DO
655 nrmu = dnrm2( nru, z( irowu, icolz+i ),
656 $ 2 )
657 CALL dscal( nru, one/nrmu,
658 $ z( irowu,icolz+i ), 2 )
659 END IF
660 END DO
661 DO i = 0, min( nsl-1, nrv-1 )
662 nrmv = dnrm2( nrv, z( irowv, icolz+i ), 2 )
663 IF( nrmv.EQ.zero ) THEN
664 info = n*2 + 1
665 RETURN
666 END IF
667 CALL dscal( nrv, -one/nrmv,
668 $ z( irowv,icolz+i ), 2 )
669 IF( nrmv.NE.one .AND.
670 $ abs( nrmv-ortol )*sqrt2.GT.one )
671 $ THEN
672 DO j = 0, i-1
673 zjtji = -ddot( nrv, z( irowv,
674 $ icolz+j ),
675 $ 2, z( irowv, icolz+i ), 2 )
676 CALL daxpy( nru, zjtji,
677 $ z( irowv, icolz+j ), 2,
678 $ z( irowv, icolz+i ), 2 )
679 END DO
680 nrmv = dnrm2( nrv, z( irowv, icolz+i ),
681 $ 2 )
682 CALL dscal( nrv, one/nrmv,
683 $ z( irowv,icolz+i ), 2 )
684 END IF
685 END DO
686 IF( vutgk.EQ.zero .AND.
687 $ idptr.LT.idend .AND.
688 $ mod(ntgk,2).GT.0 ) THEN
689*
690* D=0 in the middle of the active submatrix (one
691* eigenvalue is equal to zero): save the corresponding
692* eigenvector for later use (when bottom of the
693* active submatrix is reached).
694*
695 split = .true.
696 z( irowz:irowz+ntgk-1,n+1 ) =
697 $ z( irowz:irowz+ntgk-1,ns+nsl )
698 z( irowz:irowz+ntgk-1,ns+nsl ) =
699 $ zero
700 END IF
701 END IF !** WANTZ **!
702*
703 nsl = min( nsl, nru )
704 sveq0 = .false.
705*
706* Absolute values of the eigenvalues of TGK.
707*
708 DO i = 0, nsl-1
709 s( isbeg+i ) = abs( s( isbeg+i ) )
710 END DO
711*
712* Update pointers for TGK, S and Z.
713*
714 isbeg = isbeg + nsl
715 irowz = irowz + ntgk
716 icolz = icolz + nsl
717 irowu = irowz
718 irowv = irowz + 1
719 isplt = idptr + 1
720 ns = ns + nsl
721 nru = 0
722 nrv = 0
723 END IF !** NTGK.GT.0 **!
724 IF( irowz.LT.n*2 .AND. wantz ) THEN
725 z( 1:irowz-1, icolz ) = zero
726 END IF
727 END DO !** IDPTR loop **!
728 IF( split .AND. wantz ) THEN
729*
730* Bring back eigenvector corresponding
731* to eigenvalue equal to zero.
732*
733 z( idbeg:idend-ntgk+1,isbeg-1 ) =
734 $ z( idbeg:idend-ntgk+1,isbeg-1 ) +
735 $ z( idbeg:idend-ntgk+1,n+1 )
736 z( idbeg:idend-ntgk+1,n+1 ) = 0
737 END IF
738 irowv = irowv - 1
739 irowu = irowu + 1
740 idbeg = ieptr + 1
741 sveq0 = .false.
742 split = .false.
743 END IF !** Check for split in E **!
744 END DO !** IEPTR loop **!
745*
746* Sort the singular values into decreasing order (insertion sort on
747* singular values, but only one transposition per singular vector)
748*
749 DO i = 1, ns-1
750 k = 1
751 smin = s( 1 )
752 DO j = 2, ns + 1 - i
753 IF( s( j ).LE.smin ) THEN
754 k = j
755 smin = s( j )
756 END IF
757 END DO
758 IF( k.NE.ns+1-i ) THEN
759 s( k ) = s( ns+1-i )
760 s( ns+1-i ) = smin
761 IF( wantz ) CALL dswap( n*2, z( 1,k ), 1, z( 1,ns+1-i ),
762 $ 1 )
763 END IF
764 END DO
765*
766* If RANGE=I, check for singular values/vectors to be discarded.
767*
768 IF( indsv ) THEN
769 k = iu - il + 1
770 IF( k.LT.ns ) THEN
771 s( k+1:ns ) = zero
772 IF( wantz ) z( 1:n*2,k+1:ns ) = zero
773 ns = k
774 END IF
775 END IF
776*
777* Reorder Z: U = Z( 1:N,1:NS ), V = Z( N+1:N*2,1:NS ).
778* If B is a lower diagonal, swap U and V.
779*
780 IF( wantz ) THEN
781 DO i = 1, ns
782 CALL dcopy( n*2, z( 1,i ), 1, work, 1 )
783 IF( lower ) THEN
784 CALL dcopy( n, work( 2 ), 2, z( n+1,i ), 1 )
785 CALL dcopy( n, work( 1 ), 2, z( 1 ,i ), 1 )
786 ELSE
787 CALL dcopy( n, work( 2 ), 2, z( 1 ,i ), 1 )
788 CALL dcopy( n, work( 1 ), 2, z( n+1,i ), 1 )
789 END IF
790 END DO
791 END IF
792*
793 RETURN
794*
795* End of DBDSVDX
796*
797 END
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
daxpy
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
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
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
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:108
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
dstevx
subroutine dstevx(jobz, range, n, d, e, vl, vu, il, iu, abstol, m, w, z, ldz, work, iwork, ifail, info)
DSTEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matrice...
Definition dstevx.f:226
dswap
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82