contribution from Zlatko Drmac

Note: I still need to work on merging [C/Z]GEJSV, but there is much more work
on these two files.  We will see when this can be done.

6 files changed, 415 insertions, 370 deletions

diff --git a/SRC/cgesvj.f b/SRC/cgesvj.f
index 7fd001a7..d94113fa 100644
--- a/SRC/cgesvj.f
+++ b/SRC/cgesvj.f
@@ -36,15 +36,15 @@
 *>
 *> \verbatim
 *>
-* CGESVJ computes the singular value decomposition (SVD) of a complex
-* M-by-N matrix A, where M >= N. The SVD of A is written as
-*                                    [++]   [xx]   [x0]   [xx]
-*              A = U * SIGMA * V^*,  [++] = [xx] * [ox] * [xx]
-*                                    [++]   [xx]
-* where SIGMA is an N-by-N diagonal matrix, U is an M-by-N orthonormal
-* matrix, and 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 right singular vectors of A, respectively.
+*> CGESVJ computes the singular value decomposition (SVD) of a complex
+*> M-by-N matrix A, where M >= N. The SVD of A is written as
+*>                                    [++]   [xx]   [x0]   [xx]
+*>              A = U * SIGMA * V^*,  [++] = [xx] * [ox] * [xx]
+*>                                    [++]   [xx]
+*> where SIGMA is an N-by-N diagonal matrix, U is an M-by-N orthonormal
+*> matrix, and 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 right singular vectors of A, respectively.
 *> \endverbatim
 *
 *  Arguments:
@@ -64,7 +64,7 @@
 *>          JOBU is CHARACTER*1
 *>          Specifies whether to compute the left singular vectors
 *>          (columns of U):
-*>          = 'U': The left singular vectors corresponding to the nonzero
+*>          = 'U' or 'F': The left singular vectors corresponding to the nonzero
 *>                 singular values are computed and returned in the leading
 *>                 columns of A. See more details in the description of A.
 *>                 The default numerical orthogonality threshold is set to
@@ -88,7 +88,7 @@
 *>          JOBV is CHARACTER*1
 *>          Specifies whether to compute the right singular vectors, that
 *>          is, the matrix V:
-*>          = 'V' : the matrix V is computed and returned in the array V
+*>          = 'V' or 'J': the matrix V is computed and returned in the array V
 *>          = 'A' : the Jacobi rotations are applied to the MV-by-N
 *>                  array V. In other words, the right singular vector
 *>                  matrix V is not computed explicitly; instead it is
@@ -206,19 +206,22 @@
 *>
 *> \param[in,out] CWORK
 *> \verbatim
-*>          CWORK is COMPLEX array, dimension M+N.
-*>          Used as work space.
+*>          CWORK is COMPLEX array, dimension max(1,LWORK).
+*>          Used as workspace.
+*>          If on entry LWORK .EQ. -1, then a workspace query is assumed and
+*>          no computation is done; CWORK(1) is set to the minial (and optimal)
+*>          length of CWORK.
 *> \endverbatim
 *>
 *> \param[in] LWORK
 *> \verbatim
-*>          LWORK is INTEGER
+*>          LWORK is INTEGER.
 *>          Length of CWORK, LWORK >= M+N.
 *> \endverbatim
 *>
 *> \param[in,out] RWORK
 *> \verbatim
-*>          RWORK is REAL array, dimension max(6,M+N).
+*>          RWORK is REAL array, dimension max(6,LRWORK).
 *>          On entry,
 *>          If JOBU .EQ. 'C' :
 *>          RWORK(1) = CTOL, where CTOL defines the threshold for convergence.
@@ -244,11 +247,14 @@
 *>          RWORK(6) = the largest absolute value over all sines of the
 *>                    Jacobi rotation angles in the last sweep. It can be
 *>                    useful for a post festum analysis.
+*>         If on entry LRWORK .EQ. -1, then a workspace query is assumed and
+*>         no computation is done; RWORK(1) is set to the minial (and optimal)
+*>         length of RWORK.
 *> \endverbatim
 *>
 *> \param[in] LRWORK
 *> \verbatim
-*>          LRWORK is INTEGER
+*>         LRWORK is INTEGER
 *>         Length of RWORK, LRWORK >= MAX(6,N).
 *> \endverbatim
 *>
@@ -261,7 +267,7 @@
 *>                (NSWEEP=30) of sweeps. The output may still be useful.
 *>                See the description of RWORK.
 *> \endverbatim
-*
+*>
 *  Authors:
 *  ========
 *
@@ -277,6 +283,8 @@
 *> \par Further Details:
 *  =====================
 *>
+*> \verbatim
+*>
 *> The orthogonal N-by-N matrix V is obtained as a product of Jacobi plane
 *> rotations. In the case of underflow of the tangent of the Jacobi angle, a
 *> modified Jacobi transformation of Drmac [3] is used. Pivot strategy uses
@@ -294,12 +302,19 @@
 *> number interval ( UNDERFLOW , OVERFLOW ). In extreme cases, even
 *> denormalized singular values can be computed with the corresponding
 *> gradual loss of accurate digits.
-*>
-*> \par Contributors:
+*> \endverbatim
+*
+*> \par Contributor:
 *  ==================
 *>
-*> Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
+*> \verbatim
+*>
+*>  ============
 *>
+*>  Zlatko Drmac (Zagreb, Croatia)
+*>
+*> \endverbatim
+*
 *> \par References:
 *  ================
 *>
@@ -319,13 +334,17 @@
 *> [6] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV,
 *>    QSVD, (H,K)-SVD computations.
 *>    Department of Mathematics, University of Zagreb, 2008, 2015.
-*>
-*> \par Bugs, Examples and Comments:
+*> \endverbatim
+*
+*> \par Bugs, examples and comments:
 *  =================================
 *>
-*> Please report all bugs and send interesting test examples and comments to
-*> drmac@math.hr. Thank you.
-*
+*> \verbatim
+*>  ===========================
+*>  Please report all bugs and send interesting test examples and comments to
+*>  drmac@math.hr. Thank you.
+*> \endverbatim
+*>
 *  =====================================================================
       SUBROUTINE CGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
      $                   LDV, CWORK, LWORK, RWORK, LRWORK, INFO )
@@ -358,24 +377,23 @@
 *     .. Local Scalars ..
       COMPLEX AAPQ, OMPQ
       REAL    AAPP, AAPP0, AAPQ1, AAQQ, APOAQ, AQOAP, BIG,
-     $        BIGTHETA, CS, CTOL, EPSLN, LARGE, MXAAPQ,
+     $        BIGTHETA, CS, CTOL, EPSLN, MXAAPQ,
      $        MXSINJ, ROOTBIG, ROOTEPS, ROOTSFMIN, ROOTTOL,
      $        SKL, SFMIN, SMALL, SN, T, TEMP1, THETA, THSIGN, TOL
       INTEGER BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1,
      $        ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, N2, N34,
      $        N4, NBL, NOTROT, p, PSKIPPED, q, ROWSKIP, SWBAND
-      LOGICAL APPLV, GOSCALE, LOWER, LSVEC, NOSCALE, ROTOK,
+      LOGICAL APPLV, GOSCALE, LOWER, LQUERY, LSVEC, NOSCALE, ROTOK,
      $        RSVEC, UCTOL, UPPER
 *     ..
 *     ..
 *     .. Intrinsic Functions ..
-      INTRINSIC ABS, AMAX1, AMIN1, CONJG, FLOAT, MIN0, MAX0,
-     $          SIGN, SQRT
+      INTRINSIC ABS, MAX, MIN, CONJG, REAL, SIGN, SQRT
 *     ..
 *     .. External Functions ..
 *     ..
 *     from BLAS
-      REAL                      SCNRM2
+      REAL               SCNRM2
       COMPLEX            CDOTC
       EXTERNAL           CDOTC, SCNRM2
       INTEGER            ISAMAX
@@ -398,13 +416,14 @@
 *
 *     Test the input arguments
 *
-      LSVEC = LSAME( JOBU, 'U' )
+      LSVEC = LSAME( JOBU, 'U' ) .OR. LSAME( JOBU, 'F' )
       UCTOL = LSAME( JOBU, 'C' )
-      RSVEC = LSAME( JOBV, 'V' )
+      RSVEC = LSAME( JOBV, 'V' ) .OR. LSAME( JOBV, 'J' )
       APPLV = LSAME( JOBV, 'A' )
       UPPER = LSAME( JOBA, 'U' )
       LOWER = LSAME( JOBA, 'L' )
 *
+      LQUERY = ( LWORK .EQ. -1 ) .OR. ( LRWORK .EQ. -1 )
       IF( .NOT.( UPPER .OR. LOWER .OR. LSAME( JOBA, 'G' ) ) ) THEN
          INFO = -1
       ELSE IF( .NOT.( LSVEC .OR. UCTOL .OR. LSAME( JOBU, 'N' ) ) ) THEN
@@ -424,9 +443,9 @@
          INFO = -11
       ELSE IF( UCTOL .AND. ( RWORK( 1 ).LE.ONE ) ) THEN
          INFO = -12
-      ELSE IF( LWORK.LT.( M+N ) ) THEN
+      ELSE IF( LWORK.LT.( M+N ) .AND. ( .NOT.LQUERY ) ) THEN
          INFO = -13
-      ELSE IF( LRWORK.LT.MAX0( N, 6 ) ) THEN
+      ELSE IF( LRWORK.LT.MAX( N, 6 ) .AND. ( .NOT.LQUERY ) ) THEN
          INFO = -15
       ELSE
          INFO = 0
@@ -436,6 +455,10 @@
       IF( INFO.NE.0 ) THEN
          CALL XERBLA( 'CGESVJ', -INFO )
          RETURN
+      ELSE IF ( LQUERY ) THEN
+         CWORK(1) = M + N
+         RWORK(1) = MAX( N, 6 )
+         RETURN
       END IF
 *
 * #:) Quick return for void matrix
@@ -455,9 +478,9 @@
       ELSE
 *        ... default
          IF( LSVEC .OR. RSVEC .OR. APPLV ) THEN
-            CTOL = SQRT( FLOAT( M ) )
+            CTOL = SQRT( REAL( M ) )
          ELSE
-            CTOL = FLOAT( M )
+            CTOL = REAL( M )
          END IF
       END IF
 *     ... and the machine dependent parameters are
@@ -468,16 +491,16 @@
       SFMIN = SLAMCH( 'SafeMinimum' )
       ROOTSFMIN = SQRT( SFMIN )
       SMALL = SFMIN / EPSLN
-      BIG = SLAMCH( 'Overflow' )
-*     BIG         = ONE    / SFMIN
+*      BIG = SLAMCH( 'Overflow' )
+      BIG     = ONE  / SFMIN
       ROOTBIG = ONE / ROOTSFMIN
-      LARGE = BIG / SQRT( FLOAT( M*N ) )
+*     LARGE = BIG / SQRT( REAL( M*N ) )
       BIGTHETA = ONE / ROOTEPS
 *
       TOL = CTOL*EPSLN
       ROOTTOL = SQRT( TOL )
 *
-      IF( FLOAT( M )*EPSLN.GE.ONE ) THEN
+      IF( REAL( M )*EPSLN.GE.ONE ) THEN
          INFO = -4
          CALL XERBLA( 'CGESVJ', -INFO )
          RETURN
@@ -502,7 +525,7 @@
 *     SQRT(N)*max_i SVA(i) does not overflow. If INFinite entries
 *     in A are detected, the procedure returns with INFO=-6.
 *
-      SKL = ONE / SQRT( FLOAT( M )*FLOAT( N ) )
+      SKL = ONE / SQRT( REAL( M )*REAL( N ) )
       NOSCALE = .TRUE.
       GOSCALE = .TRUE.
 *
@@ -592,8 +615,8 @@
       AAPP = ZERO
       AAQQ = BIG
       DO 4781 p = 1, N
-         IF( SVA( p ).NE.ZERO )AAQQ = AMIN1( AAQQ, SVA( p ) )
-         AAPP = AMAX1( AAPP, SVA( p ) )
+         IF( SVA( p ).NE.ZERO )AAQQ = MIN( AAQQ, SVA( p ) )
+         AAPP = MAX( AAPP, SVA( p ) )
  4781 CONTINUE
 *
 * #:) Quick return for zero matrix
@@ -631,22 +654,22 @@
 *     avoid underflows/overflows in computing Jacobi rotations.
 *
       SN = SQRT( SFMIN / EPSLN )
-      TEMP1 = SQRT( BIG / FLOAT( N ) )
+      TEMP1 = SQRT( BIG / REAL( N ) )
       IF( ( AAPP.LE.SN ) .OR. ( AAQQ.GE.TEMP1 ) .OR.
      $    ( ( SN.LE.AAQQ ) .AND. ( AAPP.LE.TEMP1 ) ) ) THEN
-         TEMP1 = AMIN1( BIG, TEMP1 / AAPP )
+         TEMP1 = MIN( BIG, TEMP1 / AAPP )
 *         AAQQ  = AAQQ*TEMP1
 *         AAPP  = AAPP*TEMP1
       ELSE IF( ( AAQQ.LE.SN ) .AND. ( AAPP.LE.TEMP1 ) ) THEN
-         TEMP1 = AMIN1( SN / AAQQ, BIG / ( AAPP*SQRT( FLOAT( N ) ) ) )
+         TEMP1 = MIN( SN / AAQQ, BIG / ( AAPP*SQRT( REAL( N ) ) ) )
 *         AAQQ  = AAQQ*TEMP1
 *         AAPP  = AAPP*TEMP1
       ELSE IF( ( AAQQ.GE.SN ) .AND. ( AAPP.GE.TEMP1 ) ) THEN
-         TEMP1 = AMAX1( SN / AAQQ, TEMP1 / AAPP )
+         TEMP1 = MAX( SN / AAQQ, TEMP1 / AAPP )
 *         AAQQ  = AAQQ*TEMP1
 *         AAPP  = AAPP*TEMP1
       ELSE IF( ( AAQQ.LE.SN ) .AND. ( AAPP.GE.TEMP1 ) ) THEN
-         TEMP1 = AMIN1( SN / AAQQ, BIG / ( SQRT( FLOAT( N ) )*AAPP ) )
+         TEMP1 = MIN( SN / AAQQ, BIG / ( SQRT( REAL( N ) )*AAPP ) )
 *         AAQQ  = AAQQ*TEMP1
 *         AAPP  = AAPP*TEMP1
       ELSE
@@ -683,7 +706,7 @@
 *     The boundaries are determined dynamically, based on the number of
 *     pivots above a threshold.
 *
-      KBL = MIN0( 8, N )
+      KBL = MIN( 8, N )
 *[TP] KBL is a tuning parameter that defines the tile size in the
 *     tiling of the p-q loops of pivot pairs. In general, an optimal
 *     value of KBL depends on the matrix dimensions and on the
@@ -695,7 +718,7 @@
       BLSKIP = KBL**2
 *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
 *
-      ROWSKIP = MIN0( 5, KBL )
+      ROWSKIP = MIN( 5, KBL )
 *[TP] ROWSKIP is a tuning parameter.
 *
       LKAHEAD = 1
@@ -706,7 +729,7 @@
 *     invokes cubic convergence. Big part of this cycle is done inside
 *     canonical subspaces of dimensions less than M.
 *
-      IF( ( LOWER .OR. UPPER ) .AND. ( N.GT.MAX0( 64, 4*KBL ) ) ) THEN
+      IF( ( LOWER .OR. UPPER ) .AND. ( N.GT.MAX( 64, 4*KBL ) ) ) THEN
 *[TP] The number of partition levels and the actual partition are
 *     tuning parameters.
          N4 = N / 4
@@ -804,11 +827,11 @@
 *
             igl = ( ibr-1 )*KBL + 1
 *
-            DO 1002 ir1 = 0, MIN0( LKAHEAD, NBL-ibr )
+            DO 1002 ir1 = 0, MIN( LKAHEAD, NBL-ibr )
 *
                igl = igl + ir1*KBL
 *
-               DO 2001 p = igl, MIN0( igl+KBL-1, N-1 )
+               DO 2001 p = igl, MIN( igl+KBL-1, N-1 )
 *
 *     .. de Rijk's pivoting
 *
@@ -857,7 +880,7 @@
 *
                      PSKIPPED = 0
 *
-                     DO 2002 q = p + 1, MIN0( igl+KBL-1, N )
+                     DO 2002 q = p + 1, MIN( igl+KBL-1, N )
 *
                         AAQQ = SVA( q )
 *
@@ -881,7 +904,7 @@
                               ROTOK = AAPP.LE.( AAQQ / SMALL )
                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
                                  AAPQ = ( CDOTC( M, A( 1, p ), 1,
-     $                                    A( 1, q ), 1 ) / AAQQ ) / AAPP
+     $                                    A( 1, q ), 1 ) / AAPP ) / AAQQ
                               ELSE
                                  CALL CCOPY( M, A( 1, q ), 1,
      $                                        CWORK(N+1), 1 )
@@ -895,11 +918,12 @@
 *
 *                           AAPQ = AAPQ * CONJG( CWORK(p) ) * CWORK(q)
                            AAPQ1  = -ABS(AAPQ)
-                           MXAAPQ = AMAX1( MXAAPQ, -AAPQ1 )
+                           MXAAPQ = MAX( MXAAPQ, -AAPQ1 )
 *
 *        TO rotate or NOT to rotate, THAT is the question ...
 *
                            IF( ABS( AAPQ1 ).GT.TOL ) THEN
+                               OMPQ = AAPQ / ABS(AAPQ)
 *
 *           .. rotate
 *[RTD]      ROTATED = ROTATED + ONE
@@ -912,7 +936,6 @@
 *
                               IF( ROTOK ) THEN
 *
-                                 OMPQ = AAPQ / ABS(AAPQ)
                                  AQOAP = AAQQ / AAPP
                                  APOAQ = AAPP / AAQQ
                                  THETA = -HALF*ABS( AQOAP-APOAQ )/AAPQ1
@@ -929,11 +952,11 @@
      $                                  V(1,q), 1, CS, CONJG(OMPQ)*T )
                                     END IF
 
-                                    SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
+                                    SVA( q ) = AAQQ*SQRT( MAX( ZERO,
      $                                          ONE+T*APOAQ*AAPQ1 ) )
-                                    AAPP = AAPP*SQRT( AMAX1( ZERO,
+                                    AAPP = AAPP*SQRT( MAX( ZERO,
      $                                          ONE-T*AQOAP*AAPQ1 ) )
-                                    MXSINJ = AMAX1( MXSINJ, ABS( T ) )
+                                    MXSINJ = MAX( MXSINJ, ABS( T ) )
 *
                                  ELSE
 *
@@ -945,10 +968,10 @@
                                     CS = SQRT( ONE / ( ONE+T*T ) )
                                     SN = T*CS
 *
-                                    MXSINJ = AMAX1( MXSINJ, ABS( SN ) )
-                                    SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
+                                    MXSINJ = MAX( MXSINJ, ABS( SN ) )
+                                    SVA( q ) = AAQQ*SQRT( MAX( ZERO,
      $                                          ONE+T*APOAQ*AAPQ1 ) )
-                                    AAPP = AAPP*SQRT( AMAX1( ZERO,
+                                    AAPP = AAPP*SQRT( MAX( ZERO,
      $                                      ONE-T*AQOAP*AAPQ1 ) )
 *
                                     CALL CROT( M, A(1,p), 1, A(1,q), 1,
@@ -973,9 +996,9 @@
      $                                       A( 1, q ), 1 )
                                  CALL CLASCL( 'G', 0, 0, ONE, AAQQ, M,
      $                                        1, A( 1, q ), LDA, IERR )
-                                 SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
+                                 SVA( q ) = AAQQ*SQRT( MAX( ZERO,
      $                                      ONE-AAPQ1*AAPQ1 ) )
-                                 MXSINJ = AMAX1( MXSINJ, SFMIN )
+                                 MXSINJ = MAX( MXSINJ, SFMIN )
                               END IF
 *           END IF ROTOK THEN ... ELSE
 *
@@ -1039,7 +1062,7 @@
                   ELSE
                      SVA( p ) = AAPP
                      IF( ( ir1.EQ.0 ) .AND. ( AAPP.EQ.ZERO ) )
-     $                   NOTROT = NOTROT + MIN0( igl+KBL-1, N ) - p
+     $                   NOTROT = NOTROT + MIN( igl+KBL-1, N ) - p
                   END IF
 *
  2001          CONTINUE
@@ -1059,14 +1082,14 @@
 *        doing the block at ( ibr, jbc )
 *
                IJBLSK = 0
-               DO 2100 p = igl, MIN0( igl+KBL-1, N )
+               DO 2100 p = igl, MIN( igl+KBL-1, N )
 *
                   AAPP = SVA( p )
                   IF( AAPP.GT.ZERO ) THEN
 *
                      PSKIPPED = 0
 *
-                     DO 2200 q = jgl, MIN0( jgl+KBL-1, N )
+                     DO 2200 q = jgl, MIN( jgl+KBL-1, N )
 *
                         AAQQ = SVA( q )
                         IF( AAQQ.GT.ZERO ) THEN
@@ -1102,7 +1125,8 @@
                               END IF
                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
                                  AAPQ = ( CDOTC( M, A( 1, p ), 1,
-     $                                   A( 1, q ), 1 ) / AAQQ ) / AAPP
+     $                                 A( 1, q ), 1 ) / MAX(AAQQ,AAPP) )
+     $                                                / MIN(AAQQ,AAPP)
                               ELSE
                                  CALL CCOPY( M, A( 1, q ), 1,
      $                                       CWORK(N+1), 1 )
@@ -1116,11 +1140,12 @@
 *
 *                           AAPQ = AAPQ * CONJG(CWORK(p))*CWORK(q)
                            AAPQ1  = -ABS(AAPQ)
-                           MXAAPQ = AMAX1( MXAAPQ, -AAPQ1 )
+                           MXAAPQ = MAX( MXAAPQ, -AAPQ1 )
 *
 *        TO rotate or NOT to rotate, THAT is the question ...
 *
                            IF( ABS( AAPQ1 ).GT.TOL ) THEN
+                              OMPQ = AAPQ / ABS(AAPQ)
                               NOTROT = 0
 *[RTD]      ROTATED  = ROTATED + 1
                               PSKIPPED = 0
@@ -1128,7 +1153,6 @@
 *
                               IF( ROTOK ) THEN
 *
-                                 OMPQ = AAPQ / ABS(AAPQ)
                                  AQOAP = AAQQ / AAPP
                                  APOAQ = AAPP / AAQQ
                                  THETA = -HALF*ABS( AQOAP-APOAQ )/ AAPQ1
@@ -1143,11 +1167,11 @@
                                         CALL CROT( MVL, V(1,p), 1,
      $                                  V(1,q), 1, CS, CONJG(OMPQ)*T )
                                     END IF
-                                    SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
+                                    SVA( q ) = AAQQ*SQRT( MAX( ZERO,
      $                                         ONE+T*APOAQ*AAPQ1 ) )
-                                    AAPP = AAPP*SQRT( AMAX1( ZERO,
+                                    AAPP = AAPP*SQRT( MAX( ZERO,
      $                                     ONE-T*AQOAP*AAPQ1 ) )
-                                    MXSINJ = AMAX1( MXSINJ, ABS( T ) )
+                                    MXSINJ = MAX( MXSINJ, ABS( T ) )
                                  ELSE
 *
 *                 .. choose correct signum for THETA and rotate
@@ -1158,10 +1182,10 @@
      $                                  SQRT( ONE+THETA*THETA ) )
                                     CS = SQRT( ONE / ( ONE+T*T ) )
                                     SN = T*CS
-                                    MXSINJ = AMAX1( MXSINJ, ABS( SN ) )
-                                    SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
+                                    MXSINJ = MAX( MXSINJ, ABS( SN ) )
+                                    SVA( q ) = AAQQ*SQRT( MAX( ZERO,
      $                                         ONE+T*APOAQ*AAPQ1 ) )
-                                    AAPP = AAPP*SQRT( AMAX1( ZERO,
+                                    AAPP = AAPP*SQRT( MAX( ZERO,
      $                                         ONE-T*AQOAP*AAPQ1 ) )
 *
                                     CALL CROT( M, A(1,p), 1, A(1,q), 1,
@@ -1189,9 +1213,9 @@
                                     CALL CLASCL( 'G', 0, 0, ONE, AAQQ,
      $                                           M, 1, A( 1, q ), LDA,
      $                                           IERR )
-                                    SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
+                                    SVA( q ) = AAQQ*SQRT( MAX( ZERO,
      $                                         ONE-AAPQ1*AAPQ1 ) )
-                                    MXSINJ = AMAX1( MXSINJ, SFMIN )
+                                    MXSINJ = MAX( MXSINJ, SFMIN )
                                ELSE
                                    CALL CCOPY( M, A( 1, q ), 1,
      $                                          CWORK(N+1), 1 )
@@ -1206,9 +1230,9 @@
                                     CALL CLASCL( 'G', 0, 0, ONE, AAPP,
      $                                           M, 1, A( 1, p ), LDA,
      $                                           IERR )
-                                    SVA( p ) = AAPP*SQRT( AMAX1( ZERO,
+                                    SVA( p ) = AAPP*SQRT( MAX( ZERO,
      $                                         ONE-AAPQ1*AAPQ1 ) )
-                                    MXSINJ = AMAX1( MXSINJ, SFMIN )
+                                    MXSINJ = MAX( MXSINJ, SFMIN )
                                END IF
                               END IF
 *           END IF ROTOK THEN ... ELSE
@@ -1276,7 +1300,7 @@
                   ELSE
 *
                      IF( AAPP.EQ.ZERO )NOTROT = NOTROT +
-     $                   MIN0( jgl+KBL-1, N ) - jgl + 1
+     $                   MIN( jgl+KBL-1, N ) - jgl + 1
                      IF( AAPP.LT.ZERO )NOTROT = 0
 *
                   END IF
@@ -1287,7 +1311,7 @@
 *     end of the jbc-loop
  2011       CONTINUE
 *2011 bailed out of the jbc-loop
-            DO 2012 p = igl, MIN0( igl+KBL-1, N )
+            DO 2012 p = igl, MIN( igl+KBL-1, N )
                SVA( p ) = ABS( SVA( p ) )
  2012       CONTINUE
 ***
@@ -1310,8 +1334,8 @@
          IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.
      $       ( ISWROT.LE.N ) ) )SWBAND = i
 *
-         IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.SQRT( FLOAT( N ) )*
-     $       TOL ) .AND. ( FLOAT( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
+         IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.SQRT( REAL( N ) )*
+     $       TOL ) .AND. ( REAL( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
             GO TO 1994
          END IF
 *
@@ -1359,8 +1383,9 @@
 *     Normalize the left singular vectors.
 *
       IF( LSVEC .OR. UCTOL ) THEN
-         DO 1998 p = 1, N2
-            CALL CSSCAL( M, ONE / SVA( p ), A( 1, p ), 1 )
+         DO 1998 p = 1, N4
+*           CALL CSSCAL( M, ONE / SVA( p ), A( 1, p ), 1 )
+            CALL CLASCL( 'G',0,0, SVA(p), ONE, M, 1, A(1,p), M, IERR )
  1998    CONTINUE
       END IF
 *
@@ -1388,15 +1413,15 @@
 *     then some of the singular values may overflow or underflow and
 *     the spectrum is given in this factored representation.
 *
-      RWORK( 2 ) = FLOAT( N4 )
+      RWORK( 2 ) = REAL( N4 )
 *     N4 is the number of computed nonzero singular values of A.
 *
-      RWORK( 3 ) = FLOAT( N2 )
+      RWORK( 3 ) = REAL( N2 )
 *     N2 is the number of singular values of A greater than SFMIN.
 *     If N2<N, SVA(N2:N) contains ZEROS and/or denormalized numbers
 *     that may carry some information.
 *
-      RWORK( 4 ) = FLOAT( i )
+      RWORK( 4 ) = REAL( i )
 *     i is the index of the last sweep before declaring convergence.
 *
       RWORK( 5 ) = MXAAPQ
@@ -1412,3 +1437,4 @@
 *     .. END OF CGESVJ
 *     ..
       END
+*
diff --git a/SRC/cgsvj0.f b/SRC/cgsvj0.f
index b00da9aa..631ef4ad 100644
--- a/SRC/cgsvj0.f
+++ b/SRC/cgsvj0.f
@@ -203,10 +203,10 @@
 *> CGSVJ0 is used just to enable CGESVJ to call a simplified version of
 *> itself to work on a submatrix of the original matrix.
 *>
-*> \par Contributors:
+*> \par Contributor:
 *  ==================
 *>
-*> Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
+*> Zlatko Drmac (Zagreb, Croatia)
 *>
 *> \par Bugs, Examples and Comments:
 *  =================================
@@ -255,7 +255,7 @@
 *     ..
 *     ..
 *     .. Intrinsic Functions ..
-      INTRINSIC ABS, AMAX1, CONJG, FLOAT, MIN0, SIGN, SQRT
+      INTRINSIC ABS, MAX, CONJG, REAL, MIN, SIGN, SQRT
 *     ..
 *     .. External Functions ..
       REAL               SCNRM2
@@ -338,7 +338,7 @@
 *     The boundaries are determined dynamically, based on the number of
 *     pivots above a threshold.
 *
-      KBL = MIN0( 8, N )
+      KBL = MIN( 8, N )
 *[TP] KBL is a tuning parameter that defines the tile size in the
 *     tiling of the p-q loops of pivot pairs. In general, an optimal
 *     value of KBL depends on the matrix dimensions and on the
@@ -350,7 +350,7 @@
       BLSKIP = KBL**2
 *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
 *
-      ROWSKIP = MIN0( 5, KBL )
+      ROWSKIP = MIN( 5, KBL )
 *[TP] ROWSKIP is a tuning parameter.
 *
       LKAHEAD = 1
@@ -384,11 +384,11 @@
 *
             igl = ( ibr-1 )*KBL + 1
 *
-            DO 1002 ir1 = 0, MIN0( LKAHEAD, NBL-ibr )
+            DO 1002 ir1 = 0, MIN( LKAHEAD, NBL-ibr )
 *
                igl = igl + ir1*KBL
 *
-               DO 2001 p = igl, MIN0( igl+KBL-1, N-1 )
+               DO 2001 p = igl, MIN( igl+KBL-1, N-1 )
 *
 *     .. de Rijk's pivoting
 *
@@ -437,7 +437,7 @@
 *
                      PSKIPPED = 0
 *
-                     DO 2002 q = p + 1, MIN0( igl+KBL-1, N )
+                     DO 2002 q = p + 1, MIN( igl+KBL-1, N )
 *
                         AAQQ = SVA( q )
 *
@@ -461,7 +461,7 @@
                               ROTOK = AAPP.LE.( AAQQ / SMALL )
                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
                                  AAPQ = ( CDOTC( M, A( 1, p ), 1,
-     $                                    A( 1, q ), 1 ) / AAQQ ) / AAPP
+     $                                    A( 1, q ), 1 ) / AAPP ) / AAQQ
                               ELSE
                                  CALL CCOPY( M, A( 1, q ), 1,
      $                                        WORK, 1 )
@@ -473,14 +473,14 @@
                               END IF
                            END IF
 *
-                           OMPQ = AAPQ / ABS(AAPQ)
 *                           AAPQ = AAPQ * CONJG( CWORK(p) ) * CWORK(q)
                            AAPQ1  = -ABS(AAPQ)
-                           MXAAPQ = AMAX1( MXAAPQ, -AAPQ1 )
+                           MXAAPQ = MAX( MXAAPQ, -AAPQ1 )
 *
 *        TO rotate or NOT to rotate, THAT is the question ...
 *
                            IF( ABS( AAPQ1 ).GT.TOL ) THEN
+                              OMPQ = AAPQ / ABS(AAPQ)
 *
 *           .. rotate
 *[RTD]      ROTATED = ROTATED + ONE
@@ -509,11 +509,11 @@
      $                                  V(1,q), 1, CS, CONJG(OMPQ)*T )
                                     END IF
 
-                                    SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
+                                    SVA( q ) = AAQQ*SQRT( MAX( ZERO,
      $                                          ONE+T*APOAQ*AAPQ1 ) )
-                                    AAPP = AAPP*SQRT( AMAX1( ZERO,
+                                    AAPP = AAPP*SQRT( MAX( ZERO,
      $                                          ONE-T*AQOAP*AAPQ1 ) )
-                                    MXSINJ = AMAX1( MXSINJ, ABS( T ) )
+                                    MXSINJ = MAX( MXSINJ, ABS( T ) )
 *
                                  ELSE
 *
@@ -525,10 +525,10 @@
                                     CS = SQRT( ONE / ( ONE+T*T ) )
                                     SN = T*CS
 *
-                                    MXSINJ = AMAX1( MXSINJ, ABS( SN ) )
-                                    SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
+                                    MXSINJ = MAX( MXSINJ, ABS( SN ) )
+                                    SVA( q ) = AAQQ*SQRT( MAX( ZERO,
      $                                          ONE+T*APOAQ*AAPQ1 ) )
-                                    AAPP = AAPP*SQRT( AMAX1( ZERO,
+                                    AAPP = AAPP*SQRT( MAX( ZERO,
      $                                      ONE-T*AQOAP*AAPQ1 ) )
 *
                                     CALL CROT( M, A(1,p), 1, A(1,q), 1,
@@ -553,9 +553,9 @@
      $                                       A( 1, q ), 1 )
                                  CALL CLASCL( 'G', 0, 0, ONE, AAQQ, M,
      $                                        1, A( 1, q ), LDA, IERR )
-                                 SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
+                                 SVA( q ) = AAQQ*SQRT( MAX( ZERO,
      $                                      ONE-AAPQ1*AAPQ1 ) )
-                                 MXSINJ = AMAX1( MXSINJ, SFMIN )
+                                 MXSINJ = MAX( MXSINJ, SFMIN )
                               END IF
 *           END IF ROTOK THEN ... ELSE
 *
@@ -619,7 +619,7 @@
                   ELSE
                      SVA( p ) = AAPP
                      IF( ( ir1.EQ.0 ) .AND. ( AAPP.EQ.ZERO ) )
-     $                   NOTROT = NOTROT + MIN0( igl+KBL-1, N ) - p
+     $                   NOTROT = NOTROT + MIN( igl+KBL-1, N ) - p
                   END IF
 *
  2001          CONTINUE
@@ -639,14 +639,14 @@
 *        doing the block at ( ibr, jbc )
 *
                IJBLSK = 0
-               DO 2100 p = igl, MIN0( igl+KBL-1, N )
+               DO 2100 p = igl, MIN( igl+KBL-1, N )
 *
                   AAPP = SVA( p )
                   IF( AAPP.GT.ZERO ) THEN
 *
                      PSKIPPED = 0
 *
-                     DO 2200 q = jgl, MIN0( jgl+KBL-1, N )
+                     DO 2200 q = jgl, MIN( jgl+KBL-1, N )
 *
                         AAQQ = SVA( q )
                         IF( AAQQ.GT.ZERO ) THEN
@@ -682,7 +682,8 @@
                               END IF
                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
                                  AAPQ = ( CDOTC( M, A( 1, p ), 1,
-     $                                   A( 1, q ), 1 ) / AAQQ ) / AAPP
+     $                                 A( 1, q ), 1 ) / MAX(AAQQ,AAPP) )
+     $                                               / MIN(AAQQ,AAPP)
                               ELSE
                                  CALL CCOPY( M, A( 1, q ), 1,
      $                                       WORK, 1 )
@@ -694,14 +695,14 @@
                               END IF
                            END IF
 *
-                           OMPQ = AAPQ / ABS(AAPQ)
 *                           AAPQ = AAPQ * CONJG(CWORK(p))*CWORK(q)
                            AAPQ1  = -ABS(AAPQ)
-                           MXAAPQ = AMAX1( MXAAPQ, -AAPQ1 )
+                           MXAAPQ = MAX( MXAAPQ, -AAPQ1 )
 *
 *        TO rotate or NOT to rotate, THAT is the question ...
 *
                            IF( ABS( AAPQ1 ).GT.TOL ) THEN
+                              OMPQ = AAPQ / ABS(AAPQ)
                               NOTROT = 0
 *[RTD]      ROTATED  = ROTATED + 1
                               PSKIPPED = 0
@@ -723,11 +724,11 @@
                                         CALL CROT( MVL, V(1,p), 1,
      $                                  V(1,q), 1, CS, CONJG(OMPQ)*T )
                                     END IF
-                                    SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
+                                    SVA( q ) = AAQQ*SQRT( MAX( ZERO,
      $                                         ONE+T*APOAQ*AAPQ1 ) )
-                                    AAPP = AAPP*SQRT( AMAX1( ZERO,
+                                    AAPP = AAPP*SQRT( MAX( ZERO,
      $                                     ONE-T*AQOAP*AAPQ1 ) )
-                                    MXSINJ = AMAX1( MXSINJ, ABS( T ) )
+                                    MXSINJ = MAX( MXSINJ, ABS( T ) )
                                  ELSE
 *
 *                 .. choose correct signum for THETA and rotate
@@ -738,10 +739,10 @@
      $                                  SQRT( ONE+THETA*THETA ) )
                                     CS = SQRT( ONE / ( ONE+T*T ) )
                                     SN = T*CS
-                                    MXSINJ = AMAX1( MXSINJ, ABS( SN ) )
-                                    SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
+                                    MXSINJ = MAX( MXSINJ, ABS( SN ) )
+                                    SVA( q ) = AAQQ*SQRT( MAX( ZERO,
      $                                         ONE+T*APOAQ*AAPQ1 ) )
-                                    AAPP = AAPP*SQRT( AMAX1( ZERO,
+                                    AAPP = AAPP*SQRT( MAX( ZERO,
      $                                         ONE-T*AQOAP*AAPQ1 ) )
 *
                                     CALL CROT( M, A(1,p), 1, A(1,q), 1,
@@ -769,9 +770,9 @@
                                     CALL CLASCL( 'G', 0, 0, ONE, AAQQ,
      $                                           M, 1, A( 1, q ), LDA,
      $                                           IERR )
-                                    SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
+                                    SVA( q ) = AAQQ*SQRT( MAX( ZERO,
      $                                         ONE-AAPQ1*AAPQ1 ) )
-                                    MXSINJ = AMAX1( MXSINJ, SFMIN )
+                                    MXSINJ = MAX( MXSINJ, SFMIN )
                                ELSE
                                    CALL CCOPY( M, A( 1, q ), 1,
      $                                          WORK, 1 )
@@ -786,9 +787,9 @@
                                     CALL CLASCL( 'G', 0, 0, ONE, AAPP,
      $                                           M, 1, A( 1, p ), LDA,
      $                                           IERR )
-                                    SVA( p ) = AAPP*SQRT( AMAX1( ZERO,
+                                    SVA( p ) = AAPP*SQRT( MAX( ZERO,
      $                                         ONE-AAPQ1*AAPQ1 ) )
-                                    MXSINJ = AMAX1( MXSINJ, SFMIN )
+                                    MXSINJ = MAX( MXSINJ, SFMIN )
                                END IF
                               END IF
 *           END IF ROTOK THEN ... ELSE
@@ -856,7 +857,7 @@
                   ELSE
 *
                      IF( AAPP.EQ.ZERO )NOTROT = NOTROT +
-     $                   MIN0( jgl+KBL-1, N ) - jgl + 1
+     $                   MIN( jgl+KBL-1, N ) - jgl + 1
                      IF( AAPP.LT.ZERO )NOTROT = 0
 *
                   END IF
@@ -867,7 +868,7 @@
 *     end of the jbc-loop
  2011       CONTINUE
 *2011 bailed out of the jbc-loop
-            DO 2012 p = igl, MIN0( igl+KBL-1, N )
+            DO 2012 p = igl, MIN( igl+KBL-1, N )
                SVA( p ) = ABS( SVA( p ) )
  2012       CONTINUE
 ***
@@ -890,8 +891,8 @@
          IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.
      $       ( ISWROT.LE.N ) ) )SWBAND = i
 *
-         IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.SQRT( FLOAT( N ) )*
-     $       TOL ) .AND. ( FLOAT( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
+         IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.SQRT( REAL( N ) )*
+     $       TOL ) .AND. ( REAL( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
             GO TO 1994
          END IF
 *
diff --git a/SRC/cgsvj1.f b/SRC/cgsvj1.f
index d36df3fa..e649b824 100644
--- a/SRC/cgsvj1.f
+++ b/SRC/cgsvj1.f
@@ -27,8 +27,8 @@
 *       CHARACTER*1        JOBV
 *       ..
 *       .. Array Arguments ..
-*       COMPLEX        A( LDA, * ), D( N ), V( LDV, * ), WORK( LWORK )
-*       REAL           SVA( N )
+*       COMPLEX            A( LDA, * ), D( N ), V( LDV, * ), WORK( LWORK )
+*       REAL               SVA( N )
 *       ..
 *
 *
@@ -227,10 +227,10 @@
 *
 *> \ingroup complexOTHERcomputational
 *
-*> \par Contributors:
+*> \par Contributor:
 *  ==================
 *>
-*> Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
+*> Zlatko Drmac (Zagreb, Croatia)
 *
 *  =====================================================================
       SUBROUTINE CGSVJ1( JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV,
@@ -260,7 +260,7 @@
 *     .. Local Scalars ..
       COMPLEX            AAPQ, OMPQ
       REAL               AAPP, AAPP0, AAPQ1, AAQQ, APOAQ, AQOAP, BIG,
-     $                   BIGTHETA, CS, LARGE, MXAAPQ, MXSINJ, ROOTBIG,
+     $                   BIGTHETA, CS, MXAAPQ, MXSINJ, ROOTBIG,
      $                   ROOTEPS, ROOTSFMIN, ROOTTOL, SMALL, SN, T,
      $                   TEMP1, THETA, THSIGN
       INTEGER            BLSKIP, EMPTSW, i, ibr, igl, IERR, IJBLSK,
@@ -270,7 +270,7 @@
 *     ..
 *     ..
 *     .. Intrinsic Functions ..
-      INTRINSIC          ABS, AMAX1, CONJG, FLOAT, MIN0, SIGN, SQRT
+      INTRINSIC          ABS, MAX, CONJG, REAL, MIN, SIGN, SQRT
 *     ..
 *     .. External Functions ..
       REAL               SCNRM2
@@ -334,7 +334,7 @@
       SMALL = SFMIN / EPS
       BIG = ONE / SFMIN
       ROOTBIG = ONE / ROOTSFMIN
-      LARGE = BIG / SQRT( FLOAT( M*N ) )
+*     LARGE = BIG / SQRT( REAL( M*N ) )
       BIGTHETA = ONE / ROOTEPS
       ROOTTOL = SQRT( TOL )
 *
@@ -347,7 +347,7 @@
 *
 *     .. Row-cyclic pivot strategy with de Rijk's pivoting ..
 *
-      KBL = MIN0( 8, N )
+      KBL = MIN( 8, N )
       NBLR = N1 / KBL
       IF( ( NBLR*KBL ).NE.N1 )NBLR = NBLR + 1
 
@@ -358,7 +358,7 @@
       BLSKIP = ( KBL**2 ) + 1
 *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
 
-      ROWSKIP = MIN0( 5, KBL )
+      ROWSKIP = MIN( 5, KBL )
 *[TP] ROWSKIP is a tuning parameter.
       SWBAND = 0
 *[TP] SWBAND is a tuning parameter. It is meaningful and effective
@@ -408,14 +408,14 @@
 *        doing the block at ( ibr, jbc )
 *
                IJBLSK = 0
-               DO 2100 p = igl, MIN0( igl+KBL-1, N1 )
+               DO 2100 p = igl, MIN( igl+KBL-1, N1 )
 *
                   AAPP = SVA( p )
                   IF( AAPP.GT.ZERO ) THEN
 *
                      PSKIPPED = 0
 *
-                     DO 2200 q = jgl, MIN0( jgl+KBL-1, N )
+                     DO 2200 q = jgl, MIN( jgl+KBL-1, N )
 *
                         AAQQ = SVA( q )
                         IF( AAQQ.GT.ZERO ) THEN
@@ -451,7 +451,8 @@
                               END IF
                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
                                  AAPQ = ( CDOTC( M, A( 1, p ), 1,
-     $                                   A( 1, q ), 1 ) / AAQQ ) / AAPP
+     $                                 A( 1, q ), 1 ) / MAX(AAQQ,AAPP) )
+     $                                               / MIN(AAQQ,AAPP)
                               ELSE
                                  CALL CCOPY( M, A( 1, q ), 1,
      $                                       WORK, 1 )
@@ -463,14 +464,14 @@
                               END IF
                            END IF
 *
-                           OMPQ = AAPQ / ABS(AAPQ)
 *                           AAPQ = AAPQ * CONJG(CWORK(p))*CWORK(q)
                            AAPQ1  = -ABS(AAPQ)
-                           MXAAPQ = AMAX1( MXAAPQ, -AAPQ1 )
+                           MXAAPQ = MAX( MXAAPQ, -AAPQ1 )
 *
 *        TO rotate or NOT to rotate, THAT is the question ...
 *
                            IF( ABS( AAPQ1 ).GT.TOL ) THEN
+                              OMPQ = AAPQ / ABS(AAPQ)
                               NOTROT = 0
 *[RTD]      ROTATED  = ROTATED + 1
                               PSKIPPED = 0
@@ -492,11 +493,11 @@
                                         CALL CROT( MVL, V(1,p), 1,
      $                                  V(1,q), 1, CS, CONJG(OMPQ)*T )
                                     END IF
-                                    SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
+                                    SVA( q ) = AAQQ*SQRT( MAX( ZERO,
      $                                         ONE+T*APOAQ*AAPQ1 ) )
-                                    AAPP = AAPP*SQRT( AMAX1( ZERO,
+                                    AAPP = AAPP*SQRT( MAX( ZERO,
      $                                     ONE-T*AQOAP*AAPQ1 ) )
-                                    MXSINJ = AMAX1( MXSINJ, ABS( T ) )
+                                    MXSINJ = MAX( MXSINJ, ABS( T ) )
                                  ELSE
 *
 *                 .. choose correct signum for THETA and rotate
@@ -507,10 +508,10 @@
      $                                  SQRT( ONE+THETA*THETA ) )
                                     CS = SQRT( ONE / ( ONE+T*T ) )
                                     SN = T*CS
-                                    MXSINJ = AMAX1( MXSINJ, ABS( SN ) )
-                                    SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
+                                    MXSINJ = MAX( MXSINJ, ABS( SN ) )
+                                    SVA( q ) = AAQQ*SQRT( MAX( ZERO,
      $                                         ONE+T*APOAQ*AAPQ1 ) )
-                                    AAPP = AAPP*SQRT( AMAX1( ZERO,
+                                    AAPP = AAPP*SQRT( MAX( ZERO,
      $                                         ONE-T*AQOAP*AAPQ1 ) )
 *
                                     CALL CROT( M, A(1,p), 1, A(1,q), 1,
@@ -538,9 +539,9 @@
                                     CALL CLASCL( 'G', 0, 0, ONE, AAQQ,
      $                                           M, 1, A( 1, q ), LDA,
      $                                           IERR )
-                                    SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
+                                    SVA( q ) = AAQQ*SQRT( MAX( ZERO,
      $                                         ONE-AAPQ1*AAPQ1 ) )
-                                    MXSINJ = AMAX1( MXSINJ, SFMIN )
+                                    MXSINJ = MAX( MXSINJ, SFMIN )
                                ELSE
                                    CALL CCOPY( M, A( 1, q ), 1,
      $                                          WORK, 1 )
@@ -555,9 +556,9 @@
                                     CALL CLASCL( 'G', 0, 0, ONE, AAPP,
      $                                           M, 1, A( 1, p ), LDA,
      $                                           IERR )
-                                    SVA( p ) = AAPP*SQRT( AMAX1( ZERO,
+                                    SVA( p ) = AAPP*SQRT( MAX( ZERO,
      $                                         ONE-AAPQ1*AAPQ1 ) )
-                                    MXSINJ = AMAX1( MXSINJ, SFMIN )
+                                    MXSINJ = MAX( MXSINJ, SFMIN )
                                END IF
                               END IF
 *           END IF ROTOK THEN ... ELSE
@@ -625,7 +626,7 @@
                   ELSE
 *
                      IF( AAPP.EQ.ZERO )NOTROT = NOTROT +
-     $                   MIN0( jgl+KBL-1, N ) - jgl + 1
+     $                   MIN( jgl+KBL-1, N ) - jgl + 1
                      IF( AAPP.LT.ZERO )NOTROT = 0
 *
                   END IF
@@ -636,7 +637,7 @@
 *     end of the jbc-loop
  2011       CONTINUE
 *2011 bailed out of the jbc-loop
-            DO 2012 p = igl, MIN0( igl+KBL-1, N )
+            DO 2012 p = igl, MIN( igl+KBL-1, N )
                SVA( p ) = ABS( SVA( p ) )
  2012       CONTINUE
 ***
@@ -659,8 +660,8 @@
          IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.
      $       ( ISWROT.LE.N ) ) )SWBAND = i
 *
-         IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.SQRT( FLOAT( N ) )*
-     $       TOL ) .AND. ( FLOAT( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
+         IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.SQRT( REAL( N ) )*
+     $       TOL ) .AND. ( REAL( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
             GO TO 1994
          END IF
 *
diff --git a/SRC/zgesvj.f b/SRC/zgesvj.f
index 8ac489ba..2a43f1f3 100644
--- a/SRC/zgesvj.f
+++ b/SRC/zgesvj.f
@@ -64,11 +64,11 @@
 *>          JOBU is CHARACTER*1
 *>          Specifies whether to compute the left singular vectors
 *>          (columns of U):
-*>          = 'U': The left singular vectors corresponding to the nonzero
+*>          = 'U' or 'F': The left singular vectors corresponding to the nonzero
 *>                 singular values are computed and returned in the leading
 *>                 columns of A. See more details in the description of A.
 *>                 The default numerical orthogonality threshold is set to
-*>                 approximately TOL=CTOL*EPS, CTOL=DSQRT(M), EPS=DLAMCH('E').
+*>                 approximately TOL=CTOL*EPS, CTOL=SQRT(M), EPS=DLAMCH('E').
 *>          = 'C': Analogous to JOBU='U', except that user can control the
 *>                 level of numerical orthogonality of the computed left
 *>                 singular vectors. TOL can be set to TOL = CTOL*EPS, where
@@ -88,10 +88,10 @@
 *>          JOBV is CHARACTER*1
 *>          Specifies whether to compute the right singular vectors, that
 *>          is, the matrix V:
-*>          = 'V' : the matrix V is computed and returned in the array V
+*>          = 'V' or 'J': the matrix V is computed and returned in the array V
 *>          = 'A' : the Jacobi rotations are applied to the MV-by-N
 *>                  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
 *>                  first MV rows of V.
 *>          = 'N' : the matrix V is not computed and the array V is not
@@ -206,8 +206,11 @@
 *>
 *> \param[in,out] CWORK
 *> \verbatim
-*>          CWORK is COMPLEX*16 array, dimension M+N.
-*>          Used as work space.
+*>          CWORK is COMPLEX*16 array, dimension max(1,LWORK).
+*>          Used as workspace.
+*>          If on entry LWORK .EQ. -1, then a workspace query is assumed and
+*>          no computation is done; CWORK(1) is set to the minial (and optimal)
+*>          length of CWORK.
 *> \endverbatim
 *>
 *> \param[in] LWORK
@@ -218,7 +221,7 @@
 *>
 *> \param[in,out] RWORK
 *> \verbatim
-*>          RWORK is DOUBLE PRECISION array, dimension max(6,M+N).
+*>          RWORK is DOUBLE PRECISION array, dimension max(6,LRWORK).
 *>          On entry,
 *>          If JOBU .EQ. 'C' :
 *>          RWORK(1) = CTOL, where CTOL defines the threshold for convergence.
@@ -244,11 +247,14 @@
 *>          RWORK(6) = the largest absolute value over all sines of the
 *>                    Jacobi rotation angles in the last sweep. It can be
 *>                    useful for a post festum analysis.
+*>         If on entry LRWORK .EQ. -1, then a workspace query is assumed and
+*>         no computation is done; RWORK(1) is set to the minial (and optimal)
+*>         length of RWORK.
 *> \endverbatim
 *>
 *> \param[in] LRWORK
 *> \verbatim
-*>          LRWORK is INTEGER
+*>         LRWORK is INTEGER
 *>         Length of RWORK, LRWORK >= MAX(6,N).
 *> \endverbatim
 *>
@@ -272,7 +278,7 @@
 *
 *> \date June 2016
 *
-*> \ingroup doubleGEcomputational
+*> \ingroup complex16GEcomputational
 *
 *> \par Further Details:
 *  =====================
@@ -298,14 +304,15 @@
 *> gradual loss of accurate digits.
 *> \endverbatim
 *
-*> \par Contributors:
+*> \par Contributor:
 *  ==================
 *>
 *> \verbatim
 *>
 *>  ============
 *>
-*>  Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
+*>  Zlatko Drmac (Zagreb, Croatia)
+*>
 *> \endverbatim
 *
 *> \par References:
@@ -329,8 +336,8 @@
 *>    Department of Mathematics, University of Zagreb, 2008, 2015.
 *> \endverbatim
 *
-*>  \par Bugs, examples and comments:
-*   =================================
+*> \par Bugs, examples and comments:
+*  =================================
 *>
 *> \verbatim
 *>  ===========================
@@ -369,20 +376,19 @@
 *     ..
 *     .. Local Scalars ..
       COMPLEX*16 AAPQ, OMPQ
-      DOUBLE PRECISION    AAPP, AAPP0, AAPQ1, AAQQ, APOAQ, AQOAP, BIG,
-     $        BIGTHETA, CS, CTOL, EPSLN, LARGE, MXAAPQ,
-     $        MXSINJ, ROOTBIG, ROOTEPS, ROOTSFMIN, ROOTTOL,
-     $        SKL, SFMIN, SMALL, SN, T, TEMP1, THETA, THSIGN, TOL
+      DOUBLE PRECISION AAPP, AAPP0, AAPQ1, AAQQ, APOAQ, AQOAP, BIG,
+     $       BIGTHETA, CS, CTOL, EPSLN, MXAAPQ,
+     $       MXSINJ, ROOTBIG, ROOTEPS, ROOTSFMIN, ROOTTOL,
+     $       SKL, SFMIN, SMALL, SN, T, TEMP1, THETA, THSIGN, TOL
       INTEGER BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1,
      $        ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, N2, N34,
      $        N4, NBL, NOTROT, p, PSKIPPED, q, ROWSKIP, SWBAND
-      LOGICAL APPLV, GOSCALE, LOWER, LSVEC, NOSCALE, ROTOK,
+      LOGICAL APPLV, GOSCALE, LOWER, LQUERY, LSVEC, NOSCALE, ROTOK,
      $        RSVEC, UCTOL, UPPER
 *     ..
 *     ..
 *     .. Intrinsic Functions ..
-      INTRINSIC ABS, DMAX1, DMIN1, DCONJG, DBLE, MIN0, MAX0,
-     $          DSIGN, DSQRT
+      INTRINSIC ABS, MAX, MIN, CONJG, DBLE, SIGN, SQRT
 *     ..
 *     .. External Functions ..
 *     ..
@@ -410,13 +416,14 @@
 *
 *     Test the input arguments
 *
-      LSVEC = LSAME( JOBU, 'U' )
+      LSVEC = LSAME( JOBU, 'U' ) .OR. LSAME( JOBU, 'F' )
       UCTOL = LSAME( JOBU, 'C' )
-      RSVEC = LSAME( JOBV, 'V' )
+      RSVEC = LSAME( JOBV, 'V' ) .OR. LSAME( JOBV, 'J' )
       APPLV = LSAME( JOBV, 'A' )
       UPPER = LSAME( JOBA, 'U' )
       LOWER = LSAME( JOBA, 'L' )
 *
+      LQUERY = ( LWORK .EQ. -1 ) .OR. ( LRWORK .EQ. -1 )
       IF( .NOT.( UPPER .OR. LOWER .OR. LSAME( JOBA, 'G' ) ) ) THEN
          INFO = -1
       ELSE IF( .NOT.( LSVEC .OR. UCTOL .OR. LSAME( JOBU, 'N' ) ) ) THEN
@@ -436,9 +443,9 @@
          INFO = -11
       ELSE IF( UCTOL .AND. ( RWORK( 1 ).LE.ONE ) ) THEN
          INFO = -12
-      ELSE IF( LWORK.LT.( M+N ) ) THEN
+      ELSE IF( ( LWORK.LT.( M+N ) ) .AND. ( .NOT.LQUERY ) ) THEN
          INFO = -13
-      ELSE IF( LRWORK.LT.MAX0( N, 6 ) ) THEN
+      ELSE IF( ( LRWORK.LT.MAX( N, 6 ) ) .AND. ( .NOT.LQUERY ) ) THEN
          INFO = -15
       ELSE
          INFO = 0
@@ -448,6 +455,10 @@
       IF( INFO.NE.0 ) THEN
          CALL XERBLA( 'ZGESVJ', -INFO )
          RETURN
+      ELSE IF ( LQUERY ) THEN
+         CWORK(1) = M + N
+         RWORK(1) = MAX( N, 6 )
+         RETURN
       END IF
 *
 * #:) Quick return for void matrix
@@ -467,27 +478,27 @@
       ELSE
 *        ... default
          IF( LSVEC .OR. RSVEC .OR. APPLV ) THEN
-            CTOL = DSQRT( DBLE( M ) )
+            CTOL = SQRT( DBLE( M ) )
          ELSE
             CTOL = DBLE( M )
          END IF
       END IF
 *     ... and the machine dependent parameters are
-*[!]  (Make sure that DLAMCH() works properly on the target machine.)
+*[!]  (Make sure that SLAMCH() works properly on the target machine.)
 *
       EPSLN = DLAMCH( 'Epsilon' )
-      ROOTEPS = DSQRT( EPSLN )
+      ROOTEPS = SQRT( EPSLN )
       SFMIN = DLAMCH( 'SafeMinimum' )
-      ROOTSFMIN = DSQRT( SFMIN )
+      ROOTSFMIN = SQRT( SFMIN )
       SMALL = SFMIN / EPSLN
       BIG = DLAMCH( 'Overflow' )
 *     BIG         = ONE    / SFMIN
       ROOTBIG = ONE / ROOTSFMIN
-      LARGE = BIG / DSQRT( DBLE( M*N ) )
+*      LARGE = BIG / SQRT( DBLE( M*N ) )
       BIGTHETA = ONE / ROOTEPS
 *
       TOL = CTOL*EPSLN
-      ROOTTOL = DSQRT( TOL )
+      ROOTTOL = SQRT( TOL )
 *
       IF( DBLE( M )*EPSLN.GE.ONE ) THEN
          INFO = -4
@@ -514,7 +525,7 @@
 *     SQRT(N)*max_i SVA(i) does not overflow. If INFinite entries
 *     in A are detected, the procedure returns with INFO=-6.
 *
-      SKL = ONE / DSQRT( DBLE( M )*DBLE( N ) )
+      SKL = ONE / SQRT( DBLE( M )*DBLE( N ) )
       NOSCALE = .TRUE.
       GOSCALE = .TRUE.
 *
@@ -529,7 +540,7 @@
                CALL XERBLA( 'ZGESVJ', -INFO )
                RETURN
             END IF
-            AAQQ = DSQRT( AAQQ )
+            AAQQ = SQRT( AAQQ )
             IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN
                SVA( p ) = AAPP*AAQQ
             ELSE
@@ -554,7 +565,7 @@
                CALL XERBLA( 'ZGESVJ', -INFO )
                RETURN
             END IF
-            AAQQ = DSQRT( AAQQ )
+            AAQQ = SQRT( AAQQ )
             IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN
                SVA( p ) = AAPP*AAQQ
             ELSE
@@ -579,7 +590,7 @@
                CALL XERBLA( 'ZGESVJ', -INFO )
                RETURN
             END IF
-            AAQQ = DSQRT( AAQQ )
+            AAQQ = SQRT( AAQQ )
             IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN
                SVA( p ) = AAPP*AAQQ
             ELSE
@@ -604,8 +615,8 @@
       AAPP = ZERO
       AAQQ = BIG
       DO 4781 p = 1, N
-         IF( SVA( p ).NE.ZERO )AAQQ = DMIN1( AAQQ, SVA( p ) )
-         AAPP = DMAX1( AAPP, SVA( p ) )
+         IF( SVA( p ).NE.ZERO )AAQQ = MIN( AAQQ, SVA( p ) )
+         AAPP = MAX( AAPP, SVA( p ) )
  4781 CONTINUE
 *
 * #:) Quick return for zero matrix
@@ -642,23 +653,23 @@
 *     Protect small singular values from underflow, and try to
 *     avoid underflows/overflows in computing Jacobi rotations.
 *
-      SN = DSQRT( SFMIN / EPSLN )
-      TEMP1 = DSQRT( BIG / DBLE( N ) )
+      SN = SQRT( SFMIN / EPSLN )
+      TEMP1 = SQRT( BIG / DBLE( N ) )
       IF( ( AAPP.LE.SN ) .OR. ( AAQQ.GE.TEMP1 ) .OR.
      $    ( ( SN.LE.AAQQ ) .AND. ( AAPP.LE.TEMP1 ) ) ) THEN
-         TEMP1 = DMIN1( BIG, TEMP1 / AAPP )
+         TEMP1 = MIN( BIG, TEMP1 / AAPP )
 *         AAQQ  = AAQQ*TEMP1
 *         AAPP  = AAPP*TEMP1
       ELSE IF( ( AAQQ.LE.SN ) .AND. ( AAPP.LE.TEMP1 ) ) THEN
-         TEMP1 = DMIN1( SN / AAQQ, BIG / (AAPP*DSQRT( DBLE(N)) ) )
+         TEMP1 = MIN( SN / AAQQ, BIG / (AAPP*SQRT( DBLE(N)) ) )
 *         AAQQ  = AAQQ*TEMP1
 *         AAPP  = AAPP*TEMP1
       ELSE IF( ( AAQQ.GE.SN ) .AND. ( AAPP.GE.TEMP1 ) ) THEN
-         TEMP1 = DMAX1( SN / AAQQ, TEMP1 / AAPP )
+         TEMP1 = MAX( SN / AAQQ, TEMP1 / AAPP )
 *         AAQQ  = AAQQ*TEMP1
 *         AAPP  = AAPP*TEMP1
       ELSE IF( ( AAQQ.LE.SN ) .AND. ( AAPP.GE.TEMP1 ) ) THEN
-         TEMP1 = DMIN1( SN / AAQQ, BIG / ( DSQRT( DBLE( N ) )*AAPP ) )
+         TEMP1 = MIN( SN / AAQQ, BIG / ( SQRT( DBLE( N ) )*AAPP ) )
 *         AAQQ  = AAQQ*TEMP1
 *         AAPP  = AAPP*TEMP1
       ELSE
@@ -695,7 +706,7 @@
 *     The boundaries are determined dynamically, based on the number of
 *     pivots above a threshold.
 *
-      KBL = MIN0( 8, N )
+      KBL = MIN( 8, N )
 *[TP] KBL is a tuning parameter that defines the tile size in the
 *     tiling of the p-q loops of pivot pairs. In general, an optimal
 *     value of KBL depends on the matrix dimensions and on the
@@ -707,7 +718,7 @@
       BLSKIP = KBL**2
 *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
 *
-      ROWSKIP = MIN0( 5, KBL )
+      ROWSKIP = MIN( 5, KBL )
 *[TP] ROWSKIP is a tuning parameter.
 *
       LKAHEAD = 1
@@ -718,7 +729,7 @@
 *     invokes cubic convergence. Big part of this cycle is done inside
 *     canonical subspaces of dimensions less than M.
 *
-      IF( ( LOWER .OR. UPPER ) .AND. ( N.GT.MAX0( 64, 4*KBL ) ) ) THEN
+      IF( ( LOWER .OR. UPPER ) .AND. ( N.GT.MAX( 64, 4*KBL ) ) ) THEN
 *[TP] The number of partition levels and the actual partition are
 *     tuning parameters.
          N4 = N / 4
@@ -816,11 +827,11 @@
 *
             igl = ( ibr-1 )*KBL + 1
 *
-            DO 1002 ir1 = 0, MIN0( LKAHEAD, NBL-ibr )
+            DO 1002 ir1 = 0, MIN( LKAHEAD, NBL-ibr )
 *
                igl = igl + ir1*KBL
 *
-               DO 2001 p = igl, MIN0( igl+KBL-1, N-1 )
+               DO 2001 p = igl, MIN( igl+KBL-1, N-1 )
 *
 *     .. de Rijk's pivoting
 *
@@ -858,7 +869,7 @@
                         TEMP1 = ZERO
                         AAPP = ONE
                         CALL ZLASSQ( M, A( 1, p ), 1, TEMP1, AAPP )
-                        SVA( p ) = TEMP1*DSQRT( AAPP )
+                        SVA( p ) = TEMP1*SQRT( AAPP )
                      END IF
                      AAPP = SVA( p )
                   ELSE
@@ -869,7 +880,7 @@
 *
                      PSKIPPED = 0
 *
-                     DO 2002 q = p + 1, MIN0( igl+KBL-1, N )
+                     DO 2002 q = p + 1, MIN( igl+KBL-1, N )
 *
                         AAQQ = SVA( q )
 *
@@ -893,7 +904,7 @@
                               ROTOK = AAPP.LE.( AAQQ / SMALL )
                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
                                  AAPQ = ( ZDOTC( M, A( 1, p ), 1,
-     $                                    A( 1, q ), 1 ) / AAQQ ) / AAPP
+     $                                    A( 1, q ), 1 ) / AAPP ) / AAQQ
                               ELSE
                                  CALL ZCOPY( M, A( 1, q ), 1,
      $                                        CWORK(N+1), 1 )
@@ -905,13 +916,15 @@
                               END IF
                            END IF
 *
-*                           AAPQ = AAPQ * DCONJG( CWORK(p) ) * CWORK(q)
+
+*                           AAPQ = AAPQ * CONJG( CWORK(p) ) * CWORK(q)
                            AAPQ1  = -ABS(AAPQ)
-                           MXAAPQ = DMAX1( MXAAPQ, -AAPQ1 )
+                           MXAAPQ = MAX( MXAAPQ, -AAPQ1 )
 *
 *        TO rotate or NOT to rotate, THAT is the question ...
 *
                            IF( ABS( AAPQ1 ).GT.TOL ) THEN
+                           OMPQ = AAPQ / ABS(AAPQ)
 *
 *           .. rotate
 *[RTD]      ROTATED = ROTATED + ONE
@@ -924,8 +937,7 @@
 *
                               IF( ROTOK ) THEN
 *
-                            	OMPQ = AAPQ / ABS(AAPQ)
-                                AQOAP = AAQQ / AAPP
+                                 AQOAP = AAQQ / AAPP
                                  APOAQ = AAPP / AAQQ
                                  THETA = -HALF*ABS( AQOAP-APOAQ )/AAPQ1
 *
@@ -935,39 +947,39 @@
                                     CS = ONE
 
                                     CALL ZROT( M, A(1,p), 1, A(1,q), 1,
-     $                                          CS, DCONJG(OMPQ)*T )
+     $                                          CS, CONJG(OMPQ)*T )
                                     IF ( RSVEC ) THEN
                                         CALL ZROT( MVL, V(1,p), 1,
-     $                                  V(1,q), 1, CS, DCONJG(OMPQ)*T )
+     $                                  V(1,q), 1, CS, CONJG(OMPQ)*T )
                                     END IF
 
-                                    SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
+                                    SVA( q ) = AAQQ*SQRT( MAX( ZERO,
      $                                          ONE+T*APOAQ*AAPQ1 ) )
-                                    AAPP = AAPP*DSQRT( DMAX1( ZERO,
+                                    AAPP = AAPP*SQRT( MAX( ZERO,
      $                                          ONE-T*AQOAP*AAPQ1 ) )
-                                    MXSINJ = DMAX1( MXSINJ, ABS( T ) )
+                                    MXSINJ = MAX( MXSINJ, ABS( T ) )
 *
                                  ELSE
 *
 *                 .. choose correct signum for THETA and rotate
 *
-                                    THSIGN = -DSIGN( ONE, AAPQ1 )
+                                    THSIGN = -SIGN( ONE, AAPQ1 )
                                     T = ONE / ( THETA+THSIGN*
-     $                                   DSQRT( ONE+THETA*THETA ) )
-                                    CS = DSQRT( ONE / ( ONE+T*T ) )
+     $                                   SQRT( ONE+THETA*THETA ) )
+                                    CS = SQRT( ONE / ( ONE+T*T ) )
                                     SN = T*CS
 *
-                                    MXSINJ = DMAX1( MXSINJ, ABS( SN ) )
-                                    SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
+                                    MXSINJ = MAX( MXSINJ, ABS( SN ) )
+                                    SVA( q ) = AAQQ*SQRT( MAX( ZERO,
      $                                          ONE+T*APOAQ*AAPQ1 ) )
-                                    AAPP = AAPP*DSQRT( DMAX1( ZERO,
+                                    AAPP = AAPP*SQRT( MAX( ZERO,
      $                                      ONE-T*AQOAP*AAPQ1 ) )
 *
                                     CALL ZROT( M, A(1,p), 1, A(1,q), 1,
-     $                                          CS, DCONJG(OMPQ)*SN )
+     $                                          CS, CONJG(OMPQ)*SN )
                                     IF ( RSVEC ) THEN
                                         CALL ZROT( MVL, V(1,p), 1,
-     $                                  V(1,q), 1, CS, DCONJG(OMPQ)*SN )
+     $                                  V(1,q), 1, CS, CONJG(OMPQ)*SN )
                                     END IF
                                  END IF
                                  CWORK(p) = -CWORK(q) * OMPQ
@@ -985,9 +997,9 @@
      $                                       A( 1, q ), 1 )
                                  CALL ZLASCL( 'G', 0, 0, ONE, AAQQ, M,
      $                                        1, A( 1, q ), LDA, IERR )
-                                 SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
+                                 SVA( q ) = AAQQ*SQRT( MAX( ZERO,
      $                                      ONE-AAPQ1*AAPQ1 ) )
-                                 MXSINJ = DMAX1( MXSINJ, SFMIN )
+                                 MXSINJ = MAX( MXSINJ, SFMIN )
                               END IF
 *           END IF ROTOK THEN ... ELSE
 *
@@ -1004,7 +1016,7 @@
                                     AAQQ = ONE
                                     CALL ZLASSQ( M, A( 1, q ), 1, T,
      $                                           AAQQ )
-                                    SVA( q ) = T*DSQRT( AAQQ )
+                                    SVA( q ) = T*SQRT( AAQQ )
                                  END IF
                               END IF
                               IF( ( AAPP / AAPP0 ).LE.ROOTEPS ) THEN
@@ -1016,7 +1028,7 @@
                                     AAPP = ONE
                                     CALL ZLASSQ( M, A( 1, p ), 1, T,
      $                                           AAPP )
-                                    AAPP = T*DSQRT( AAPP )
+                                    AAPP = T*SQRT( AAPP )
                                  END IF
                                  SVA( p ) = AAPP
                               END IF
@@ -1051,7 +1063,7 @@
                   ELSE
                      SVA( p ) = AAPP
                      IF( ( ir1.EQ.0 ) .AND. ( AAPP.EQ.ZERO ) )
-     $                   NOTROT = NOTROT + MIN0( igl+KBL-1, N ) - p
+     $                   NOTROT = NOTROT + MIN( igl+KBL-1, N ) - p
                   END IF
 *
  2001          CONTINUE
@@ -1071,14 +1083,14 @@
 *        doing the block at ( ibr, jbc )
 *
                IJBLSK = 0
-               DO 2100 p = igl, MIN0( igl+KBL-1, N )
+               DO 2100 p = igl, MIN( igl+KBL-1, N )
 *
                   AAPP = SVA( p )
                   IF( AAPP.GT.ZERO ) THEN
 *
                      PSKIPPED = 0
 *
-                     DO 2200 q = jgl, MIN0( jgl+KBL-1, N )
+                     DO 2200 q = jgl, MIN( jgl+KBL-1, N )
 *
                         AAQQ = SVA( q )
                         IF( AAQQ.GT.ZERO ) THEN
@@ -1114,7 +1126,8 @@
                               END IF
                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
                                  AAPQ = ( ZDOTC( M, A( 1, p ), 1,
-     $                                   A( 1, q ), 1 ) / AAQQ ) / AAPP
+     $                                 A( 1, q ), 1 ) / MAX(AAQQ,AAPP) )
+     $                                                / MIN(AAQQ,AAPP)
                               ELSE
                                  CALL ZCOPY( M, A( 1, q ), 1,
      $                                       CWORK(N+1), 1 )
@@ -1126,13 +1139,15 @@
                               END IF
                            END IF
 *
-*                           AAPQ = AAPQ * DCONJG(CWORK(p))*CWORK(q)
+
+*                           AAPQ = AAPQ * CONJG(CWORK(p))*CWORK(q)
                            AAPQ1  = -ABS(AAPQ)
-                           MXAAPQ = DMAX1( MXAAPQ, -AAPQ1 )
+                           MXAAPQ = MAX( MXAAPQ, -AAPQ1 )
 *
 *        TO rotate or NOT to rotate, THAT is the question ...
 *
                            IF( ABS( AAPQ1 ).GT.TOL ) THEN
+                              OMPQ = AAPQ / ABS(AAPQ)
                               NOTROT = 0
 *[RTD]      ROTATED  = ROTATED + 1
                               PSKIPPED = 0
@@ -1140,7 +1155,6 @@
 *
                               IF( ROTOK ) THEN
 *
-                           		 OMPQ = AAPQ / ABS(AAPQ)
                                  AQOAP = AAQQ / AAPP
                                  APOAQ = AAPP / AAQQ
                                  THETA = -HALF*ABS( AQOAP-APOAQ )/ AAPQ1
@@ -1150,37 +1164,37 @@
                                     T  = HALF / THETA
                                     CS = ONE
                                     CALL ZROT( M, A(1,p), 1, A(1,q), 1,
-     $                                          CS, DCONJG(OMPQ)*T )
+     $                                          CS, CONJG(OMPQ)*T )
                                     IF( RSVEC ) THEN
                                         CALL ZROT( MVL, V(1,p), 1,
-     $                                  V(1,q), 1, CS, DCONJG(OMPQ)*T )
+     $                                  V(1,q), 1, CS, CONJG(OMPQ)*T )
                                     END IF
-                                    SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
+                                    SVA( q ) = AAQQ*SQRT( MAX( ZERO,
      $                                         ONE+T*APOAQ*AAPQ1 ) )
-                                    AAPP = AAPP*DSQRT( DMAX1( ZERO,
+                                    AAPP = AAPP*SQRT( MAX( ZERO,
      $                                     ONE-T*AQOAP*AAPQ1 ) )
-                                    MXSINJ = DMAX1( MXSINJ, ABS( T ) )
+                                    MXSINJ = MAX( MXSINJ, ABS( T ) )
                                  ELSE
 *
 *                 .. choose correct signum for THETA and rotate
 *
-                                    THSIGN = -DSIGN( ONE, AAPQ1 )
+                                    THSIGN = -SIGN( ONE, AAPQ1 )
                                     IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN
                                     T = ONE / ( THETA+THSIGN*
-     $                                  DSQRT( ONE+THETA*THETA ) )
-                                    CS = DSQRT( ONE / ( ONE+T*T ) )
+     $                                  SQRT( ONE+THETA*THETA ) )
+                                    CS = SQRT( ONE / ( ONE+T*T ) )
                                     SN = T*CS
-                                    MXSINJ = DMAX1( MXSINJ, ABS( SN ) )
-                                    SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
+                                    MXSINJ = MAX( MXSINJ, ABS( SN ) )
+                                    SVA( q ) = AAQQ*SQRT( MAX( ZERO,
      $                                         ONE+T*APOAQ*AAPQ1 ) )
-                                    AAPP = AAPP*DSQRT( DMAX1( ZERO,
+                                    AAPP = AAPP*SQRT( MAX( ZERO,
      $                                         ONE-T*AQOAP*AAPQ1 ) )
 *
                                     CALL ZROT( M, A(1,p), 1, A(1,q), 1,
-     $                                          CS, DCONJG(OMPQ)*SN )
+     $                                          CS, CONJG(OMPQ)*SN )
                                     IF( RSVEC ) THEN
                                         CALL ZROT( MVL, V(1,p), 1,
-     $                                  V(1,q), 1, CS, DCONJG(OMPQ)*SN )
+     $                                  V(1,q), 1, CS, CONJG(OMPQ)*SN )
                                     END IF
                                  END IF
                                  CWORK(p) = -CWORK(q) * OMPQ
@@ -1201,9 +1215,9 @@
                                     CALL ZLASCL( 'G', 0, 0, ONE, AAQQ,
      $                                           M, 1, A( 1, q ), LDA,
      $                                           IERR )
-                                    SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
+                                    SVA( q ) = AAQQ*SQRT( MAX( ZERO,
      $                                         ONE-AAPQ1*AAPQ1 ) )
-                                    MXSINJ = DMAX1( MXSINJ, SFMIN )
+                                    MXSINJ = MAX( MXSINJ, SFMIN )
                                ELSE
                                    CALL ZCOPY( M, A( 1, q ), 1,
      $                                          CWORK(N+1), 1 )
@@ -1213,14 +1227,14 @@
                                     CALL ZLASCL( 'G', 0, 0, AAPP, ONE,
      $                                           M, 1, A( 1, p ), LDA,
      $                                           IERR )
-                                    CALL ZAXPY( M, -DCONJG(AAPQ),
+                                    CALL ZAXPY( M, -CONJG(AAPQ),
      $                                   CWORK(N+1), 1, A( 1, p ), 1 )
                                     CALL ZLASCL( 'G', 0, 0, ONE, AAPP,
      $                                           M, 1, A( 1, p ), LDA,
      $                                           IERR )
-                                    SVA( p ) = AAPP*DSQRT( DMAX1( ZERO,
+                                    SVA( p ) = AAPP*SQRT( MAX( ZERO,
      $                                         ONE-AAPQ1*AAPQ1 ) )
-                                    MXSINJ = DMAX1( MXSINJ, SFMIN )
+                                    MXSINJ = MAX( MXSINJ, SFMIN )
                                END IF
                               END IF
 *           END IF ROTOK THEN ... ELSE
@@ -1237,7 +1251,7 @@
                                     AAQQ = ONE
                                     CALL ZLASSQ( M, A( 1, q ), 1, T,
      $                                           AAQQ )
-                                    SVA( q ) = T*DSQRT( AAQQ )
+                                    SVA( q ) = T*SQRT( AAQQ )
                                  END IF
                               END IF
                               IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN
@@ -1249,7 +1263,7 @@
                                     AAPP = ONE
                                     CALL ZLASSQ( M, A( 1, p ), 1, T,
      $                                           AAPP )
-                                    AAPP = T*DSQRT( AAPP )
+                                    AAPP = T*SQRT( AAPP )
                                  END IF
                                  SVA( p ) = AAPP
                               END IF
@@ -1288,7 +1302,7 @@
                   ELSE
 *
                      IF( AAPP.EQ.ZERO )NOTROT = NOTROT +
-     $                   MIN0( jgl+KBL-1, N ) - jgl + 1
+     $                   MIN( jgl+KBL-1, N ) - jgl + 1
                      IF( AAPP.LT.ZERO )NOTROT = 0
 *
                   END IF
@@ -1299,7 +1313,7 @@
 *     end of the jbc-loop
  2011       CONTINUE
 *2011 bailed out of the jbc-loop
-            DO 2012 p = igl, MIN0( igl+KBL-1, N )
+            DO 2012 p = igl, MIN( igl+KBL-1, N )
                SVA( p ) = ABS( SVA( p ) )
  2012       CONTINUE
 ***
@@ -1314,7 +1328,7 @@
             T = ZERO
             AAPP = ONE
             CALL ZLASSQ( M, A( 1, N ), 1, T, AAPP )
-            SVA( N ) = T*DSQRT( AAPP )
+            SVA( N ) = T*SQRT( AAPP )
          END IF
 *
 *     Additional steering devices
@@ -1322,7 +1336,7 @@
          IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.
      $       ( ISWROT.LE.N ) ) )SWBAND = i
 *
-         IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.DSQRT( DBLE( N ) )*
+         IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.SQRT( DBLE( N ) )*
      $       TOL ) .AND. ( DBLE( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
             GO TO 1994
          END IF
@@ -1371,8 +1385,9 @@
 *     Normalize the left singular vectors.
 *
       IF( LSVEC .OR. UCTOL ) THEN
-         DO 1998 p = 1, N2
-            CALL ZDSCAL( M, ONE / SVA( p ), A( 1, p ), 1 )
+         DO 1998 p = 1, N4
+*            CALL ZDSCAL( M, ONE / SVA( p ), A( 1, p ), 1 )
+            CALL ZLASCL( 'G',0,0, SVA(p), ONE, M, 1, A(1,p), M, IERR )
  1998    CONTINUE
       END IF
 *
@@ -1390,7 +1405,7 @@
      $    .OR. ( ( SKL.LT.ONE ) .AND. ( SVA( MAX( N2, 1 ) ) .GT.
      $    ( SFMIN / SKL ) ) ) ) THEN
          DO 2400 p = 1, N
-            SVA( P ) = SKL*SVA( P )
+            SVA( p ) = SKL*SVA( p )
  2400    CONTINUE
          SKL = ONE
       END IF
diff --git a/SRC/zgsvj0.f b/SRC/zgsvj0.f
index 8eb57436..a22af86c 100644
--- a/SRC/zgsvj0.f
+++ b/SRC/zgsvj0.f
@@ -203,12 +203,12 @@
 *> ZGSVJ0 is used just to enable ZGESVJ to call a simplified version of
 *> itself to work on a submatrix of the original matrix.
 *>
-*> Contributors:
+*> Contributor:
 * =============
 *>
-*> Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
+*> Zlatko Drmac (Zagreb, Croatia)
 *>
-*> Bugs, Examples and Comments:
+*> \par Bugs, Examples and Comments:
 * ============================
 *>
 *> Please report all bugs and send interesting test examples and comments to
@@ -255,7 +255,7 @@
 *     ..
 *     ..
 *     .. Intrinsic Functions ..
-      INTRINSIC ABS, DMAX1, DCONJG, DBLE, MIN0, DSIGN, DSQRT
+      INTRINSIC ABS, MAX, CONJG, DBLE, MIN, SIGN, SQRT
 *     ..
 *     .. External Functions ..
       DOUBLE PRECISION   DZNRM2
@@ -314,13 +314,13 @@
       END IF
       RSVEC = RSVEC .OR. APPLV
 
-      ROOTEPS = DSQRT( EPS )
-      ROOTSFMIN = DSQRT( SFMIN )
+      ROOTEPS = SQRT( EPS )
+      ROOTSFMIN = SQRT( SFMIN )
       SMALL = SFMIN / EPS
       BIG = ONE / SFMIN
       ROOTBIG = ONE / ROOTSFMIN
       BIGTHETA = ONE / ROOTEPS
-      ROOTTOL = DSQRT( TOL )
+      ROOTTOL = SQRT( TOL )
 *
 *     .. Row-cyclic Jacobi SVD algorithm with column pivoting ..
 *
@@ -338,7 +338,7 @@
 *     The boundaries are determined dynamically, based on the number of
 *     pivots above a threshold.
 *
-      KBL = MIN0( 8, N )
+      KBL = MIN( 8, N )
 *[TP] KBL is a tuning parameter that defines the tile size in the
 *     tiling of the p-q loops of pivot pairs. In general, an optimal
 *     value of KBL depends on the matrix dimensions and on the
@@ -350,7 +350,7 @@
       BLSKIP = KBL**2
 *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
 *
-      ROWSKIP = MIN0( 5, KBL )
+      ROWSKIP = MIN( 5, KBL )
 *[TP] ROWSKIP is a tuning parameter.
 *
       LKAHEAD = 1
@@ -384,11 +384,11 @@
 *
             igl = ( ibr-1 )*KBL + 1
 *
-            DO 1002 ir1 = 0, MIN0( LKAHEAD, NBL-ibr )
+            DO 1002 ir1 = 0, MIN( LKAHEAD, NBL-ibr )
 *
                igl = igl + ir1*KBL
 *
-               DO 2001 p = igl, MIN0( igl+KBL-1, N-1 )
+               DO 2001 p = igl, MIN( igl+KBL-1, N-1 )
 *
 *     .. de Rijk's pivoting
 *
@@ -426,7 +426,7 @@
                         TEMP1 = ZERO
                         AAPP = ONE
                         CALL ZLASSQ( M, A( 1, p ), 1, TEMP1, AAPP )
-                        SVA( p ) = TEMP1*DSQRT( AAPP )
+                        SVA( p ) = TEMP1*SQRT( AAPP )
                      END IF
                      AAPP = SVA( p )
                   ELSE
@@ -437,7 +437,7 @@
 *
                      PSKIPPED = 0
 *
-                     DO 2002 q = p + 1, MIN0( igl+KBL-1, N )
+                     DO 2002 q = p + 1, MIN( igl+KBL-1, N )
 *
                         AAQQ = SVA( q )
 *
@@ -461,7 +461,7 @@
                               ROTOK = AAPP.LE.( AAQQ / SMALL )
                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
                                  AAPQ = ( ZDOTC( M, A( 1, p ), 1,
-     $                                    A( 1, q ), 1 ) / AAQQ ) / AAPP
+     $                                    A( 1, q ), 1 ) / AAPP ) / AAQQ
                               ELSE
                                  CALL ZCOPY( M, A( 1, q ), 1,
      $                                        WORK, 1 )
@@ -473,14 +473,14 @@
                               END IF
                            END IF
 *
-                           OMPQ = AAPQ / ABS(AAPQ)
-*                           AAPQ = AAPQ * DCONJG( CWORK(p) ) * CWORK(q)
+*                           AAPQ = AAPQ * CONJG( CWORK(p) ) * CWORK(q)
                            AAPQ1  = -ABS(AAPQ)
-                           MXAAPQ = DMAX1( MXAAPQ, -AAPQ1 )
+                           MXAAPQ = MAX( MXAAPQ, -AAPQ1 )
 *
 *        TO rotate or NOT to rotate, THAT is the question ...
 *
                            IF( ABS( AAPQ1 ).GT.TOL ) THEN
+                              OMPQ = AAPQ / ABS(AAPQ)
 *
 *           .. rotate
 *[RTD]      ROTATED = ROTATED + ONE
@@ -503,39 +503,39 @@
                                     CS = ONE
 
                                     CALL ZROT( M, A(1,p), 1, A(1,q), 1,
-     $                                          CS, DCONJG(OMPQ)*T )
+     $                                          CS, CONJG(OMPQ)*T )
                                     IF ( RSVEC ) THEN
                                         CALL ZROT( MVL, V(1,p), 1,
-     $                                  V(1,q), 1, CS, DCONJG(OMPQ)*T )
+     $                                  V(1,q), 1, CS, CONJG(OMPQ)*T )
                                     END IF
 
-                                    SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
+                                    SVA( q ) = AAQQ*SQRT( MAX( ZERO,
      $                                          ONE+T*APOAQ*AAPQ1 ) )
-                                    AAPP = AAPP*DSQRT( DMAX1( ZERO,
+                                    AAPP = AAPP*SQRT( MAX( ZERO,
      $                                          ONE-T*AQOAP*AAPQ1 ) )
-                                    MXSINJ = DMAX1( MXSINJ, ABS( T ) )
+                                    MXSINJ = MAX( MXSINJ, ABS( T ) )
 *
                                  ELSE
 *
 *                 .. choose correct signum for THETA and rotate
 *
-                                    THSIGN = -DSIGN( ONE, AAPQ1 )
+                                    THSIGN = -SIGN( ONE, AAPQ1 )
                                     T = ONE / ( THETA+THSIGN*
-     $                                   DSQRT( ONE+THETA*THETA ) )
-                                    CS = DSQRT( ONE / ( ONE+T*T ) )
+     $                                   SQRT( ONE+THETA*THETA ) )
+                                    CS = SQRT( ONE / ( ONE+T*T ) )
                                     SN = T*CS
 *
-                                    MXSINJ = DMAX1( MXSINJ, ABS( SN ) )
-                                    SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
+                                    MXSINJ = MAX( MXSINJ, ABS( SN ) )
+                                    SVA( q ) = AAQQ*SQRT( MAX( ZERO,
      $                                          ONE+T*APOAQ*AAPQ1 ) )
-                                    AAPP = AAPP*DSQRT( DMAX1( ZERO,
+                                    AAPP = AAPP*SQRT( MAX( ZERO,
      $                                      ONE-T*AQOAP*AAPQ1 ) )
 *
                                     CALL ZROT( M, A(1,p), 1, A(1,q), 1,
-     $                                          CS, DCONJG(OMPQ)*SN )
+     $                                          CS, CONJG(OMPQ)*SN )
                                     IF ( RSVEC ) THEN
                                         CALL ZROT( MVL, V(1,p), 1,
-     $                                  V(1,q), 1, CS, DCONJG(OMPQ)*SN )
+     $                                  V(1,q), 1, CS, CONJG(OMPQ)*SN )
                                     END IF
                                  END IF
                                  D(p) = -D(q) * OMPQ
@@ -553,9 +553,9 @@
      $                                       A( 1, q ), 1 )
                                  CALL ZLASCL( 'G', 0, 0, ONE, AAQQ, M,
      $                                        1, A( 1, q ), LDA, IERR )
-                                 SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
+                                 SVA( q ) = AAQQ*SQRT( MAX( ZERO,
      $                                      ONE-AAPQ1*AAPQ1 ) )
-                                 MXSINJ = DMAX1( MXSINJ, SFMIN )
+                                 MXSINJ = MAX( MXSINJ, SFMIN )
                               END IF
 *           END IF ROTOK THEN ... ELSE
 *
@@ -572,7 +572,7 @@
                                     AAQQ = ONE
                                     CALL ZLASSQ( M, A( 1, q ), 1, T,
      $                                           AAQQ )
-                                    SVA( q ) = T*DSQRT( AAQQ )
+                                    SVA( q ) = T*SQRT( AAQQ )
                                  END IF
                               END IF
                               IF( ( AAPP / AAPP0 ).LE.ROOTEPS ) THEN
@@ -584,7 +584,7 @@
                                     AAPP = ONE
                                     CALL ZLASSQ( M, A( 1, p ), 1, T,
      $                                           AAPP )
-                                    AAPP = T*DSQRT( AAPP )
+                                    AAPP = T*SQRT( AAPP )
                                  END IF
                                  SVA( p ) = AAPP
                               END IF
@@ -619,7 +619,7 @@
                   ELSE
                      SVA( p ) = AAPP
                      IF( ( ir1.EQ.0 ) .AND. ( AAPP.EQ.ZERO ) )
-     $                   NOTROT = NOTROT + MIN0( igl+KBL-1, N ) - p
+     $                   NOTROT = NOTROT + MIN( igl+KBL-1, N ) - p
                   END IF
 *
  2001          CONTINUE
@@ -639,14 +639,14 @@
 *        doing the block at ( ibr, jbc )
 *
                IJBLSK = 0
-               DO 2100 p = igl, MIN0( igl+KBL-1, N )
+               DO 2100 p = igl, MIN( igl+KBL-1, N )
 *
                   AAPP = SVA( p )
                   IF( AAPP.GT.ZERO ) THEN
 *
                      PSKIPPED = 0
 *
-                     DO 2200 q = jgl, MIN0( jgl+KBL-1, N )
+                     DO 2200 q = jgl, MIN( jgl+KBL-1, N )
 *
                         AAQQ = SVA( q )
                         IF( AAQQ.GT.ZERO ) THEN
@@ -682,7 +682,8 @@
                               END IF
                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
                                  AAPQ = ( ZDOTC( M, A( 1, p ), 1,
-     $                                   A( 1, q ), 1 ) / AAQQ ) / AAPP
+     $                                 A( 1, q ), 1 ) / MAX(AAQQ,AAPP) )
+     $                                               / MIN(AAQQ,AAPP)
                               ELSE
                                  CALL ZCOPY( M, A( 1, q ), 1,
      $                                       WORK, 1 )
@@ -694,14 +695,14 @@
                               END IF
                            END IF
 *
-                           OMPQ = AAPQ / ABS(AAPQ)
-*                           AAPQ = AAPQ * DCONJG(CWORK(p))*CWORK(q)
+*                           AAPQ = AAPQ * CONJG(CWORK(p))*CWORK(q)
                            AAPQ1  = -ABS(AAPQ)
-                           MXAAPQ = DMAX1( MXAAPQ, -AAPQ1 )
+                           MXAAPQ = MAX( MXAAPQ, -AAPQ1 )
 *
 *        TO rotate or NOT to rotate, THAT is the question ...
 *
                            IF( ABS( AAPQ1 ).GT.TOL ) THEN
+                              OMPQ = AAPQ / ABS(AAPQ)
                               NOTROT = 0
 *[RTD]      ROTATED  = ROTATED + 1
                               PSKIPPED = 0
@@ -718,37 +719,37 @@
                                     T  = HALF / THETA
                                     CS = ONE
                                     CALL ZROT( M, A(1,p), 1, A(1,q), 1,
-     $                                          CS, DCONJG(OMPQ)*T )
+     $                                          CS, CONJG(OMPQ)*T )
                                     IF( RSVEC ) THEN
                                         CALL ZROT( MVL, V(1,p), 1,
-     $                                  V(1,q), 1, CS, DCONJG(OMPQ)*T )
+     $                                  V(1,q), 1, CS, CONJG(OMPQ)*T )
                                     END IF
-                                    SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
+                                    SVA( q ) = AAQQ*SQRT( MAX( ZERO,
      $                                         ONE+T*APOAQ*AAPQ1 ) )
-                                    AAPP = AAPP*DSQRT( DMAX1( ZERO,
+                                    AAPP = AAPP*SQRT( MAX( ZERO,
      $                                     ONE-T*AQOAP*AAPQ1 ) )
-                                    MXSINJ = DMAX1( MXSINJ, ABS( T ) )
+                                    MXSINJ = MAX( MXSINJ, ABS( T ) )
                                  ELSE
 *
 *                 .. choose correct signum for THETA and rotate
 *
-                                    THSIGN = -DSIGN( ONE, AAPQ1 )
+                                    THSIGN = -SIGN( ONE, AAPQ1 )
                                     IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN
                                     T = ONE / ( THETA+THSIGN*
-     $                                  DSQRT( ONE+THETA*THETA ) )
-                                    CS = DSQRT( ONE / ( ONE+T*T ) )
+     $                                  SQRT( ONE+THETA*THETA ) )
+                                    CS = SQRT( ONE / ( ONE+T*T ) )
                                     SN = T*CS
-                                    MXSINJ = DMAX1( MXSINJ, ABS( SN ) )
-                                    SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
+                                    MXSINJ = MAX( MXSINJ, ABS( SN ) )
+                                    SVA( q ) = AAQQ*SQRT( MAX( ZERO,
      $                                         ONE+T*APOAQ*AAPQ1 ) )
-                                    AAPP = AAPP*DSQRT( DMAX1( ZERO,
+                                    AAPP = AAPP*SQRT( MAX( ZERO,
      $                                         ONE-T*AQOAP*AAPQ1 ) )
 *
                                     CALL ZROT( M, A(1,p), 1, A(1,q), 1,
-     $                                          CS, DCONJG(OMPQ)*SN )
+     $                                          CS, CONJG(OMPQ)*SN )
                                     IF( RSVEC ) THEN
                                         CALL ZROT( MVL, V(1,p), 1,
-     $                                  V(1,q), 1, CS, DCONJG(OMPQ)*SN )
+     $                                  V(1,q), 1, CS, CONJG(OMPQ)*SN )
                                     END IF
                                  END IF
                                  D(p) = -D(q) * OMPQ
@@ -769,9 +770,9 @@
                                     CALL ZLASCL( 'G', 0, 0, ONE, AAQQ,
      $                                           M, 1, A( 1, q ), LDA,
      $                                           IERR )
-                                    SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
+                                    SVA( q ) = AAQQ*SQRT( MAX( ZERO,
      $                                         ONE-AAPQ1*AAPQ1 ) )
-                                    MXSINJ = DMAX1( MXSINJ, SFMIN )
+                                    MXSINJ = MAX( MXSINJ, SFMIN )
                                ELSE
                                    CALL ZCOPY( M, A( 1, q ), 1,
      $                                          WORK, 1 )
@@ -781,14 +782,14 @@
                                     CALL ZLASCL( 'G', 0, 0, AAPP, ONE,
      $                                           M, 1, A( 1, p ), LDA,
      $                                           IERR )
-                                    CALL ZAXPY( M, -DCONJG(AAPQ),
+                                    CALL ZAXPY( M, -CONJG(AAPQ),
      $                                   WORK, 1, A( 1, p ), 1 )
                                     CALL ZLASCL( 'G', 0, 0, ONE, AAPP,
      $                                           M, 1, A( 1, p ), LDA,
      $                                           IERR )
-                                    SVA( p ) = AAPP*DSQRT( DMAX1( ZERO,
+                                    SVA( p ) = AAPP*SQRT( MAX( ZERO,
      $                                         ONE-AAPQ1*AAPQ1 ) )
-                                    MXSINJ = DMAX1( MXSINJ, SFMIN )
+                                    MXSINJ = MAX( MXSINJ, SFMIN )
                                END IF
                               END IF
 *           END IF ROTOK THEN ... ELSE
@@ -805,7 +806,7 @@
                                     AAQQ = ONE
                                     CALL ZLASSQ( M, A( 1, q ), 1, T,
      $                                           AAQQ )
-                                    SVA( q ) = T*DSQRT( AAQQ )
+                                    SVA( q ) = T*SQRT( AAQQ )
                                  END IF
                               END IF
                               IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN
@@ -817,7 +818,7 @@
                                     AAPP = ONE
                                     CALL ZLASSQ( M, A( 1, p ), 1, T,
      $                                           AAPP )
-                                    AAPP = T*DSQRT( AAPP )
+                                    AAPP = T*SQRT( AAPP )
                                  END IF
                                  SVA( p ) = AAPP
                               END IF
@@ -856,7 +857,7 @@
                   ELSE
 *
                      IF( AAPP.EQ.ZERO )NOTROT = NOTROT +
-     $                   MIN0( jgl+KBL-1, N ) - jgl + 1
+     $                   MIN( jgl+KBL-1, N ) - jgl + 1
                      IF( AAPP.LT.ZERO )NOTROT = 0
 *
                   END IF
@@ -867,7 +868,7 @@
 *     end of the jbc-loop
  2011       CONTINUE
 *2011 bailed out of the jbc-loop
-            DO 2012 p = igl, MIN0( igl+KBL-1, N )
+            DO 2012 p = igl, MIN( igl+KBL-1, N )
                SVA( p ) = ABS( SVA( p ) )
  2012       CONTINUE
 ***
@@ -882,7 +883,7 @@
             T = ZERO
             AAPP = ONE
             CALL ZLASSQ( M, A( 1, N ), 1, T, AAPP )
-            SVA( N ) = T*DSQRT( AAPP )
+            SVA( N ) = T*SQRT( AAPP )
          END IF
 *
 *     Additional steering devices
@@ -890,7 +891,7 @@
          IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.
      $       ( ISWROT.LE.N ) ) )SWBAND = i
 *
-         IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.DSQRT( DBLE( N ) )*
+         IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.SQRT( DBLE( N ) )*
      $       TOL ) .AND. ( DBLE( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
             GO TO 1994
          END IF
@@ -910,7 +911,7 @@
 *
       INFO = 0
 * #:) INFO = 0 confirms successful iterations.
- 1995  CONTINUE
+ 1995 CONTINUE
 *
 *     Sort the vector SVA() of column norms.
       DO 5991 p = 1, N - 1
diff --git a/SRC/zgsvj1.f b/SRC/zgsvj1.f
index 89ce3d01..9c764c89 100644
--- a/SRC/zgsvj1.f
+++ b/SRC/zgsvj1.f
@@ -27,8 +27,8 @@
 *       CHARACTER*1        JOBV
 *       ..
 *       .. Array Arguments ..
-*       COMPLEX*16     A( LDA, * ), D( N ), V( LDV, * ), WORK( LWORK )
-*       DOUBLE PRECISION           SVA( N )
+*       COMPLEX*16         A( LDA, * ), D( N ), V( LDV, * ), WORK( LWORK )
+*       DOUBLE PRECISION   SVA( N )
 *       ..
 *
 *
@@ -227,10 +227,10 @@
 *
 *> \ingroup complex16OTHERcomputational
 *
-*> \par Contributors:
+*> \par Contributor:
 *  ==================
 *>
-*> Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
+*> Zlatko Drmac (Zagreb, Croatia)
 *
 *  =====================================================================
       SUBROUTINE ZGSVJ1( JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV,
@@ -261,7 +261,7 @@
 *     .. Local Scalars ..
       COMPLEX*16         AAPQ, OMPQ
       DOUBLE PRECISION   AAPP, AAPP0, AAPQ1, AAQQ, APOAQ, AQOAP, BIG,
-     $                   BIGTHETA, CS, LARGE, MXAAPQ, MXSINJ, ROOTBIG,
+     $                   BIGTHETA, CS, MXAAPQ, MXSINJ, ROOTBIG,
      $                   ROOTEPS, ROOTSFMIN, ROOTTOL, SMALL, SN, T,
      $                   TEMP1, THETA, THSIGN
       INTEGER            BLSKIP, EMPTSW, i, ibr, igl, IERR, IJBLSK,
@@ -271,7 +271,7 @@
 *     ..
 *     ..
 *     .. Intrinsic Functions ..
-      INTRINSIC          ABS, DCONJG, DMAX1, DBLE, MIN0, DSIGN, DSQRT
+      INTRINSIC          ABS, CONJG, MAX, DBLE, MIN, SIGN, SQRT
 *     ..
 *     .. External Functions ..
       DOUBLE PRECISION   DZNRM2
@@ -330,14 +330,14 @@
       END IF
       RSVEC = RSVEC .OR. APPLV
 
-      ROOTEPS = DSQRT( EPS )
-      ROOTSFMIN = DSQRT( SFMIN )
+      ROOTEPS = SQRT( EPS )
+      ROOTSFMIN = SQRT( SFMIN )
       SMALL = SFMIN / EPS
       BIG = ONE / SFMIN
       ROOTBIG = ONE / ROOTSFMIN
-      LARGE = BIG / DSQRT( DBLE( M*N ) )
+*     LARGE = BIG / SQRT( DBLE( M*N ) )
       BIGTHETA = ONE / ROOTEPS
-      ROOTTOL = DSQRT( TOL )
+      ROOTTOL = SQRT( TOL )
 *
 *     .. Initialize the right singular vector matrix ..
 *
@@ -348,7 +348,7 @@
 *
 *     .. Row-cyclic pivot strategy with de Rijk's pivoting ..
 *
-      KBL = MIN0( 8, N )
+      KBL = MIN( 8, N )
       NBLR = N1 / KBL
       IF( ( NBLR*KBL ).NE.N1 )NBLR = NBLR + 1
 
@@ -359,7 +359,7 @@
       BLSKIP = ( KBL**2 ) + 1
 *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
 
-      ROWSKIP = MIN0( 5, KBL )
+      ROWSKIP = MIN( 5, KBL )
 *[TP] ROWSKIP is a tuning parameter.
       SWBAND = 0
 *[TP] SWBAND is a tuning parameter. It is meaningful and effective
@@ -409,14 +409,14 @@
 *        doing the block at ( ibr, jbc )
 *
                IJBLSK = 0
-               DO 2100 p = igl, MIN0( igl+KBL-1, N1 )
+               DO 2100 p = igl, MIN( igl+KBL-1, N1 )
 *
                   AAPP = SVA( p )
                   IF( AAPP.GT.ZERO ) THEN
 *
                      PSKIPPED = 0
 *
-                     DO 2200 q = jgl, MIN0( jgl+KBL-1, N )
+                     DO 2200 q = jgl, MIN( jgl+KBL-1, N )
 *
                         AAQQ = SVA( q )
                         IF( AAQQ.GT.ZERO ) THEN
@@ -452,7 +452,8 @@
                               END IF
                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
                                  AAPQ = ( ZDOTC( M, A( 1, p ), 1,
-     $                                   A( 1, q ), 1 ) / AAQQ ) / AAPP
+     $                                 A( 1, q ), 1 ) / MAX(AAQQ,AAPP) )
+     $                                               / MIN(AAQQ,AAPP)
                               ELSE
                                  CALL ZCOPY( M, A( 1, q ), 1,
      $                                       WORK, 1 )
@@ -464,14 +465,14 @@
                               END IF
                            END IF
 *
-                           OMPQ = AAPQ / ABS(AAPQ)
-*                           AAPQ = AAPQ * DCONJG(CWORK(p))*CWORK(q)
+*                           AAPQ = AAPQ * CONJG(CWORK(p))*CWORK(q)
                            AAPQ1  = -ABS(AAPQ)
-                           MXAAPQ = DMAX1( MXAAPQ, -AAPQ1 )
+                           MXAAPQ = MAX( MXAAPQ, -AAPQ1 )
 *
 *        TO rotate or NOT to rotate, THAT is the question ...
 *
                            IF( ABS( AAPQ1 ).GT.TOL ) THEN
+                              OMPQ = AAPQ / ABS(AAPQ)
                               NOTROT = 0
 *[RTD]      ROTATED  = ROTATED + 1
                               PSKIPPED = 0
@@ -488,37 +489,37 @@
                                     T  = HALF / THETA
                                     CS = ONE
                                     CALL ZROT( M, A(1,p), 1, A(1,q), 1,
-     $                                          CS, DCONJG(OMPQ)*T )
+     $                                          CS, CONJG(OMPQ)*T )
                                     IF( RSVEC ) THEN
                                         CALL ZROT( MVL, V(1,p), 1,
-     $                                  V(1,q), 1, CS, DCONJG(OMPQ)*T )
+     $                                  V(1,q), 1, CS, CONJG(OMPQ)*T )
                                     END IF
-                                    SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
+                                    SVA( q ) = AAQQ*SQRT( MAX( ZERO,
      $                                         ONE+T*APOAQ*AAPQ1 ) )
-                                    AAPP = AAPP*DSQRT( DMAX1( ZERO,
+                                    AAPP = AAPP*SQRT( MAX( ZERO,
      $                                     ONE-T*AQOAP*AAPQ1 ) )
-                                    MXSINJ = DMAX1( MXSINJ, ABS( T ) )
+                                    MXSINJ = MAX( MXSINJ, ABS( T ) )
                                  ELSE
 *
 *                 .. choose correct signum for THETA and rotate
 *
-                                    THSIGN = -DSIGN( ONE, AAPQ1 )
+                                    THSIGN = -SIGN( ONE, AAPQ1 )
                                     IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN
                                     T = ONE / ( THETA+THSIGN*
-     $                                  DSQRT( ONE+THETA*THETA ) )
-                                    CS = DSQRT( ONE / ( ONE+T*T ) )
+     $                                  SQRT( ONE+THETA*THETA ) )
+                                    CS = SQRT( ONE / ( ONE+T*T ) )
                                     SN = T*CS
-                                    MXSINJ = DMAX1( MXSINJ, ABS( SN ) )
-                                    SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
+                                    MXSINJ = MAX( MXSINJ, ABS( SN ) )
+                                    SVA( q ) = AAQQ*SQRT( MAX( ZERO,
      $                                         ONE+T*APOAQ*AAPQ1 ) )
-                                    AAPP = AAPP*DSQRT( DMAX1( ZERO,
+                                    AAPP = AAPP*SQRT( MAX( ZERO,
      $                                         ONE-T*AQOAP*AAPQ1 ) )
 *
                                     CALL ZROT( M, A(1,p), 1, A(1,q), 1,
-     $                                          CS, DCONJG(OMPQ)*SN )
+     $                                          CS, CONJG(OMPQ)*SN )
                                     IF( RSVEC ) THEN
                                         CALL ZROT( MVL, V(1,p), 1,
-     $                                  V(1,q), 1, CS, DCONJG(OMPQ)*SN )
+     $                                  V(1,q), 1, CS, CONJG(OMPQ)*SN )
                                     END IF
                                  END IF
                                  D(p) = -D(q) * OMPQ
@@ -539,9 +540,9 @@
                                     CALL ZLASCL( 'G', 0, 0, ONE, AAQQ,
      $                                           M, 1, A( 1, q ), LDA,
      $                                           IERR )
-                                    SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
+                                    SVA( q ) = AAQQ*SQRT( MAX( ZERO,
      $                                         ONE-AAPQ1*AAPQ1 ) )
-                                    MXSINJ = DMAX1( MXSINJ, SFMIN )
+                                    MXSINJ = MAX( MXSINJ, SFMIN )
                                ELSE
                                    CALL ZCOPY( M, A( 1, q ), 1,
      $                                          WORK, 1 )
@@ -551,14 +552,14 @@
                                     CALL ZLASCL( 'G', 0, 0, AAPP, ONE,
      $                                           M, 1, A( 1, p ), LDA,
      $                                           IERR )
-                                    CALL ZAXPY( M, -DCONJG(AAPQ),
+                                    CALL ZAXPY( M, -CONJG(AAPQ),
      $                                   WORK, 1, A( 1, p ), 1 )
                                     CALL ZLASCL( 'G', 0, 0, ONE, AAPP,
      $                                           M, 1, A( 1, p ), LDA,
      $                                           IERR )
-                                    SVA( p ) = AAPP*DSQRT( DMAX1( ZERO,
+                                    SVA( p ) = AAPP*SQRT( MAX( ZERO,
      $                                         ONE-AAPQ1*AAPQ1 ) )
-                                    MXSINJ = DMAX1( MXSINJ, SFMIN )
+                                    MXSINJ = MAX( MXSINJ, SFMIN )
                                END IF
                               END IF
 *           END IF ROTOK THEN ... ELSE
@@ -575,7 +576,7 @@
                                     AAQQ = ONE
                                     CALL ZLASSQ( M, A( 1, q ), 1, T,
      $                                           AAQQ )
-                                    SVA( q ) = T*DSQRT( AAQQ )
+                                    SVA( q ) = T*SQRT( AAQQ )
                                  END IF
                               END IF
                               IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN
@@ -587,7 +588,7 @@
                                     AAPP = ONE
                                     CALL ZLASSQ( M, A( 1, p ), 1, T,
      $                                           AAPP )
-                                    AAPP = T*DSQRT( AAPP )
+                                    AAPP = T*SQRT( AAPP )
                                  END IF
                                  SVA( p ) = AAPP
                               END IF
@@ -626,7 +627,7 @@
                   ELSE
 *
                      IF( AAPP.EQ.ZERO )NOTROT = NOTROT +
-     $                   MIN0( jgl+KBL-1, N ) - jgl + 1
+     $                   MIN( jgl+KBL-1, N ) - jgl + 1
                      IF( AAPP.LT.ZERO )NOTROT = 0
 *
                   END IF
@@ -637,7 +638,7 @@
 *     end of the jbc-loop
  2011       CONTINUE
 *2011 bailed out of the jbc-loop
-            DO 2012 p = igl, MIN0( igl+KBL-1, N )
+            DO 2012 p = igl, MIN( igl+KBL-1, N )
                SVA( p ) = ABS( SVA( p ) )
  2012       CONTINUE
 ***
@@ -652,7 +653,7 @@
             T = ZERO
             AAPP = ONE
             CALL ZLASSQ( M, A( 1, N ), 1, T, AAPP )
-            SVA( N ) = T*DSQRT( AAPP )
+            SVA( N ) = T*SQRT( AAPP )
          END IF
 *
 *     Additional steering devices
@@ -660,7 +661,7 @@
          IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.
      $       ( ISWROT.LE.N ) ) )SWBAND = i
 *
-         IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.DSQRT( DBLE( N ) )*
+         IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.SQRT( DBLE( N ) )*
      $       TOL ) .AND. ( DBLE( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
             GO TO 1994
          END IF