Thanks to Lawrence Mulholland from NAG for improving the newly released complex

Jacobi SVD, and thanks to Zlatko Drmac for following so fast and improving the
code.  See: http://icl.cs.utk.edu/lapack-forum/viewtopic.php?t=4844
for more information about all this. Julien.

2 files changed, 46 insertions, 46 deletions

diff --git a/SRC/cgejsv.f b/SRC/cgejsv.f
index 4b3e9056..c003b5f5 100644
--- a/SRC/cgejsv.f
+++ b/SRC/cgejsv.f
@@ -556,13 +556,13 @@
 *     ..
 *     .. External Functions ..
       REAL      SLAMCH, SCNRM2
-      INTEGER   ISAMAX
+      INTEGER   ISAMAX, ICAMAX
       LOGICAL   LSAME
-      EXTERNAL  ISAMAX, LSAME, SLAMCH, SCNRM2
+      EXTERNAL  ISAMAX, ICAMAX, LSAME, SLAMCH, SCNRM2
 *     ..
 *     .. External Subroutines ..
       EXTERNAL  CCOPY,  CGELQF, CGEQP3, CGEQRF, CLACPY, CLASCL,
-     $          CLASET, CLASSQ, CLASWP, CUNGQR, CUNMLQ,
+     $          SLASCL, CLASET, CLASSQ, SLASSQ, CLASWP, CUNGQR, CUNMLQ,
      $          CUNMQR, CPOCON, SSCAL, CSSCAL, CSWAP,  CTRSM,  XERBLA
 *
       EXTERNAL  CGESVJ
@@ -803,7 +803,7 @@
  1950       CONTINUE
          ELSE
             DO 1904 p = 1, M
-               RWORK(M+N+p) = SCALEM*ABS( A(p,ISAMAX(N,A(p,1),LDA)) )
+               RWORK(M+N+p) = SCALEM*ABS( A(p,ICAMAX(N,A(p,1),LDA)) )
                AATMAX = AMAX1( AATMAX, RWORK(M+N+p) )
                AATMIN = AMIN1( AATMIN, RWORK(M+N+p) )
  1904       CONTINUE
@@ -824,7 +824,7 @@
 *
          XSC   = ZERO
          TEMP1 = ONE
-         CALL CLASSQ( N, SVA, 1, XSC, TEMP1 )
+         CALL SLASSQ( N, SVA, 1, XSC, TEMP1 )
          TEMP1 = ONE / TEMP1
 *
          ENTRA = ZERO
@@ -903,7 +903,7 @@
       BIG1   = SQRT( BIG )
       TEMP1  = SQRT( BIG / FLOAT(N) )
 *
-      CALL CLASCL( 'G', 0, 0, AAPP, TEMP1, N, 1, SVA, N, IERR )
+      CALL SLASCL( 'G', 0, 0, AAPP, TEMP1, N, 1, SVA, N, IERR )
       IF ( AAQQ .GT. (AAPP * SFMIN) ) THEN
           AAQQ = ( AAQQ / AAPP ) * TEMP1
       ELSE
@@ -1221,7 +1221,7 @@
                CALL CCOPY( NR-p+1, V(p,p), LDV, V(p,p), 1 )
                CALL CLACGV( NR-p+1, V(p,p), 1 ) 
  8998       CONTINUE
-            CALL CLASET( 'Upper', NR-1, NR-1, ZERO, ZERO, V(1,2), LDV )
+            CALL CLASET('Upper', NR-1, NR-1, CZERO, CZERO, V(1,2), LDV)
 *
             CALL CGESVJ( 'Lower', 'U','N', NR, NR, V,LDV, SVA, NR, U,
      $                  LDU, CWORK(N+1), LWORK-N, RWORK, LRWORK, INFO )
@@ -1517,9 +1517,9 @@
                   CALL CTRSM('L','U','C','N',NR,NR,CONE,CWORK(2*N+1),
      $                 N,V,LDV)
                   IF ( NR .LT. N ) THEN
-                   CALL CLASET('A',N-NR,NR,ZERO,CZERO,V(NR+1,1),LDV)
-                   CALL CLASET('A',NR,N-NR,ZERO,CZERO,V(1,NR+1),LDV)
-                   CALL CLASET('A',N-NR,N-NR,ZERO,CONE,V(NR+1,NR+1),LDV)
+                  CALL CLASET('A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV)
+                  CALL CLASET('A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV)
+                  CALL CLASET('A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV)
                   END IF
                   CALL CUNMQR('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)
@@ -1775,9 +1775,9 @@
          NUMRANK = NINT(RWORK(2))
 
          IF ( NR .LT. N ) THEN
-            CALL CLASET( 'A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV )
-            CALL CLASET( 'A',NR,N-NR,ZERO,ZERO,V(1,NR+1),LDV )
-            CALL CLASET( 'A',N-NR,N-NR,ZERO,ONE,V(NR+1,NR+1),LDV )
+            CALL CLASET( 'A',N-NR,NR,CZERO,CZERO,V(NR+1,1),LDV )
+            CALL CLASET( 'A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV )
+            CALL CLASET( 'A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV )
          END IF
 
          CALL CUNMQR( 'L','N',N,N,NR,CWORK(2*N+1),N,CWORK(N+1),
@@ -1832,7 +1832,7 @@
 *     Undo scaling, if necessary (and possible)
 *
       IF ( USCAL2 .LE. (BIG/SVA(1))*USCAL1 ) THEN
-         CALL CLASCL( 'G', 0, 0, USCAL1, USCAL2, NR, 1, SVA, N, IERR )
+         CALL SLASCL( 'G', 0, 0, USCAL1, USCAL2, NR, 1, SVA, N, IERR )
          USCAL1 = ONE
          USCAL2 = ONE
       END IF
diff --git a/SRC/zgejsv.f b/SRC/zgejsv.f
index 62274f3a..644e67a3 100644
--- a/SRC/zgejsv.f
+++ b/SRC/zgejsv.f
@@ -554,18 +554,18 @@
      $        NOSCAL, ROWPIV, RSVEC,  TRANSP
 *     ..
 *     .. Intrinsic Functions ..
-      INTRINSIC ABS,  DCMPLX, DCONJG, DLOG, DMAX1, DMIN1, DFLOAT,
+      INTRINSIC ABS,  DCMPLX, DCONJG, DLOG, DMAX1, DMIN1, DBLE,
      $          MAX0, MIN0, NINT,  DSQRT
 *     ..
 *     .. External Functions ..
       DOUBLE PRECISION      DLAMCH, DZNRM2
-      INTEGER   IDAMAX
+      INTEGER   IDAMAX, IZAMAX
       LOGICAL   LSAME
-      EXTERNAL  IDAMAX, LSAME, DLAMCH, DZNRM2
+      EXTERNAL  IDAMAX, IZAMAX, LSAME, DLAMCH, DZNRM2
 *     ..
 *     .. External Subroutines ..
-      EXTERNAL  ZCOPY,  ZGELQF, ZGEQP3, ZGEQRF, ZLACPY, ZLASCL,
-     $          ZLASET, ZLASSQ, ZLASWP, ZUNGQR, ZUNMLQ,
+      EXTERNAL  DLASSQ, ZCOPY,  ZGELQF, ZGEQP3, ZGEQRF, ZLACPY, ZLASCL,
+     $          DLASCL, ZLASET, ZLASSQ, ZLASWP, ZUNGQR, ZUNMLQ,
      $          ZUNMQR, ZPOCON, DSCAL, ZDSCAL, ZSWAP,  ZTRSM,  XERBLA
 *
       EXTERNAL  ZGESVJ
@@ -665,7 +665,7 @@
 *     overflow. It is possible that this scaling pushes the smallest
 *     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.
       GOSCAL  = .TRUE.
       DO 1874 p = 1, N
@@ -807,7 +807,7 @@
  1950       CONTINUE
          ELSE
             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) )
                AATMIN = DMIN1( AATMIN, RWORK(M+N+p) )
  1904       CONTINUE
@@ -828,7 +828,7 @@
 *
          XSC   = ZERO
          TEMP1 = ONE
-         CALL ZLASSQ( N, SVA, 1, XSC, TEMP1 )
+         CALL DLASSQ( N, SVA, 1, XSC, TEMP1 )
          TEMP1 = ONE / TEMP1
 *
          ENTRA = ZERO
@@ -836,7 +836,7 @@
             BIG1  = ( ( SVA(p) / XSC )**2 ) * TEMP1
             IF ( BIG1 .NE. ZERO ) ENTRA = ENTRA + BIG1 * DLOG(BIG1)
  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.
 *        It is derived from the diagonal of  A^* * A.  Do the same with the
@@ -849,7 +849,7 @@
             BIG1 = ( ( RWORK(p) / XSC )**2 ) * TEMP1
             IF ( BIG1 .NE. ZERO ) ENTRAT = ENTRAT + BIG1 * DLOG(BIG1)
  1114    CONTINUE
-         ENTRAT = - ENTRAT / DLOG(DFLOAT(M))
+         ENTRAT = - ENTRAT / DLOG(DBLE(M))
 *
 *        Analyze the entropies and decide A or A^*. Smaller entropy
 *        usually means better input for the algorithm.
@@ -905,9 +905,9 @@
 *     one should use ZGESVJ instead of ZGEJSV.
 *
       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
           AAQQ = ( AAQQ / AAPP ) * TEMP1
       ELSE
@@ -1009,7 +1009,7 @@
 *        sigma_i < N*EPSLN*||A|| are flushed to zero. This is an
 *        agressive enforcement of lower numerical rank by introducing 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
             IF ( ABS(A(p,p)) .GE. (TEMP1*ABS(A(1,1))) ) THEN
                NR = NR + 1
@@ -1056,7 +1056,7 @@
             TEMP1  = ABS(A(p,p)) / SVA(IWORK(p))
             MAXPRJ = DMIN1( MAXPRJ, TEMP1 )
  3051    CONTINUE
-         IF ( MAXPRJ**2 .GE. ONE - DFLOAT(N)*EPSLN ) ALMORT = .TRUE.
+         IF ( MAXPRJ**2 .GE. ONE - DBLE(N)*EPSLN ) ALMORT = .TRUE.
       END IF
 *
 *
@@ -1136,7 +1136,7 @@
 *
             IF ( L2PERT ) THEN
 *              XSC = SQRT(SMALL)
-               XSC = EPSLN / DFLOAT(N)
+               XSC = EPSLN / DBLE(N)
                DO 4947 q = 1, NR
                   CTEMP = DCMPLX(XSC*ABS(A(q,q)),ZERO)
                   DO 4949 p = 1, N
@@ -1168,7 +1168,7 @@
 *           to drown denormals
             IF ( L2PERT ) THEN
 *              XSC = SQRT(SMALL)
-               XSC = EPSLN / DFLOAT(N)
+               XSC = EPSLN / DBLE(N)
                DO 1947 q = 1, NR
                   CTEMP = DCMPLX(XSC*ABS(A(q,q)),ZERO)
                   DO 1949 p = 1, NR
@@ -1226,7 +1226,7 @@
                CALL ZCOPY( NR-p+1, V(p,p), LDV, V(p,p), 1 )
                CALL ZLACGV( NR-p+1, V(p,p), 1 ) 
  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,
      $                  LDU, CWORK(N+1), LWORK-N, RWORK, LRWORK, INFO )
@@ -1363,10 +1363,10 @@
             CONDR1 = ONE / DSQRT(TEMP1)
 *           .. here need a second oppinion on the condition number
 *           .. then assume worst case scenario
-*           R1 is OK for inverse <=> CONDR1 .LT. DFLOAT(N)
-*           more conservative    <=> CONDR1 .LT. SQRT(DFLOAT(N))
+*           R1 is OK for inverse <=> CONDR1 .LT. DBLE(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.
 *
             IF ( CONDR1 .LT. COND_OK ) THEN
@@ -1521,9 +1521,9 @@
                   CALL ZTRSM('L','U','C','N',NR,NR,CONE,CWORK(2*N+1),
      $                 N,V,LDV)
                   IF ( NR .LT. N ) THEN
-                   CALL ZLASET('A',N-NR,NR,ZERO,CZERO,V(NR+1,1),LDV)
-                   CALL ZLASET('A',NR,N-NR,ZERO,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,NR,CZERO,CZERO,V(NR+1,1),LDV)
+                  CALL ZLASET('A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV)
+                  CALL ZLASET('A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV)
                   END IF
                   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)
@@ -1605,7 +1605,7 @@
 *           first QRF. Also, scale the columns to make them unit in
 *           Euclidean norm. This applies to all cases.
 *
-            TEMP1 = DSQRT(DFLOAT(N)) * EPSLN
+            TEMP1 = DSQRT(DBLE(N)) * EPSLN
             DO 1972 q = 1, N
                DO 972 p = 1, N
                   CWORK(2*N+N*NR+NR+IWORK(p)) = V(p,q)
@@ -1635,7 +1635,7 @@
      $           LDU, CWORK(N+1), LWORK-N, IERR )
 
 *           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
                XSC = ONE / DZNRM2( M, U(1,p), 1 )
                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
@@ -1684,7 +1684,7 @@
             DO 6972 p = 1, N
                CALL ZCOPY( N, CWORK(N+p), N, V(IWORK(p),1), LDV )
  6972       CONTINUE
-            TEMP1 = DSQRT(DFLOAT(N))*EPSLN
+            TEMP1 = DSQRT(DBLE(N))*EPSLN
             DO 6971 p = 1, N
                XSC = ONE / DZNRM2( N, V(1,p), 1 )
                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
@@ -1702,7 +1702,7 @@
             END IF
             CALL ZUNMQR( 'Left', 'No Tr', M, N1, N, A, LDA, CWORK, U,
      $           LDU, CWORK(N+1), LWORK-N, IERR )
-            TEMP1 = DSQRT(DFLOAT(M))*EPSLN
+            TEMP1 = DSQRT(DBLE(M))*EPSLN
             DO 6973 p = 1, N1
                XSC = ONE / DZNRM2( M, U(1,p), 1 )
                IF ( (XSC .LT. (ONE-TEMP1)) .OR. (XSC .GT. (ONE+TEMP1)) )
@@ -1779,9 +1779,9 @@
          NUMRANK = NINT(RWORK(2))
 
          IF ( NR .LT. N ) THEN
-            CALL ZLASET( 'A',N-NR,NR,ZERO,ZERO,V(NR+1,1),LDV )
-            CALL ZLASET( 'A',NR,N-NR,ZERO,ZERO,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,NR,CZERO,CZERO,V(NR+1,1),LDV )
+            CALL ZLASET( 'A',NR,N-NR,CZERO,CZERO,V(1,NR+1),LDV )
+            CALL ZLASET( 'A',N-NR,N-NR,CZERO,CONE,V(NR+1,NR+1),LDV )
          END IF
 
          CALL ZUNMQR( 'L','N',N,N,NR,CWORK(2*N+1),N,CWORK(N+1),
@@ -1791,7 +1791,7 @@
 *           first QRF. Also, scale the columns to make them unit in
 *           Euclidean norm. This applies to all cases.
 *
-            TEMP1 = DSQRT(DFLOAT(N)) * EPSLN
+            TEMP1 = DSQRT(DBLE(N)) * EPSLN
             DO 7972 q = 1, N
                DO 8972 p = 1, N
                   CWORK(2*N+N*NR+NR+IWORK(p)) = V(p,q)
@@ -1836,7 +1836,7 @@
 *     Undo scaling, if necessary (and possible)
 *
       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
          USCAL2 = ONE
       END IF