[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.1 and 1.2

version 1.1, 2015/11/26 11:44:21 version 1.2, 2016/08/27 15:27:12
Line 27 Line 27
 *     INTEGER     INFO, LDA, LDU, LDV, LWORK, M, N  *     INTEGER     INFO, LDA, LDU, LDV, LWORK, M, N
 *     ..  *     ..
 *     .. Array Arguments ..  *     .. Array Arguments ..
 *     DOUBLE COMPLEX     A( LDA, * ),  U( LDU, * ), V( LDV, * ), CWORK( LWORK )  *     COMPLEX*16     A( LDA, * ),  U( LDU, * ), V( LDV, * ), CWORK( LWORK )
 *     DOUBLE PRECISION   SVA( N ), RWORK( LRWORK )        *     DOUBLE PRECISION   SVA( N ), RWORK( LRWORK )      
 *     INTEGER     IWORK( * )  *     INTEGER     IWORK( * )
 *     CHARACTER*1 JOBA, JOBP, JOBR, JOBT, JOBU, JOBV  *     CHARACTER*1 JOBA, JOBP, JOBR, JOBT, JOBU, JOBV
Line 39 Line 39
 *>  *>
 *> \verbatim  *> \verbatim
 *>  *>
 * ZGEJSV computes the singular value decomposition (SVD) of a real M-by-N  *> ZGEJSV computes the singular value decomposition (SVD) of a complex M-by-N
 * matrix [A], where M >= N. The SVD of [A] is written as  *> matrix [A], where M >= N. The SVD of [A] is written as
 *  *>
 *              [A] = [U] * [SIGMA] * [V]^*,  *>              [A] = [U] * [SIGMA] * [V]^*,
 *  *>
 * where [SIGMA] is an N-by-N (M-by-N) matrix which is zero except for its N  *> where [SIGMA] is an N-by-N (M-by-N) matrix which is zero except for its N
 * diagonal elements, [U] is an M-by-N (or M-by-M) orthonormal matrix, and  *> diagonal elements, [U] is an M-by-N (or M-by-M) unitary matrix, and
 * [V] is an N-by-N orthogonal matrix. The diagonal elements of [SIGMA] are  *> [V] is an N-by-N unitary matrix. The diagonal elements of [SIGMA] are
 * the singular values of [A]. The columns of [U] and [V] are the left and  *> the singular values of [A]. The columns of [U] and [V] are the left and
 * the right singular vectors of [A], respectively. The matrices [U] and [V]  *> the right singular vectors of [A], respectively. The matrices [U] and [V]
 * are computed and stored in the arrays U and V, respectively. The diagonal  *> are computed and stored in the arrays U and V, respectively. The diagonal
 * of [SIGMA] is computed and stored in the array SVA.  *> of [SIGMA] is computed and stored in the array SVA.
 *  *> \endverbatim
 *  Arguments:  *>
 *  ==========  *>  Arguments:
   *>  ==========
 *>  *>
 *> \param[in] JOBA  *> \param[in] JOBA
 *> \verbatim  *> \verbatim
Line 193 Line 194
 *>  *>
 *> \param[in,out] A  *> \param[in,out] A
 *> \verbatim  *> \verbatim
 *>          A is DOUBLE COMPLEX array, dimension (LDA,N)  *>          A is COMPLEX*16 array, dimension (LDA,N)
 *>          On entry, the M-by-N matrix A.  *>          On entry, the M-by-N matrix A.
 *> \endverbatim  *> \endverbatim
 *>  *>
Line 221 Line 222
 *>  *>
 *> \param[out] U  *> \param[out] U
 *> \verbatim  *> \verbatim
 *>          U is DOUBLE COMPLEX array, dimension ( LDU, N )  *>          U is COMPLEX*16 array, dimension ( LDU, N )
 *>          If JOBU = 'U', then U contains on exit the M-by-N matrix of  *>          If JOBU = 'U', then U contains on exit the M-by-N matrix of
 *>                         the left singular vectors.  *>                         the left singular vectors.
 *>          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
Line 234 Line 235
 *>                         copied back to the V array. This 'W' option is just  *>                         copied back to the V array. This 'W' option is just
 *>                         a reminder to the caller that in this case U is  *>                         a reminder to the caller that in this case U is
 *>                         reserved as workspace of length N*N.  *>                         reserved as workspace of length N*N.
 *>          If JOBU = 'N'  U is not referenced.  *>          If JOBU = 'N'  U is not referenced, unless JOBT='T'.
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \param[in] LDU  *> \param[in] LDU
Line 246 Line 247
 *>  *>
 *> \param[out] V  *> \param[out] V
 *> \verbatim  *> \verbatim
 *>          V is DOUBLE COMPLEX 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.EQ.'U' AND JOBT.EQ.'T' AND M.EQ.N),
Line 256 Line 257
 *>                         copied back to the U array. This 'W' option is just  *>                         copied back to the U array. This 'W' option is just
 *>                         a reminder to the caller that in this case V is  *>                         a reminder to the caller that in this case V is
 *>                         reserved as workspace of length N*N.  *>                         reserved as workspace of length N*N.
 *>          If JOBV = 'N'  V is not referenced.  *>          If JOBV = 'N'  V is not referenced, unless JOBT='T'.
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \param[in] LDV  *> \param[in] LDV
Line 268 Line 269
 *>  *>
 *> \param[out] CWORK  *> \param[out] CWORK
 *> \verbatim  *> \verbatim
 *> CWORK (workspace)  *>          CWORK is COMPLEX*16 array, dimension at least LWORK.     
 *>          CWORK is DOUBLE COMPLEX array, dimension at least LWORK.       
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \param[in] LWORK  *> \param[in] LWORK
Line 279 Line 279
 *>          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.EQ.'N', JOBV.EQ.'N' ) and
 *>            1.1 .. no scaled condition estimate required (JOBE.EQ.'N'):  *>            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
 *>               is LWORK >= N + (N+1)*NB. Here NB is the optimal  *>               is LWORK >= N + (N+1)*NB. Here NB is the optimal
Line 293 Line 293
 *>               is LWORK >= max(N+(N+1)*NB, N*N+3*N).  *>               is LWORK >= max(N+(N+1)*NB, N*N+3*N).
 *>               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 >= max(N+LWORK(ZGEQP3),N+LWORK(ZGEQRF), 
 *>                                                     N+N*N+LWORK(CPOCON)).  *>                                                     N+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.EQ.'V'),
 *>             (JOBU.EQ.'N')  *>             (JOBU.EQ.'N')
 *>            -> the minimal requirement is LWORK >= 3*N.  *>            -> the minimal requirement is LWORK >= 3*N.
 *>            -> For optimal performance, LWORK >= max(N+(N+1)*NB, 3*N,2*N+N*NB),  *>            -> For optimal performance, LWORK >= max(N+(N+1)*NB, 3*N,2*N+N*NB),
 *>               where NB is the optimal block size for ZGEQP3, ZGEQRF, ZGELQ,  *>               where NB is the optimal block size for ZGEQP3, ZGEQRF, ZGELQF,
 *>               CUNMLQ. In general, the optimal length LWORK is computed as  *>               ZUNMLQ. In general, the optimal length LWORK is computed as
 *>               LWORK >= max(N+LWORK(ZGEQP3), N+LWORK(CPOCON), N+LWORK(ZGESVJ),  *>               LWORK >= max(N+LWORK(ZGEQP3), N+LWORK(ZPOCON), N+LWORK(ZGESVJ),
 *>                       N+LWORK(ZGELQF), 2*N+LWORK(ZGEQRF), N+LWORK(CUNMLQ)).  *>                       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
 *>            -> 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),  *>               if JOBU.EQ.'U' :: LWORK >= max(3*N, N+(N+1)*NB, 2*N+N*NB),
 *>               where NB is the optimal block size for ZGEQP3, ZGEQRF, CUNMQR.  *>               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(CPOCON),  *>               LWORK >= max(N+LWORK(ZGEQP3),N+LWORK(ZPOCON),
 *>                        2*N+LWORK(ZGEQRF), N+LWORK(CUNMQR)).   *>                        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.EQ.'U' or JOBU.EQ.'F') and 
 *>            4.1. if JOBV.EQ.'V'    *>            4.1. if JOBV.EQ.'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.EQ.'J' the minimal requirement is 
 *>               LWORK >= 4*N+N*N.  *>               LWORK >= 4*N+N*N.
 *>            In both cases, the allocated CWORK can accomodate blocked runs  *>            In both cases, the allocated CWORK can accommodate blocked runs
 *>            of ZGEQP3, ZGEQRF, ZGELQF, SUNMQR, CUNMLQ.  *>            of ZGEQP3, ZGEQRF, ZGELQF, ZUNMQR, ZUNMLQ.
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \param[out] RWORK  *> \param[out] RWORK
Line 433 Line 433
 *> \author Univ. of Colorado Denver   *> \author Univ. of Colorado Denver 
 *> \author NAG Ltd.   *> \author NAG Ltd. 
 *  *
 *> \date November 2015  *> \date June 2016
 *  *
 *> \ingroup complex16GEsing  *> \ingroup complex16GEsing
 *  *
Line 491 Line 491
 *>  *>
 *> \verbatim  *> \verbatim
 *>  *>
 * [1] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I.  *> [1] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I.
 *     SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342.  *>     SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342.
 *     LAPACK Working note 169.  *>     LAPACK Working note 169.
 * [2] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II.  *> [2] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II.
 *     SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362.  *>     SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362.
 *     LAPACK Working note 170.  *>     LAPACK Working note 170.
 * [3] Z. Drmac and Z. Bujanovic: On the failure of rank-revealing QR  *> [3] Z. Drmac and Z. Bujanovic: On the failure of rank-revealing QR
 *     factorization software - a case study.  *>     factorization software - a case study.
 *     ACM Trans. math. Softw. Vol. 35, No 2 (2008), pp. 1-28.  *>     ACM Trans. Math. Softw. Vol. 35, No 2 (2008), pp. 1-28.
 *     LAPACK Working note 176.  *>     LAPACK Working note 176.
 * [4] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV,  *> [4] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV,
 *     QSVD, (H,K)-SVD computations.  *>     QSVD, (H,K)-SVD computations.
 *     Department of Mathematics, University of Zagreb, 2008.  *>     Department of Mathematics, University of Zagreb, 2008.
 *> \endverbatim  *> \endverbatim
 *  *
 *>  \par Bugs, examples and comments:  *>  \par Bugs, examples and comments:
Line 517 Line 517
      $                   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.6.0) --  *  -- LAPACK computational routine (version 3.6.1) --
 *  -- 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  *     June 2016
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       IMPLICIT    NONE        IMPLICIT    NONE
       INTEGER     INFO, LDA, LDU, LDV, LWORK, LRWORK, M, N        INTEGER     INFO, LDA, LDU, LDV, LWORK, LRWORK, M, N
 *     ..  *     ..
 *     .. Array Arguments ..  *     .. Array Arguments ..
       DOUBLE COMPLEX   A( LDA, * ), U( LDU, * ), V( LDV, * ),         COMPLEX*16       A( LDA, * ), U( LDU, * ), V( LDV, * ), 
      $                 CWORK( LWORK )       $                 CWORK( LWORK )
       DOUBLE PRECISION SVA( N ), RWORK( * )        DOUBLE PRECISION SVA( N ), RWORK( * )
       INTEGER          IWORK( * )        INTEGER          IWORK( * )
Line 539 Line 539
 *     .. Local Parameters ..  *     .. Local Parameters ..
       DOUBLE PRECISION ZERO,         ONE        DOUBLE PRECISION ZERO,         ONE
       PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )        PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
       DOUBLE COMPLEX            CZERO,       CONE        COMPLEX*16                CZERO,       CONE
       PARAMETER  ( CZERO = ( 0.0D0, 0.0D0 ), CONE = ( 1.0D0, 0.0D0 ) )        PARAMETER  ( CZERO = ( 0.0D0, 0.0D0 ), CONE = ( 1.0D0, 0.0D0 ) )
 *     ..  *     ..
 *     .. Local Scalars ..  *     .. Local Scalars ..
       DOUBLE COMPLEX   CTEMP        COMPLEX*16       CTEMP
       DOUBLE PRECISION AAPP,    AAQQ,   AATMAX, AATMIN, BIG,    BIG1,           DOUBLE PRECISION AAPP,    AAQQ,   AATMAX, AATMIN, BIG,    BIG1,   
      $                 COND_OK, CONDR1, CONDR2, ENTRA,  ENTRAT, EPSLN,         $                 COND_OK, CONDR1, CONDR2, ENTRA,  ENTRAT, EPSLN,  
      $                 MAXPRJ,  SCALEM, SCONDA, SFMIN,  SMALL,  TEMP1,         $                 MAXPRJ,  SCALEM, SCONDA, SFMIN,  SMALL,  TEMP1,  
Line 554 Line 554
      $        NOSCAL, ROWPIV, RSVEC,  TRANSP       $        NOSCAL, ROWPIV, RSVEC,  TRANSP
 *     ..  *     ..
 *     .. Intrinsic Functions ..  *     .. Intrinsic Functions ..
       INTRINSIC ABS,  DCMPLX, DCONJG, DLOG, DMAX1, DMIN1, DFLOAT,        INTRINSIC ABS,  DCMPLX, DCONJG, DLOG, DMAX1, DMIN1, DBLE,
      $          MAX0, MIN0, NINT,  DSQRT       $          MAX0, MIN0, NINT,  DSQRT
 *     ..  *     ..
 *     .. External Functions ..  *     .. External Functions ..
       DOUBLE PRECISION      DLAMCH, DZNRM2        DOUBLE PRECISION      DLAMCH, DZNRM2
       INTEGER   IDAMAX        INTEGER   IDAMAX, IZAMAX
       LOGICAL   LSAME        LOGICAL   LSAME
       EXTERNAL  IDAMAX, LSAME, DLAMCH, DZNRM2        EXTERNAL  IDAMAX, IZAMAX, LSAME, DLAMCH, DZNRM2
 *     ..  *     ..
 *     .. External Subroutines ..  *     .. External Subroutines ..
       EXTERNAL  ZCOPY,  ZGELQF, ZGEQP3, ZGEQRF, ZLACPY, ZLASCL,        EXTERNAL  DLASSQ, ZCOPY,  ZGELQF, ZGEQP3, ZGEQRF, ZLACPY, ZLASCL,
      $          ZLASET, ZLASSQ, ZLASWP, ZUNGQR, ZUNMLQ,       $          DLASCL, ZLASET, ZLASSQ, ZLASWP, ZUNGQR, ZUNMLQ,
      $          ZUNMQR, ZPOCON, DSCAL, ZDSCAL, ZSWAP,  ZTRSM,  XERBLA       $          ZUNMQR, ZPOCON, DSCAL, ZDSCAL, ZSWAP,  ZTRSM,  XERBLA
 *  *
       EXTERNAL  ZGESVJ        EXTERNAL  ZGESVJ
Line 640 Line 640
 *  *
 *     Quick return for void matrix (Y3K safe)  *     Quick return for void matrix (Y3K safe)
 * #:)  * #:)
       IF ( ( M .EQ. 0 ) .OR. ( N .EQ. 0 ) ) RETURN        IF ( ( M .EQ. 0 ) .OR. ( N .EQ. 0 ) ) THEN
            IWORK(1:3) = 0
            RWORK(1:7) = 0
            RETURN
         ENDIF
 *  *
 *     Determine whether the matrix U should be M x N or M x M  *     Determine whether the matrix U should be M x N or M x M
 *  *
Line 665 Line 669
 *     overflow. It is possible that this scaling pushes the smallest  *     overflow. It is possible that this scaling pushes the smallest
 *     column norm left from the underflow threshold (extreme case).  *     column norm left from the underflow threshold (extreme case).
 *  *
       SCALEM  = ONE / DSQRT(DFLOAT(M)*DFLOAT(N))        SCALEM  = ONE / DSQRT(DBLE(M)*DBLE(N))
       NOSCAL  = .TRUE.        NOSCAL  = .TRUE.
       GOSCAL  = .TRUE.        GOSCAL  = .TRUE.
       DO 1874 p = 1, N        DO 1874 p = 1, N
Line 807 Line 811
  1950       CONTINUE   1950       CONTINUE
          ELSE           ELSE
             DO 1904 p = 1, M              DO 1904 p = 1, M
                RWORK(M+N+p) = SCALEM*ABS( A(p,IDAMAX(N,A(p,1),LDA)) )                 RWORK(M+N+p) = SCALEM*ABS( A(p,IZAMAX(N,A(p,1),LDA)) )
                AATMAX = DMAX1( AATMAX, RWORK(M+N+p) )                 AATMAX = DMAX1( AATMAX, RWORK(M+N+p) )
                AATMIN = DMIN1( AATMIN, RWORK(M+N+p) )                 AATMIN = DMIN1( AATMIN, RWORK(M+N+p) )
  1904       CONTINUE   1904       CONTINUE
Line 828 Line 832
 *  *
          XSC   = ZERO           XSC   = ZERO
          TEMP1 = ONE           TEMP1 = ONE
          CALL ZLASSQ( N, SVA, 1, XSC, TEMP1 )           CALL DLASSQ( N, SVA, 1, XSC, TEMP1 )
          TEMP1 = ONE / TEMP1           TEMP1 = ONE / TEMP1
 *  *
          ENTRA = ZERO           ENTRA = ZERO
Line 836 Line 840
             BIG1  = ( ( SVA(p) / XSC )**2 ) * TEMP1              BIG1  = ( ( SVA(p) / XSC )**2 ) * TEMP1
             IF ( BIG1 .NE. ZERO ) ENTRA = ENTRA + BIG1 * DLOG(BIG1)              IF ( BIG1 .NE. ZERO ) ENTRA = ENTRA + BIG1 * DLOG(BIG1)
  1113    CONTINUE   1113    CONTINUE
          ENTRA = - ENTRA / DLOG(DFLOAT(N))           ENTRA = - ENTRA / DLOG(DBLE(N))
 *  *
 *        Now, SVA().^2/Trace(A^* * A) is a point in the probability simplex.  *        Now, SVA().^2/Trace(A^* * A) is a point in the probability simplex.
 *        It is derived from the diagonal of  A^* * A.  Do the same with the  *        It is derived from the diagonal of  A^* * A.  Do the same with the
Line 849 Line 853
             BIG1 = ( ( RWORK(p) / XSC )**2 ) * TEMP1              BIG1 = ( ( RWORK(p) / XSC )**2 ) * TEMP1
             IF ( BIG1 .NE. ZERO ) ENTRAT = ENTRAT + BIG1 * DLOG(BIG1)              IF ( BIG1 .NE. ZERO ) ENTRAT = ENTRAT + BIG1 * DLOG(BIG1)
  1114    CONTINUE   1114    CONTINUE
          ENTRAT = - ENTRAT / DLOG(DFLOAT(M))           ENTRAT = - ENTRAT / DLOG(DBLE(M))
 *  *
 *        Analyze the entropies and decide A or A^*. Smaller entropy  *        Analyze the entropies and decide A or A^*. Smaller entropy
 *        usually means better input for the algorithm.  *        usually means better input for the algorithm.
Line 905 Line 909
 *     one should use ZGESVJ instead of ZGEJSV.  *     one should use ZGESVJ instead of ZGEJSV.
 *  *
       BIG1   = DSQRT( BIG )        BIG1   = DSQRT( BIG )
       TEMP1  = DSQRT( BIG / DFLOAT(N) )        TEMP1  = DSQRT( BIG / DBLE(N) )
 *  *
       CALL ZLASCL( 'G', 0, 0, AAPP, TEMP1, N, 1, SVA, N, IERR )        CALL DLASCL( 'G', 0, 0, AAPP, TEMP1, N, 1, SVA, N, IERR )
       IF ( AAQQ .GT. (AAPP * SFMIN) ) THEN        IF ( AAQQ .GT. (AAPP * SFMIN) ) THEN
           AAQQ = ( AAQQ / AAPP ) * TEMP1            AAQQ = ( AAQQ / AAPP ) * TEMP1
       ELSE        ELSE
Line 1009 Line 1013
 *        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  *        agressive 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(DFLOAT(N))*EPSLN           TEMP1 = DSQRT(DBLE(N))*EPSLN
          DO 3001 p = 2, N           DO 3001 p = 2, N
             IF ( ABS(A(p,p)) .GE. (TEMP1*ABS(A(1,1))) ) THEN              IF ( ABS(A(p,p)) .GE. (TEMP1*ABS(A(1,1))) ) THEN
                NR = NR + 1                 NR = NR + 1
Line 1056 Line 1060
             TEMP1  = ABS(A(p,p)) / SVA(IWORK(p))              TEMP1  = ABS(A(p,p)) / SVA(IWORK(p))
             MAXPRJ = DMIN1( MAXPRJ, TEMP1 )              MAXPRJ = DMIN1( MAXPRJ, TEMP1 )
  3051    CONTINUE   3051    CONTINUE
          IF ( MAXPRJ**2 .GE. ONE - DFLOAT(N)*EPSLN ) ALMORT = .TRUE.           IF ( MAXPRJ**2 .GE. ONE - DBLE(N)*EPSLN ) ALMORT = .TRUE.
       END IF        END IF
 *  *
 *  *
Line 1136 Line 1140
 *  *
             IF ( L2PERT ) THEN              IF ( L2PERT ) THEN
 *              XSC = SQRT(SMALL)  *              XSC = SQRT(SMALL)
                XSC = EPSLN / DFLOAT(N)                 XSC = EPSLN / DBLE(N)
                DO 4947 q = 1, NR                 DO 4947 q = 1, NR
                   CTEMP = DCMPLX(XSC*ABS(A(q,q)),ZERO)                    CTEMP = DCMPLX(XSC*ABS(A(q,q)),ZERO)
                   DO 4949 p = 1, N                    DO 4949 p = 1, N
Line 1168 Line 1172
 *           to drown denormals  *           to drown denormals
             IF ( L2PERT ) THEN              IF ( L2PERT ) THEN
 *              XSC = SQRT(SMALL)  *              XSC = SQRT(SMALL)
                XSC = EPSLN / DFLOAT(N)                 XSC = EPSLN / DBLE(N)
                DO 1947 q = 1, NR                 DO 1947 q = 1, NR
                   CTEMP = DCMPLX(XSC*ABS(A(q,q)),ZERO)                    CTEMP = DCMPLX(XSC*ABS(A(q,q)),ZERO)
                   DO 1949 p = 1, NR                    DO 1949 p = 1, NR
Line 1226 Line 1230
                CALL ZCOPY( NR-p+1, V(p,p), LDV, V(p,p), 1 )                 CALL ZCOPY( NR-p+1, V(p,p), LDV, V(p,p), 1 )
                CALL ZLACGV( NR-p+1, V(p,p), 1 )                  CALL ZLACGV( NR-p+1, V(p,p), 1 ) 
  8998       CONTINUE   8998       CONTINUE
             CALL ZLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )              CALL ZLASET('Upper', NR-1, NR-1, CZERO, CZERO, V(1,2), LDV)
 *  *
             CALL ZGESVJ( 'Lower', 'U','N', NR, NR, V,LDV, SVA, NR, U,              CALL ZGESVJ( 'Lower', 'U','N', NR, NR, V,LDV, SVA, NR, U,
      $                  LDU, CWORK(N+1), LWORK-N, RWORK, LRWORK, INFO )       $                  LDU, CWORK(N+1), LWORK-N, RWORK, LRWORK, INFO )
Line 1363 Line 1367
             CONDR1 = ONE / DSQRT(TEMP1)              CONDR1 = ONE / DSQRT(TEMP1)
 *           .. here need a second oppinion on the condition number  *           .. here need a second oppinion on the condition number
 *           .. then assume worst case scenario  *           .. then assume worst case scenario
 *           R1 is OK for inverse <=> CONDR1 .LT. DFLOAT(N)  *           R1 is OK for inverse <=> CONDR1 .LT. DBLE(N)
 *           more conservative    <=> CONDR1 .LT. SQRT(DFLOAT(N))  *           more conservative    <=> CONDR1 .LT. SQRT(DBLE(N))
 *  *
             COND_OK = DSQRT(DSQRT(DFLOAT(NR)))              COND_OK = DSQRT(DSQRT(DBLE(NR)))
 *[TP]       COND_OK is a tuning parameter.  *[TP]       COND_OK is a tuning parameter.
 *  *
             IF ( CONDR1 .LT. COND_OK ) THEN              IF ( CONDR1 .LT. COND_OK ) THEN
Line 1521 Line 1525
                   CALL ZTRSM('L','U','C','N',NR,NR,CONE,CWORK(2*N+1),                    CALL ZTRSM('L','U','C','N',NR,NR,CONE,CWORK(2*N+1),
      $                 N,V,LDV)       $                 N,V,LDV)
                   IF ( NR .LT. N ) THEN                    IF ( NR .LT. N ) THEN
                    CALL ZLASET('A',N-NR,NR,ZERO,CZERO,V(NR+1,1),LDV)                    CALL ZLASET('A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV)
                    CALL ZLASET('A',NR,N-NR,ZERO,CZERO,V(1,NR+1),LDV)                    CALL ZLASET('A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV)
                    CALL ZLASET('A',N-NR,N-NR,ZERO,CONE,V(NR+1,NR+1),LDV)                    CALL ZLASET('A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV)
                   END IF                    END IF
                   CALL ZUNMQR('L','N',N,N,NR,CWORK(2*N+1),N,CWORK(N+1),                    CALL ZUNMQR('L','N',N,N,NR,CWORK(2*N+1),N,CWORK(N+1),
      $                V,LDV,CWORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR)       $                V,LDV,CWORK(2*N+N*NR+NR+1),LWORK-2*N-N*NR-NR,IERR)
Line 1605 Line 1609
 *           first QRF. Also, scale the columns to make them unit in  *           first QRF. Also, scale the columns to make them unit in
 *           Euclidean norm. This applies to all cases.  *           Euclidean norm. This applies to all cases.
 *  *
             TEMP1 = DSQRT(DFLOAT(N)) * EPSLN              TEMP1 = DSQRT(DBLE(N)) * EPSLN
             DO 1972 q = 1, N              DO 1972 q = 1, N
                DO 972 p = 1, N                 DO 972 p = 1, N
                   CWORK(2*N+N*NR+NR+IWORK(p)) = V(p,q)                    CWORK(2*N+N*NR+NR+IWORK(p)) = V(p,q)
Line 1635 Line 1639
      $           LDU, CWORK(N+1), LWORK-N, IERR )       $           LDU, CWORK(N+1), LWORK-N, IERR )
   
 *           The columns of U are normalized. The cost is O(M*N) flops.  *           The columns of U are normalized. The cost is O(M*N) flops.
             TEMP1 = DSQRT(DFLOAT(M)) * EPSLN              TEMP1 = DSQRT(DBLE(M)) * EPSLN
             DO 1973 p = 1, NR              DO 1973 p = 1, NR
                XSC = ONE / DZNRM2( M, U(1,p), 1 )                 XSC = ONE / DZNRM2( M, U(1,p), 1 )
                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )                 IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
Line 1684 Line 1688
             DO 6972 p = 1, N              DO 6972 p = 1, N
                CALL ZCOPY( N, CWORK(N+p), N, V(IWORK(p),1), LDV )                 CALL ZCOPY( N, CWORK(N+p), N, V(IWORK(p),1), LDV )
  6972       CONTINUE   6972       CONTINUE
             TEMP1 = DSQRT(DFLOAT(N))*EPSLN              TEMP1 = DSQRT(DBLE(N))*EPSLN
             DO 6971 p = 1, N              DO 6971 p = 1, N
                XSC = ONE / DZNRM2( N, V(1,p), 1 )                 XSC = ONE / DZNRM2( N, V(1,p), 1 )
                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )                 IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
Line 1702 Line 1706
             END IF              END IF
             CALL ZUNMQR( 'Left', 'No Tr', M, N1, N, A, LDA, CWORK, U,              CALL ZUNMQR( 'Left', 'No Tr', M, N1, N, A, LDA, CWORK, U,
      $           LDU, CWORK(N+1), LWORK-N, IERR )       $           LDU, CWORK(N+1), LWORK-N, IERR )
             TEMP1 = DSQRT(DFLOAT(M))*EPSLN              TEMP1 = DSQRT(DBLE(M))*EPSLN
             DO 6973 p = 1, N1              DO 6973 p = 1, N1
                XSC = ONE / DZNRM2( M, U(1,p), 1 )                 XSC = ONE / DZNRM2( M, U(1,p), 1 )
                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )                 IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
Line 1779 Line 1783
          NUMRANK = NINT(RWORK(2))           NUMRANK = NINT(RWORK(2))
   
          IF ( NR .LT. N ) THEN           IF ( NR .LT. N ) THEN
             CALL ZLASET( 'A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV )              CALL ZLASET( 'A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV )
             CALL ZLASET( 'A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV )              CALL ZLASET( 'A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV )
             CALL ZLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV )              CALL ZLASET( 'A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV )
          END IF           END IF
   
          CALL ZUNMQR( 'L','N',N,N,NR,CWORK(2*N+1),N,CWORK(N+1),           CALL ZUNMQR( 'L','N',N,N,NR,CWORK(2*N+1),N,CWORK(N+1),
Line 1791 Line 1795
 *           first QRF. Also, scale the columns to make them unit in  *           first QRF. Also, scale the columns to make them unit in
 *           Euclidean norm. This applies to all cases.  *           Euclidean norm. This applies to all cases.
 *  *
             TEMP1 = DSQRT(DFLOAT(N)) * EPSLN              TEMP1 = DSQRT(DBLE(N)) * EPSLN
             DO 7972 q = 1, N              DO 7972 q = 1, N
                DO 8972 p = 1, N                 DO 8972 p = 1, N
                   CWORK(2*N+N*NR+NR+IWORK(p)) = V(p,q)                    CWORK(2*N+N*NR+NR+IWORK(p)) = V(p,q)
Line 1836 Line 1840
 *     Undo scaling, if necessary (and possible)  *     Undo scaling, if necessary (and possible)
 *  *
       IF ( USCAL2 .LE. (BIG/SVA(1))*USCAL1 ) THEN        IF ( USCAL2 .LE. (BIG/SVA(1))*USCAL1 ) THEN
          CALL ZLASCL( 'G', 0, 0, USCAL1, USCAL2, NR, 1, SVA, N, IERR )           CALL DLASCL( 'G', 0, 0, USCAL1, USCAL2, NR, 1, SVA, N, IERR )
          USCAL1 = ONE           USCAL1 = ONE
          USCAL2 = ONE           USCAL2 = ONE
       END IF        END IF
Removed from v.1.1  
changed lines
  Added in v.1.2

CVSweb interface <joel.bertrand@systella.fr>