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

Diff for /rpl/lapack/lapack/dgejsv.f between versions 1.19 and 1.20

version 1.19, 2018/05/29 07:17:51 version 1.20, 2020/05/21 21:45:56
Line 82 Line 82
 *>             desirable, then this option is advisable. The input matrix A  *>             desirable, then this option is advisable. The input matrix A
 *>             is preprocessed with QR factorization with FULL (row and  *>             is preprocessed with QR factorization with FULL (row and
 *>             column) pivoting.  *>             column) pivoting.
 *>       = 'G'  Computation as with 'F' with an additional estimate of the  *>       = 'G': Computation as with 'F' with an additional estimate of the
 *>             condition number of B, where A=D*B. If A has heavily weighted  *>             condition number of B, where A=D*B. If A has heavily weighted
 *>             rows, then using this condition number gives too pessimistic  *>             rows, then using this condition number gives too pessimistic
 *>             error bound.  *>             error bound.
Line 133 Line 133
 *>        specified range. If A .NE. 0 is scaled so that the largest singular  *>        specified range. If A .NE. 0 is scaled so that the largest singular
 *>        value of c*A is around DSQRT(BIG), BIG=SLAMCH('O'), then JOBR issues  *>        value of c*A is around DSQRT(BIG), BIG=SLAMCH('O'), then JOBR issues
 *>        the licence to kill columns of A whose norm in c*A is less than  *>        the licence to kill columns of A whose norm in c*A is less than
 *>        DSQRT(SFMIN) (for JOBR.EQ.'R'), or less than SMALL=SFMIN/EPSLN,  *>        DSQRT(SFMIN) (for JOBR = 'R'), or less than SMALL=SFMIN/EPSLN,
 *>        where SFMIN=SLAMCH('S'), EPSLN=SLAMCH('E').  *>        where SFMIN=SLAMCH('S'), EPSLN=SLAMCH('E').
 *>       = 'N': Do not kill small columns of c*A. This option assumes that  *>       = 'N': Do not kill small columns of c*A. This option assumes that
 *>             BLAS and QR factorizations and triangular solvers are  *>             BLAS and QR factorizations and triangular solvers are
Line 230 Line 230
 *>          If JOBU = 'F', then U contains on exit the M-by-M matrix of  *>          If JOBU = 'F', then U contains on exit the M-by-M matrix of
 *>                         the left singular vectors, including an ONB  *>                         the left singular vectors, including an ONB
 *>                         of the orthogonal complement of the Range(A).  *>                         of the orthogonal complement of the Range(A).
 *>          If JOBU = 'W'  .AND. (JOBV.EQ.'V' .AND. JOBT.EQ.'T' .AND. M.EQ.N),  *>          If JOBU = 'W'  .AND. (JOBV = 'V' .AND. JOBT = 'T' .AND. M = N),
 *>                         then U is used as workspace if the procedure  *>                         then U is used as workspace if the procedure
 *>                         replaces A with A^t. In that case, [V] is computed  *>                         replaces A with A^t. In that case, [V] is computed
 *>                         in U as left singular vectors of A^t and then  *>                         in U as left singular vectors of A^t and then
Line 252 Line 252
 *>          V is DOUBLE PRECISION array, dimension ( LDV, N )  *>          V is DOUBLE PRECISION array, dimension ( LDV, N )
 *>          If JOBV = 'V', 'J' then V contains on exit the N-by-N matrix of  *>          If JOBV = 'V', 'J' then V contains on exit the N-by-N matrix of
 *>                         the right singular vectors;  *>                         the right singular vectors;
 *>          If JOBV = 'W', AND (JOBU.EQ.'U' AND JOBT.EQ.'T' AND M.EQ.N),  *>          If JOBV = 'W', AND (JOBU = 'U' AND JOBT = 'T' AND M = N),
 *>                         then V is used as workspace if the pprocedure  *>                         then V is used as workspace if the pprocedure
 *>                         replaces A with A^t. In that case, [U] is computed  *>                         replaces A with A^t. In that case, [U] is computed
 *>                         in V as right singular vectors of A^t and then  *>                         in V as right singular vectors of A^t and then
Line 272 Line 272
 *> \param[out] WORK  *> \param[out] WORK
 *> \verbatim  *> \verbatim
 *>          WORK is DOUBLE PRECISION array, dimension (LWORK)  *>          WORK is DOUBLE PRECISION array, dimension (LWORK)
 *>          On exit, if N.GT.0 .AND. M.GT.0 (else not referenced),  *>          On exit, if N > 0 .AND. M > 0 (else not referenced),
 *>          WORK(1) = SCALE = WORK(2) / WORK(1) is the scaling factor such  *>          WORK(1) = SCALE = WORK(2) / WORK(1) is the scaling factor such
 *>                    that SCALE*SVA(1:N) are the computed singular values  *>                    that SCALE*SVA(1:N) are the computed singular values
 *>                    of A. (See the description of SVA().)  *>                    of A. (See the description of SVA().)
 *>          WORK(2) = See the description of WORK(1).  *>          WORK(2) = See the description of WORK(1).
 *>          WORK(3) = SCONDA is an estimate for the condition number of  *>          WORK(3) = SCONDA is an estimate for the condition number of
 *>                    column equilibrated A. (If JOBA .EQ. 'E' or 'G')  *>                    column equilibrated A. (If JOBA = 'E' or 'G')
 *>                    SCONDA is an estimate of DSQRT(||(R^t * R)^(-1)||_1).  *>                    SCONDA is an estimate of DSQRT(||(R^t * R)^(-1)||_1).
 *>                    It is computed using DPOCON. It holds  *>                    It is computed using DPOCON. It holds
 *>                    N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA  *>                    N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA
Line 297 Line 297
 *>                    triangular factor in the first QR factorization.  *>                    triangular factor in the first QR factorization.
 *>          WORK(5) = an estimate of the scaled condition number of the  *>          WORK(5) = an estimate of the scaled condition number of the
 *>                    triangular factor in the second QR factorization.  *>                    triangular factor in the second QR factorization.
 *>          The following two parameters are computed if JOBT .EQ. 'T'.  *>          The following two parameters are computed if JOBT = 'T'.
 *>          They are provided for a developer/implementer who is familiar  *>          They are provided for a developer/implementer who is familiar
 *>          with the details of the method.  *>          with the details of the method.
 *>  *>
Line 313 Line 313
 *>          Length of WORK to confirm proper allocation of work space.  *>          Length of WORK to confirm proper allocation of work space.
 *>          LWORK depends on the job:  *>          LWORK depends on the job:
 *>  *>
 *>          If only SIGMA is needed ( JOBU.EQ.'N', JOBV.EQ.'N' ) and  *>          If only SIGMA is needed (JOBU = 'N', JOBV = 'N') and
 *>            -> .. no scaled condition estimate required (JOBE.EQ.'N'):  *>            -> .. no scaled condition estimate required (JOBE = 'N'):
 *>               LWORK >= max(2*M+N,4*N+1,7). This is the minimal requirement.  *>               LWORK >= max(2*M+N,4*N+1,7). This is the minimal requirement.
 *>               ->> For optimal performance (blocked code) the optimal value  *>               ->> For optimal performance (blocked code) the optimal value
 *>               is LWORK >= max(2*M+N,3*N+(N+1)*NB,7). Here NB is the optimal  *>               is LWORK >= max(2*M+N,3*N+(N+1)*NB,7). Here NB is the optimal
Line 330 Line 330
 *>               LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DGEQRF),  *>               LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DGEQRF),
 *>                                                     N+N*N+LWORK(DPOCON),7).  *>                                                     N+N*N+LWORK(DPOCON),7).
 *>  *>
 *>          If SIGMA and the right singular vectors are needed (JOBV.EQ.'V'),  *>          If SIGMA and the right singular vectors are needed (JOBV = 'V'),
 *>            -> the minimal requirement is LWORK >= max(2*M+N,4*N+1,7).  *>            -> the minimal requirement is LWORK >= max(2*M+N,4*N+1,7).
 *>            -> For optimal performance, LWORK >= max(2*M+N,3*N+(N+1)*NB,7),  *>            -> For optimal performance, LWORK >= max(2*M+N,3*N+(N+1)*NB,7),
 *>               where NB is the optimal block size for DGEQP3, DGEQRF, DGELQF,  *>               where NB is the optimal block size for DGEQP3, DGEQRF, DGELQF,
Line 341 Line 341
 *>          If SIGMA and the left singular vectors are needed  *>          If SIGMA and the left singular vectors are needed
 *>            -> the minimal requirement is LWORK >= max(2*M+N,4*N+1,7).  *>            -> the minimal requirement is LWORK >= max(2*M+N,4*N+1,7).
 *>            -> For optimal performance:  *>            -> For optimal performance:
 *>               if JOBU.EQ.'U' :: LWORK >= max(2*M+N,3*N+(N+1)*NB,7),  *>               if JOBU = 'U' :: LWORK >= max(2*M+N,3*N+(N+1)*NB,7),
 *>               if JOBU.EQ.'F' :: LWORK >= max(2*M+N,3*N+(N+1)*NB,N+M*NB,7),  *>               if JOBU = 'F' :: LWORK >= max(2*M+N,3*N+(N+1)*NB,N+M*NB,7),
 *>               where NB is the optimal block size for DGEQP3, DGEQRF, DORMQR.  *>               where NB is the optimal block size for DGEQP3, DGEQRF, DORMQR.
 *>               In general, the optimal length LWORK is computed as  *>               In general, the optimal length LWORK is computed as
 *>               LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DPOCON),  *>               LWORK >= max(2*M+N,N+LWORK(DGEQP3),N+LWORK(DPOCON),
 *>                        2*N+LWORK(DGEQRF), N+LWORK(DORMQR)).  *>                        2*N+LWORK(DGEQRF), N+LWORK(DORMQR)).
 *>               Here LWORK(DORMQR) equals N*NB (for JOBU.EQ.'U') or  *>               Here LWORK(DORMQR) equals N*NB (for JOBU = 'U') or
 *>               M*NB (for JOBU.EQ.'F').  *>               M*NB (for JOBU = 'F').
 *>  *>
 *>          If the full SVD is needed: (JOBU.EQ.'U' or JOBU.EQ.'F') and  *>          If the full SVD is needed: (JOBU = 'U' or JOBU = 'F') and
 *>            -> if JOBV.EQ.'V'  *>            -> if JOBV = 'V'
 *>               the minimal requirement is LWORK >= max(2*M+N,6*N+2*N*N).  *>               the minimal requirement is LWORK >= max(2*M+N,6*N+2*N*N).
 *>            -> if JOBV.EQ.'J' the minimal requirement is  *>            -> if JOBV = 'J' the minimal requirement is
 *>               LWORK >= max(2*M+N, 4*N+N*N,2*N+N*N+6).  *>               LWORK >= max(2*M+N, 4*N+N*N,2*N+N*N+6).
 *>            -> For optimal performance, LWORK should be additionally  *>            -> For optimal performance, LWORK should be additionally
 *>               larger than N+M*NB, where NB is the optimal block size  *>               larger than N+M*NB, where NB is the optimal block size
Line 369 Line 369
 *>                     of JOBA and JOBR.  *>                     of JOBA and JOBR.
 *>          IWORK(2) = the number of the computed nonzero singular values  *>          IWORK(2) = the number of the computed nonzero singular values
 *>          IWORK(3) = if nonzero, a warning message:  *>          IWORK(3) = if nonzero, a warning message:
 *>                     If IWORK(3).EQ.1 then some of the column norms of A  *>                     If IWORK(3) = 1 then some of the column norms of A
 *>                     were denormalized floats. The requested high accuracy  *>                     were denormalized floats. The requested high accuracy
 *>                     is not warranted by the data.  *>                     is not warranted by the data.
 *> \endverbatim  *> \endverbatim
Line 377 Line 377
 *> \param[out] INFO  *> \param[out] INFO
 *> \verbatim  *> \verbatim
 *>          INFO is INTEGER  *>          INFO is INTEGER
 *>           < 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 :  successful exit;  *>           = 0:  successful exit;
 *>           > 0 :  DGEJSV  did not converge in the maximal allowed number  *>           > 0:  DGEJSV  did not converge in the maximal allowed number
 *>                  of sweeps. The computed values may be inaccurate.  *>                 of sweeps. The computed values may be inaccurate.
 *> \endverbatim  *> \endverbatim
 *  *
 *  Authors:  *  Authors:
Line 953 Line 953
       IF ( L2ABER ) THEN        IF ( L2ABER ) THEN
 *        Standard absolute error bound suffices. All sigma_i with  *        Standard absolute error bound suffices. All sigma_i with
 *        sigma_i < N*EPSLN*||A|| are flushed to zero. This is an  *        sigma_i < N*EPSLN*||A|| are flushed to zero. This is an
 *        agressive enforcement of lower numerical rank by introducing a  *        aggressive enforcement of lower numerical rank by introducing a
 *        backward error of the order of N*EPSLN*||A||.  *        backward error of the order of N*EPSLN*||A||.
          TEMP1 = DSQRT(DBLE(N))*EPSLN           TEMP1 = DSQRT(DBLE(N))*EPSLN
          DO 3001 p = 2, N           DO 3001 p = 2, N
Line 965 Line 965
  3001    CONTINUE   3001    CONTINUE
  3002    CONTINUE   3002    CONTINUE
       ELSE IF ( L2RANK ) THEN        ELSE IF ( L2RANK ) THEN
 *        .. similarly as above, only slightly more gentle (less agressive).  *        .. similarly as above, only slightly more gentle (less aggressive).
 *        Sudden drop on the diagonal of R1 is used as the criterion for  *        Sudden drop on the diagonal of R1 is used as the criterion for
 *        close-to-rank-deficient.  *        close-to-rank-deficient.
          TEMP1 = DSQRT(SFMIN)           TEMP1 = DSQRT(SFMIN)
Line 1294 Line 1294
             CALL DPOCON('Lower',NR,WORK(2*N+1),NR,ONE,TEMP1,              CALL DPOCON('Lower',NR,WORK(2*N+1),NR,ONE,TEMP1,
      $                   WORK(2*N+NR*NR+1),IWORK(M+2*N+1),IERR)       $                   WORK(2*N+NR*NR+1),IWORK(M+2*N+1),IERR)
             CONDR1 = ONE / DSQRT(TEMP1)              CONDR1 = ONE / DSQRT(TEMP1)
 *           .. here need a second oppinion on the condition number  *           .. here need a second opinion on the condition number
 *           .. then assume worst case scenario  *           .. then assume worst case scenario
 *           R1 is OK for inverse <=> CONDR1 .LT. DBLE(N)  *           R1 is OK for inverse <=> CONDR1 .LT. DBLE(N)
 *           more conservative    <=> CONDR1 .LT. DSQRT(DBLE(N))  *           more conservative    <=> CONDR1 .LT. DSQRT(DBLE(N))
Line 1335 Line 1335
             ELSE              ELSE
 *  *
 *              .. ill-conditioned case: second QRF with pivoting  *              .. ill-conditioned case: second QRF with pivoting
 *              Note that windowed pivoting would be equaly good  *              Note that windowed pivoting would be equally good
 *              numerically, and more run-time efficient. So, in  *              numerically, and more run-time efficient. So, in
 *              an optimal implementation, the next call to DGEQP3  *              an optimal implementation, the next call to DGEQP3
 *              should be replaced with eg. CALL SGEQPX (ACM TOMS #782)  *              should be replaced with eg. CALL SGEQPX (ACM TOMS #782)
Line 1388 Line 1388
 *  *
                IF ( CONDR2 .GE. COND_OK ) THEN                 IF ( CONDR2 .GE. COND_OK ) THEN
 *                 .. save the Householder vectors used for Q3  *                 .. save the Householder vectors used for Q3
 *                 (this overwrittes the copy of R2, as it will not be  *                 (this overwrites the copy of R2, as it will not be
 *                 needed in this branch, but it does not overwritte the  *                 needed in this branch, but it does not overwritte the
 *                 Huseholder vectors of Q2.).  *                 Huseholder vectors of Q2.).
                   CALL DLACPY( 'U', NR, NR, V, LDV, WORK(2*N+1), N )                    CALL DLACPY( 'U', NR, NR, V, LDV, WORK(2*N+1), N )
Line 1638 Line 1638
 *  *
 *        This branch deploys a preconditioned Jacobi SVD with explicitly  *        This branch deploys a preconditioned Jacobi SVD with explicitly
 *        accumulated rotations. It is included as optional, mainly for  *        accumulated rotations. It is included as optional, mainly for
 *        experimental purposes. It does perfom well, and can also be used.  *        experimental purposes. It does perform well, and can also be used.
 *        In this implementation, this branch will be automatically activated  *        In this implementation, this branch will be automatically activated
 *        if the  condition number sigma_max(A) / sigma_min(A) is predicted  *        if the  condition number sigma_max(A) / sigma_min(A) is predicted
 *        to be greater than the overflow threshold. This is because the  *        to be greater than the overflow threshold. This is because the
Removed from v.1.19  
changed lines
  Added in v.1.20

CVSweb interface <joel.bertrand@systella.fr>