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