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

Diff for /rpl/lapack/lapack/zgejsv.f between versions 1.6 and 1.9

version 1.6, 2018/05/29 06:55:22 version 1.9, 2023/08/07 08:39:17
Line 80 Line 80
 *>              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=B*D. If A has heavily weighted  *>              condition number of B, where A=B*D. 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.
 *>       = 'A': Small singular values are not well determined by the data   *>       = 'A': Small singular values are not well determined by the data 
 *>              and are considered as noisy; the matrix is treated as  *>              and are considered as noisy; the matrix is treated as
 *>              numerically rank defficient. The error in the computed  *>              numerically rank deficient. The error in the computed
 *>              singular values is bounded by f(m,n)*epsilon*||A||.  *>              singular values is bounded by f(m,n)*epsilon*||A||.
 *>              The computed SVD A = U * S * V^* restores A up to  *>              The computed SVD A = U * S * V^* restores A up to
 *>              f(m,n)*epsilon*||A||.  *>              f(m,n)*epsilon*||A||.
Line 117 Line 117
 *>       = 'V': N columns of V are returned in the array V; Jacobi rotations  *>       = 'V': N columns of V are returned in the array V; Jacobi rotations
 *>              are not explicitly accumulated.  *>              are not explicitly accumulated.
 *>       = 'J': N columns of V are returned in the array V, but they are  *>       = 'J': N columns of V are returned in the array V, but they are
 *>              computed as the product of Jacobi rotations, if JOBT .EQ. 'N'.  *>              computed as the product of Jacobi rotations, if JOBT = 'N'.
 *>       = 'W': V may be used as workspace of length N*N. See the description  *>       = 'W': V may be used as workspace of length N*N. See the description
 *>              of V.  *>              of V.
 *>       = 'N': V is not computed.  *>       = 'N': V is not computed.
Line 131 Line 131
 *>         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 SQRT(BIG), BIG=DLAMCH('O'), then JOBR issues  *>         value of c*A is around SQRT(BIG), BIG=DLAMCH('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
 *>         SQRT(SFMIN) (for JOBR.EQ.'R'), or less than SMALL=SFMIN/EPSLN,  *>         SQRT(SFMIN) (for JOBR = 'R'), or less than SMALL=SFMIN/EPSLN,
 *>         where SFMIN=DLAMCH('S'), EPSLN=DLAMCH('E').  *>         where SFMIN=DLAMCH('S'), EPSLN=DLAMCH('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 229 Line 229
 *>          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^*. In that case, [V] is computed  *>                         replaces A with A^*. In that case, [V] is computed
 *>                         in U as left singular vectors of A^* and then  *>                         in U as left singular vectors of A^* and then
Line 251 Line 251
 *>          V is COMPLEX*16 array, dimension ( LDV, N )  *>          V is COMPLEX*16 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^*. In that case, [U] is computed  *>                         replaces A with A^*. In that case, [U] is computed
 *>                         in V as right singular vectors of A^* and then  *>                         in V as right singular vectors of A^* and then
Line 282 Line 282
 *>          Length of CWORK to confirm proper allocation of workspace.  *>          Length of CWORK to confirm proper allocation of workspace.
 *>          LWORK depends on the job:  *>          LWORK depends on the job:
 *>  *>
 *>          1. If only SIGMA is needed ( JOBU.EQ.'N', JOBV.EQ.'N' ) and  *>          1. If only SIGMA is needed ( JOBU = 'N', JOBV = 'N' ) and
 *>            1.1 .. no scaled condition estimate required (JOBA.NE.'E'.AND.JOBA.NE.'G'):  *>            1.1 .. no scaled condition estimate required (JOBA.NE.'E'.AND.JOBA.NE.'G'):
 *>               LWORK >= 2*N+1. This is the minimal requirement.  *>               LWORK >= 2*N+1. This is the minimal requirement.
 *>               ->> For optimal performance (blocked code) the optimal value  *>               ->> For optimal performance (blocked code) the optimal value
Line 298 Line 298
 *>               In general, the optimal length LWORK is computed as  *>               In general, the optimal length LWORK is computed as
 *>               LWORK >= max(N+LWORK(ZGEQP3),N+LWORK(ZGEQRF), LWORK(ZGESVJ),  *>               LWORK >= max(N+LWORK(ZGEQP3),N+LWORK(ZGEQRF), LWORK(ZGESVJ),
 *>                            N*N+LWORK(ZPOCON)).  *>                            N*N+LWORK(ZPOCON)).
 *>          2. If SIGMA and the right singular vectors are needed (JOBV.EQ.'V'),  *>          2. If SIGMA and the right singular vectors are needed (JOBV = 'V'),
 *>             (JOBU.EQ.'N')  *>             (JOBU = 'N')
 *>            2.1   .. no scaled condition estimate requested (JOBE.EQ.'N'):      *>            2.1   .. no scaled condition estimate requested (JOBE = 'N'):    
 *>            -> the minimal requirement is LWORK >= 3*N.  *>            -> the minimal requirement is LWORK >= 3*N.
 *>            -> For optimal performance,   *>            -> For optimal performance, 
 *>               LWORK >= max(N+(N+1)*NB, 2*N+N*NB)=2*N+N*NB,  *>               LWORK >= max(N+(N+1)*NB, 2*N+N*NB)=2*N+N*NB,
Line 318 Line 318
 *>               LWORK >= max(N+LWORK(ZGEQP3), LWORK(ZPOCON), N+LWORK(ZGESVJ),  *>               LWORK >= max(N+LWORK(ZGEQP3), LWORK(ZPOCON), N+LWORK(ZGESVJ),
 *>                       N+LWORK(ZGELQF), 2*N+LWORK(ZGEQRF), N+LWORK(ZUNMLQ)).     *>                       N+LWORK(ZGELQF), 2*N+LWORK(ZGEQRF), N+LWORK(ZUNMLQ)).   
 *>          3. If SIGMA and the left singular vectors are needed  *>          3. If SIGMA and the left singular vectors are needed
 *>            3.1  .. no scaled condition estimate requested (JOBE.EQ.'N'):  *>            3.1  .. no scaled condition estimate requested (JOBE = 'N'):
 *>            -> the minimal requirement is LWORK >= 3*N.  *>            -> the minimal requirement is LWORK >= 3*N.
 *>            -> For optimal performance:  *>            -> For optimal performance:
 *>               if JOBU.EQ.'U' :: LWORK >= max(3*N, N+(N+1)*NB, 2*N+N*NB)=2*N+N*NB,  *>               if JOBU = 'U' :: LWORK >= max(3*N, N+(N+1)*NB, 2*N+N*NB)=2*N+N*NB,
 *>               where NB is the optimal block size for ZGEQP3, ZGEQRF, ZUNMQR.  *>               where NB is the optimal block size for ZGEQP3, ZGEQRF, ZUNMQR.
 *>               In general, the optimal length LWORK is computed as  *>               In general, the optimal length LWORK is computed as
 *>               LWORK >= max(N+LWORK(ZGEQP3), 2*N+LWORK(ZGEQRF), N+LWORK(ZUNMQR)).   *>               LWORK >= max(N+LWORK(ZGEQP3), 2*N+LWORK(ZGEQRF), N+LWORK(ZUNMQR)). 
Line 329 Line 329
 *>               required (JOBA='E', or 'G').  *>               required (JOBA='E', or 'G').
 *>            -> the minimal requirement is LWORK >= 3*N.  *>            -> the minimal requirement is LWORK >= 3*N.
 *>            -> For optimal performance:  *>            -> For optimal performance:
 *>               if JOBU.EQ.'U' :: LWORK >= max(3*N, N+(N+1)*NB, 2*N+N*NB)=2*N+N*NB,  *>               if JOBU = 'U' :: LWORK >= max(3*N, N+(N+1)*NB, 2*N+N*NB)=2*N+N*NB,
 *>               where NB is the optimal block size for ZGEQP3, ZGEQRF, ZUNMQR.  *>               where NB is the optimal block size for ZGEQP3, ZGEQRF, ZUNMQR.
 *>               In general, the optimal length LWORK is computed as  *>               In general, the optimal length LWORK is computed as
 *>               LWORK >= max(N+LWORK(ZGEQP3),N+LWORK(ZPOCON),  *>               LWORK >= max(N+LWORK(ZGEQP3),N+LWORK(ZPOCON),
 *>                        2*N+LWORK(ZGEQRF), N+LWORK(ZUNMQR)).  *>                        2*N+LWORK(ZGEQRF), N+LWORK(ZUNMQR)).
 *>          4. If the full SVD is needed: (JOBU.EQ.'U' or JOBU.EQ.'F') and   *>          4. If the full SVD is needed: (JOBU = 'U' or JOBU = 'F') and 
 *>            4.1. if JOBV.EQ.'V'    *>            4.1. if JOBV = 'V'  
 *>               the minimal requirement is LWORK >= 5*N+2*N*N.   *>               the minimal requirement is LWORK >= 5*N+2*N*N. 
 *>            4.2. if JOBV.EQ.'J' the minimal requirement is   *>            4.2. if JOBV = 'J' the minimal requirement is 
 *>               LWORK >= 4*N+N*N.  *>               LWORK >= 4*N+N*N.
 *>            In both cases, the allocated CWORK can accommodate blocked runs  *>            In both cases, the allocated CWORK can accommodate blocked runs
 *>            of ZGEQP3, ZGEQRF, ZGELQF, SUNMQR, ZUNMLQ.  *>            of ZGEQP3, ZGEQRF, ZGELQF, SUNMQR, ZUNMLQ.
Line 356 Line 356
 *>                    of A. (See the description of SVA().)  *>                    of A. (See the description of SVA().)
 *>          RWORK(2) = See the description of RWORK(1).  *>          RWORK(2) = See the description of RWORK(1).
 *>          RWORK(3) = SCONDA is an estimate for the condition number of  *>          RWORK(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 SQRT(||(R^* * R)^(-1)||_1).  *>                    SCONDA is an estimate of SQRT(||(R^* * R)^(-1)||_1).
 *>                    It is computed using SPOCON. It holds  *>                    It is computed using ZPOCON. It holds
 *>                    N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA  *>                    N^(-1/4) * SCONDA <= ||R^(-1)||_2 <= N^(1/4) * SCONDA
 *>                    where R is the triangular factor from the QRF of A.  *>                    where R is the triangular factor from the QRF of A.
 *>                    However, if R is truncated and the numerical rank is  *>                    However, if R is truncated and the numerical rank is
Line 367 Line 367
 *>                    singular values might be lost.  *>                    singular values might be lost.
 *>  *>
 *>          If full SVD is needed, the following two condition numbers are  *>          If full SVD is needed, the following two condition numbers are
 *>          useful for the analysis of the algorithm. They are provied for  *>          useful for the analysis of the algorithm. They are provided for
 *>          a developer/implementer who is familiar with the details of  *>          a developer/implementer who is familiar with the details of
 *>          the method.  *>          the method.
 *>  *>
Line 375 Line 375
 *>                    triangular factor in the first QR factorization.  *>                    triangular factor in the first QR factorization.
 *>          RWORK(5) = an estimate of the scaled condition number of the  *>          RWORK(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.
 *>          RWORK(6) = the entropy of A^* * A :: this is the Shannon entropy  *>          RWORK(6) = the entropy of A^* * A :: this is the Shannon entropy
Line 456 Line 456
 *>                     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.
 *>          IWORK(4) = 1 or -1. If IWORK(4) .EQ. 1, then the procedure used A^* to  *>          IWORK(4) = 1 or -1. If IWORK(4) = 1, then the procedure used A^* to
 *>                     do the job as specified by the JOB parameters.  *>                     do the job as specified by the JOB parameters.
 *>          If the call to ZGEJSV is a workspace query (indicated by LWORK .EQ. -1 or  *>          If the call to ZGEJSV is a workspace query (indicated by LWORK = -1 or
 *>          LRWORK .EQ. -1), then on exit IWORK(1) contains the required length of   *>          LRWORK = -1), then on exit IWORK(1) contains the required length of 
 *>          IWORK for the job parameters used in the call.  *>          IWORK for the job parameters used in the call.
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \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 :  ZGEJSV  did not converge in the maximal allowed number  *>           > 0:  ZGEJSV  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 483 Line 483
 *> \author Univ. of Colorado Denver  *> \author Univ. of Colorado Denver
 *> \author NAG Ltd.  *> \author NAG Ltd.
 *  *
 *> \date June 2016  
 *  
 *> \ingroup complex16GEsing  *> \ingroup complex16GEsing
 *  *
 *> \par Further Details:  *> \par Further Details:
Line 569 Line 567
      $                   M, N, A, LDA, SVA, U, LDU, V, LDV,       $                   M, N, A, LDA, SVA, U, LDU, V, LDV,
      $                   CWORK, LWORK, RWORK, LRWORK, IWORK, INFO )       $                   CWORK, LWORK, RWORK, LRWORK, IWORK, INFO )
 *  *
 *  -- LAPACK computational routine (version 3.7.1) --  *  -- 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..--
 *     June 2017  
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       IMPLICIT    NONE        IMPLICIT    NONE
Line 704 Line 701
           LWSVDJ  = MAX( 2 * N, 1 )                     LWSVDJ  = MAX( 2 * N, 1 )         
           LWSVDJV = MAX( 2 * N, 1 )            LWSVDJV = MAX( 2 * N, 1 )
 *         .. minimal REAL workspace length for ZGEQP3, ZPOCON, ZGESVJ  *         .. minimal REAL workspace length for ZGEQP3, ZPOCON, ZGESVJ
           LRWQP3  = N             LRWQP3  = 2 * N 
           LRWCON  = N             LRWCON  = N 
           LRWSVDJ = N             LRWSVDJ = N 
           IF ( LQUERY ) THEN             IF ( LQUERY ) THEN 
               CALL ZGEQP3( M, N, A, LDA, IWORK, CDUMMY, CDUMMY, -1,                 CALL ZGEQP3( M, N, A, LDA, IWORK, CDUMMY, CDUMMY, -1, 
      $             RDUMMY, IERR )       $             RDUMMY, IERR )
               LWRK_ZGEQP3 = CDUMMY(1)                LWRK_ZGEQP3 = INT( CDUMMY(1) )
               CALL ZGEQRF( N, N, A, LDA, CDUMMY, CDUMMY,-1, IERR )                CALL ZGEQRF( N, N, A, LDA, CDUMMY, CDUMMY,-1, IERR )
               LWRK_ZGEQRF = CDUMMY(1)                LWRK_ZGEQRF = INT( CDUMMY(1) )
               CALL ZGELQF( N, N, A, LDA, CDUMMY, CDUMMY,-1, IERR )                CALL ZGELQF( N, N, A, LDA, CDUMMY, CDUMMY,-1, IERR )
               LWRK_ZGELQF = CDUMMY(1)                             LWRK_ZGELQF = INT( CDUMMY(1) )
           END IF            END IF
           MINWRK  = 2            MINWRK  = 2
           OPTWRK  = 2            OPTWRK  = 2
Line 730 Line 727
               IF ( LQUERY ) THEN                 IF ( LQUERY ) THEN 
                   CALL ZGESVJ( 'L', 'N', 'N', N, N, A, LDA, SVA, N, V,                     CALL ZGESVJ( 'L', 'N', 'N', N, N, A, LDA, SVA, N, V, 
      $                 LDV, CDUMMY, -1, RDUMMY, -1, IERR )       $                 LDV, CDUMMY, -1, RDUMMY, -1, IERR )
                   LWRK_ZGESVJ = CDUMMY(1)                    LWRK_ZGESVJ = INT( CDUMMY(1) )
                   IF ( ERREST ) THEN                     IF ( ERREST ) THEN 
                       OPTWRK = MAX( N+LWRK_ZGEQP3, N**2+LWCON,                         OPTWRK = MAX( N+LWRK_ZGEQP3, N**2+LWCON, 
      $                              N+LWRK_ZGEQRF, LWRK_ZGESVJ )       $                              N+LWRK_ZGEQRF, LWRK_ZGESVJ )
Line 766 Line 763
              IF ( LQUERY ) THEN               IF ( LQUERY ) THEN
                  CALL ZGESVJ( 'L', 'U', 'N', N,N, U, LDU, SVA, N, A,                   CALL ZGESVJ( 'L', 'U', 'N', N,N, U, LDU, SVA, N, A,
      $                LDA, CDUMMY, -1, RDUMMY, -1, IERR )       $                LDA, CDUMMY, -1, RDUMMY, -1, IERR )
                  LWRK_ZGESVJ = CDUMMY(1)                   LWRK_ZGESVJ = INT( CDUMMY(1) )
                  CALL ZUNMLQ( 'L', 'C', N, N, N, A, LDA, CDUMMY,                   CALL ZUNMLQ( 'L', 'C', N, N, N, A, LDA, CDUMMY,
      $                V, LDV, CDUMMY, -1, IERR )       $                V, LDV, CDUMMY, -1, IERR )
                  LWRK_ZUNMLQ = CDUMMY(1)                                   LWRK_ZUNMLQ = INT( CDUMMY(1) )
                  IF ( ERREST ) THEN                    IF ( ERREST ) THEN 
                  OPTWRK = MAX( N+LWRK_ZGEQP3, LWCON, LWRK_ZGESVJ,                    OPTWRK = MAX( N+LWRK_ZGEQP3, LWCON, LWRK_ZGESVJ, 
      $                         N+LWRK_ZGELQF, 2*N+LWRK_ZGEQRF,       $                         N+LWRK_ZGELQF, 2*N+LWRK_ZGEQRF,
Line 805 Line 802
              IF ( LQUERY ) THEN               IF ( LQUERY ) THEN
                  CALL ZGESVJ( 'L', 'U', 'N', N,N, U, LDU, SVA, N, A,                   CALL ZGESVJ( 'L', 'U', 'N', N,N, U, LDU, SVA, N, A,
      $                LDA, CDUMMY, -1, RDUMMY, -1, IERR )       $                LDA, CDUMMY, -1, RDUMMY, -1, IERR )
                  LWRK_ZGESVJ = CDUMMY(1)                   LWRK_ZGESVJ = INT( CDUMMY(1) )
                  CALL ZUNMQR( 'L', 'N', M, N, N, A, LDA, CDUMMY, U,                   CALL ZUNMQR( 'L', 'N', M, N, N, A, LDA, CDUMMY, U,
      $               LDU, CDUMMY, -1, IERR )       $               LDU, CDUMMY, -1, IERR )
                  LWRK_ZUNMQRM = CDUMMY(1)                   LWRK_ZUNMQRM = INT( CDUMMY(1) )
                  IF ( ERREST ) THEN                   IF ( ERREST ) THEN
                  OPTWRK = N + MAX( LWRK_ZGEQP3, LWCON, N+LWRK_ZGEQRF,                   OPTWRK = N + MAX( LWRK_ZGEQP3, LWCON, N+LWRK_ZGEQRF,
      $                             LWRK_ZGESVJ, LWRK_ZUNMQRM )       $                             LWRK_ZGESVJ, LWRK_ZUNMQRM )
Line 867 Line 864
              IF ( LQUERY ) THEN               IF ( LQUERY ) THEN
                  CALL ZUNMQR( 'L', 'N', M, N, N, A, LDA, CDUMMY, U,                   CALL ZUNMQR( 'L', 'N', M, N, N, A, LDA, CDUMMY, U,
      $                LDU, CDUMMY, -1, IERR )       $                LDU, CDUMMY, -1, IERR )
                  LWRK_ZUNMQRM = CDUMMY(1)                   LWRK_ZUNMQRM = INT( CDUMMY(1) )
                  CALL ZUNMQR( 'L', 'N', N, N, N, A, LDA, CDUMMY, U,                   CALL ZUNMQR( 'L', 'N', N, N, N, A, LDA, CDUMMY, U,
      $                LDU, CDUMMY, -1, IERR )       $                LDU, CDUMMY, -1, IERR )
                  LWRK_ZUNMQR = CDUMMY(1)                   LWRK_ZUNMQR = INT( CDUMMY(1) )
                  IF ( .NOT. JRACC ) THEN                   IF ( .NOT. JRACC ) THEN
                      CALL ZGEQP3( N,N, A, LDA, IWORK, CDUMMY,CDUMMY, -1,                       CALL ZGEQP3( N,N, A, LDA, IWORK, CDUMMY,CDUMMY, -1,
      $                    RDUMMY, IERR )       $                    RDUMMY, IERR )
                      LWRK_ZGEQP3N = CDUMMY(1)                       LWRK_ZGEQP3N = INT( CDUMMY(1) )
                      CALL ZGESVJ( 'L', 'U', 'N', N, N, U, LDU, SVA,                       CALL ZGESVJ( 'L', 'U', 'N', N, N, U, LDU, SVA,
      $                    N, V, LDV, CDUMMY, -1, RDUMMY, -1, IERR )       $                    N, V, LDV, CDUMMY, -1, RDUMMY, -1, IERR )
                      LWRK_ZGESVJ = CDUMMY(1)                       LWRK_ZGESVJ = INT( CDUMMY(1) )
                      CALL ZGESVJ( 'U', 'U', 'N', N, N, U, LDU, SVA,                       CALL ZGESVJ( 'U', 'U', 'N', N, N, U, LDU, SVA,
      $                    N, V, LDV, CDUMMY, -1, RDUMMY, -1, IERR )       $                    N, V, LDV, CDUMMY, -1, RDUMMY, -1, IERR )
                      LWRK_ZGESVJU = CDUMMY(1)                       LWRK_ZGESVJU = INT( CDUMMY(1) )
                      CALL ZGESVJ( 'L', 'U', 'V', N, N, U, LDU, SVA,                       CALL ZGESVJ( 'L', 'U', 'V', N, N, U, LDU, SVA,
      $                    N, V, LDV, CDUMMY, -1, RDUMMY, -1, IERR )       $                    N, V, LDV, CDUMMY, -1, RDUMMY, -1, IERR )
                      LWRK_ZGESVJV = CDUMMY(1)                       LWRK_ZGESVJV = INT( CDUMMY(1) )
                      CALL ZUNMLQ( 'L', 'C', N, N, N, A, LDA, CDUMMY,                       CALL ZUNMLQ( 'L', 'C', N, N, N, A, LDA, CDUMMY,
      $                    V, LDV, CDUMMY, -1, IERR )       $                    V, LDV, CDUMMY, -1, IERR )
                      LWRK_ZUNMLQ = CDUMMY(1)                       LWRK_ZUNMLQ = INT( CDUMMY(1) )
                      IF ( ERREST ) THEN                        IF ( ERREST ) THEN 
                        OPTWRK = MAX( N+LWRK_ZGEQP3, N+LWCON,                          OPTWRK = MAX( N+LWRK_ZGEQP3, N+LWCON, 
      $                          2*N+N**2+LWCON, 2*N+LWRK_ZGEQRF,        $                          2*N+N**2+LWCON, 2*N+LWRK_ZGEQRF, 
Line 915 Line 912
                  ELSE                   ELSE
                      CALL ZGESVJ( 'L', 'U', 'V', N, N, U, LDU, SVA,                       CALL ZGESVJ( 'L', 'U', 'V', N, N, U, LDU, SVA,
      $                    N, V, LDV, CDUMMY, -1, RDUMMY, -1, IERR )       $                    N, V, LDV, CDUMMY, -1, RDUMMY, -1, IERR )
                      LWRK_ZGESVJV = CDUMMY(1)                       LWRK_ZGESVJV = INT( CDUMMY(1) )
                      CALL ZUNMQR( 'L', 'N', N, N, N, CDUMMY, N, CDUMMY,                       CALL ZUNMQR( 'L', 'N', N, N, N, CDUMMY, N, CDUMMY,
      $                    V, LDV, CDUMMY, -1, IERR )       $                    V, LDV, CDUMMY, -1, IERR )
                      LWRK_ZUNMQR = CDUMMY(1)                       LWRK_ZUNMQR = INT( CDUMMY(1) )
                      CALL ZUNMQR( 'L', 'N', M, N, N, A, LDA, CDUMMY, U,                       CALL ZUNMQR( 'L', 'N', M, N, N, A, LDA, CDUMMY, U,
      $                    LDU, CDUMMY, -1, IERR )       $                    LDU, CDUMMY, -1, IERR )
                      LWRK_ZUNMQRM = CDUMMY(1)                          LWRK_ZUNMQRM = INT( CDUMMY(1) )
                      IF ( ERREST ) THEN                        IF ( ERREST ) THEN 
                         OPTWRK = MAX( N+LWRK_ZGEQP3, N+LWCON,                             OPTWRK = MAX( N+LWRK_ZGEQP3, N+LWCON,   
      $                           2*N+LWRK_ZGEQRF, 2*N+N**2,         $                           2*N+LWRK_ZGEQRF, 2*N+N**2,  
Line 942 Line 939
              END IF                END IF 
           END IF            END IF
           MINWRK = MAX( 2, MINWRK )            MINWRK = MAX( 2, MINWRK )
           OPTWRK = MAX( 2, OPTWRK )            OPTWRK = MAX( MINWRK, OPTWRK )
           IF ( LWORK  .LT. MINWRK  .AND. (.NOT.LQUERY) ) INFO = - 17            IF ( LWORK  .LT. MINWRK  .AND. (.NOT.LQUERY) ) INFO = - 17
           IF ( LRWORK .LT. MINRWRK .AND. (.NOT.LQUERY) ) INFO = - 19               IF ( LRWORK .LT. MINRWRK .AND. (.NOT.LQUERY) ) INFO = - 19   
       END IF        END IF
Line 1316 Line 1313
 *     (eg speed by replacing global with restricted window pivoting, such  *     (eg speed by replacing global with restricted window pivoting, such
 *     as in xGEQPX from TOMS # 782). Good results will be obtained using  *     as in xGEQPX from TOMS # 782). Good results will be obtained using
 *     xGEQPX with properly (!) chosen numerical parameters.  *     xGEQPX with properly (!) chosen numerical parameters.
 *     Any improvement of ZGEQP3 improves overal performance of ZGEJSV.  *     Any improvement of ZGEQP3 improves overall performance of ZGEJSV.
 *  *
 *     A * P1 = Q1 * [ R1^* 0]^*:  *     A * P1 = Q1 * [ R1^* 0]^*:
       DO 1963 p = 1, N        DO 1963 p = 1, N
Line 1338 Line 1335
       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 = SQRT(DBLE(N))*EPSLN           TEMP1 = SQRT(DBLE(N))*EPSLN
          DO 3001 p = 2, N           DO 3001 p = 2, N
Line 1350 Line 1347
  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 = SQRT(SFMIN)           TEMP1 = SQRT(SFMIN)
Line 1720 Line 1717
             CALL ZPOCON('L',NR,CWORK(2*N+1),NR,ONE,TEMP1,              CALL ZPOCON('L',NR,CWORK(2*N+1),NR,ONE,TEMP1,
      $                   CWORK(2*N+NR*NR+1),RWORK,IERR)       $                   CWORK(2*N+NR*NR+1),RWORK,IERR)
             CONDR1 = ONE / SQRT(TEMP1)              CONDR1 = ONE / SQRT(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. SQRT(DBLE(N))  *           more conservative    <=> CONDR1 .LT. SQRT(DBLE(N))
Line 1765 Line 1762
             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 ZGEQP3  *              an optimal implementation, the next call to ZGEQP3
 *              should be replaced with eg. CALL ZGEQPX (ACM TOMS #782)  *              should be replaced with eg. CALL ZGEQPX (ACM TOMS #782)
Line 1823 Line 1820
 *  *
                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 ZLACPY( 'U', NR, NR, V, LDV, CWORK(2*N+1), N )                    CALL ZLACPY( 'U', NR, NR, V, LDV, CWORK(2*N+1), N )
Line 2079 Line 2076
 *  *
 *        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.6  
changed lines
  Added in v.1.9

CVSweb interface <joel.bertrand@systella.fr>