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

Diff for /rpl/lapack/lapack/dgesvj.f between versions 1.14 and 1.21

version 1.14, 2015/11/26 11:44:15 version 1.21, 2023/08/07 08:38:50
Line 2 Line 2
 *  *
 *  =========== DOCUMENTATION ===========  *  =========== DOCUMENTATION ===========
 *  *
 * Online html documentation available at   * Online html documentation available at
 *            http://www.netlib.org/lapack/explore-html/   *            http://www.netlib.org/lapack/explore-html/
 *  *
 *> \htmlonly  *> \htmlonly
 *> Download DGESVJ + dependencies   *> Download DGESVJ + dependencies
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesvj.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesvj.f">
 *> [TGZ]</a>   *> [TGZ]</a>
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesvj.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesvj.f">
 *> [ZIP]</a>   *> [ZIP]</a>
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesvj.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesvj.f">
 *> [TXT]</a>  *> [TXT]</a>
 *> \endhtmlonly   *> \endhtmlonly
 *  *
 *  Definition:  *  Definition:
 *  ===========  *  ===========
 *  *
 *       SUBROUTINE DGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,  *       SUBROUTINE DGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
 *                          LDV, WORK, LWORK, INFO )  *                          LDV, WORK, LWORK, INFO )
 *   *
 *       .. Scalar Arguments ..  *       .. Scalar Arguments ..
 *       INTEGER            INFO, LDA, LDV, LWORK, M, MV, N  *       INTEGER            INFO, LDA, LDV, LWORK, M, MV, N
 *       CHARACTER*1        JOBA, JOBU, JOBV  *       CHARACTER*1        JOBA, JOBU, JOBV
Line 29 Line 29
 *       DOUBLE PRECISION   A( LDA, * ), SVA( N ), V( LDV, * ),  *       DOUBLE PRECISION   A( LDA, * ), SVA( N ), V( LDV, * ),
 *      $                   WORK( LWORK )  *      $                   WORK( LWORK )
 *       ..  *       ..
 *    *
 *  *
 *> \par Purpose:  *> \par Purpose:
 *  =============  *  =============
Line 54 Line 54
 *  *
 *> \param[in] JOBA  *> \param[in] JOBA
 *> \verbatim  *> \verbatim
 *>          JOBA is CHARACTER* 1  *>          JOBA is CHARACTER*1
 *>          Specifies the structure of A.  *>          Specifies the structure of A.
 *>          = 'L': The input matrix A is lower triangular;  *>          = 'L': The input matrix A is lower triangular;
 *>          = 'U': The input matrix A is upper triangular;  *>          = 'U': The input matrix A is upper triangular;
Line 90 Line 90
 *>          JOBV is CHARACTER*1  *>          JOBV is CHARACTER*1
 *>          Specifies whether to compute the right singular vectors, that  *>          Specifies whether to compute the right singular vectors, that
 *>          is, the matrix V:  *>          is, the matrix V:
 *>          = 'V' : the matrix V is computed and returned in the array V  *>          = 'V':  the matrix V is computed and returned in the array V
 *>          = 'A' : the Jacobi rotations are applied to the MV-by-N  *>          = 'A':  the Jacobi rotations are applied to the MV-by-N
 *>                  array V. In other words, the right singular vector  *>                  array V. In other words, the right singular vector
 *>                  matrix V is not computed explicitly, instead it is  *>                  matrix V is not computed explicitly, instead it is
 *>                  applied to an MV-by-N matrix initially stored in the  *>                  applied to an MV-by-N matrix initially stored in the
 *>                  first MV rows of V.  *>                  first MV rows of V.
 *>          = 'N' : the matrix V is not computed and the array V is not  *>          = 'N':  the matrix V is not computed and the array V is not
 *>                  referenced  *>                  referenced
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \param[in] M  *> \param[in] M
 *> \verbatim  *> \verbatim
 *>          M is INTEGER  *>          M is INTEGER
 *>          The number of rows of the input matrix A. 1/DLAMCH('E') > M >= 0.    *>          The number of rows of the input matrix A. 1/DLAMCH('E') > M >= 0.
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \param[in] N  *> \param[in] N
Line 118 Line 118
 *>          A is DOUBLE PRECISION array, dimension (LDA,N)  *>          A is DOUBLE PRECISION array, dimension (LDA,N)
 *>          On entry, the M-by-N matrix A.  *>          On entry, the M-by-N matrix A.
 *>          On exit :  *>          On exit :
 *>          If JOBU .EQ. 'U' .OR. JOBU .EQ. 'C' :  *>          If JOBU = 'U' .OR. JOBU = 'C' :
 *>                 If INFO .EQ. 0 :  *>                 If INFO = 0 :
 *>                 RANKA orthonormal columns of U are returned in the  *>                 RANKA orthonormal columns of U are returned in the
 *>                 leading RANKA columns of the array A. Here RANKA <= N  *>                 leading RANKA columns of the array A. Here RANKA <= N
 *>                 is the number of computed singular values of A that are  *>                 is the number of computed singular values of A that are
Line 129 Line 129
 *>                 in the array WORK as RANKA=NINT(WORK(2)). Also see the  *>                 in the array WORK as RANKA=NINT(WORK(2)). Also see the
 *>                 descriptions of SVA and WORK. The computed columns of U  *>                 descriptions of SVA and WORK. The computed columns of U
 *>                 are mutually numerically orthogonal up to approximately  *>                 are mutually numerically orthogonal up to approximately
 *>                 TOL=DSQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU.EQ.'C'),  *>                 TOL=DSQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU = 'C'),
 *>                 see the description of JOBU.  *>                 see the description of JOBU.
 *>                 If INFO .GT. 0 :  *>                 If INFO > 0 :
 *>                 the procedure DGESVJ did not converge in the given number  *>                 the procedure DGESVJ did not converge in the given number
 *>                 of iterations (sweeps). In that case, the computed  *>                 of iterations (sweeps). In that case, the computed
 *>                 columns of U may not be orthogonal up to TOL. The output  *>                 columns of U may not be orthogonal up to TOL. The output
Line 140 Line 140
 *>                 input matrix A in the sense that the residual  *>                 input matrix A in the sense that the residual
 *>                 ||A-SCALE*U*SIGMA*V^T||_2 / ||A||_2 is small.  *>                 ||A-SCALE*U*SIGMA*V^T||_2 / ||A||_2 is small.
 *>  *>
 *>          If JOBU .EQ. 'N' :  *>          If JOBU = 'N' :
 *>                 If INFO .EQ. 0 :  *>                 If INFO = 0 :
 *>                 Note that the left singular vectors are 'for free' in the  *>                 Note that the left singular vectors are 'for free' in the
 *>                 one-sided Jacobi SVD algorithm. However, if only the  *>                 one-sided Jacobi SVD algorithm. However, if only the
 *>                 singular values are needed, the level of numerical  *>                 singular values are needed, the level of numerical
Line 150 Line 150
 *>                 numerically orthogonal up to approximately M*EPS. Thus,  *>                 numerically orthogonal up to approximately M*EPS. Thus,
 *>                 on exit, A contains the columns of U scaled with the  *>                 on exit, A contains the columns of U scaled with the
 *>                 corresponding singular values.  *>                 corresponding singular values.
 *>                 If INFO .GT. 0 :  *>                 If INFO > 0 :
 *>                 the procedure DGESVJ did not converge in the given number  *>                 the procedure DGESVJ did not converge in the given number
 *>                 of iterations (sweeps).  *>                 of iterations (sweeps).
 *> \endverbatim  *> \endverbatim
Line 165 Line 165
 *> \verbatim  *> \verbatim
 *>          SVA is DOUBLE PRECISION array, dimension (N)  *>          SVA is DOUBLE PRECISION array, dimension (N)
 *>          On exit :  *>          On exit :
 *>          If INFO .EQ. 0 :  *>          If INFO = 0 :
 *>          depending on the value SCALE = WORK(1), we have:  *>          depending on the value SCALE = WORK(1), we have:
 *>                 If SCALE .EQ. ONE :  *>                 If SCALE = ONE :
 *>                 SVA(1:N) contains the computed singular values of A.  *>                 SVA(1:N) contains the computed singular values of A.
 *>                 During the computation SVA contains the Euclidean column  *>                 During the computation SVA contains the Euclidean column
 *>                 norms of the iterated matrices in the array A.  *>                 norms of the iterated matrices in the array A.
Line 175 Line 175
 *>                 The singular values of A are SCALE*SVA(1:N), and this  *>                 The singular values of A are SCALE*SVA(1:N), and this
 *>                 factored representation is due to the fact that some of the  *>                 factored representation is due to the fact that some of the
 *>                 singular values of A might underflow or overflow.  *>                 singular values of A might underflow or overflow.
 *>          If INFO .GT. 0 :  *>          If INFO > 0 :
 *>          the procedure DGESVJ did not converge in the given number of  *>          the procedure DGESVJ did not converge in the given number of
 *>          iterations (sweeps) and SCALE*SVA(1:N) may not be accurate.  *>          iterations (sweeps) and SCALE*SVA(1:N) may not be accurate.
 *> \endverbatim  *> \endverbatim
Line 183 Line 183
 *> \param[in] MV  *> \param[in] MV
 *> \verbatim  *> \verbatim
 *>          MV is INTEGER  *>          MV is INTEGER
 *>          If JOBV .EQ. 'A', then the product of Jacobi rotations in DGESVJ  *>          If JOBV = 'A', then the product of Jacobi rotations in DGESVJ
 *>          is applied to the first MV rows of V. See the description of JOBV.  *>          is applied to the first MV rows of V. See the description of JOBV.
 *> \endverbatim  *> \endverbatim
 *>  *>
Line 201 Line 201
 *> \param[in] LDV  *> \param[in] LDV
 *> \verbatim  *> \verbatim
 *>          LDV is INTEGER  *>          LDV is INTEGER
 *>          The leading dimension of the array V, LDV .GE. 1.  *>          The leading dimension of the array V, LDV >= 1.
 *>          If JOBV .EQ. 'V', then LDV .GE. max(1,N).  *>          If JOBV = 'V', then LDV >= max(1,N).
 *>          If JOBV .EQ. 'A', then LDV .GE. max(1,MV) .  *>          If JOBV = 'A', then LDV >= max(1,MV) .
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \param[in,out] WORK  *> \param[in,out] WORK
 *> \verbatim  *> \verbatim
 *>          WORK is DOUBLE PRECISION array, dimension max(4,M+N).  *>          WORK is DOUBLE PRECISION array, dimension (LWORK)
 *>          On entry :  *>          On entry :
 *>          If JOBU .EQ. 'C' :  *>          If JOBU = 'C' :
 *>          WORK(1) = CTOL, where CTOL defines the threshold for convergence.  *>          WORK(1) = CTOL, where CTOL defines the threshold for convergence.
 *>                    The process stops if all columns of A are mutually  *>                    The process stops if all columns of A are mutually
 *>                    orthogonal up to CTOL*EPS, EPS=DLAMCH('E').  *>                    orthogonal up to CTOL*EPS, EPS=DLAMCH('E').
Line 230 Line 230
 *>          WORK(5) = max_{i.NE.j} |COS(A(:,i),A(:,j))| in the last sweep.  *>          WORK(5) = max_{i.NE.j} |COS(A(:,i),A(:,j))| in the last sweep.
 *>                    This is useful information in cases when DGESVJ did  *>                    This is useful information in cases when DGESVJ did
 *>                    not converge, as it can be used to estimate whether  *>                    not converge, as it can be used to estimate whether
 *>                    the output is stil useful and for post festum analysis.  *>                    the output is still useful and for post festum analysis.
 *>          WORK(6) = the largest absolute value over all sines of the  *>          WORK(6) = the largest absolute value over all sines of the
 *>                    Jacobi rotation angles in the last sweep. It can be  *>                    Jacobi rotation angles in the last sweep. It can be
 *>                    useful for a post festum analysis.  *>                    useful for a post festum analysis.
Line 245 Line 245
 *> \param[out] INFO  *> \param[out] INFO
 *> \verbatim  *> \verbatim
 *>          INFO is INTEGER  *>          INFO is INTEGER
 *>          = 0 : successful exit.  *>          = 0:  successful exit.
 *>          < 0 : if INFO = -i, then the i-th argument had an illegal value  *>          < 0:  if INFO = -i, then the i-th argument had an illegal value
 *>          > 0 : DGESVJ did not converge in the maximal allowed number (30)  *>          > 0:  DGESVJ did not converge in the maximal allowed number (30)
 *>                of sweeps. The output may still be useful. See the  *>                of sweeps. The output may still be useful. See the
 *>                description of WORK.  *>                description of WORK.
 *> \endverbatim  *> \endverbatim
Line 255 Line 255
 *  Authors:  *  Authors:
 *  ========  *  ========
 *  *
 *> \author Univ. of Tennessee   *> \author Univ. of Tennessee
 *> \author Univ. of California Berkeley   *> \author Univ. of California Berkeley
 *> \author Univ. of Colorado Denver   *> \author Univ. of Colorado Denver
 *> \author NAG Ltd.   *> \author NAG Ltd.
 *  
 *> \date November 2015  
 *  *
 *> \ingroup doubleGEcomputational  *> \ingroup doubleGEcomputational
 *  *
Line 281 Line 279
 *>  spectral condition number. The best performance of this Jacobi SVD  *>  spectral condition number. The best performance of this Jacobi SVD
 *>  procedure is achieved if used in an  accelerated version of Drmac and  *>  procedure is achieved if used in an  accelerated version of Drmac and
 *>  Veselic [5,6], and it is the kernel routine in the SIGMA library [7].  *>  Veselic [5,6], and it is the kernel routine in the SIGMA library [7].
 *>  Some tunning parameters (marked with [TP]) are available for the  *>  Some tuning parameters (marked with [TP]) are available for the
 *>  implementer.  *>  implementer.
 *>  The computational range for the nonzero singular values is the  machine  *>  The computational range for the nonzero singular values is the  machine
 *>  number interval ( UNDERFLOW , OVERFLOW ). In extreme cases, even  *>  number interval ( UNDERFLOW , OVERFLOW ). In extreme cases, even
Line 337 Line 335
       SUBROUTINE DGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,        SUBROUTINE DGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
      $                   LDV, WORK, LWORK, INFO )       $                   LDV, WORK, LWORK, INFO )
 *  *
 *  -- LAPACK computational routine (version 3.6.0) --  *  -- 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..--
 *     November 2015  
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       INTEGER            INFO, LDA, LDV, LWORK, M, MV, N        INTEGER            INFO, LDA, LDV, LWORK, M, MV, N
Line 1261 Line 1258
                                     MXSINJ = MAX( MXSINJ, DABS( SN ) )                                      MXSINJ = MAX( MXSINJ, DABS( SN ) )
                                     SVA( q ) = AAQQ*DSQRT( MAX( ZERO,                                      SVA( q ) = AAQQ*DSQRT( MAX( ZERO,
      $                                         ONE+T*APOAQ*AAPQ ) )       $                                         ONE+T*APOAQ*AAPQ ) )
                                     AAPP = AAPP*DSQRT( MAX( ZERO,                                       AAPP = AAPP*DSQRT( MAX( ZERO,
      $                                     ONE-T*AQOAP*AAPQ ) )       $                                     ONE-T*AQOAP*AAPQ ) )
 *  *
                                     APOAQ = WORK( p ) / WORK( q )                                      APOAQ = WORK( p ) / WORK( q )
Removed from v.1.14  
changed lines
  Added in v.1.21

CVSweb interface <joel.bertrand@systella.fr>