[BACK]Return to dbdsqr.f CVS log [TXT] [DIR] Up to [local] / rpl / lapack / lapack

Diff for /rpl/lapack/lapack/dbdsqr.f between versions 1.4 and 1.19

version 1.4, 2010/08/06 15:32:22 version 1.19, 2023/08/07 08:38:47
Line 1 Line 1
   *> \brief \b DBDSQR
   *
   *  =========== DOCUMENTATION ===========
   *
   * Online html documentation available at
   *            http://www.netlib.org/lapack/explore-html/
   *
   *> \htmlonly
   *> Download DBDSQR + dependencies
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dbdsqr.f">
   *> [TGZ]</a>
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dbdsqr.f">
   *> [ZIP]</a>
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dbdsqr.f">
   *> [TXT]</a>
   *> \endhtmlonly
   *
   *  Definition:
   *  ===========
   *
   *       SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
   *                          LDU, C, LDC, WORK, INFO )
   *
   *       .. Scalar Arguments ..
   *       CHARACTER          UPLO
   *       INTEGER            INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU
   *       ..
   *       .. Array Arguments ..
   *       DOUBLE PRECISION   C( LDC, * ), D( * ), E( * ), U( LDU, * ),
   *      $                   VT( LDVT, * ), WORK( * )
   *       ..
   *
   *
   *> \par Purpose:
   *  =============
   *>
   *> \verbatim
   *>
   *> DBDSQR computes the singular values and, optionally, the right and/or
   *> left singular vectors from the singular value decomposition (SVD) of
   *> a real N-by-N (upper or lower) bidiagonal matrix B using the implicit
   *> zero-shift QR algorithm.  The SVD of B has the form
   *>
   *>    B = Q * S * P**T
   *>
   *> where S is the diagonal matrix of singular values, Q is an orthogonal
   *> matrix of left singular vectors, and P is an orthogonal matrix of
   *> right singular vectors.  If left singular vectors are requested, this
   *> subroutine actually returns U*Q instead of Q, and, if right singular
   *> vectors are requested, this subroutine returns P**T*VT instead of
   *> P**T, for given real input matrices U and VT.  When U and VT are the
   *> orthogonal matrices that reduce a general matrix A to bidiagonal
   *> form:  A = U*B*VT, as computed by DGEBRD, then
   *>
   *>    A = (U*Q) * S * (P**T*VT)
   *>
   *> is the SVD of A.  Optionally, the subroutine may also compute Q**T*C
   *> for a given real input matrix C.
   *>
   *> See "Computing  Small Singular Values of Bidiagonal Matrices With
   *> Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,
   *> LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,
   *> no. 5, pp. 873-912, Sept 1990) and
   *> "Accurate singular values and differential qd algorithms," by
   *> B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics
   *> Department, University of California at Berkeley, July 1992
   *> for a detailed description of the algorithm.
   *> \endverbatim
   *
   *  Arguments:
   *  ==========
   *
   *> \param[in] UPLO
   *> \verbatim
   *>          UPLO is CHARACTER*1
   *>          = 'U':  B is upper bidiagonal;
   *>          = 'L':  B is lower bidiagonal.
   *> \endverbatim
   *>
   *> \param[in] N
   *> \verbatim
   *>          N is INTEGER
   *>          The order of the matrix B.  N >= 0.
   *> \endverbatim
   *>
   *> \param[in] NCVT
   *> \verbatim
   *>          NCVT is INTEGER
   *>          The number of columns of the matrix VT. NCVT >= 0.
   *> \endverbatim
   *>
   *> \param[in] NRU
   *> \verbatim
   *>          NRU is INTEGER
   *>          The number of rows of the matrix U. NRU >= 0.
   *> \endverbatim
   *>
   *> \param[in] NCC
   *> \verbatim
   *>          NCC is INTEGER
   *>          The number of columns of the matrix C. NCC >= 0.
   *> \endverbatim
   *>
   *> \param[in,out] D
   *> \verbatim
   *>          D is DOUBLE PRECISION array, dimension (N)
   *>          On entry, the n diagonal elements of the bidiagonal matrix B.
   *>          On exit, if INFO=0, the singular values of B in decreasing
   *>          order.
   *> \endverbatim
   *>
   *> \param[in,out] E
   *> \verbatim
   *>          E is DOUBLE PRECISION array, dimension (N-1)
   *>          On entry, the N-1 offdiagonal elements of the bidiagonal
   *>          matrix B.
   *>          On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E
   *>          will contain the diagonal and superdiagonal elements of a
   *>          bidiagonal matrix orthogonally equivalent to the one given
   *>          as input.
   *> \endverbatim
   *>
   *> \param[in,out] VT
   *> \verbatim
   *>          VT is DOUBLE PRECISION array, dimension (LDVT, NCVT)
   *>          On entry, an N-by-NCVT matrix VT.
   *>          On exit, VT is overwritten by P**T * VT.
   *>          Not referenced if NCVT = 0.
   *> \endverbatim
   *>
   *> \param[in] LDVT
   *> \verbatim
   *>          LDVT is INTEGER
   *>          The leading dimension of the array VT.
   *>          LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.
   *> \endverbatim
   *>
   *> \param[in,out] U
   *> \verbatim
   *>          U is DOUBLE PRECISION array, dimension (LDU, N)
   *>          On entry, an NRU-by-N matrix U.
   *>          On exit, U is overwritten by U * Q.
   *>          Not referenced if NRU = 0.
   *> \endverbatim
   *>
   *> \param[in] LDU
   *> \verbatim
   *>          LDU is INTEGER
   *>          The leading dimension of the array U.  LDU >= max(1,NRU).
   *> \endverbatim
   *>
   *> \param[in,out] C
   *> \verbatim
   *>          C is DOUBLE PRECISION array, dimension (LDC, NCC)
   *>          On entry, an N-by-NCC matrix C.
   *>          On exit, C is overwritten by Q**T * C.
   *>          Not referenced if NCC = 0.
   *> \endverbatim
   *>
   *> \param[in] LDC
   *> \verbatim
   *>          LDC is INTEGER
   *>          The leading dimension of the array C.
   *>          LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.
   *> \endverbatim
   *>
   *> \param[out] WORK
   *> \verbatim
   *>          WORK is DOUBLE PRECISION array, dimension (4*(N-1))
   *> \endverbatim
   *>
   *> \param[out] INFO
   *> \verbatim
   *>          INFO is INTEGER
   *>          = 0:  successful exit
   *>          < 0:  If INFO = -i, the i-th argument had an illegal value
   *>          > 0:
   *>             if NCVT = NRU = NCC = 0,
   *>                = 1, a split was marked by a positive value in E
   *>                = 2, current block of Z not diagonalized after 30*N
   *>                     iterations (in inner while loop)
   *>                = 3, termination criterion of outer while loop not met
   *>                     (program created more than N unreduced blocks)
   *>             else NCVT = NRU = NCC = 0,
   *>                   the algorithm did not converge; D and E contain the
   *>                   elements of a bidiagonal matrix which is orthogonally
   *>                   similar to the input matrix B;  if INFO = i, i
   *>                   elements of E have not converged to zero.
   *> \endverbatim
   *
   *> \par Internal Parameters:
   *  =========================
   *>
   *> \verbatim
   *>  TOLMUL  DOUBLE PRECISION, default = max(10,min(100,EPS**(-1/8)))
   *>          TOLMUL controls the convergence criterion of the QR loop.
   *>          If it is positive, TOLMUL*EPS is the desired relative
   *>             precision in the computed singular values.
   *>          If it is negative, abs(TOLMUL*EPS*sigma_max) is the
   *>             desired absolute accuracy in the computed singular
   *>             values (corresponds to relative accuracy
   *>             abs(TOLMUL*EPS) in the largest singular value.
   *>          abs(TOLMUL) should be between 1 and 1/EPS, and preferably
   *>             between 10 (for fast convergence) and .1/EPS
   *>             (for there to be some accuracy in the results).
   *>          Default is to lose at either one eighth or 2 of the
   *>             available decimal digits in each computed singular value
   *>             (whichever is smaller).
   *>
   *>  MAXITR  INTEGER, default = 6
   *>          MAXITR controls the maximum number of passes of the
   *>          algorithm through its inner loop. The algorithms stops
   *>          (and so fails to converge) if the number of passes
   *>          through the inner loop exceeds MAXITR*N**2.
   *>
   *> \endverbatim
   *
   *> \par Note:
   *  ===========
   *>
   *> \verbatim
   *>  Bug report from Cezary Dendek.
   *>  On March 23rd 2017, the INTEGER variable MAXIT = MAXITR*N**2 is
   *>  removed since it can overflow pretty easily (for N larger or equal
   *>  than 18,919). We instead use MAXITDIVN = MAXITR*N.
   *> \endverbatim
   *
   *  Authors:
   *  ========
   *
   *> \author Univ. of Tennessee
   *> \author Univ. of California Berkeley
   *> \author Univ. of Colorado Denver
   *> \author NAG Ltd.
   *
   *> \ingroup auxOTHERcomputational
   *
   *  =====================================================================
       SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,        SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U,
      $                   LDU, C, LDC, WORK, INFO )       $                   LDU, C, LDC, WORK, INFO )
 *  *
 *  -- LAPACK routine (version 3.2) --  *  -- LAPACK computational routine --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --  *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--  *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 *     January 2007  
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       CHARACTER          UPLO        CHARACTER          UPLO
Line 15 Line 252
      $                   VT( LDVT, * ), WORK( * )       $                   VT( LDVT, * ), WORK( * )
 *     ..  *     ..
 *  *
 *  Purpose  
 *  =======  
 *  
 *  DBDSQR computes the singular values and, optionally, the right and/or  
 *  left singular vectors from the singular value decomposition (SVD) of  
 *  a real N-by-N (upper or lower) bidiagonal matrix B using the implicit  
 *  zero-shift QR algorithm.  The SVD of B has the form  
 *   
 *     B = Q * S * P**T  
 *   
 *  where S is the diagonal matrix of singular values, Q is an orthogonal  
 *  matrix of left singular vectors, and P is an orthogonal matrix of  
 *  right singular vectors.  If left singular vectors are requested, this  
 *  subroutine actually returns U*Q instead of Q, and, if right singular  
 *  vectors are requested, this subroutine returns P**T*VT instead of  
 *  P**T, for given real input matrices U and VT.  When U and VT are the  
 *  orthogonal matrices that reduce a general matrix A to bidiagonal  
 *  form:  A = U*B*VT, as computed by DGEBRD, then  
 *  
 *     A = (U*Q) * S * (P**T*VT)  
 *  
 *  is the SVD of A.  Optionally, the subroutine may also compute Q**T*C  
 *  for a given real input matrix C.  
 *  
 *  See "Computing  Small Singular Values of Bidiagonal Matrices With  
 *  Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan,  
 *  LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11,  
 *  no. 5, pp. 873-912, Sept 1990) and  
 *  "Accurate singular values and differential qd algorithms," by  
 *  B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics  
 *  Department, University of California at Berkeley, July 1992  
 *  for a detailed description of the algorithm.  
 *  
 *  Arguments  
 *  =========  
 *  
 *  UPLO    (input) CHARACTER*1  
 *          = 'U':  B is upper bidiagonal;  
 *          = 'L':  B is lower bidiagonal.  
 *  
 *  N       (input) INTEGER  
 *          The order of the matrix B.  N >= 0.  
 *  
 *  NCVT    (input) INTEGER  
 *          The number of columns of the matrix VT. NCVT >= 0.  
 *  
 *  NRU     (input) INTEGER  
 *          The number of rows of the matrix U. NRU >= 0.  
 *  
 *  NCC     (input) INTEGER  
 *          The number of columns of the matrix C. NCC >= 0.  
 *  
 *  D       (input/output) DOUBLE PRECISION array, dimension (N)  
 *          On entry, the n diagonal elements of the bidiagonal matrix B.  
 *          On exit, if INFO=0, the singular values of B in decreasing  
 *          order.  
 *  
 *  E       (input/output) DOUBLE PRECISION array, dimension (N-1)  
 *          On entry, the N-1 offdiagonal elements of the bidiagonal  
 *          matrix B.   
 *          On exit, if INFO = 0, E is destroyed; if INFO > 0, D and E  
 *          will contain the diagonal and superdiagonal elements of a  
 *          bidiagonal matrix orthogonally equivalent to the one given  
 *          as input.  
 *  
 *  VT      (input/output) DOUBLE PRECISION array, dimension (LDVT, NCVT)  
 *          On entry, an N-by-NCVT matrix VT.  
 *          On exit, VT is overwritten by P**T * VT.  
 *          Not referenced if NCVT = 0.  
 *  
 *  LDVT    (input) INTEGER  
 *          The leading dimension of the array VT.  
 *          LDVT >= max(1,N) if NCVT > 0; LDVT >= 1 if NCVT = 0.  
 *  
 *  U       (input/output) DOUBLE PRECISION array, dimension (LDU, N)  
 *          On entry, an NRU-by-N matrix U.  
 *          On exit, U is overwritten by U * Q.  
 *          Not referenced if NRU = 0.  
 *  
 *  LDU     (input) INTEGER  
 *          The leading dimension of the array U.  LDU >= max(1,NRU).  
 *  
 *  C       (input/output) DOUBLE PRECISION array, dimension (LDC, NCC)  
 *          On entry, an N-by-NCC matrix C.  
 *          On exit, C is overwritten by Q**T * C.  
 *          Not referenced if NCC = 0.  
 *  
 *  LDC     (input) INTEGER  
 *          The leading dimension of the array C.  
 *          LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0.  
 *  
 *  WORK    (workspace) DOUBLE PRECISION array, dimension (4*N)  
 *  
 *  INFO    (output) INTEGER  
 *          = 0:  successful exit  
 *          < 0:  If INFO = -i, the i-th argument had an illegal value  
 *          > 0:  
 *             if NCVT = NRU = NCC = 0,  
 *                = 1, a split was marked by a positive value in E  
 *                = 2, current block of Z not diagonalized after 30*N  
 *                     iterations (in inner while loop)  
 *                = 3, termination criterion of outer while loop not met   
 *                     (program created more than N unreduced blocks)  
 *             else NCVT = NRU = NCC = 0,  
 *                   the algorithm did not converge; D and E contain the  
 *                   elements of a bidiagonal matrix which is orthogonally  
 *                   similar to the input matrix B;  if INFO = i, i  
 *                   elements of E have not converged to zero.  
 *  
 *  Internal Parameters  
 *  ===================  
 *  
 *  TOLMUL  DOUBLE PRECISION, default = max(10,min(100,EPS**(-1/8)))  
 *          TOLMUL controls the convergence criterion of the QR loop.  
 *          If it is positive, TOLMUL*EPS is the desired relative  
 *             precision in the computed singular values.  
 *          If it is negative, abs(TOLMUL*EPS*sigma_max) is the  
 *             desired absolute accuracy in the computed singular  
 *             values (corresponds to relative accuracy  
 *             abs(TOLMUL*EPS) in the largest singular value.  
 *          abs(TOLMUL) should be between 1 and 1/EPS, and preferably  
 *             between 10 (for fast convergence) and .1/EPS  
 *             (for there to be some accuracy in the results).  
 *          Default is to lose at either one eighth or 2 of the  
 *             available decimal digits in each computed singular value  
 *             (whichever is smaller).  
 *  
 *  MAXITR  INTEGER, default = 6  
 *          MAXITR controls the maximum number of passes of the  
 *          algorithm through its inner loop. The algorithms stops  
 *          (and so fails to converge) if the number of passes  
 *          through the inner loop exceeds MAXITR*N**2.  
 *  
 *  =====================================================================  *  =====================================================================
 *  *
 *     .. Parameters ..  *     .. Parameters ..
Line 170 Line 274
 *     ..  *     ..
 *     .. Local Scalars ..  *     .. Local Scalars ..
       LOGICAL            LOWER, ROTATE        LOGICAL            LOWER, ROTATE
       INTEGER            I, IDIR, ISUB, ITER, J, LL, LLL, M, MAXIT, NM1,        INTEGER            I, IDIR, ISUB, ITER, ITERDIVN, J, LL, LLL, M,
      $                   NM12, NM13, OLDLL, OLDM       $                   MAXITDIVN, NM1, NM12, NM13, OLDLL, OLDM
       DOUBLE PRECISION   ABSE, ABSS, COSL, COSR, CS, EPS, F, G, H, MU,        DOUBLE PRECISION   ABSE, ABSS, COSL, COSR, CS, EPS, F, G, H, MU,
      $                   OLDCS, OLDSN, R, SHIFT, SIGMN, SIGMX, SINL,       $                   OLDCS, OLDSN, R, SHIFT, SIGMN, SIGMX, SINL,
      $                   SINR, SLL, SMAX, SMIN, SMINL, SMINOA,       $                   SINR, SLL, SMAX, SMIN, SMINL, SMINOA,
Line 231 Line 335
 *  *
       IF( .NOT.ROTATE ) THEN        IF( .NOT.ROTATE ) THEN
          CALL DLASQ1( N, D, E, WORK, INFO )           CALL DLASQ1( N, D, E, WORK, INFO )
          RETURN  *
   *     If INFO equals 2, dqds didn't finish, try to finish
   *
            IF( INFO .NE. 2 ) RETURN
            INFO = 0
       END IF        END IF
 *  *
       NM1 = N - 1        NM1 = N - 1
Line 300 Line 408
    40    CONTINUE     40    CONTINUE
    50    CONTINUE     50    CONTINUE
          SMINOA = SMINOA / SQRT( DBLE( N ) )           SMINOA = SMINOA / SQRT( DBLE( N ) )
          THRESH = MAX( TOL*SMINOA, MAXITR*N*N*UNFL )           THRESH = MAX( TOL*SMINOA, MAXITR*(N*(N*UNFL)) )
       ELSE        ELSE
 *  *
 *        Absolute accuracy desired  *        Absolute accuracy desired
 *  *
          THRESH = MAX( ABS( TOL )*SMAX, MAXITR*N*N*UNFL )           THRESH = MAX( ABS( TOL )*SMAX, MAXITR*(N*(N*UNFL)) )
       END IF        END IF
 *  *
 *     Prepare for main iteration loop for the singular values  *     Prepare for main iteration loop for the singular values
 *     (MAXIT is the maximum number of passes through the inner  *     (MAXIT is the maximum number of passes through the inner
 *     loop permitted before nonconvergence signalled.)  *     loop permitted before nonconvergence signalled.)
 *  *
       MAXIT = MAXITR*N*N        MAXITDIVN = MAXITR*N
       ITER = 0        ITERDIVN = 0
         ITER = -1
       OLDLL = -1        OLDLL = -1
       OLDM = -1        OLDM = -1
 *  *
Line 329 Line 438
 *  *
       IF( M.LE.1 )        IF( M.LE.1 )
      $   GO TO 160       $   GO TO 160
       IF( ITER.GT.MAXIT )  *
      $   GO TO 200        IF( ITER.GE.N ) THEN
            ITER = ITER - N
            ITERDIVN = ITERDIVN + 1
            IF( ITERDIVN.GE.MAXITDIVN )
        $      GO TO 200
         END IF
 *  *
 *     Find diagonal block of matrix to work on  *     Find diagonal block of matrix to work on
 *  *
Removed from v.1.4  
changed lines
  Added in v.1.19

CVSweb interface <joel.bertrand@systella.fr>