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

Diff for /rpl/lapack/lapack/dbdsdc.f between versions 1.2 and 1.17

version 1.2, 2010/04/21 13:45:12 version 1.17, 2017/06/17 10:53:46
Line 1 Line 1
   *> \brief \b DBDSDC
   *
   *  =========== DOCUMENTATION ===========
   *
   * Online html documentation available at
   *            http://www.netlib.org/lapack/explore-html/
   *
   *> \htmlonly
   *> Download DBDSDC + dependencies
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dbdsdc.f">
   *> [TGZ]</a>
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dbdsdc.f">
   *> [ZIP]</a>
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dbdsdc.f">
   *> [TXT]</a>
   *> \endhtmlonly
   *
   *  Definition:
   *  ===========
   *
   *       SUBROUTINE DBDSDC( UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ,
   *                          WORK, IWORK, INFO )
   *
   *       .. Scalar Arguments ..
   *       CHARACTER          COMPQ, UPLO
   *       INTEGER            INFO, LDU, LDVT, N
   *       ..
   *       .. Array Arguments ..
   *       INTEGER            IQ( * ), IWORK( * )
   *       DOUBLE PRECISION   D( * ), E( * ), Q( * ), U( LDU, * ),
   *      $                   VT( LDVT, * ), WORK( * )
   *       ..
   *
   *
   *> \par Purpose:
   *  =============
   *>
   *> \verbatim
   *>
   *> DBDSDC computes the singular value decomposition (SVD) of a real
   *> N-by-N (upper or lower) bidiagonal matrix B:  B = U * S * VT,
   *> using a divide and conquer method, 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. DBDSDC can be used to compute all singular values,
   *> and optionally, singular vectors or singular vectors in compact form.
   *>
   *> This code makes very mild assumptions about floating point
   *> arithmetic. It will work on machines with a guard digit in
   *> add/subtract, or on those binary machines without guard digits
   *> which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
   *> It could conceivably fail on hexadecimal or decimal machines
   *> without guard digits, but we know of none.  See DLASD3 for details.
   *>
   *> The code currently calls DLASDQ if singular values only are desired.
   *> However, it can be slightly modified to compute singular values
   *> using the divide and conquer method.
   *> \endverbatim
   *
   *  Arguments:
   *  ==========
   *
   *> \param[in] UPLO
   *> \verbatim
   *>          UPLO is CHARACTER*1
   *>          = 'U':  B is upper bidiagonal.
   *>          = 'L':  B is lower bidiagonal.
   *> \endverbatim
   *>
   *> \param[in] COMPQ
   *> \verbatim
   *>          COMPQ is CHARACTER*1
   *>          Specifies whether singular vectors are to be computed
   *>          as follows:
   *>          = 'N':  Compute singular values only;
   *>          = 'P':  Compute singular values and compute singular
   *>                  vectors in compact form;
   *>          = 'I':  Compute singular values and singular vectors.
   *> \endverbatim
   *>
   *> \param[in] N
   *> \verbatim
   *>          N is INTEGER
   *>          The order of the matrix B.  N >= 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.
   *> \endverbatim
   *>
   *> \param[in,out] E
   *> \verbatim
   *>          E is DOUBLE PRECISION array, dimension (N-1)
   *>          On entry, the elements of E contain the offdiagonal
   *>          elements of the bidiagonal matrix whose SVD is desired.
   *>          On exit, E has been destroyed.
   *> \endverbatim
   *>
   *> \param[out] U
   *> \verbatim
   *>          U is DOUBLE PRECISION array, dimension (LDU,N)
   *>          If  COMPQ = 'I', then:
   *>             On exit, if INFO = 0, U contains the left singular vectors
   *>             of the bidiagonal matrix.
   *>          For other values of COMPQ, U is not referenced.
   *> \endverbatim
   *>
   *> \param[in] LDU
   *> \verbatim
   *>          LDU is INTEGER
   *>          The leading dimension of the array U.  LDU >= 1.
   *>          If singular vectors are desired, then LDU >= max( 1, N ).
   *> \endverbatim
   *>
   *> \param[out] VT
   *> \verbatim
   *>          VT is DOUBLE PRECISION array, dimension (LDVT,N)
   *>          If  COMPQ = 'I', then:
   *>             On exit, if INFO = 0, VT**T contains the right singular
   *>             vectors of the bidiagonal matrix.
   *>          For other values of COMPQ, VT is not referenced.
   *> \endverbatim
   *>
   *> \param[in] LDVT
   *> \verbatim
   *>          LDVT is INTEGER
   *>          The leading dimension of the array VT.  LDVT >= 1.
   *>          If singular vectors are desired, then LDVT >= max( 1, N ).
   *> \endverbatim
   *>
   *> \param[out] Q
   *> \verbatim
   *>          Q is DOUBLE PRECISION array, dimension (LDQ)
   *>          If  COMPQ = 'P', then:
   *>             On exit, if INFO = 0, Q and IQ contain the left
   *>             and right singular vectors in a compact form,
   *>             requiring O(N log N) space instead of 2*N**2.
   *>             In particular, Q contains all the DOUBLE PRECISION data in
   *>             LDQ >= N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1))))
   *>             words of memory, where SMLSIZ is returned by ILAENV and
   *>             is equal to the maximum size of the subproblems at the
   *>             bottom of the computation tree (usually about 25).
   *>          For other values of COMPQ, Q is not referenced.
   *> \endverbatim
   *>
   *> \param[out] IQ
   *> \verbatim
   *>          IQ is INTEGER array, dimension (LDIQ)
   *>          If  COMPQ = 'P', then:
   *>             On exit, if INFO = 0, Q and IQ contain the left
   *>             and right singular vectors in a compact form,
   *>             requiring O(N log N) space instead of 2*N**2.
   *>             In particular, IQ contains all INTEGER data in
   *>             LDIQ >= N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1))))
   *>             words of memory, where SMLSIZ is returned by ILAENV and
   *>             is equal to the maximum size of the subproblems at the
   *>             bottom of the computation tree (usually about 25).
   *>          For other values of COMPQ, IQ is not referenced.
   *> \endverbatim
   *>
   *> \param[out] WORK
   *> \verbatim
   *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
   *>          If COMPQ = 'N' then LWORK >= (4 * N).
   *>          If COMPQ = 'P' then LWORK >= (6 * N).
   *>          If COMPQ = 'I' then LWORK >= (3 * N**2 + 4 * N).
   *> \endverbatim
   *>
   *> \param[out] IWORK
   *> \verbatim
   *>          IWORK is INTEGER array, dimension (8*N)
   *> \endverbatim
   *>
   *> \param[out] INFO
   *> \verbatim
   *>          INFO is INTEGER
   *>          = 0:  successful exit.
   *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
   *>          > 0:  The algorithm failed to compute a singular value.
   *>                The update process of divide and conquer failed.
   *> \endverbatim
   *
   *  Authors:
   *  ========
   *
   *> \author Univ. of Tennessee
   *> \author Univ. of California Berkeley
   *> \author Univ. of Colorado Denver
   *> \author NAG Ltd.
   *
   *> \date June 2016
   *
   *> \ingroup auxOTHERcomputational
   *
   *> \par Contributors:
   *  ==================
   *>
   *>     Ming Gu and Huan Ren, Computer Science Division, University of
   *>     California at Berkeley, USA
   *>
   *  =====================================================================
       SUBROUTINE DBDSDC( UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ,        SUBROUTINE DBDSDC( UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ,
      $                   WORK, IWORK, INFO )       $                   WORK, IWORK, INFO )
 *  *
 *  -- LAPACK routine (version 3.2) --  *  -- LAPACK computational routine (version 3.7.0) --
 *  -- 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..--
 *     November 2006  *     June 2016
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       CHARACTER          COMPQ, UPLO        CHARACTER          COMPQ, UPLO
Line 16 Line 220
      $                   VT( LDVT, * ), WORK( * )       $                   VT( LDVT, * ), WORK( * )
 *     ..  *     ..
 *  *
 *  Purpose  
 *  =======  
 *  
 *  DBDSDC computes the singular value decomposition (SVD) of a real  
 *  N-by-N (upper or lower) bidiagonal matrix B:  B = U * S * VT,  
 *  using a divide and conquer method, 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. DBDSDC can be used to compute all singular values,  
 *  and optionally, singular vectors or singular vectors in compact form.  
 *  
 *  This code makes very mild assumptions about floating point  
 *  arithmetic. It will work on machines with a guard digit in  
 *  add/subtract, or on those binary machines without guard digits  
 *  which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.  
 *  It could conceivably fail on hexadecimal or decimal machines  
 *  without guard digits, but we know of none.  See DLASD3 for details.  
 *  
 *  The code currently calls DLASDQ if singular values only are desired.  
 *  However, it can be slightly modified to compute singular values  
 *  using the divide and conquer method.  
 *  
 *  Arguments  
 *  =========  
 *  
 *  UPLO    (input) CHARACTER*1  
 *          = 'U':  B is upper bidiagonal.  
 *          = 'L':  B is lower bidiagonal.  
 *  
 *  COMPQ   (input) CHARACTER*1  
 *          Specifies whether singular vectors are to be computed  
 *          as follows:  
 *          = 'N':  Compute singular values only;  
 *          = 'P':  Compute singular values and compute singular  
 *                  vectors in compact form;  
 *          = 'I':  Compute singular values and singular vectors.  
 *  
 *  N       (input) INTEGER  
 *          The order of the matrix B.  N >= 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.  
 *  
 *  E       (input/output) DOUBLE PRECISION array, dimension (N-1)  
 *          On entry, the elements of E contain the offdiagonal  
 *          elements of the bidiagonal matrix whose SVD is desired.  
 *          On exit, E has been destroyed.  
 *  
 *  U       (output) DOUBLE PRECISION array, dimension (LDU,N)  
 *          If  COMPQ = 'I', then:  
 *             On exit, if INFO = 0, U contains the left singular vectors  
 *             of the bidiagonal matrix.  
 *          For other values of COMPQ, U is not referenced.  
 *  
 *  LDU     (input) INTEGER  
 *          The leading dimension of the array U.  LDU >= 1.  
 *          If singular vectors are desired, then LDU >= max( 1, N ).  
 *  
 *  VT      (output) DOUBLE PRECISION array, dimension (LDVT,N)  
 *          If  COMPQ = 'I', then:  
 *             On exit, if INFO = 0, VT' contains the right singular  
 *             vectors of the bidiagonal matrix.  
 *          For other values of COMPQ, VT is not referenced.  
 *  
 *  LDVT    (input) INTEGER  
 *          The leading dimension of the array VT.  LDVT >= 1.  
 *          If singular vectors are desired, then LDVT >= max( 1, N ).  
 *  
 *  Q       (output) DOUBLE PRECISION array, dimension (LDQ)  
 *          If  COMPQ = 'P', then:  
 *             On exit, if INFO = 0, Q and IQ contain the left  
 *             and right singular vectors in a compact form,  
 *             requiring O(N log N) space instead of 2*N**2.  
 *             In particular, Q contains all the DOUBLE PRECISION data in  
 *             LDQ >= N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1))))  
 *             words of memory, where SMLSIZ is returned by ILAENV and  
 *             is equal to the maximum size of the subproblems at the  
 *             bottom of the computation tree (usually about 25).  
 *          For other values of COMPQ, Q is not referenced.  
 *  
 *  IQ      (output) INTEGER array, dimension (LDIQ)  
 *          If  COMPQ = 'P', then:  
 *             On exit, if INFO = 0, Q and IQ contain the left  
 *             and right singular vectors in a compact form,  
 *             requiring O(N log N) space instead of 2*N**2.  
 *             In particular, IQ contains all INTEGER data in  
 *             LDIQ >= N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1))))  
 *             words of memory, where SMLSIZ is returned by ILAENV and  
 *             is equal to the maximum size of the subproblems at the  
 *             bottom of the computation tree (usually about 25).  
 *          For other values of COMPQ, IQ is not referenced.  
 *  
 *  WORK    (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK))  
 *          If COMPQ = 'N' then LWORK >= (4 * N).  
 *          If COMPQ = 'P' then LWORK >= (6 * N).  
 *          If COMPQ = 'I' then LWORK >= (3 * N**2 + 4 * N).  
 *  
 *  IWORK   (workspace) INTEGER array, dimension (8*N)  
 *  
 *  INFO    (output) INTEGER  
 *          = 0:  successful exit.  
 *          < 0:  if INFO = -i, the i-th argument had an illegal value.  
 *          > 0:  The algorithm failed to compute an singular value.  
 *                The update process of divide and conquer failed.  
 *  
 *  Further Details  
 *  ===============  
 *  
 *  Based on contributions by  
 *     Ming Gu and Huan Ren, Computer Science Division, University of  
 *     California at Berkeley, USA  
 *  
 *  =====================================================================  *  =====================================================================
 *  Changed dimension statement in comment describing E from (N) to  *  Changed dimension statement in comment describing E from (N) to
 *  (N-1).  Sven, 17 Feb 05.  *  (N-1).  Sven, 17 Feb 05.
Line 220 Line 311
       WSTART = 1        WSTART = 1
       QSTART = 3        QSTART = 3
       IF( ICOMPQ.EQ.1 ) THEN        IF( ICOMPQ.EQ.1 ) THEN
          CALL DCOPY( N, D, 1, Q( 1 ), 1 )           CALL DCOPY( N,   D, 1, Q( 1 ),   1 )
          CALL DCOPY( N-1, E, 1, Q( N+1 ), 1 )           CALL DCOPY( N-1, E, 1, Q( N+1 ), 1 )
       END IF        END IF
       IF( IUPLO.EQ.2 ) THEN        IF( IUPLO.EQ.2 ) THEN
Line 244 Line 335
 *     If ICOMPQ = 0, use DLASDQ to compute the singular values.  *     If ICOMPQ = 0, use DLASDQ to compute the singular values.
 *  *
       IF( ICOMPQ.EQ.0 ) THEN        IF( ICOMPQ.EQ.0 ) THEN
   *        Ignore WSTART, instead using WORK( 1 ), since the two vectors
   *        for CS and -SN above are added only if ICOMPQ == 2,
   *        and adding them exceeds documented WORK size of 4*n.
          CALL DLASDQ( 'U', 0, N, 0, 0, 0, D, E, VT, LDVT, U, LDU, U,           CALL DLASDQ( 'U', 0, N, 0, 0, 0, D, E, VT, LDVT, U, LDU, U,
      $                LDU, WORK( WSTART ), INFO )       $                LDU, WORK( 1 ), INFO )
          GO TO 40           GO TO 40
       END IF        END IF
 *  *
Line 287 Line 381
       CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, IERR )        CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, IERR )
       CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, IERR )        CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, IERR )
 *  *
       EPS = DLAMCH( 'Epsilon' )        EPS = (0.9D+0)*DLAMCH( 'Epsilon' )
 *  *
       MLVL = INT( LOG( DBLE( N ) / DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1        MLVL = INT( LOG( DBLE( N ) / DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1
       SMLSZP = SMLSIZ + 1        SMLSZP = SMLSIZ + 1
Line 321 Line 415
       DO 30 I = 1, NM1        DO 30 I = 1, NM1
          IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN           IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN
 *  *
 *        Subproblem found. First determine its size and then  *           Subproblem found. First determine its size and then
 *        apply divide and conquer on it.  *           apply divide and conquer on it.
 *  *
             IF( I.LT.NM1 ) THEN              IF( I.LT.NM1 ) THEN
 *  *
 *        A subproblem with E(I) small for I < NM1.  *              A subproblem with E(I) small for I < NM1.
 *  *
                NSIZE = I - START + 1                 NSIZE = I - START + 1
             ELSE IF( ABS( E( I ) ).GE.EPS ) THEN              ELSE IF( ABS( E( I ) ).GE.EPS ) THEN
 *  *
 *        A subproblem with E(NM1) not too small but I = NM1.  *              A subproblem with E(NM1) not too small but I = NM1.
 *  *
                NSIZE = N - START + 1                 NSIZE = N - START + 1
             ELSE              ELSE
 *  *
 *        A subproblem with E(NM1) small. This implies an  *              A subproblem with E(NM1) small. This implies an
 *        1-by-1 subproblem at D(N). Solve this 1-by-1 problem  *              1-by-1 subproblem at D(N). Solve this 1-by-1 problem
 *        first.  *              first.
 *  *
                NSIZE = I - START + 1                 NSIZE = I - START + 1
                IF( ICOMPQ.EQ.2 ) THEN                 IF( ICOMPQ.EQ.2 ) THEN
Line 368 Line 462
      $                      Q( START+( IC+QSTART-2 )*N ),       $                      Q( START+( IC+QSTART-2 )*N ),
      $                      Q( START+( IS+QSTART-2 )*N ),       $                      Q( START+( IS+QSTART-2 )*N ),
      $                      WORK( WSTART ), IWORK, INFO )       $                      WORK( WSTART ), IWORK, INFO )
                IF( INFO.NE.0 ) THEN              END IF
                   RETURN              IF( INFO.NE.0 ) THEN
                END IF                 RETURN
             END IF              END IF
             START = I + 1              START = I + 1
          END IF           END IF
Removed from v.1.2  
changed lines
  Added in v.1.17

CVSweb interface <joel.bertrand@systella.fr>