LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine dbdsvdx ( character  UPLO,
character  JOBZ,
character  RANGE,
integer  N,
double precision, dimension( * )  D,
double precision, dimension( * )  E,
double precision  VL,
double precision  VU,
integer  IL,
integer  IU,
integer  NS,
double precision, dimension( * )  S,
double precision, dimension( ldz, * )  Z,
integer  LDZ,
double precision, dimension( * )  WORK,
integer, dimension( * )  IWORK,
integer  INFO 
)

DBDSVDX

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

Purpose:
  DBDSVDX 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 ], DBDSVDX 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 DGESVD/DGESDD) or b) compute s,v and reorder 
  the values (and corresponding vectors). DBDSVDX implements a) by 
  calling DSTEVX (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 DOUBLE PRECISION array, dimension (N)
          The n diagonal elements of the bidiagonal matrix B.
[in]E
          E is DOUBLE PRECISION 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 DOUBLE PRECISION
          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 DOUBLE PRECISION
          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 DOUBLE PRECISION array, dimension (N)
          The first NS elements contain the selected singular values in
          ascending order.
[out]Z
          Z is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DSTEVX. The indices of the eigenvectors
                   (as returned by DSTEVX) 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.
Date
June 2016

Definition at line 228 of file dbdsvdx.f.

228 *
229 * -- LAPACK driver routine (version 3.6.1) --
230 * -- LAPACK is a software package provided by Univ. of Tennessee, --
231 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
232 * November 2016
233 *
234 * .. Scalar Arguments ..
235  CHARACTER jobz, range, uplo
236  INTEGER il, info, iu, ldz, n, ns
237  DOUBLE PRECISION vl, vu
238 * ..
239 * .. Array Arguments ..
240  INTEGER iwork( * )
241  DOUBLE PRECISION d( * ), e( * ), s( * ), work( * ),
242  $ z( ldz, * )
243 * ..
244 *
245 * =====================================================================
246 *
247 * .. Parameters ..
248  DOUBLE PRECISION zero, one, ten, hndrd, meigth
249  parameter ( zero = 0.0d0, one = 1.0d0, ten = 10.0d0,
250  $ hndrd = 100.0d0, meigth = -0.1250d0 )
251  DOUBLE PRECISION fudge
252  parameter ( fudge = 2.0d0 )
253 * ..
254 * .. Local Scalars ..
255  CHARACTER rngvx
256  LOGICAL allsv, indsv, lower, split, sveq0, valsv, wantz
257  INTEGER i, icolz, idbeg, idend, idtgk, idptr, ieptr,
258  $ ietgk, iifail, iiwork, iltgk, irowu, irowv,
259  $ irowz, isbeg, isplt, itemp, iutgk, j, k,
260  $ ntgk, nru, nrv, nsl
261  DOUBLE PRECISION abstol, eps, emin, mu, nrmu, nrmv, ortol, smax,
262  $ smin, sqrt2, thresh, tol, ulp,
263  $ vltgk, vutgk, zjtji
264 * ..
265 * .. External Functions ..
266  LOGICAL lsame
267  INTEGER idamax
268  DOUBLE PRECISION ddot, dlamch, dnrm2
269  EXTERNAL idamax, lsame, daxpy, ddot, dlamch, dnrm2
270 * ..
271 * .. External Subroutines ..
272  EXTERNAL dcopy, dlaset, dscal, dswap
273 * ..
274 * .. Intrinsic Functions ..
275  INTRINSIC abs, dble, sign, sqrt
276 * ..
277 * .. Executable Statements ..
278 *
279 * Test the input parameters.
280 *
281  allsv = lsame( range, 'A' )
282  valsv = lsame( range, 'V' )
283  indsv = lsame( range, 'I' )
284  wantz = lsame( jobz, 'V' )
285  lower = lsame( uplo, 'L' )
286 *
287  info = 0
288  IF( .NOT.lsame( uplo, 'U' ) .AND. .NOT.lower ) THEN
289  info = -1
290  ELSE IF( .NOT.( wantz .OR. lsame( jobz, 'N' ) ) ) THEN
291  info = -2
292  ELSE IF( .NOT.( allsv .OR. valsv .OR. indsv ) ) THEN
293  info = -3
294  ELSE IF( n.LT.0 ) THEN
295  info = -4
296  ELSE IF( n.GT.0 ) THEN
297  IF( valsv ) THEN
298  IF( vl.LT.zero ) THEN
299  info = -7
300  ELSE IF( vu.LE.vl ) THEN
301  info = -8
302  END IF
303  ELSE IF( indsv ) THEN
304  IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
305  info = -9
306  ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
307  info = -10
308  END IF
309  END IF
310  END IF
311  IF( info.EQ.0 ) THEN
312  IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n*2 ) ) info = -14
313  END IF
314 *
315  IF( info.NE.0 ) THEN
316  CALL xerbla( 'DBDSVDX', -info )
317  RETURN
318  END IF
319 *
320 * Quick return if possible (N.LE.1)
321 *
322  ns = 0
323  IF( n.EQ.0 ) RETURN
324 *
325  IF( n.EQ.1 ) THEN
326  IF( allsv .OR. indsv ) THEN
327  ns = 1
328  s( 1 ) = abs( d( 1 ) )
329  ELSE
330  IF( vl.LT.abs( d( 1 ) ) .AND. vu.GE.abs( d( 1 ) ) ) THEN
331  ns = 1
332  s( 1 ) = abs( d( 1 ) )
333  END IF
334  END IF
335  IF( wantz ) THEN
336  z( 1, 1 ) = sign( one, d( 1 ) )
337  z( 2, 1 ) = one
338  END IF
339  RETURN
340  END IF
341 *
342  abstol = 2*dlamch( 'Safe Minimum' )
343  ulp = dlamch( 'Precision' )
344  eps = dlamch( 'Epsilon' )
345  sqrt2 = sqrt( 2.0d0 )
346  ortol = sqrt( ulp )
347 *
348 * Criterion for splitting is taken from DBDSQR when singular
349 * values are computed to relative accuracy TOL. (See J. Demmel and
350 * W. Kahan, Accurate singular values of bidiagonal matrices, SIAM
351 * J. Sci. and Stat. Comput., 11:873–912, 1990.)
352 *
353  tol = max( ten, min( hndrd, eps**meigth ) )*eps
354 *
355 * Compute approximate maximum, minimum singular values.
356 *
357  i = idamax( n, d, 1 )
358  smax = abs( d( i ) )
359  i = idamax( n-1, e, 1 )
360  smax = max( smax, abs( e( i ) ) )
361 *
362 * Compute threshold for neglecting D's and E's.
363 *
364  smin = abs( d( 1 ) )
365  IF( smin.NE.zero ) THEN
366  mu = smin
367  DO i = 2, n
368  mu = abs( d( i ) )*( mu / ( mu+abs( e( i-1 ) ) ) )
369  smin = min( smin, mu )
370  IF( smin.EQ.zero ) EXIT
371  END DO
372  END IF
373  smin = smin / sqrt( dble( n ) )
374  thresh = tol*smin
375 *
376 * Check for zeros in D and E (splits), i.e. submatrices.
377 *
378  DO i = 1, n-1
379  IF( abs( d( i ) ).LE.thresh ) d( i ) = zero
380  IF( abs( e( i ) ).LE.thresh ) e( i ) = zero
381  END DO
382  IF( abs( d( n ) ).LE.thresh ) d( n ) = zero
383 *
384 * Pointers for arrays used by DSTEVX.
385 *
386  idtgk = 1
387  ietgk = idtgk + n*2
388  itemp = ietgk + n*2
389  iifail = 1
390  iiwork = iifail + n*2
391 *
392 * Set RNGVX, which corresponds to RANGE for DSTEVX in TGK mode.
393 * VL,VU or IL,IU are redefined to conform to implementation a)
394 * described in the leading comments.
395 *
396  iltgk = 0
397  iutgk = 0
398  vltgk = zero
399  vutgk = zero
400 *
401  IF( allsv ) THEN
402 *
403 * All singular values will be found. We aim at -s (see
404 * leading comments) with RNGVX = 'I'. IL and IU are set
405 * later (as ILTGK and IUTGK) according to the dimension
406 * of the active submatrix.
407 *
408  rngvx = 'I'
409  IF( wantz ) CALL dlaset( 'F', n*2, n+1, zero, zero, z, ldz )
410  ELSE IF( valsv ) THEN
411 *
412 * Find singular values in a half-open interval. We aim
413 * at -s (see leading comments) and we swap VL and VU
414 * (as VUTGK and VLTGK), changing their signs.
415 *
416  rngvx = 'V'
417  vltgk = -vu
418  vutgk = -vl
419  work( idtgk:idtgk+2*n-1 ) = zero
420  CALL dcopy( n, d, 1, work( ietgk ), 2 )
421  CALL dcopy( n-1, e, 1, work( ietgk+1 ), 2 )
422  CALL dstevx( 'N', 'V', n*2, work( idtgk ), work( ietgk ),
423  $ vltgk, vutgk, iltgk, iltgk, abstol, ns, s,
424  $ z, ldz, work( itemp ), iwork( iiwork ),
425  $ iwork( iifail ), info )
426  IF( ns.EQ.0 ) THEN
427  RETURN
428  ELSE
429  IF( wantz ) CALL dlaset( 'F', n*2, ns, zero, zero, z, ldz )
430  END IF
431  ELSE IF( indsv ) THEN
432 *
433 * Find the IL-th through the IU-th singular values. We aim
434 * at -s (see leading comments) and indices are mapped into
435 * values, therefore mimicking DSTEBZ, where
436 *
437 * GL = GL - FUDGE*TNORM*ULP*N - FUDGE*TWO*PIVMIN
438 * GU = GU + FUDGE*TNORM*ULP*N + FUDGE*PIVMIN
439 *
440  iltgk = il
441  iutgk = iu
442  rngvx = 'V'
443  work( idtgk:idtgk+2*n-1 ) = zero
444  CALL dcopy( n, d, 1, work( ietgk ), 2 )
445  CALL dcopy( n-1, e, 1, work( ietgk+1 ), 2 )
446  CALL dstevx( 'N', 'I', n*2, work( idtgk ), work( ietgk ),
447  $ vltgk, vltgk, iltgk, iltgk, abstol, ns, s,
448  $ z, ldz, work( itemp ), iwork( iiwork ),
449  $ iwork( iifail ), info )
450  vltgk = s( 1 ) - fudge*smax*ulp*n
451  work( idtgk:idtgk+2*n-1 ) = zero
452  CALL dcopy( n, d, 1, work( ietgk ), 2 )
453  CALL dcopy( n-1, e, 1, work( ietgk+1 ), 2 )
454  CALL dstevx( 'N', 'I', n*2, work( idtgk ), work( ietgk ),
455  $ vutgk, vutgk, iutgk, iutgk, abstol, ns, s,
456  $ z, ldz, work( itemp ), iwork( iiwork ),
457  $ iwork( iifail ), info )
458  vutgk = s( 1 ) + fudge*smax*ulp*n
459  vutgk = min( vutgk, zero )
460 *
461 * If VLTGK=VUTGK, DSTEVX returns an error message,
462 * so if needed we change VUTGK slightly.
463 *
464  IF( vltgk.EQ.vutgk ) vltgk = vltgk - tol
465 *
466  IF( wantz ) CALL dlaset( 'F', n*2, iu-il+1, zero, zero, z, ldz)
467  END IF
468 *
469 * Initialize variables and pointers for S, Z, and WORK.
470 *
471 * NRU, NRV: number of rows in U and V for the active submatrix
472 * IDBEG, ISBEG: offsets for the entries of D and S
473 * IROWZ, ICOLZ: offsets for the rows and columns of Z
474 * IROWU, IROWV: offsets for the rows of U and V
475 *
476  ns = 0
477  nru = 0
478  nrv = 0
479  idbeg = 1
480  isbeg = 1
481  irowz = 1
482  icolz = 1
483  irowu = 2
484  irowv = 1
485  split = .false.
486  sveq0 = .false.
487 *
488 * Form the tridiagonal TGK matrix.
489 *
490  s( 1:n ) = zero
491  work( ietgk+2*n-1 ) = zero
492  work( idtgk:idtgk+2*n-1 ) = zero
493  CALL dcopy( n, d, 1, work( ietgk ), 2 )
494  CALL dcopy( n-1, e, 1, work( ietgk+1 ), 2 )
495 *
496 *
497 * Check for splits in two levels, outer level
498 * in E and inner level in D.
499 *
500  DO ieptr = 2, n*2, 2
501  IF( work( ietgk+ieptr-1 ).EQ.zero ) THEN
502 *
503 * Split in E (this piece of B is square) or bottom
504 * of the (input bidiagonal) matrix.
505 *
506  isplt = idbeg
507  idend = ieptr - 1
508  DO idptr = idbeg, idend, 2
509  IF( work( ietgk+idptr-1 ).EQ.zero ) THEN
510 *
511 * Split in D (rectangular submatrix). Set the number
512 * of rows in U and V (NRU and NRV) accordingly.
513 *
514  IF( idptr.EQ.idbeg ) THEN
515 *
516 * D=0 at the top.
517 *
518  sveq0 = .true.
519  IF( idbeg.EQ.idend) THEN
520  nru = 1
521  nrv = 1
522  END IF
523  ELSE IF( idptr.EQ.idend ) THEN
524 *
525 * D=0 at the bottom.
526 *
527  sveq0 = .true.
528  nru = (idend-isplt)/2 + 1
529  nrv = nru
530  IF( isplt.NE.idbeg ) THEN
531  nru = nru + 1
532  END IF
533  ELSE
534  IF( isplt.EQ.idbeg ) THEN
535 *
536 * Split: top rectangular submatrix.
537 *
538  nru = (idptr-idbeg)/2
539  nrv = nru + 1
540  ELSE
541 *
542 * Split: middle square submatrix.
543 *
544  nru = (idptr-isplt)/2 + 1
545  nrv = nru
546  END IF
547  END IF
548  ELSE IF( idptr.EQ.idend ) THEN
549 *
550 * Last entry of D in the active submatrix.
551 *
552  IF( isplt.EQ.idbeg ) THEN
553 *
554 * No split (trivial case).
555 *
556  nru = (idend-idbeg)/2 + 1
557  nrv = nru
558  ELSE
559 *
560 * Split: bottom rectangular submatrix.
561 *
562  nrv = (idend-isplt)/2 + 1
563  nru = nrv + 1
564  END IF
565  END IF
566 *
567  ntgk = nru + nrv
568 *
569  IF( ntgk.GT.0 ) THEN
570 *
571 * Compute eigenvalues/vectors of the active
572 * submatrix according to RANGE:
573 * if RANGE='A' (ALLSV) then RNGVX = 'I'
574 * if RANGE='V' (VALSV) then RNGVX = 'V'
575 * if RANGE='I' (INDSV) then RNGVX = 'V'
576 *
577  iltgk = 1
578  iutgk = ntgk / 2
579  IF( allsv .OR. vutgk.EQ.zero ) THEN
580  IF( sveq0 .OR.
581  $ smin.LT.eps .OR.
582  $ mod(ntgk,2).GT.0 ) THEN
583 * Special case: eigenvalue equal to zero or very
584 * small, additional eigenvector is needed.
585  iutgk = iutgk + 1
586  END IF
587  END IF
588 *
589 * Workspace needed by DSTEVX:
590 * WORK( ITEMP: ): 2*5*NTGK
591 * IWORK( 1: ): 2*6*NTGK
592 *
593  CALL dstevx( jobz, rngvx, ntgk, 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, icolz+j ),
649  $ 2, z( irowu, icolz+i ), 2 )
650  CALL daxpy( nru, zjtji,
651  $ z( irowu, icolz+j ), 2,
652  $ z( irowu, icolz+i ), 2 )
653  END DO
654  nrmu = dnrm2( nru, z( irowu, icolz+i ), 2 )
655  CALL dscal( nru, one/nrmu,
656  $ z( irowu,icolz+i ), 2 )
657  END IF
658  END DO
659  DO i = 0, min( nsl-1, nrv-1 )
660  nrmv = dnrm2( nrv, z( irowv, icolz+i ), 2 )
661  IF( nrmv.EQ.zero ) THEN
662  info = n*2 + 1
663  RETURN
664  END IF
665  CALL dscal( nrv, -one/nrmv,
666  $ z( irowv,icolz+i ), 2 )
667  IF( nrmv.NE.one .AND.
668  $ abs( nrmv-ortol )*sqrt2.GT.one )
669  $ THEN
670  DO j = 0, i-1
671  zjtji = -ddot( nrv, z( irowv, icolz+j ),
672  $ 2, z( irowv, icolz+i ), 2 )
673  CALL daxpy( nru, zjtji,
674  $ z( irowv, icolz+j ), 2,
675  $ z( irowv, icolz+i ), 2 )
676  END DO
677  nrmv = dnrm2( nrv, z( irowv, icolz+i ), 2 )
678  CALL dscal( nrv, one/nrmv,
679  $ z( irowv,icolz+i ), 2 )
680  END IF
681  END DO
682  IF( vutgk.EQ.zero .AND.
683  $ idptr.LT.idend .AND.
684  $ mod(ntgk,2).GT.0 ) THEN
685 *
686 * D=0 in the middle of the active submatrix (one
687 * eigenvalue is equal to zero): save the corresponding
688 * eigenvector for later use (when bottom of the
689 * active submatrix is reached).
690 *
691  split = .true.
692  z( irowz:irowz+ntgk-1,n+1 ) =
693  $ z( irowz:irowz+ntgk-1,ns+nsl )
694  z( irowz:irowz+ntgk-1,ns+nsl ) =
695  $ zero
696  END IF
697  END IF !** WANTZ **!
698 *
699  nsl = min( nsl, nru )
700  sveq0 = .false.
701 *
702 * Absolute values of the eigenvalues of TGK.
703 *
704  DO i = 0, nsl-1
705  s( isbeg+i ) = abs( s( isbeg+i ) )
706  END DO
707 *
708 * Update pointers for TGK, S and Z.
709 *
710  isbeg = isbeg + nsl
711  irowz = irowz + ntgk
712  icolz = icolz + nsl
713  irowu = irowz
714  irowv = irowz + 1
715  isplt = idptr + 1
716  ns = ns + nsl
717  nru = 0
718  nrv = 0
719  END IF !** NTGK.GT.0 **!
720  IF( irowz.LT.n*2 .AND. wantz ) THEN
721  z( 1:irowz-1, icolz ) = zero
722  END IF
723  END DO !** IDPTR loop **!
724  IF( split .AND. wantz ) THEN
725 *
726 * Bring back eigenvector corresponding
727 * to eigenvalue equal to zero.
728 *
729  z( idbeg:idend-ntgk+1,isbeg-1 ) =
730  $ z( idbeg:idend-ntgk+1,isbeg-1 ) +
731  $ z( idbeg:idend-ntgk+1,n+1 )
732  z( idbeg:idend-ntgk+1,n+1 ) = 0
733  END IF
734  irowv = irowv - 1
735  irowu = irowu + 1
736  idbeg = ieptr + 1
737  sveq0 = .false.
738  split = .false.
739  END IF !** Check for split in E **!
740  END DO !** IEPTR loop **!
741 *
742 * Sort the singular values into decreasing order (insertion sort on
743 * singular values, but only one transposition per singular vector)
744 *
745  DO i = 1, ns-1
746  k = 1
747  smin = s( 1 )
748  DO j = 2, ns + 1 - i
749  IF( s( j ).LE.smin ) THEN
750  k = j
751  smin = s( j )
752  END IF
753  END DO
754  IF( k.NE.ns+1-i ) THEN
755  s( k ) = s( ns+1-i )
756  s( ns+1-i ) = smin
757  IF( wantz ) CALL dswap( n*2, z( 1,k ), 1, z( 1,ns+1-i ), 1 )
758  END IF
759  END DO
760 *
761 * If RANGE=I, check for singular values/vectors to be discarded.
762 *
763  IF( indsv ) THEN
764  k = iu - il + 1
765  IF( k.LT.ns ) THEN
766  s( k+1:ns ) = zero
767  IF( wantz ) z( 1:n*2,k+1:ns ) = zero
768  ns = k
769  END IF
770  END IF
771 *
772 * Reorder Z: U = Z( 1:N,1:NS ), V = Z( N+1:N*2,1:NS ).
773 * If B is a lower diagonal, swap U and V.
774 *
775  IF( wantz ) THEN
776  DO i = 1, ns
777  CALL dcopy( n*2, z( 1,i ), 1, work, 1 )
778  IF( lower ) THEN
779  CALL dcopy( n, work( 2 ), 2, z( n+1,i ), 1 )
780  CALL dcopy( n, work( 1 ), 2, z( 1 ,i ), 1 )
781  ELSE
782  CALL dcopy( n, work( 2 ), 2, z( 1 ,i ), 1 )
783  CALL dcopy( n, work( 1 ), 2, z( n+1,i ), 1 )
784  END IF
785  END DO
786  END IF
787 *
788  RETURN
789 *
790 * End of DBDSVDX
791 *
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:112
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
integer function idamax(N, DX, INCX)
IDAMAX
Definition: idamax.f:53
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
double precision function ddot(N, DX, INCX, DY, INCY)
DDOT
Definition: ddot.f:53
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:54
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:53
double precision function dnrm2(N, X, INCX)
DNRM2
Definition: dnrm2.f:56
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55
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 matric...
Definition: dstevx.f:229
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function: