Return to dlasq1.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de lapack.

    1: *> \brief \b DLASQ1 computes the singular values of a real square bidiagonal matrix. Used by sbdsqr.
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at
    6: *            http://www.netlib.org/lapack/explore-html/
    7: *
    8: *> \htmlonly
    9: *> Download DLASQ1 + dependencies
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasq1.f">
   11: *> [TGZ]</a>
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasq1.f">
   13: *> [ZIP]</a>
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasq1.f">
   15: *> [TXT]</a>
   16: *> \endhtmlonly
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE DLASQ1( N, D, E, WORK, INFO )
   22: *
   23: *       .. Scalar Arguments ..
   24: *       INTEGER            INFO, N
   25: *       ..
   26: *       .. Array Arguments ..
   27: *       DOUBLE PRECISION   D( * ), E( * ), WORK( * )
   28: *       ..
   29: *
   30: *
   31: *> \par Purpose:
   32: *  =============
   33: *>
   34: *> \verbatim
   35: *>
   36: *> DLASQ1 computes the singular values of a real N-by-N bidiagonal
   37: *> matrix with diagonal D and off-diagonal E. The singular values
   38: *> are computed to high relative accuracy, in the absence of
   39: *> denormalization, underflow and overflow. The algorithm was first
   40: *> presented in
   41: *>
   42: *> "Accurate singular values and differential qd algorithms" by K. V.
   43: *> Fernando and B. N. Parlett, Numer. Math., Vol-67, No. 2, pp. 191-230,
   44: *> 1994,
   45: *>
   46: *> and the present implementation is described in "An implementation of
   47: *> the dqds Algorithm (Positive Case)", LAPACK Working Note.
   48: *> \endverbatim
   49: *
   50: *  Arguments:
   51: *  ==========
   52: *
   53: *> \param[in] N
   54: *> \verbatim
   55: *>          N is INTEGER
   56: *>        The number of rows and columns in the matrix. N >= 0.
   57: *> \endverbatim
   58: *>
   59: *> \param[in,out] D
   60: *> \verbatim
   61: *>          D is DOUBLE PRECISION array, dimension (N)
   62: *>        On entry, D contains the diagonal elements of the
   63: *>        bidiagonal matrix whose SVD is desired. On normal exit,
   64: *>        D contains the singular values in decreasing order.
   65: *> \endverbatim
   66: *>
   67: *> \param[in,out] E
   68: *> \verbatim
   69: *>          E is DOUBLE PRECISION array, dimension (N)
   70: *>        On entry, elements E(1:N-1) contain the off-diagonal elements
   71: *>        of the bidiagonal matrix whose SVD is desired.
   72: *>        On exit, E is overwritten.
   73: *> \endverbatim
   74: *>
   75: *> \param[out] WORK
   76: *> \verbatim
   77: *>          WORK is DOUBLE PRECISION array, dimension (4*N)
   78: *> \endverbatim
   79: *>
   80: *> \param[out] INFO
   81: *> \verbatim
   82: *>          INFO is INTEGER
   83: *>        = 0: successful exit
   84: *>        < 0: if INFO = -i, the i-th argument had an illegal value
   85: *>        > 0: the algorithm failed
   86: *>             = 1, a split was marked by a positive value in E
   87: *>             = 2, current block of Z not diagonalized after 100*N
   88: *>                  iterations (in inner while loop)  On exit D and E
   89: *>                  represent a matrix with the same singular values
   90: *>                  which the calling subroutine could use to finish the
   91: *>                  computation, or even feed back into DLASQ1
   92: *>             = 3, termination criterion of outer while loop not met
   93: *>                  (program created more than N unreduced blocks)
   94: *> \endverbatim
   95: *
   96: *  Authors:
   97: *  ========
   98: *
   99: *> \author Univ. of Tennessee
  100: *> \author Univ. of California Berkeley
  101: *> \author Univ. of Colorado Denver
  102: *> \author NAG Ltd.
  103: *
  104: *> \date December 2016
  105: *
  106: *> \ingroup auxOTHERcomputational
  107: *
  108: *  =====================================================================
  109:       SUBROUTINE DLASQ1( N, D, E, WORK, INFO )
  110: *
  111: *  -- LAPACK computational routine (version 3.7.0) --
  112: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  113: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  114: *     December 2016
  115: *
  116: *     .. Scalar Arguments ..
  117:       INTEGER            INFO, N
  118: *     ..
  119: *     .. Array Arguments ..
  120:       DOUBLE PRECISION   D( * ), E( * ), WORK( * )
  121: *     ..
  122: *
  123: *  =====================================================================
  124: *
  125: *     .. Parameters ..
  126:       DOUBLE PRECISION   ZERO
  127:       PARAMETER          ( ZERO = 0.0D0 )
  128: *     ..
  129: *     .. Local Scalars ..
  130:       INTEGER            I, IINFO
  131:       DOUBLE PRECISION   EPS, SCALE, SAFMIN, SIGMN, SIGMX
  132: *     ..
  133: *     .. External Subroutines ..
  134:       EXTERNAL           DCOPY, DLAS2, DLASCL, DLASQ2, DLASRT, XERBLA
  135: *     ..
  136: *     .. External Functions ..
  137:       DOUBLE PRECISION   DLAMCH
  138:       EXTERNAL           DLAMCH
  139: *     ..
  140: *     .. Intrinsic Functions ..
  141:       INTRINSIC          ABS, MAX, SQRT
  142: *     ..
  143: *     .. Executable Statements ..
  144: *
  145:       INFO = 0
  146:       IF( N.LT.0 ) THEN
  147:          INFO = -1
  148:          CALL XERBLA( 'DLASQ1', -INFO )
  149:          RETURN
  150:       ELSE IF( N.EQ.0 ) THEN
  151:          RETURN
  152:       ELSE IF( N.EQ.1 ) THEN
  153:          D( 1 ) = ABS( D( 1 ) )
  154:          RETURN
  155:       ELSE IF( N.EQ.2 ) THEN
  156:          CALL DLAS2( D( 1 ), E( 1 ), D( 2 ), SIGMN, SIGMX )
  157:          D( 1 ) = SIGMX
  158:          D( 2 ) = SIGMN
  159:          RETURN
  160:       END IF
  161: *
  162: *     Estimate the largest singular value.
  163: *
  164:       SIGMX = ZERO
  165:       DO 10 I = 1, N - 1
  166:          D( I ) = ABS( D( I ) )
  167:          SIGMX = MAX( SIGMX, ABS( E( I ) ) )
  168:    10 CONTINUE
  169:       D( N ) = ABS( D( N ) )
  170: *
  171: *     Early return if SIGMX is zero (matrix is already diagonal).
  172: *
  173:       IF( SIGMX.EQ.ZERO ) THEN
  174:          CALL DLASRT( 'D', N, D, IINFO )
  175:          RETURN
  176:       END IF
  177: *
  178:       DO 20 I = 1, N
  179:          SIGMX = MAX( SIGMX, D( I ) )
  180:    20 CONTINUE
  181: *
  182: *     Copy D and E into WORK (in the Z format) and scale (squaring the
  183: *     input data makes scaling by a power of the radix pointless).
  184: *
  185:       EPS = DLAMCH( 'Precision' )
  186:       SAFMIN = DLAMCH( 'Safe minimum' )
  187:       SCALE = SQRT( EPS / SAFMIN )
  188:       CALL DCOPY( N, D, 1, WORK( 1 ), 2 )
  189:       CALL DCOPY( N-1, E, 1, WORK( 2 ), 2 )
  190:       CALL DLASCL( 'G', 0, 0, SIGMX, SCALE, 2*N-1, 1, WORK, 2*N-1,
  191:      $             IINFO )
  192: *
  193: *     Compute the q's and e's.
  194: *
  195:       DO 30 I = 1, 2*N - 1
  196:          WORK( I ) = WORK( I )**2
  197:    30 CONTINUE
  198:       WORK( 2*N ) = ZERO
  199: *
  200:       CALL DLASQ2( N, WORK, INFO )
  201: *
  202:       IF( INFO.EQ.0 ) THEN
  203:          DO 40 I = 1, N
  204:             D( I ) = SQRT( WORK( I ) )
  205:    40    CONTINUE
  206:          CALL DLASCL( 'G', 0, 0, SCALE, SIGMX, N, 1, D, N, IINFO )
  207:       ELSE IF( INFO.EQ.2 ) THEN
  208: *
  209: *     Maximum number of iterations exceeded.  Move data from WORK
  210: *     into D and E so the calling subroutine can try to finish
  211: *
  212:          DO I = 1, N
  213:             D( I ) = SQRT( WORK( 2*I-1 ) )
  214:             E( I ) = SQRT( WORK( 2*I ) )
  215:          END DO
  216:          CALL DLASCL( 'G', 0, 0, SCALE, SIGMX, N, 1, D, N, IINFO )
  217:          CALL DLASCL( 'G', 0, 0, SCALE, SIGMX, N, 1, E, N, IINFO )
  218:       END IF
  219: *
  220:       RETURN
  221: *
  222: *     End of DLASQ1
  223: *
  224:       END