LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ sbdsvdx()

subroutine sbdsvdx ( character  UPLO,
character  JOBZ,
character  RANGE,
integer  N,
real, dimension( * )  D,
real, dimension( * )  E,
real  VL,
real  VU,
integer  IL,
integer  IU,
integer  NS,
real, dimension( * )  S,
real, dimension( ldz, * )  Z,
integer  LDZ,
real, dimension( * )  WORK,
integer, dimension( * )  IWORK,
integer  INFO 
)

SBDSVDX

Download SBDSVDX + dependencies [TGZ] [ZIP] [TXT]

Purpose:
  SBDSVDX computes the singular value decomposition (SVD) of a real
  N-by-N (upper or lower) bidiagonal matrix B, B = U * S * VT,
  where S is a diagonal matrix with non-negative diagonal elements
  (the singular values of B), and U and VT are orthogonal matrices
  of left and right singular vectors, respectively.

  Given an upper bidiagonal B with diagonal D = [ d_1 d_2 ... d_N ]
  and superdiagonal E = [ e_1 e_2 ... e_N-1 ], SBDSVDX computes the
  singular value decompositon of B through the eigenvalues and
  eigenvectors of the N*2-by-N*2 tridiagonal matrix

        |  0  d_1                |
        | d_1  0  e_1            |
  TGK = |     e_1  0  d_2        |
        |         d_2  .   .     |
        |              .   .   . |

  If (s,u,v) is a singular triplet of B with ||u|| = ||v|| = 1, then
  (+/-s,q), ||q|| = 1, are eigenpairs of TGK, with q = P * ( u' +/-v' ) /
  sqrt(2) = ( v_1 u_1 v_2 u_2 ... v_n u_n ) / sqrt(2), and
  P = [ e_{n+1} e_{1} e_{n+2} e_{2} ... ].

  Given a TGK matrix, one can either a) compute -s,-v and change signs
  so that the singular values (and corresponding vectors) are already in
  descending order (as in SGESVD/SGESDD) or b) compute s,v and reorder
  the values (and corresponding vectors). SBDSVDX implements a) by
  calling SSTEVX (bisection plus inverse iteration, to be replaced
  with a version of the Multiple Relative Robust Representation
  algorithm. (See P. Willems and B. Lang, A framework for the MR^3
  algorithm: theory and implementation, SIAM J. Sci. Comput.,
  35:740-766, 2013.)
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          = 'U':  B is upper bidiagonal;
          = 'L':  B is lower bidiagonal.
[in]JOBZ
          JOBZ is CHARACTER*1
          = 'N':  Compute singular values only;
          = 'V':  Compute singular values and singular vectors.
[in]RANGE
          RANGE is CHARACTER*1
          = 'A': all singular values will be found.
          = 'V': all singular values in the half-open interval [VL,VU)
                 will be found.
          = 'I': the IL-th through IU-th singular values will be found.
[in]N
          N is INTEGER
          The order of the bidiagonal matrix.  N >= 0.
[in]D
          D is REAL array, dimension (N)
          The n diagonal elements of the bidiagonal matrix B.
[in]E
          E is REAL array, dimension (max(1,N-1))
          The (n-1) superdiagonal elements of the bidiagonal matrix
          B in elements 1 to N-1.
[in]VL
         VL is REAL
          If RANGE='V', the lower bound of the interval to
          be searched for singular values. VU > VL.
          Not referenced if RANGE = 'A' or 'I'.
[in]VU
         VU is REAL
          If RANGE='V', the upper bound of the interval to
          be searched for singular values. VU > VL.
          Not referenced if RANGE = 'A' or 'I'.
[in]IL
          IL is INTEGER
          If RANGE='I', the index of the
          smallest singular value to be returned.
          1 <= IL <= IU <= min(M,N), if min(M,N) > 0.
          Not referenced if RANGE = 'A' or 'V'.
[in]IU
          IU is INTEGER
          If RANGE='I', the index of the
          largest singular value to be returned.
          1 <= IL <= IU <= min(M,N), if min(M,N) > 0.
          Not referenced if RANGE = 'A' or 'V'.
[out]NS
          NS is INTEGER
          The total number of singular values found.  0 <= NS <= N.
          If RANGE = 'A', NS = N, and if RANGE = 'I', NS = IU-IL+1.
[out]S
          S is REAL array, dimension (N)
          The first NS elements contain the selected singular values in
          ascending order.
[out]Z
          Z is REAL array, dimension (2*N,K)
          If JOBZ = 'V', then if INFO = 0 the first NS columns of Z
          contain the singular vectors of the matrix B corresponding to
          the selected singular values, with U in rows 1 to N and V
          in rows N+1 to N*2, i.e.
          Z = [ U ]
              [ V ]
          If JOBZ = 'N', then Z is not referenced.
          Note: The user must ensure that at least K = NS+1 columns are
          supplied in the array Z; if RANGE = 'V', the exact value of
          NS is not known in advance and an upper bound must be used.
[in]LDZ
          LDZ is INTEGER
          The leading dimension of the array Z. LDZ >= 1, and if
          JOBZ = 'V', LDZ >= max(2,N*2).
[out]WORK
          WORK is REAL array, dimension (14*N)
[out]IWORK
          IWORK is INTEGER array, dimension (12*N)
          If JOBZ = 'V', then if INFO = 0, the first NS elements of
          IWORK are zero. If INFO > 0, then IWORK contains the indices
          of the eigenvectors that failed to converge in DSTEVX.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
          > 0:  if INFO = i, then i eigenvectors failed to converge
                   in SSTEVX. The indices of the eigenvectors
                   (as returned by SSTEVX) are stored in the
                   array IWORK.
                if INFO = N*2 + 1, an internal error occurred.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 224 of file sbdsvdx.f.

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 *
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
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:71
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
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
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
real function sdot(N, SX, INCX, SY, INCY)
SDOT
Definition: sdot.f:82
real(wp) function snrm2(n, x, incx)
SNRM2
Definition: snrm2.f90:89
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
Definition: saxpy.f:89
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: