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

Diff for /rpl/lapack/lapack/dgesvj.f between versions 1.1 and 1.11

version 1.1, 2010/08/07 13:21:03 version 1.11, 2012/12/14 12:30:20
Line 1 Line 1
       SUBROUTINE DGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,  *> \brief \b DGESVJ
      +                   LDV, WORK, LWORK, INFO )  *
   *  =========== DOCUMENTATION ===========
 *  *
 *  -- LAPACK routine (version 3.2.2)                                    --  * Online html documentation available at 
   *            http://www.netlib.org/lapack/explore-html/ 
 *  *
 *  -- Contributed by Zlatko Drmac of the University of Zagreb and     --  *> \htmlonly
 *  -- Kresimir Veselic of the Fernuniversitaet Hagen                  --  *> Download DGESVJ + dependencies 
 *  -- June 2010                                                       --  *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesvj.f"> 
   *> [TGZ]</a> 
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesvj.f"> 
   *> [ZIP]</a> 
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesvj.f"> 
   *> [TXT]</a>
   *> \endhtmlonly 
   *
   *  Definition:
   *  ===========
   *
   *       SUBROUTINE DGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
   *                          LDV, WORK, LWORK, INFO )
   * 
   *       .. Scalar Arguments ..
   *       INTEGER            INFO, LDA, LDV, LWORK, M, MV, N
   *       CHARACTER*1        JOBA, JOBU, JOBV
   *       ..
   *       .. Array Arguments ..
   *       DOUBLE PRECISION   A( LDA, * ), SVA( N ), V( LDV, * ),
   *      $                   WORK( LWORK )
   *       ..
   *  
   *
   *> \par Purpose:
   *  =============
   *>
   *> \verbatim
   *>
   *> DGESVJ computes the singular value decomposition (SVD) of a real
   *> M-by-N matrix A, where M >= N. The SVD of A is written as
   *>                                    [++]   [xx]   [x0]   [xx]
   *>              A = U * SIGMA * V^t,  [++] = [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 orthogonal 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:
   *  ==========
   *
   *> \param[in] JOBA
   *> \verbatim
   *>          JOBA is CHARACTER* 1
   *>          Specifies the structure of A.
   *>          = 'L': The input matrix A is lower triangular;
   *>          = 'U': The input matrix A is upper triangular;
   *>          = 'G': The input matrix A is general M-by-N matrix, M >= N.
   *> \endverbatim
   *>
   *> \param[in] JOBU
   *> \verbatim
   *>          JOBU is CHARACTER*1
   *>          Specifies whether to compute the left singular vectors
   *>          (columns of U):
   *>          = 'U': 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').
   *>          = '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
   *>                 CTOL is given on input in the array WORK.
   *>                 No CTOL smaller than ONE is allowed. CTOL greater
   *>                 than 1 / EPS is meaningless. The option 'C'
   *>                 can be used if M*EPS is satisfactory orthogonality
   *>                 of the computed left singular vectors, so CTOL=M could
   *>                 save few sweeps of Jacobi rotations.
   *>                 See the descriptions of A and WORK(1).
   *>          = 'N': The matrix U is not computed. However, see the
   *>                 description of A.
   *> \endverbatim
   *>
   *> \param[in] JOBV
   *> \verbatim
   *>          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
   *>          = '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
   *>                  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
   *>                  referenced
   *> \endverbatim
   *>
   *> \param[in] M
   *> \verbatim
   *>          M is INTEGER
   *>          The number of rows of the input matrix A. 1/DLAMCH('E') > M >= 0.  
   *> \endverbatim
   *>
   *> \param[in] N
   *> \verbatim
   *>          N is INTEGER
   *>          The number of columns of the input matrix A.
   *>          M >= N >= 0.
   *> \endverbatim
   *>
   *> \param[in,out] A
   *> \verbatim
   *>          A is DOUBLE PRECISION array, dimension (LDA,N)
   *>          On entry, the M-by-N matrix A.
   *>          On exit :
   *>          If JOBU .EQ. 'U' .OR. JOBU .EQ. 'C' :
   *>                 If INFO .EQ. 0 :
   *>                 RANKA orthonormal columns of U are returned in the
   *>                 leading RANKA columns of the array A. Here RANKA <= N
   *>                 is the number of computed singular values of A that are
   *>                 above the underflow threshold DLAMCH('S'). The singular
   *>                 vectors corresponding to underflowed or zero singular
   *>                 values are not computed. The value of RANKA is returned
   *>                 in the array WORK as RANKA=NINT(WORK(2)). Also see the
   *>                 descriptions of SVA and WORK. The computed columns of U
   *>                 are mutually numerically orthogonal up to approximately
   *>                 TOL=DSQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU.EQ.'C'),
   *>                 see the description of JOBU.
   *>                 If INFO .GT. 0 :
   *>                 the procedure DGESVJ did not converge in the given number
   *>                 of iterations (sweeps). In that case, the computed
   *>                 columns of U may not be orthogonal up to TOL. The output
   *>                 U (stored in A), SIGMA (given by the computed singular
   *>                 values in SVA(1:N)) and V is still a decomposition of the
   *>                 input matrix A in the sense that the residual
   *>                 ||A-SCALE*U*SIGMA*V^T||_2 / ||A||_2 is small.
   *>
   *>          If JOBU .EQ. 'N' :
   *>                 If INFO .EQ. 0 :
   *>                 Note that the left singular vectors are 'for free' in the
   *>                 one-sided Jacobi SVD algorithm. However, if only the
   *>                 singular values are needed, the level of numerical
   *>                 orthogonality of U is not an issue and iterations are
   *>                 stopped when the columns of the iterated matrix are
   *>                 numerically orthogonal up to approximately M*EPS. Thus,
   *>                 on exit, A contains the columns of U scaled with the
   *>                 corresponding singular values.
   *>                 If INFO .GT. 0 :
   *>                 the procedure DGESVJ did not converge in the given number
   *>                 of iterations (sweeps).
   *> \endverbatim
   *>
   *> \param[in] LDA
   *> \verbatim
   *>          LDA is INTEGER
   *>          The leading dimension of the array A.  LDA >= max(1,M).
   *> \endverbatim
   *>
   *> \param[out] SVA
   *> \verbatim
   *>          SVA is DOUBLE PRECISION array, dimension (N)
   *>          On exit :
   *>          If INFO .EQ. 0 :
   *>          depending on the value SCALE = WORK(1), we have:
   *>                 If SCALE .EQ. ONE :
   *>                 SVA(1:N) contains the computed singular values of A.
   *>                 During the computation SVA contains the Euclidean column
   *>                 norms of the iterated matrices in the array A.
   *>                 If SCALE .NE. ONE :
   *>                 The singular values of A are SCALE*SVA(1:N), and this
   *>                 factored representation is due to the fact that some of the
   *>                 singular values of A might underflow or overflow.
   *>          If INFO .GT. 0 :
   *>          the procedure DGESVJ did not converge in the given number of
   *>          iterations (sweeps) and SCALE*SVA(1:N) may not be accurate.
   *> \endverbatim
   *>
   *> \param[in] MV
   *> \verbatim
   *>          MV is INTEGER
   *>          If JOBV .EQ. 'A', then the product of Jacobi rotations in DGESVJ
   *>          is applied to the first MV rows of V. See the description of JOBV.
   *> \endverbatim
   *>
   *> \param[in,out] V
   *> \verbatim
   *>          V is DOUBLE PRECISION array, dimension (LDV,N)
   *>          If JOBV = 'V', then V contains on exit the N-by-N matrix of
   *>                         the right singular vectors;
   *>          If JOBV = 'A', then V contains the product of the computed right
   *>                         singular vector matrix and the initial matrix in
   *>                         the array V.
   *>          If JOBV = 'N', then V is not referenced.
   *> \endverbatim
   *>
   *> \param[in] LDV
   *> \verbatim
   *>          LDV is INTEGER
   *>          The leading dimension of the array V, LDV .GE. 1.
   *>          If JOBV .EQ. 'V', then LDV .GE. max(1,N).
   *>          If JOBV .EQ. 'A', then LDV .GE. max(1,MV) .
   *> \endverbatim
   *>
   *> \param[in,out] WORK
   *> \verbatim
   *>          WORK is DOUBLE PRECISION array, dimension max(4,M+N).
   *>          On entry :
   *>          If JOBU .EQ. 'C' :
   *>          WORK(1) = CTOL, where CTOL defines the threshold for convergence.
   *>                    The process stops if all columns of A are mutually
   *>                    orthogonal up to CTOL*EPS, EPS=DLAMCH('E').
   *>                    It is required that CTOL >= ONE, i.e. it is not
   *>                    allowed to force the routine to obtain orthogonality
   *>                    below EPS.
   *>          On exit :
   *>          WORK(1) = SCALE is the scaling factor such that SCALE*SVA(1:N)
   *>                    are the computed singular values of A.
   *>                    (See description of SVA().)
   *>          WORK(2) = NINT(WORK(2)) is the number of the computed nonzero
   *>                    singular values.
   *>          WORK(3) = NINT(WORK(3)) is the number of the computed singular
   *>                    values that are larger than the underflow threshold.
   *>          WORK(4) = NINT(WORK(4)) is the number of sweeps of Jacobi
   *>                    rotations needed for numerical convergence.
   *>          WORK(5) = max_{i.NE.j} |COS(A(:,i),A(:,j))| in the last sweep.
   *>                    This is useful information in cases when DGESVJ did
   *>                    not converge, as it can be used to estimate whether
   *>                    the output is stil useful and for post festum analysis.
   *>          WORK(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.
   *> \endverbatim
   *>
   *> \param[in] LWORK
   *> \verbatim
   *>          LWORK is INTEGER
   *>          length of WORK, WORK >= MAX(6,M+N)
   *> \endverbatim
   *>
   *> \param[out] INFO
   *> \verbatim
   *>          INFO is INTEGER
   *>          = 0 : successful exit.
   *>          < 0 : if INFO = -i, then the i-th argument had an illegal value
   *>          > 0 : DGESVJ did not converge in the maximal allowed number (30)
   *>                of sweeps. The output may still be useful. See the
   *>                description of WORK.
   *> \endverbatim
   *
   *  Authors:
   *  ========
   *
   *> \author Univ. of Tennessee 
   *> \author Univ. of California Berkeley 
   *> \author Univ. of Colorado Denver 
   *> \author NAG Ltd. 
   *
   *> \date September 2012
   *
   *> \ingroup doubleGEcomputational
   *
   *> \par Further Details:
   *  =====================
   *>
   *> \verbatim
   *>
   *>  The orthogonal N-by-N matrix V is obtained as a product of Jacobi plane
   *>  rotations. The rotations are implemented as fast scaled rotations of
   *>  Anda and Park [1]. In the case of underflow of the Jacobi angle, a
   *>  modified Jacobi transformation of Drmac [4] is used. Pivot strategy uses
   *>  column interchanges of de Rijk [2]. The relative accuracy of the computed
   *>  singular values and the accuracy of the computed singular vectors (in
   *>  angle metric) is as guaranteed by the theory of Demmel and Veselic [3].
   *>  The condition number that determines the accuracy in the full rank case
   *>  is essentially min_{D=diag} kappa(A*D), where kappa(.) is the
   *>  spectral condition number. The best performance of this Jacobi SVD
   *>  procedure is achieved if used in an  accelerated version of Drmac and
   *>  Veselic [5,6], and it is the kernel routine in the SIGMA library [7].
   *>  Some tunning parameters (marked with [TP]) are available for the
   *>  implementer.
   *>  The computational range for the nonzero singular values is the  machine
   *>  number interval ( UNDERFLOW , OVERFLOW ). In extreme cases, even
   *>  denormalized singular values can be computed with the corresponding
   *>  gradual loss of accurate digits.
   *> \endverbatim
   *
   *> \par Contributors:
   *  ==================
   *>
   *> \verbatim
   *>
   *>  ============
   *>
   *>  Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
   *> \endverbatim
   *
   *> \par References:
   *  ================
   *>
   *> \verbatim
   *>
   *> [1] A. A. Anda and H. Park: Fast plane rotations with dynamic scaling.
   *>     SIAM J. matrix Anal. Appl., Vol. 15 (1994), pp. 162-174.
   *> [2] P. P. M. De Rijk: A one-sided Jacobi algorithm for computing the
   *>     singular value decomposition on a vector computer.
   *>     SIAM J. Sci. Stat. Comp., Vol. 10 (1998), pp. 359-371.
   *> [3] J. Demmel and K. Veselic: Jacobi method is more accurate than QR.
   *> [4] Z. Drmac: Implementation of Jacobi rotations for accurate singular
   *>     value computation in floating point arithmetic.
   *>     SIAM J. Sci. Comp., Vol. 18 (1997), pp. 1200-1222.
   *> [5] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I.
   *>     SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342.
   *>     LAPACK Working note 169.
   *> [6] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II.
   *>     SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362.
   *>     LAPACK Working note 170.
   *> [7] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV,
   *>     QSVD, (H,K)-SVD computations.
   *>     Department of Mathematics, University of Zagreb, 2008.
   *> \endverbatim
   *
   *>  \par Bugs, examples and comments:
   *   =================================
   *>
   *> \verbatim
   *>  ===========================
   *>  Please report all bugs and send interesting test examples and comments to
   *>  drmac@math.hr. Thank you.
   *> \endverbatim
   *>
   *  =====================================================================
         SUBROUTINE DGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
        $                   LDV, WORK, LWORK, INFO )
 *  *
   *  -- LAPACK computational routine (version 3.4.2) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --  *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--  *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
   *     September 2012
 *  *
 * This routine is also part of SIGMA (version 1.23, October 23. 2008.)  
 * SIGMA is a library of algorithms for highly accurate algorithms for  
 * computation of SVD, PSVD, QSVD, (H,K)-SVD, and for solution of the  
 * eigenvalue problems Hx = lambda M x, H M x = lambda x with H, M > 0.  
 *  
       IMPLICIT           NONE  
 *     ..  
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       INTEGER            INFO, LDA, LDV, LWORK, M, MV, N        INTEGER            INFO, LDA, LDV, LWORK, M, MV, N
       CHARACTER*1        JOBA, JOBU, JOBV        CHARACTER*1        JOBA, JOBU, JOBV
 *     ..  *     ..
 *     .. Array Arguments ..  *     .. Array Arguments ..
       DOUBLE PRECISION   A( LDA, * ), SVA( N ), V( LDV, * ),        DOUBLE PRECISION   A( LDA, * ), SVA( N ), V( LDV, * ),
      +                   WORK( LWORK )       $                   WORK( LWORK )
 *     ..  *     ..
 *  *
 *  Purpose  
 *  =======  
 *  
 *  DGESVJ computes the singular value decomposition (SVD) of a real  
 *  M-by-N matrix A, where M >= N. The SVD of A is written as  
 *                                     [++]   [xx]   [x0]   [xx]  
 *               A = U * SIGMA * V^t,  [++] = [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 orthogonal 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.  
 *  
 *  Further Details  
 *  ~~~~~~~~~~~~~~~  
 *  The orthogonal N-by-N matrix V is obtained as a product of Jacobi plane  
 *  rotations. The rotations are implemented as fast scaled rotations of  
 *  Anda and Park [1]. In the case of underflow of the Jacobi angle, a  
 *  modified Jacobi transformation of Drmac [4] is used. Pivot strategy uses  
 *  column interchanges of de Rijk [2]. The relative accuracy of the computed  
 *  singular values and the accuracy of the computed singular vectors (in  
 *  angle metric) is as guaranteed by the theory of Demmel and Veselic [3].  
 *  The condition number that determines the accuracy in the full rank case  
 *  is essentially min_{D=diag} kappa(A*D), where kappa(.) is the  
 *  spectral condition number. The best performance of this Jacobi SVD  
 *  procedure is achieved if used in an  accelerated version of Drmac and  
 *  Veselic [5,6], and it is the kernel routine in the SIGMA library [7].  
 *  Some tunning parameters (marked with [TP]) are available for the  
 *  implementer.  
 *  The computational range for the nonzero singular values is the  machine  
 *  number interval ( UNDERFLOW , OVERFLOW ). In extreme cases, even  
 *  denormalized singular values can be computed with the corresponding  
 *  gradual loss of accurate digits.  
 *  
 *  Contributors  
 *  ~~~~~~~~~~~~  
 *  Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)  
 *  
 *  References  
 *  ~~~~~~~~~~  
 * [1] A. A. Anda and H. Park: Fast plane rotations with dynamic scaling.  
 *     SIAM J. matrix Anal. Appl., Vol. 15 (1994), pp. 162-174.  
 * [2] P. P. M. De Rijk: A one-sided Jacobi algorithm for computing the  
 *     singular value decomposition on a vector computer.  
 *     SIAM J. Sci. Stat. Comp., Vol. 10 (1998), pp. 359-371.  
 * [3] J. Demmel and K. Veselic: Jacobi method is more accurate than QR.  
 * [4] Z. Drmac: Implementation of Jacobi rotations for accurate singular  
 *     value computation in floating point arithmetic.  
 *     SIAM J. Sci. Comp., Vol. 18 (1997), pp. 1200-1222.  
 * [5] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I.  
 *     SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342.  
 *     LAPACK Working note 169.  
 * [6] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II.  
 *     SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362.  
 *     LAPACK Working note 170.  
 * [7] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV,  
 *     QSVD, (H,K)-SVD computations.  
 *     Department of Mathematics, University of Zagreb, 2008.  
 *  
 *  Bugs, Examples and Comments  
 *  ~~~~~~~~~~~~~~~~~~~~~~~~~~~  
 *  Please report all bugs and send interesting test examples and comments to  
 *  drmac@math.hr. Thank you.  
 *  
 *  Arguments  
 *  =========  
 *  
 *  JOBA    (input) CHARACTER* 1  
 *          Specifies the structure of A.  
 *          = 'L': The input matrix A is lower triangular;  
 *          = 'U': The input matrix A is upper triangular;  
 *          = 'G': The input matrix A is general M-by-N matrix, M >= N.  
 *  
 *  JOBU    (input) CHARACTER*1  
 *          Specifies whether to compute the left singular vectors  
 *          (columns of U):  
 *          = 'U': 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').  
 *          = '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  
 *                 CTOL is given on input in the array WORK.  
 *                 No CTOL smaller than ONE is allowed. CTOL greater  
 *                 than 1 / EPS is meaningless. The option 'C'  
 *                 can be used if M*EPS is satisfactory orthogonality  
 *                 of the computed left singular vectors, so CTOL=M could  
 *                 save few sweeps of Jacobi rotations.  
 *                 See the descriptions of A and WORK(1).  
 *          = 'N': The matrix U is not computed. However, see the  
 *                 description of A.  
 *  
 *  JOBV    (input) 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  
 *          = '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  
 *                  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  
 *                  referenced  
 *  
 *  M       (input) INTEGER  
 *          The number of rows of the input matrix A.  M >= 0.  
 *  
 *  N       (input) INTEGER  
 *          The number of columns of the input matrix A.  
 *          M >= N >= 0.  
 *  
 *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)  
 *          On entry, the M-by-N matrix A.  
 *          On exit :  
 *          If JOBU .EQ. 'U' .OR. JOBU .EQ. 'C' :  
 *                 If INFO .EQ. 0 :  
 *                 RANKA orthonormal columns of U are returned in the  
 *                 leading RANKA columns of the array A. Here RANKA <= N  
 *                 is the number of computed singular values of A that are  
 *                 above the underflow threshold DLAMCH('S'). The singular  
 *                 vectors corresponding to underflowed or zero singular  
 *                 values are not computed. The value of RANKA is returned  
 *                 in the array WORK as RANKA=NINT(WORK(2)). Also see the  
 *                 descriptions of SVA and WORK. The computed columns of U  
 *                 are mutually numerically orthogonal up to approximately  
 *                 TOL=DSQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU.EQ.'C'),  
 *                 see the description of JOBU.  
 *                 If INFO .GT. 0 :  
 *                 the procedure DGESVJ did not converge in the given number  
 *                 of iterations (sweeps). In that case, the computed  
 *                 columns of U may not be orthogonal up to TOL. The output  
 *                 U (stored in A), SIGMA (given by the computed singular  
 *                 values in SVA(1:N)) and V is still a decomposition of the  
 *                 input matrix A in the sense that the residual  
 *                 ||A-SCALE*U*SIGMA*V^T||_2 / ||A||_2 is small.  
 *  
 *          If JOBU .EQ. 'N' :  
 *                 If INFO .EQ. 0 :  
 *                 Note that the left singular vectors are 'for free' in the  
 *                 one-sided Jacobi SVD algorithm. However, if only the  
 *                 singular values are needed, the level of numerical  
 *                 orthogonality of U is not an issue and iterations are  
 *                 stopped when the columns of the iterated matrix are  
 *                 numerically orthogonal up to approximately M*EPS. Thus,  
 *                 on exit, A contains the columns of U scaled with the  
 *                 corresponding singular values.  
 *                 If INFO .GT. 0 :  
 *                 the procedure DGESVJ did not converge in the given number  
 *                 of iterations (sweeps).  
 *  
 *  LDA     (input) INTEGER  
 *          The leading dimension of the array A.  LDA >= max(1,M).  
 *  
 *  SVA     (workspace/output) DOUBLE PRECISION array, dimension (N)  
 *          On exit :  
 *          If INFO .EQ. 0 :  
 *          depending on the value SCALE = WORK(1), we have:  
 *                 If SCALE .EQ. ONE :  
 *                 SVA(1:N) contains the computed singular values of A.  
 *                 During the computation SVA contains the Euclidean column  
 *                 norms of the iterated matrices in the array A.  
 *                 If SCALE .NE. ONE :  
 *                 The singular values of A are SCALE*SVA(1:N), and this  
 *                 factored representation is due to the fact that some of the  
 *                 singular values of A might underflow or overflow.  
 *          If INFO .GT. 0 :  
 *          the procedure DGESVJ did not converge in the given number of  
 *          iterations (sweeps) and SCALE*SVA(1:N) may not be accurate.  
 *  
 *  MV      (input) INTEGER  
 *          If JOBV .EQ. 'A', then the product of Jacobi rotations in DGESVJ  
 *          is applied to the first MV rows of V. See the description of JOBV.  
 *  
 *  V       (input/output) DOUBLE PRECISION array, dimension (LDV,N)  
 *          If JOBV = 'V', then V contains on exit the N-by-N matrix of  
 *                         the right singular vectors;  
 *          If JOBV = 'A', then V contains the product of the computed right  
 *                         singular vector matrix and the initial matrix in  
 *                         the array V.  
 *          If JOBV = 'N', then V is not referenced.  
 *  
 *  LDV     (input) INTEGER  
 *          The leading dimension of the array V, LDV .GE. 1.  
 *          If JOBV .EQ. 'V', then LDV .GE. max(1,N).  
 *          If JOBV .EQ. 'A', then LDV .GE. max(1,MV) .  
 *  
 *  WORK    (input/workspace/output) DOUBLE PRECISION array, dimension max(4,M+N).  
 *          On entry :  
 *          If JOBU .EQ. 'C' :  
 *          WORK(1) = CTOL, where CTOL defines the threshold for convergence.  
 *                    The process stops if all columns of A are mutually  
 *                    orthogonal up to CTOL*EPS, EPS=DLAMCH('E').  
 *                    It is required that CTOL >= ONE, i.e. it is not  
 *                    allowed to force the routine to obtain orthogonality  
 *                    below EPSILON.  
 *          On exit :  
 *          WORK(1) = SCALE is the scaling factor such that SCALE*SVA(1:N)  
 *                    are the computed singular values of A.  
 *                    (See description of SVA().)  
 *          WORK(2) = NINT(WORK(2)) is the number of the computed nonzero  
 *                    singular values.  
 *          WORK(3) = NINT(WORK(3)) is the number of the computed singular  
 *                    values that are larger than the underflow threshold.  
 *          WORK(4) = NINT(WORK(4)) is the number of sweeps of Jacobi  
 *                    rotations needed for numerical convergence.  
 *          WORK(5) = max_{i.NE.j} |COS(A(:,i),A(:,j))| in the last sweep.  
 *                    This is useful information in cases when DGESVJ did  
 *                    not converge, as it can be used to estimate whether  
 *                    the output is stil useful and for post festum analysis.  
 *          WORK(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.  
 *  
 *  LWORK   (input) INTEGER  
 *          length of WORK, WORK >= MAX(6,M+N)  
 *  
 *  INFO    (output) INTEGER  
 *          = 0 : successful exit.  
 *          < 0 : if INFO = -i, then the i-th argument had an illegal value  
 *          > 0 : DGESVJ did not converge in the maximal allowed number (30)  
 *                of sweeps. The output may still be useful. See the  
 *                description of WORK.  
 *  
 *  =====================================================================  *  =====================================================================
 *  *
 *     .. Local Parameters ..  *     .. Local Parameters ..
       DOUBLE PRECISION   ZERO, HALF, ONE, TWO        DOUBLE PRECISION   ZERO, HALF, ONE
       PARAMETER          ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0,        PARAMETER          ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0)
      +                   TWO = 2.0D0 )  
       INTEGER            NSWEEP        INTEGER            NSWEEP
       PARAMETER          ( NSWEEP = 30 )        PARAMETER          ( NSWEEP = 30 )
 *     ..  *     ..
 *     .. Local Scalars ..  *     .. Local Scalars ..
       DOUBLE PRECISION   AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG,        DOUBLE PRECISION   AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG,
      +                   BIGTHETA, CS, CTOL, EPSILON, LARGE, MXAAPQ,       $                   BIGTHETA, CS, CTOL, EPSLN, LARGE, MXAAPQ,
      +                   MXSINJ, ROOTBIG, ROOTEPS, ROOTSFMIN, ROOTTOL,       $                   MXSINJ, ROOTBIG, ROOTEPS, ROOTSFMIN, ROOTTOL,
      +                   SCALE, SFMIN, SMALL, SN, T, TEMP1, THETA,       $                   SKL, SFMIN, SMALL, SN, T, TEMP1, THETA,
      +                   THSIGN, TOL       $                   THSIGN, TOL
       INTEGER            BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1,        INTEGER            BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1,
      +                   ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, N2, N34,       $                   ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, N2, N34,
      +                   N4, NBL, NOTROT, p, PSKIPPED, q, ROWSKIP,       $                   N4, NBL, NOTROT, p, PSKIPPED, q, ROWSKIP,
      +                   SWBAND       $                   SWBAND
       LOGICAL            APPLV, GOSCALE, LOWER, LSVEC, NOSCALE, ROTOK,        LOGICAL            APPLV, GOSCALE, LOWER, LSVEC, NOSCALE, ROTOK,
      +                   RSVEC, UCTOL, UPPER       $                   RSVEC, UCTOL, UPPER
 *     ..  *     ..
 *     .. Local Arrays ..  *     .. Local Arrays ..
       DOUBLE PRECISION   FASTR( 5 )        DOUBLE PRECISION   FASTR( 5 )
Line 327 Line 424
       ELSE IF( MV.LT.0 ) THEN        ELSE IF( MV.LT.0 ) THEN
          INFO = -9           INFO = -9
       ELSE IF( ( RSVEC .AND. ( LDV.LT.N ) ) .OR.        ELSE IF( ( RSVEC .AND. ( LDV.LT.N ) ) .OR.
      +         ( APPLV .AND. ( LDV.LT.MV ) ) ) THEN       $         ( APPLV .AND. ( LDV.LT.MV ) ) ) THEN
          INFO = -11           INFO = -11
       ELSE IF( UCTOL .AND. ( WORK( 1 ).LE.ONE ) ) THEN        ELSE IF( UCTOL .AND. ( WORK( 1 ).LE.ONE ) ) THEN
          INFO = -12           INFO = -12
Line 368 Line 465
 *     ... and the machine dependent parameters are  *     ... and the machine dependent parameters are
 *[!]  (Make sure that DLAMCH() works properly on the target machine.)  *[!]  (Make sure that DLAMCH() works properly on the target machine.)
 *  *
       EPSILON = DLAMCH( 'Epsilon' )        EPSLN = DLAMCH( 'Epsilon' )
       ROOTEPS = DSQRT( EPSILON )        ROOTEPS = DSQRT( EPSLN )
       SFMIN = DLAMCH( 'SafeMinimum' )        SFMIN = DLAMCH( 'SafeMinimum' )
       ROOTSFMIN = DSQRT( SFMIN )        ROOTSFMIN = DSQRT( SFMIN )
       SMALL = SFMIN / EPSILON        SMALL = SFMIN / EPSLN
       BIG = DLAMCH( 'Overflow' )        BIG = DLAMCH( 'Overflow' )
 *     BIG         = ONE    / SFMIN  *     BIG         = ONE    / SFMIN
       ROOTBIG = ONE / ROOTSFMIN        ROOTBIG = ONE / ROOTSFMIN
       LARGE = BIG / DSQRT( DBLE( M*N ) )        LARGE = BIG / DSQRT( DBLE( M*N ) )
       BIGTHETA = ONE / ROOTEPS        BIGTHETA = ONE / ROOTEPS
 *  *
       TOL = CTOL*EPSILON        TOL = CTOL*EPSLN
       ROOTTOL = DSQRT( TOL )        ROOTTOL = DSQRT( TOL )
 *  *
       IF( DBLE( M )*EPSILON.GE.ONE ) THEN        IF( DBLE( M )*EPSLN.GE.ONE ) THEN
          INFO = -5           INFO = -4
          CALL XERBLA( 'DGESVJ', -INFO )           CALL XERBLA( 'DGESVJ', -INFO )
          RETURN           RETURN
       END IF        END IF
Line 407 Line 504
 *     DSQRT(N)*max_i SVA(i) does not overflow. If INFinite entries  *     DSQRT(N)*max_i SVA(i) does not overflow. If INFinite entries
 *     in A are detected, the procedure returns with INFO=-6.  *     in A are detected, the procedure returns with INFO=-6.
 *  *
       SCALE = ONE / DSQRT( DBLE( M )*DBLE( N ) )        SKL= ONE / DSQRT( DBLE( M )*DBLE( N ) )
       NOSCALE = .TRUE.        NOSCALE = .TRUE.
       GOSCALE = .TRUE.        GOSCALE = .TRUE.
 *  *
Line 415 Line 512
 *        the input matrix is M-by-N lower triangular (trapezoidal)  *        the input matrix is M-by-N lower triangular (trapezoidal)
          DO 1874 p = 1, N           DO 1874 p = 1, N
             AAPP = ZERO              AAPP = ZERO
             AAQQ = ZERO              AAQQ = ONE
             CALL DLASSQ( M-p+1, A( p, p ), 1, AAPP, AAQQ )              CALL DLASSQ( M-p+1, A( p, p ), 1, AAPP, AAQQ )
             IF( AAPP.GT.BIG ) THEN              IF( AAPP.GT.BIG ) THEN
                INFO = -6                 INFO = -6
Line 427 Line 524
                SVA( p ) = AAPP*AAQQ                 SVA( p ) = AAPP*AAQQ
             ELSE              ELSE
                NOSCALE = .FALSE.                 NOSCALE = .FALSE.
                SVA( p ) = AAPP*( AAQQ*SCALE )                 SVA( p ) = AAPP*( AAQQ*SKL)
                IF( GOSCALE ) THEN                 IF( GOSCALE ) THEN
                   GOSCALE = .FALSE.                    GOSCALE = .FALSE.
                   DO 1873 q = 1, p - 1                    DO 1873 q = 1, p - 1
                      SVA( q ) = SVA( q )*SCALE                       SVA( q ) = SVA( q )*SKL
  1873             CONTINUE   1873             CONTINUE
                END IF                 END IF
             END IF              END IF
Line 440 Line 537
 *        the input matrix is M-by-N upper triangular (trapezoidal)  *        the input matrix is M-by-N upper triangular (trapezoidal)
          DO 2874 p = 1, N           DO 2874 p = 1, N
             AAPP = ZERO              AAPP = ZERO
             AAQQ = ZERO              AAQQ = ONE
             CALL DLASSQ( p, A( 1, p ), 1, AAPP, AAQQ )              CALL DLASSQ( p, A( 1, p ), 1, AAPP, AAQQ )
             IF( AAPP.GT.BIG ) THEN              IF( AAPP.GT.BIG ) THEN
                INFO = -6                 INFO = -6
Line 452 Line 549
                SVA( p ) = AAPP*AAQQ                 SVA( p ) = AAPP*AAQQ
             ELSE              ELSE
                NOSCALE = .FALSE.                 NOSCALE = .FALSE.
                SVA( p ) = AAPP*( AAQQ*SCALE )                 SVA( p ) = AAPP*( AAQQ*SKL)
                IF( GOSCALE ) THEN                 IF( GOSCALE ) THEN
                   GOSCALE = .FALSE.                    GOSCALE = .FALSE.
                   DO 2873 q = 1, p - 1                    DO 2873 q = 1, p - 1
                      SVA( q ) = SVA( q )*SCALE                       SVA( q ) = SVA( q )*SKL
  2873             CONTINUE   2873             CONTINUE
                END IF                 END IF
             END IF              END IF
Line 465 Line 562
 *        the input matrix is M-by-N general dense  *        the input matrix is M-by-N general dense
          DO 3874 p = 1, N           DO 3874 p = 1, N
             AAPP = ZERO              AAPP = ZERO
             AAQQ = ZERO              AAQQ = ONE
             CALL DLASSQ( M, A( 1, p ), 1, AAPP, AAQQ )              CALL DLASSQ( M, A( 1, p ), 1, AAPP, AAQQ )
             IF( AAPP.GT.BIG ) THEN              IF( AAPP.GT.BIG ) THEN
                INFO = -6                 INFO = -6
Line 477 Line 574
                SVA( p ) = AAPP*AAQQ                 SVA( p ) = AAPP*AAQQ
             ELSE              ELSE
                NOSCALE = .FALSE.                 NOSCALE = .FALSE.
                SVA( p ) = AAPP*( AAQQ*SCALE )                 SVA( p ) = AAPP*( AAQQ*SKL)
                IF( GOSCALE ) THEN                 IF( GOSCALE ) THEN
                   GOSCALE = .FALSE.                    GOSCALE = .FALSE.
                   DO 3873 q = 1, p - 1                    DO 3873 q = 1, p - 1
                      SVA( q ) = SVA( q )*SCALE                       SVA( q ) = SVA( q )*SKL
  3873             CONTINUE   3873             CONTINUE
                END IF                 END IF
             END IF              END IF
  3874    CONTINUE   3874    CONTINUE
       END IF        END IF
 *  *
       IF( NOSCALE )SCALE = ONE        IF( NOSCALE )SKL= ONE
 *  *
 *     Move the smaller part of the spectrum from the underflow threshold  *     Move the smaller part of the spectrum from the underflow threshold
 *(!)  Start by determining the position of the nonzero entries of the  *(!)  Start by determining the position of the nonzero entries of the
Line 517 Line 614
 * #:) Quick return for one-column matrix  * #:) Quick return for one-column matrix
 *  *
       IF( N.EQ.1 ) THEN        IF( N.EQ.1 ) THEN
          IF( LSVEC )CALL DLASCL( 'G', 0, 0, SVA( 1 ), SCALE, M, 1,           IF( LSVEC )CALL DLASCL( 'G', 0, 0, SVA( 1 ), SKL, M, 1,
      +                           A( 1, 1 ), LDA, IERR )       $                           A( 1, 1 ), LDA, IERR )
          WORK( 1 ) = ONE / SCALE           WORK( 1 ) = ONE / SKL
          IF( SVA( 1 ).GE.SFMIN ) THEN           IF( SVA( 1 ).GE.SFMIN ) THEN
             WORK( 2 ) = ONE              WORK( 2 ) = ONE
          ELSE           ELSE
Line 535 Line 632
 *     Protect small singular values from underflow, and try to  *     Protect small singular values from underflow, and try to
 *     avoid underflows/overflows in computing Jacobi rotations.  *     avoid underflows/overflows in computing Jacobi rotations.
 *  *
       SN = DSQRT( SFMIN / EPSILON )        SN = DSQRT( SFMIN / EPSLN )
       TEMP1 = DSQRT( BIG / DBLE( N ) )        TEMP1 = DSQRT( BIG / DBLE( N ) )
       IF( ( AAPP.LE.SN ) .OR. ( AAQQ.GE.TEMP1 ) .OR.        IF( ( AAPP.LE.SN ) .OR. ( AAQQ.GE.TEMP1 ) .OR.
      +    ( ( SN.LE.AAQQ ) .AND. ( AAPP.LE.TEMP1 ) ) ) THEN       $    ( ( SN.LE.AAQQ ) .AND. ( AAPP.LE.TEMP1 ) ) ) THEN
          TEMP1 = DMIN1( BIG, TEMP1 / AAPP )           TEMP1 = DMIN1( BIG, TEMP1 / AAPP )
 *         AAQQ  = AAQQ*TEMP1  *         AAQQ  = AAQQ*TEMP1
 *         AAPP  = AAPP*TEMP1  *         AAPP  = AAPP*TEMP1
Line 563 Line 660
       IF( TEMP1.NE.ONE ) THEN        IF( TEMP1.NE.ONE ) THEN
          CALL DLASCL( 'G', 0, 0, ONE, TEMP1, N, 1, SVA, N, IERR )           CALL DLASCL( 'G', 0, 0, ONE, TEMP1, N, 1, SVA, N, IERR )
       END IF        END IF
       SCALE = TEMP1*SCALE        SKL= TEMP1*SKL
       IF( SCALE.NE.ONE ) THEN        IF( SKL.NE.ONE ) THEN
          CALL DLASCL( JOBA, 0, 0, ONE, SCALE, M, N, A, LDA, IERR )           CALL DLASCL( JOBA, 0, 0, ONE, SKL, M, N, A, LDA, IERR )
          SCALE = ONE / SCALE           SKL= ONE / SKL
       END IF        END IF
 *  *
 *     Row-cyclic Jacobi SVD algorithm with column pivoting  *     Row-cyclic Jacobi SVD algorithm with column pivoting
Line 638 Line 735
 *     [+ + x x]                    [x x].             [x x]  *     [+ + x x]                    [x x].             [x x]
 *  *
             CALL DGSVJ0( JOBV, M-N34, N-N34, A( N34+1, N34+1 ), LDA,              CALL DGSVJ0( JOBV, M-N34, N-N34, A( N34+1, N34+1 ), LDA,
      +                   WORK( N34+1 ), SVA( N34+1 ), MVL,       $                   WORK( N34+1 ), SVA( N34+1 ), MVL,
      +                   V( N34*q+1, N34+1 ), LDV, EPSILON, SFMIN, TOL,       $                   V( N34*q+1, N34+1 ), LDV, EPSLN, SFMIN, TOL,
      +                   2, WORK( N+1 ), LWORK-N, IERR )       $                   2, WORK( N+1 ), LWORK-N, IERR )
 *  *
             CALL DGSVJ0( JOBV, M-N2, N34-N2, A( N2+1, N2+1 ), LDA,              CALL DGSVJ0( JOBV, M-N2, N34-N2, A( N2+1, N2+1 ), LDA,
      +                   WORK( N2+1 ), SVA( N2+1 ), MVL,       $                   WORK( N2+1 ), SVA( N2+1 ), MVL,
      +                   V( N2*q+1, N2+1 ), LDV, EPSILON, SFMIN, TOL, 2,       $                   V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 2,
      +                   WORK( N+1 ), LWORK-N, IERR )       $                   WORK( N+1 ), LWORK-N, IERR )
 *  *
             CALL DGSVJ1( JOBV, M-N2, N-N2, N4, A( N2+1, N2+1 ), LDA,              CALL DGSVJ1( JOBV, M-N2, N-N2, N4, A( N2+1, N2+1 ), LDA,
      +                   WORK( N2+1 ), SVA( N2+1 ), MVL,       $                   WORK( N2+1 ), SVA( N2+1 ), MVL,
      +                   V( N2*q+1, N2+1 ), LDV, EPSILON, SFMIN, TOL, 1,       $                   V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 1,
      +                   WORK( N+1 ), LWORK-N, IERR )       $                   WORK( N+1 ), LWORK-N, IERR )
 *  *
             CALL DGSVJ0( JOBV, M-N4, N2-N4, A( N4+1, N4+1 ), LDA,              CALL DGSVJ0( JOBV, M-N4, N2-N4, A( N4+1, N4+1 ), LDA,
      +                   WORK( N4+1 ), SVA( N4+1 ), MVL,       $                   WORK( N4+1 ), SVA( N4+1 ), MVL,
      +                   V( N4*q+1, N4+1 ), LDV, EPSILON, SFMIN, TOL, 1,       $                   V( N4*q+1, N4+1 ), LDV, EPSLN, SFMIN, TOL, 1,
      +                   WORK( N+1 ), LWORK-N, IERR )       $                   WORK( N+1 ), LWORK-N, IERR )
 *  *
             CALL DGSVJ0( JOBV, M, N4, A, LDA, WORK, SVA, MVL, V, LDV,              CALL DGSVJ0( JOBV, M, N4, A, LDA, WORK, SVA, MVL, V, LDV,
      +                   EPSILON, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,       $                   EPSLN, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,
      +                   IERR )       $                   IERR )
 *  *
             CALL DGSVJ1( JOBV, M, N2, N4, A, LDA, WORK, SVA, MVL, V,              CALL DGSVJ1( JOBV, M, N2, N4, A, LDA, WORK, SVA, MVL, V,
      +                   LDV, EPSILON, SFMIN, TOL, 1, WORK( N+1 ),       $                   LDV, EPSLN, SFMIN, TOL, 1, WORK( N+1 ),
      +                   LWORK-N, IERR )       $                   LWORK-N, IERR )
 *  *
 *  *
          ELSE IF( UPPER ) THEN           ELSE IF( UPPER ) THEN
 *  *
 *  *
             CALL DGSVJ0( JOBV, N4, N4, A, LDA, WORK, SVA, MVL, V, LDV,              CALL DGSVJ0( JOBV, N4, N4, A, LDA, WORK, SVA, MVL, V, LDV,
      +                   EPSILON, SFMIN, TOL, 2, WORK( N+1 ), LWORK-N,       $                   EPSLN, SFMIN, TOL, 2, WORK( N+1 ), LWORK-N,
      +                   IERR )       $                   IERR )
 *  *
             CALL DGSVJ0( JOBV, N2, N4, A( 1, N4+1 ), LDA, WORK( N4+1 ),              CALL DGSVJ0( JOBV, N2, N4, A( 1, N4+1 ), LDA, WORK( N4+1 ),
      +                   SVA( N4+1 ), MVL, V( N4*q+1, N4+1 ), LDV,       $                   SVA( N4+1 ), MVL, V( N4*q+1, N4+1 ), LDV,
      +                   EPSILON, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,       $                   EPSLN, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,
      +                   IERR )       $                   IERR )
 *  *
             CALL DGSVJ1( JOBV, N2, N2, N4, A, LDA, WORK, SVA, MVL, V,              CALL DGSVJ1( JOBV, N2, N2, N4, A, LDA, WORK, SVA, MVL, V,
      +                   LDV, EPSILON, SFMIN, TOL, 1, WORK( N+1 ),       $                   LDV, EPSLN, SFMIN, TOL, 1, WORK( N+1 ),
      +                   LWORK-N, IERR )       $                   LWORK-N, IERR )
 *  *
             CALL DGSVJ0( JOBV, N2+N4, N4, A( 1, N2+1 ), LDA,              CALL DGSVJ0( JOBV, N2+N4, N4, A( 1, N2+1 ), LDA,
      +                   WORK( N2+1 ), SVA( N2+1 ), MVL,       $                   WORK( N2+1 ), SVA( N2+1 ), MVL,
      +                   V( N2*q+1, N2+1 ), LDV, EPSILON, SFMIN, TOL, 1,       $                   V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 1,
      +                   WORK( N+1 ), LWORK-N, IERR )       $                   WORK( N+1 ), LWORK-N, IERR )
   
          END IF           END IF
 *  *
Line 725 Line 822
                   IF( p.NE.q ) THEN                    IF( p.NE.q ) THEN
                      CALL DSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )                       CALL DSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
                      IF( RSVEC )CALL DSWAP( MVL, V( 1, p ), 1,                       IF( RSVEC )CALL DSWAP( MVL, V( 1, p ), 1,
      +                                      V( 1, q ), 1 )       $                                      V( 1, q ), 1 )
                      TEMP1 = SVA( p )                       TEMP1 = SVA( p )
                      SVA( p ) = SVA( q )                       SVA( p ) = SVA( q )
                      SVA( q ) = TEMP1                       SVA( q ) = TEMP1
Line 749 Line 846
 *        below should read "AAPP = DNRM2( M, A(1,p), 1 ) * WORK(p)".  *        below should read "AAPP = DNRM2( M, A(1,p), 1 ) * WORK(p)".
 *  *
                      IF( ( SVA( p ).LT.ROOTBIG ) .AND.                       IF( ( SVA( p ).LT.ROOTBIG ) .AND.
      +                   ( SVA( p ).GT.ROOTSFMIN ) ) THEN       $                   ( SVA( p ).GT.ROOTSFMIN ) ) THEN
                         SVA( p ) = DNRM2( M, A( 1, p ), 1 )*WORK( p )                          SVA( p ) = DNRM2( M, A( 1, p ), 1 )*WORK( p )
                      ELSE                       ELSE
                         TEMP1 = ZERO                          TEMP1 = ZERO
                         AAPP = ZERO                          AAPP = ONE
                         CALL DLASSQ( M, A( 1, p ), 1, TEMP1, AAPP )                          CALL DLASSQ( M, A( 1, p ), 1, TEMP1, AAPP )
                         SVA( p ) = TEMP1*DSQRT( AAPP )*WORK( p )                          SVA( p ) = TEMP1*DSQRT( AAPP )*WORK( p )
                      END IF                       END IF
Line 777 Line 874
                               ROTOK = ( SMALL*AAPP ).LE.AAQQ                                ROTOK = ( SMALL*AAPP ).LE.AAQQ
                               IF( AAPP.LT.( BIG / AAQQ ) ) THEN                                IF( AAPP.LT.( BIG / AAQQ ) ) THEN
                                  AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,                                   AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
      +                                  q ), 1 )*WORK( p )*WORK( q ) /       $                                  q ), 1 )*WORK( p )*WORK( q ) /
      +                                  AAQQ ) / AAPP       $                                  AAQQ ) / AAPP
                               ELSE                                ELSE
                                  CALL DCOPY( M, A( 1, p ), 1,                                   CALL DCOPY( M, A( 1, p ), 1,
      +                                       WORK( N+1 ), 1 )       $                                       WORK( N+1 ), 1 )
                                  CALL DLASCL( 'G', 0, 0, AAPP,                                   CALL DLASCL( 'G', 0, 0, AAPP,
      +                                        WORK( p ), M, 1,       $                                        WORK( p ), M, 1,
      +                                        WORK( N+1 ), LDA, IERR )       $                                        WORK( N+1 ), LDA, IERR )
                                  AAPQ = DDOT( M, WORK( N+1 ), 1,                                   AAPQ = DDOT( M, WORK( N+1 ), 1,
      +                                  A( 1, q ), 1 )*WORK( q ) / AAQQ       $                                  A( 1, q ), 1 )*WORK( q ) / AAQQ
                               END IF                                END IF
                            ELSE                             ELSE
                               ROTOK = AAPP.LE.( AAQQ / SMALL )                                ROTOK = AAPP.LE.( AAQQ / SMALL )
                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN                                IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
                                  AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,                                   AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
      +                                  q ), 1 )*WORK( p )*WORK( q ) /       $                                  q ), 1 )*WORK( p )*WORK( q ) /
      +                                  AAQQ ) / AAPP       $                                  AAQQ ) / AAPP
                               ELSE                                ELSE
                                  CALL DCOPY( M, A( 1, q ), 1,                                   CALL DCOPY( M, A( 1, q ), 1,
      +                                       WORK( N+1 ), 1 )       $                                       WORK( N+1 ), 1 )
                                  CALL DLASCL( 'G', 0, 0, AAQQ,                                   CALL DLASCL( 'G', 0, 0, AAQQ,
      +                                        WORK( q ), M, 1,       $                                        WORK( q ), M, 1,
      +                                        WORK( N+1 ), LDA, IERR )       $                                        WORK( N+1 ), LDA, IERR )
                                  AAPQ = DDOT( M, WORK( N+1 ), 1,                                   AAPQ = DDOT( M, WORK( N+1 ), 1,
      +                                  A( 1, p ), 1 )*WORK( p ) / AAPP       $                                  A( 1, p ), 1 )*WORK( p ) / AAPP
                               END IF                                END IF
                            END IF                             END IF
 *  *
Line 824 Line 921
 *  *
                                  AQOAP = AAQQ / AAPP                                   AQOAP = AAQQ / AAPP
                                  APOAQ = AAPP / AAQQ                                   APOAQ = AAPP / AAQQ
                                  THETA = -HALF*DABS( AQOAP-APOAQ ) /                                   THETA = -HALF*DABS(AQOAP-APOAQ)/AAPQ
      +                                   AAPQ  
 *  *
                                  IF( DABS( THETA ).GT.BIGTHETA ) THEN                                   IF( DABS( THETA ).GT.BIGTHETA ) THEN
 *  *
                                     T = HALF / THETA                                      T = HALF / THETA
                                     FASTR( 3 ) = T*WORK( p ) / WORK( q )                                      FASTR( 3 ) = T*WORK( p ) / WORK( q )
                                     FASTR( 4 ) = -T*WORK( q ) /                                      FASTR( 4 ) = -T*WORK( q ) /
      +                                           WORK( p )       $                                           WORK( p )
                                     CALL DROTM( M, A( 1, p ), 1,                                      CALL DROTM( M, A( 1, p ), 1,
      +                                          A( 1, q ), 1, FASTR )       $                                          A( 1, q ), 1, FASTR )
                                     IF( RSVEC )CALL DROTM( MVL,                                      IF( RSVEC )CALL DROTM( MVL,
      +                                              V( 1, p ), 1,       $                                              V( 1, p ), 1,
      +                                              V( 1, q ), 1,       $                                              V( 1, q ), 1,
      +                                              FASTR )       $                                              FASTR )
                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,                                      SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
      +                                         ONE+T*APOAQ*AAPQ ) )       $                                         ONE+T*APOAQ*AAPQ ) )
                                     AAPP = AAPP*DSQRT( ONE-T*AQOAP*                                      AAPP = AAPP*DSQRT( DMAX1( ZERO,
      +                                     AAPQ )       $                                     ONE-T*AQOAP*AAPQ ) )
                                     MXSINJ = DMAX1( MXSINJ, DABS( T ) )                                      MXSINJ = DMAX1( MXSINJ, DABS( T ) )
 *  *
                                  ELSE                                   ELSE
Line 851 Line 947
 *  *
                                     THSIGN = -DSIGN( ONE, AAPQ )                                      THSIGN = -DSIGN( ONE, AAPQ )
                                     T = ONE / ( THETA+THSIGN*                                      T = ONE / ( THETA+THSIGN*
      +                                  DSQRT( ONE+THETA*THETA ) )       $                                  DSQRT( ONE+THETA*THETA ) )
                                     CS = DSQRT( ONE / ( ONE+T*T ) )                                      CS = DSQRT( ONE / ( ONE+T*T ) )
                                     SN = T*CS                                      SN = T*CS
 *  *
                                     MXSINJ = DMAX1( MXSINJ, DABS( SN ) )                                      MXSINJ = DMAX1( MXSINJ, DABS( SN ) )
                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,                                      SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
      +                                         ONE+T*APOAQ*AAPQ ) )       $                                         ONE+T*APOAQ*AAPQ ) )
                                     AAPP = AAPP*DSQRT( DMAX1( ZERO,                                      AAPP = AAPP*DSQRT( DMAX1( ZERO,
      +                                     ONE-T*AQOAP*AAPQ ) )       $                                     ONE-T*AQOAP*AAPQ ) )
 *  *
                                     APOAQ = WORK( p ) / WORK( q )                                      APOAQ = WORK( p ) / WORK( q )
                                     AQOAP = WORK( q ) / WORK( p )                                      AQOAP = WORK( q ) / WORK( p )
Line 870 Line 966
                                           WORK( p ) = WORK( p )*CS                                            WORK( p ) = WORK( p )*CS
                                           WORK( q ) = WORK( q )*CS                                            WORK( q ) = WORK( q )*CS
                                           CALL DROTM( M, A( 1, p ), 1,                                            CALL DROTM( M, A( 1, p ), 1,
      +                                                A( 1, q ), 1,       $                                                A( 1, q ), 1,
      +                                                FASTR )       $                                                FASTR )
                                           IF( RSVEC )CALL DROTM( MVL,                                            IF( RSVEC )CALL DROTM( MVL,
      +                                        V( 1, p ), 1, V( 1, q ),       $                                        V( 1, p ), 1, V( 1, q ),
      +                                        1, FASTR )       $                                        1, FASTR )
                                        ELSE                                         ELSE
                                           CALL DAXPY( M, -T*AQOAP,                                            CALL DAXPY( M, -T*AQOAP,
      +                                                A( 1, q ), 1,       $                                                A( 1, q ), 1,
      +                                                A( 1, p ), 1 )       $                                                A( 1, p ), 1 )
                                           CALL DAXPY( M, CS*SN*APOAQ,                                            CALL DAXPY( M, CS*SN*APOAQ,
      +                                                A( 1, p ), 1,       $                                                A( 1, p ), 1,
      +                                                A( 1, q ), 1 )       $                                                A( 1, q ), 1 )
                                           WORK( p ) = WORK( p )*CS                                            WORK( p ) = WORK( p )*CS
                                           WORK( q ) = WORK( q ) / CS                                            WORK( q ) = WORK( q ) / CS
                                           IF( RSVEC ) THEN                                            IF( RSVEC ) THEN
                                              CALL DAXPY( MVL, -T*AQOAP,                                               CALL DAXPY( MVL, -T*AQOAP,
      +                                                   V( 1, q ), 1,       $                                                   V( 1, q ), 1,
      +                                                   V( 1, p ), 1 )       $                                                   V( 1, p ), 1 )
                                              CALL DAXPY( MVL,                                               CALL DAXPY( MVL,
      +                                                   CS*SN*APOAQ,       $                                                   CS*SN*APOAQ,
      +                                                   V( 1, p ), 1,       $                                                   V( 1, p ), 1,
      +                                                   V( 1, q ), 1 )       $                                                   V( 1, q ), 1 )
                                           END IF                                            END IF
                                        END IF                                         END IF
                                     ELSE                                      ELSE
                                        IF( WORK( q ).GE.ONE ) THEN                                         IF( WORK( q ).GE.ONE ) THEN
                                           CALL DAXPY( M, T*APOAQ,                                            CALL DAXPY( M, T*APOAQ,
      +                                                A( 1, p ), 1,       $                                                A( 1, p ), 1,
      +                                                A( 1, q ), 1 )       $                                                A( 1, q ), 1 )
                                           CALL DAXPY( M, -CS*SN*AQOAP,                                            CALL DAXPY( M, -CS*SN*AQOAP,
      +                                                A( 1, q ), 1,       $                                                A( 1, q ), 1,
      +                                                A( 1, p ), 1 )       $                                                A( 1, p ), 1 )
                                           WORK( p ) = WORK( p ) / CS                                            WORK( p ) = WORK( p ) / CS
                                           WORK( q ) = WORK( q )*CS                                            WORK( q ) = WORK( q )*CS
                                           IF( RSVEC ) THEN                                            IF( RSVEC ) THEN
                                              CALL DAXPY( MVL, T*APOAQ,                                               CALL DAXPY( MVL, T*APOAQ,
      +                                                   V( 1, p ), 1,       $                                                   V( 1, p ), 1,
      +                                                   V( 1, q ), 1 )       $                                                   V( 1, q ), 1 )
                                              CALL DAXPY( MVL,                                               CALL DAXPY( MVL,
      +                                                   -CS*SN*AQOAP,       $                                                   -CS*SN*AQOAP,
      +                                                   V( 1, q ), 1,       $                                                   V( 1, q ), 1,
      +                                                   V( 1, p ), 1 )       $                                                   V( 1, p ), 1 )
                                           END IF                                            END IF
                                        ELSE                                         ELSE
                                           IF( WORK( p ).GE.WORK( q ) )                                            IF( WORK( p ).GE.WORK( q ) )
      +                                        THEN       $                                        THEN
                                              CALL DAXPY( M, -T*AQOAP,                                               CALL DAXPY( M, -T*AQOAP,
      +                                                   A( 1, q ), 1,       $                                                   A( 1, q ), 1,
      +                                                   A( 1, p ), 1 )       $                                                   A( 1, p ), 1 )
                                              CALL DAXPY( M, CS*SN*APOAQ,                                               CALL DAXPY( M, CS*SN*APOAQ,
      +                                                   A( 1, p ), 1,       $                                                   A( 1, p ), 1,
      +                                                   A( 1, q ), 1 )       $                                                   A( 1, q ), 1 )
                                              WORK( p ) = WORK( p )*CS                                               WORK( p ) = WORK( p )*CS
                                              WORK( q ) = WORK( q ) / CS                                               WORK( q ) = WORK( q ) / CS
                                              IF( RSVEC ) THEN                                               IF( RSVEC ) THEN
                                                 CALL DAXPY( MVL,                                                  CALL DAXPY( MVL,
      +                                               -T*AQOAP,       $                                               -T*AQOAP,
      +                                               V( 1, q ), 1,       $                                               V( 1, q ), 1,
      +                                               V( 1, p ), 1 )       $                                               V( 1, p ), 1 )
                                                 CALL DAXPY( MVL,                                                  CALL DAXPY( MVL,
      +                                               CS*SN*APOAQ,       $                                               CS*SN*APOAQ,
      +                                               V( 1, p ), 1,       $                                               V( 1, p ), 1,
      +                                               V( 1, q ), 1 )       $                                               V( 1, q ), 1 )
                                              END IF                                               END IF
                                           ELSE                                            ELSE
                                              CALL DAXPY( M, T*APOAQ,                                               CALL DAXPY( M, T*APOAQ,
      +                                                   A( 1, p ), 1,       $                                                   A( 1, p ), 1,
      +                                                   A( 1, q ), 1 )       $                                                   A( 1, q ), 1 )
                                              CALL DAXPY( M,                                               CALL DAXPY( M,
      +                                                   -CS*SN*AQOAP,       $                                                   -CS*SN*AQOAP,
      +                                                   A( 1, q ), 1,       $                                                   A( 1, q ), 1,
      +                                                   A( 1, p ), 1 )       $                                                   A( 1, p ), 1 )
                                              WORK( p ) = WORK( p ) / CS                                               WORK( p ) = WORK( p ) / CS
                                              WORK( q ) = WORK( q )*CS                                               WORK( q ) = WORK( q )*CS
                                              IF( RSVEC ) THEN                                               IF( RSVEC ) THEN
                                                 CALL DAXPY( MVL,                                                  CALL DAXPY( MVL,
      +                                               T*APOAQ, V( 1, p ),       $                                               T*APOAQ, V( 1, p ),
      +                                               1, V( 1, q ), 1 )       $                                               1, V( 1, q ), 1 )
                                                 CALL DAXPY( MVL,                                                  CALL DAXPY( MVL,
      +                                               -CS*SN*AQOAP,       $                                               -CS*SN*AQOAP,
      +                                               V( 1, q ), 1,       $                                               V( 1, q ), 1,
      +                                               V( 1, p ), 1 )       $                                               V( 1, p ), 1 )
                                              END IF                                               END IF
                                           END IF                                            END IF
                                        END IF                                         END IF
Line 961 Line 1057
                               ELSE                                ELSE
 *              .. have to use modified Gram-Schmidt like transformation  *              .. have to use modified Gram-Schmidt like transformation
                                  CALL DCOPY( M, A( 1, p ), 1,                                   CALL DCOPY( M, A( 1, p ), 1,
      +                                       WORK( N+1 ), 1 )       $                                       WORK( N+1 ), 1 )
                                  CALL DLASCL( 'G', 0, 0, AAPP, ONE, M,                                   CALL DLASCL( 'G', 0, 0, AAPP, ONE, M,
      +                                        1, WORK( N+1 ), LDA,       $                                        1, WORK( N+1 ), LDA,
      +                                        IERR )       $                                        IERR )
                                  CALL DLASCL( 'G', 0, 0, AAQQ, ONE, M,                                   CALL DLASCL( 'G', 0, 0, AAQQ, ONE, M,
      +                                        1, A( 1, q ), LDA, IERR )       $                                        1, A( 1, q ), LDA, IERR )
                                  TEMP1 = -AAPQ*WORK( p ) / WORK( q )                                   TEMP1 = -AAPQ*WORK( p ) / WORK( q )
                                  CALL DAXPY( M, TEMP1, WORK( N+1 ), 1,                                   CALL DAXPY( M, TEMP1, WORK( N+1 ), 1,
      +                                       A( 1, q ), 1 )       $                                       A( 1, q ), 1 )
                                  CALL DLASCL( 'G', 0, 0, ONE, AAQQ, M,                                   CALL DLASCL( 'G', 0, 0, ONE, AAQQ, M,
      +                                        1, A( 1, q ), LDA, IERR )       $                                        1, A( 1, q ), LDA, IERR )
                                  SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,                                   SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
      +                                      ONE-AAPQ*AAPQ ) )       $                                      ONE-AAPQ*AAPQ ) )
                                  MXSINJ = DMAX1( MXSINJ, SFMIN )                                   MXSINJ = DMAX1( MXSINJ, SFMIN )
                               END IF                                END IF
 *           END IF ROTOK THEN ... ELSE  *           END IF ROTOK THEN ... ELSE
Line 982 Line 1078
 *           recompute SVA(q), SVA(p).  *           recompute SVA(q), SVA(p).
 *  *
                               IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )                                IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
      +                            THEN       $                            THEN
                                  IF( ( AAQQ.LT.ROOTBIG ) .AND.                                   IF( ( AAQQ.LT.ROOTBIG ) .AND.
      +                               ( AAQQ.GT.ROOTSFMIN ) ) THEN       $                               ( AAQQ.GT.ROOTSFMIN ) ) THEN
                                     SVA( q ) = DNRM2( M, A( 1, q ), 1 )*                                      SVA( q ) = DNRM2( M, A( 1, q ), 1 )*
      +                                         WORK( q )       $                                         WORK( q )
                                  ELSE                                   ELSE
                                     T = ZERO                                      T = ZERO
                                     AAQQ = ZERO                                      AAQQ = ONE
                                     CALL DLASSQ( M, A( 1, q ), 1, T,                                      CALL DLASSQ( M, A( 1, q ), 1, T,
      +                                           AAQQ )       $                                           AAQQ )
                                     SVA( q ) = T*DSQRT( AAQQ )*WORK( q )                                      SVA( q ) = T*DSQRT( AAQQ )*WORK( q )
                                  END IF                                   END IF
                               END IF                                END IF
                               IF( ( AAPP / AAPP0 ).LE.ROOTEPS ) THEN                                IF( ( AAPP / AAPP0 ).LE.ROOTEPS ) THEN
                                  IF( ( AAPP.LT.ROOTBIG ) .AND.                                   IF( ( AAPP.LT.ROOTBIG ) .AND.
      +                               ( AAPP.GT.ROOTSFMIN ) ) THEN       $                               ( AAPP.GT.ROOTSFMIN ) ) THEN
                                     AAPP = DNRM2( M, A( 1, p ), 1 )*                                      AAPP = DNRM2( M, A( 1, p ), 1 )*
      +                                     WORK( p )       $                                     WORK( p )
                                  ELSE                                   ELSE
                                     T = ZERO                                      T = ZERO
                                     AAPP = ZERO                                      AAPP = ONE
                                     CALL DLASSQ( M, A( 1, p ), 1, T,                                      CALL DLASSQ( M, A( 1, p ), 1, T,
      +                                           AAPP )       $                                           AAPP )
                                     AAPP = T*DSQRT( AAPP )*WORK( p )                                      AAPP = T*DSQRT( AAPP )*WORK( p )
                                  END IF                                   END IF
                                  SVA( p ) = AAPP                                   SVA( p ) = AAPP
Line 1023 Line 1119
                         END IF                          END IF
 *  *
                         IF( ( i.LE.SWBAND ) .AND.                          IF( ( i.LE.SWBAND ) .AND.
      +                      ( PSKIPPED.GT.ROWSKIP ) ) THEN       $                      ( PSKIPPED.GT.ROWSKIP ) ) THEN
                            IF( ir1.EQ.0 )AAPP = -AAPP                             IF( ir1.EQ.0 )AAPP = -AAPP
                            NOTROT = 0                             NOTROT = 0
                            GO TO 2103                             GO TO 2103
Line 1040 Line 1136
                   ELSE                    ELSE
                      SVA( p ) = AAPP                       SVA( p ) = AAPP
                      IF( ( ir1.EQ.0 ) .AND. ( AAPP.EQ.ZERO ) )                       IF( ( ir1.EQ.0 ) .AND. ( AAPP.EQ.ZERO ) )
      +                   NOTROT = NOTROT + MIN0( igl+KBL-1, N ) - p       $                   NOTROT = NOTROT + MIN0( igl+KBL-1, N ) - p
                   END IF                    END IF
 *  *
  2001          CONTINUE   2001          CONTINUE
Line 1085 Line 1181
                               END IF                                END IF
                               IF( AAPP.LT.( BIG / AAQQ ) ) THEN                                IF( AAPP.LT.( BIG / AAQQ ) ) THEN
                                  AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,                                   AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
      +                                  q ), 1 )*WORK( p )*WORK( q ) /       $                                  q ), 1 )*WORK( p )*WORK( q ) /
      +                                  AAQQ ) / AAPP       $                                  AAQQ ) / AAPP
                               ELSE                                ELSE
                                  CALL DCOPY( M, A( 1, p ), 1,                                   CALL DCOPY( M, A( 1, p ), 1,
      +                                       WORK( N+1 ), 1 )       $                                       WORK( N+1 ), 1 )
                                  CALL DLASCL( 'G', 0, 0, AAPP,                                   CALL DLASCL( 'G', 0, 0, AAPP,
      +                                        WORK( p ), M, 1,       $                                        WORK( p ), M, 1,
      +                                        WORK( N+1 ), LDA, IERR )       $                                        WORK( N+1 ), LDA, IERR )
                                  AAPQ = DDOT( M, WORK( N+1 ), 1,                                   AAPQ = DDOT( M, WORK( N+1 ), 1,
      +                                  A( 1, q ), 1 )*WORK( q ) / AAQQ       $                                  A( 1, q ), 1 )*WORK( q ) / AAQQ
                               END IF                                END IF
                            ELSE                             ELSE
                               IF( AAPP.GE.AAQQ ) THEN                                IF( AAPP.GE.AAQQ ) THEN
Line 1104 Line 1200
                               END IF                                END IF
                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN                                IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
                                  AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,                                   AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
      +                                  q ), 1 )*WORK( p )*WORK( q ) /       $                                  q ), 1 )*WORK( p )*WORK( q ) /
      +                                  AAQQ ) / AAPP       $                                  AAQQ ) / AAPP
                               ELSE                                ELSE
                                  CALL DCOPY( M, A( 1, q ), 1,                                   CALL DCOPY( M, A( 1, q ), 1,
      +                                       WORK( N+1 ), 1 )       $                                       WORK( N+1 ), 1 )
                                  CALL DLASCL( 'G', 0, 0, AAQQ,                                   CALL DLASCL( 'G', 0, 0, AAQQ,
      +                                        WORK( q ), M, 1,       $                                        WORK( q ), M, 1,
      +                                        WORK( N+1 ), LDA, IERR )       $                                        WORK( N+1 ), LDA, IERR )
                                  AAPQ = DDOT( M, WORK( N+1 ), 1,                                   AAPQ = DDOT( M, WORK( N+1 ), 1,
      +                                  A( 1, p ), 1 )*WORK( p ) / AAPP       $                                  A( 1, p ), 1 )*WORK( p ) / AAPP
                               END IF                                END IF
                            END IF                             END IF
 *  *
Line 1131 Line 1227
 *  *
                                  AQOAP = AAQQ / AAPP                                   AQOAP = AAQQ / AAPP
                                  APOAQ = AAPP / AAQQ                                   APOAQ = AAPP / AAQQ
                                  THETA = -HALF*DABS( AQOAP-APOAQ ) /                                   THETA = -HALF*DABS(AQOAP-APOAQ)/AAPQ
      +                                   AAPQ  
                                  IF( AAQQ.GT.AAPP0 )THETA = -THETA                                   IF( AAQQ.GT.AAPP0 )THETA = -THETA
 *  *
                                  IF( DABS( THETA ).GT.BIGTHETA ) THEN                                   IF( DABS( THETA ).GT.BIGTHETA ) THEN
                                     T = HALF / THETA                                      T = HALF / THETA
                                     FASTR( 3 ) = T*WORK( p ) / WORK( q )                                      FASTR( 3 ) = T*WORK( p ) / WORK( q )
                                     FASTR( 4 ) = -T*WORK( q ) /                                      FASTR( 4 ) = -T*WORK( q ) /
      +                                           WORK( p )       $                                           WORK( p )
                                     CALL DROTM( M, A( 1, p ), 1,                                      CALL DROTM( M, A( 1, p ), 1,
      +                                          A( 1, q ), 1, FASTR )       $                                          A( 1, q ), 1, FASTR )
                                     IF( RSVEC )CALL DROTM( MVL,                                      IF( RSVEC )CALL DROTM( MVL,
      +                                              V( 1, p ), 1,       $                                              V( 1, p ), 1,
      +                                              V( 1, q ), 1,       $                                              V( 1, q ), 1,
      +                                              FASTR )       $                                              FASTR )
                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,                                      SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
      +                                         ONE+T*APOAQ*AAPQ ) )       $                                         ONE+T*APOAQ*AAPQ ) )
                                     AAPP = AAPP*DSQRT( DMAX1( ZERO,                                      AAPP = AAPP*DSQRT( DMAX1( ZERO,
      +                                     ONE-T*AQOAP*AAPQ ) )       $                                     ONE-T*AQOAP*AAPQ ) )
                                     MXSINJ = DMAX1( MXSINJ, DABS( T ) )                                      MXSINJ = DMAX1( MXSINJ, DABS( T ) )
                                  ELSE                                   ELSE
 *  *
Line 1158 Line 1253
                                     THSIGN = -DSIGN( ONE, AAPQ )                                      THSIGN = -DSIGN( ONE, AAPQ )
                                     IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN                                      IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN
                                     T = ONE / ( THETA+THSIGN*                                      T = ONE / ( THETA+THSIGN*
      +                                  DSQRT( ONE+THETA*THETA ) )       $                                  DSQRT( ONE+THETA*THETA ) )
                                     CS = DSQRT( ONE / ( ONE+T*T ) )                                      CS = DSQRT( ONE / ( ONE+T*T ) )
                                     SN = T*CS                                      SN = T*CS
                                     MXSINJ = DMAX1( MXSINJ, DABS( SN ) )                                      MXSINJ = DMAX1( MXSINJ, DABS( SN ) )
                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,                                      SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
      +                                         ONE+T*APOAQ*AAPQ ) )       $                                         ONE+T*APOAQ*AAPQ ) )
                                     AAPP = AAPP*DSQRT( ONE-T*AQOAP*                                      AAPP = AAPP*DSQRT( DMAX1( ZERO, 
      +                                     AAPQ )       $                                     ONE-T*AQOAP*AAPQ ) )
 *  *
                                     APOAQ = WORK( p ) / WORK( q )                                      APOAQ = WORK( p ) / WORK( q )
                                     AQOAP = WORK( q ) / WORK( p )                                      AQOAP = WORK( q ) / WORK( p )
Line 1177 Line 1272
                                           WORK( p ) = WORK( p )*CS                                            WORK( p ) = WORK( p )*CS
                                           WORK( q ) = WORK( q )*CS                                            WORK( q ) = WORK( q )*CS
                                           CALL DROTM( M, A( 1, p ), 1,                                            CALL DROTM( M, A( 1, p ), 1,
      +                                                A( 1, q ), 1,       $                                                A( 1, q ), 1,
      +                                                FASTR )       $                                                FASTR )
                                           IF( RSVEC )CALL DROTM( MVL,                                            IF( RSVEC )CALL DROTM( MVL,
      +                                        V( 1, p ), 1, V( 1, q ),       $                                        V( 1, p ), 1, V( 1, q ),
      +                                        1, FASTR )       $                                        1, FASTR )
                                        ELSE                                         ELSE
                                           CALL DAXPY( M, -T*AQOAP,                                            CALL DAXPY( M, -T*AQOAP,
      +                                                A( 1, q ), 1,       $                                                A( 1, q ), 1,
      +                                                A( 1, p ), 1 )       $                                                A( 1, p ), 1 )
                                           CALL DAXPY( M, CS*SN*APOAQ,                                            CALL DAXPY( M, CS*SN*APOAQ,
      +                                                A( 1, p ), 1,       $                                                A( 1, p ), 1,
      +                                                A( 1, q ), 1 )       $                                                A( 1, q ), 1 )
                                           IF( RSVEC ) THEN                                            IF( RSVEC ) THEN
                                              CALL DAXPY( MVL, -T*AQOAP,                                               CALL DAXPY( MVL, -T*AQOAP,
      +                                                   V( 1, q ), 1,       $                                                   V( 1, q ), 1,
      +                                                   V( 1, p ), 1 )       $                                                   V( 1, p ), 1 )
                                              CALL DAXPY( MVL,                                               CALL DAXPY( MVL,
      +                                                   CS*SN*APOAQ,       $                                                   CS*SN*APOAQ,
      +                                                   V( 1, p ), 1,       $                                                   V( 1, p ), 1,
      +                                                   V( 1, q ), 1 )       $                                                   V( 1, q ), 1 )
                                           END IF                                            END IF
                                           WORK( p ) = WORK( p )*CS                                            WORK( p ) = WORK( p )*CS
                                           WORK( q ) = WORK( q ) / CS                                            WORK( q ) = WORK( q ) / CS
Line 1204 Line 1299
                                     ELSE                                      ELSE
                                        IF( WORK( q ).GE.ONE ) THEN                                         IF( WORK( q ).GE.ONE ) THEN
                                           CALL DAXPY( M, T*APOAQ,                                            CALL DAXPY( M, T*APOAQ,
      +                                                A( 1, p ), 1,       $                                                A( 1, p ), 1,
      +                                                A( 1, q ), 1 )       $                                                A( 1, q ), 1 )
                                           CALL DAXPY( M, -CS*SN*AQOAP,                                            CALL DAXPY( M, -CS*SN*AQOAP,
      +                                                A( 1, q ), 1,       $                                                A( 1, q ), 1,
      +                                                A( 1, p ), 1 )       $                                                A( 1, p ), 1 )
                                           IF( RSVEC ) THEN                                            IF( RSVEC ) THEN
                                              CALL DAXPY( MVL, T*APOAQ,                                               CALL DAXPY( MVL, T*APOAQ,
      +                                                   V( 1, p ), 1,       $                                                   V( 1, p ), 1,
      +                                                   V( 1, q ), 1 )       $                                                   V( 1, q ), 1 )
                                              CALL DAXPY( MVL,                                               CALL DAXPY( MVL,
      +                                                   -CS*SN*AQOAP,       $                                                   -CS*SN*AQOAP,
      +                                                   V( 1, q ), 1,       $                                                   V( 1, q ), 1,
      +                                                   V( 1, p ), 1 )       $                                                   V( 1, p ), 1 )
                                           END IF                                            END IF
                                           WORK( p ) = WORK( p ) / CS                                            WORK( p ) = WORK( p ) / CS
                                           WORK( q ) = WORK( q )*CS                                            WORK( q ) = WORK( q )*CS
                                        ELSE                                         ELSE
                                           IF( WORK( p ).GE.WORK( q ) )                                            IF( WORK( p ).GE.WORK( q ) )
      +                                        THEN       $                                        THEN
                                              CALL DAXPY( M, -T*AQOAP,                                               CALL DAXPY( M, -T*AQOAP,
      +                                                   A( 1, q ), 1,       $                                                   A( 1, q ), 1,
      +                                                   A( 1, p ), 1 )       $                                                   A( 1, p ), 1 )
                                              CALL DAXPY( M, CS*SN*APOAQ,                                               CALL DAXPY( M, CS*SN*APOAQ,
      +                                                   A( 1, p ), 1,       $                                                   A( 1, p ), 1,
      +                                                   A( 1, q ), 1 )       $                                                   A( 1, q ), 1 )
                                              WORK( p ) = WORK( p )*CS                                               WORK( p ) = WORK( p )*CS
                                              WORK( q ) = WORK( q ) / CS                                               WORK( q ) = WORK( q ) / CS
                                              IF( RSVEC ) THEN                                               IF( RSVEC ) THEN
                                                 CALL DAXPY( MVL,                                                  CALL DAXPY( MVL,
      +                                               -T*AQOAP,       $                                               -T*AQOAP,
      +                                               V( 1, q ), 1,       $                                               V( 1, q ), 1,
      +                                               V( 1, p ), 1 )       $                                               V( 1, p ), 1 )
                                                 CALL DAXPY( MVL,                                                  CALL DAXPY( MVL,
      +                                               CS*SN*APOAQ,       $                                               CS*SN*APOAQ,
      +                                               V( 1, p ), 1,       $                                               V( 1, p ), 1,
      +                                               V( 1, q ), 1 )       $                                               V( 1, q ), 1 )
                                              END IF                                               END IF
                                           ELSE                                            ELSE
                                              CALL DAXPY( M, T*APOAQ,                                               CALL DAXPY( M, T*APOAQ,
      +                                                   A( 1, p ), 1,       $                                                   A( 1, p ), 1,
      +                                                   A( 1, q ), 1 )       $                                                   A( 1, q ), 1 )
                                              CALL DAXPY( M,                                               CALL DAXPY( M,
      +                                                   -CS*SN*AQOAP,       $                                                   -CS*SN*AQOAP,
      +                                                   A( 1, q ), 1,       $                                                   A( 1, q ), 1,
      +                                                   A( 1, p ), 1 )       $                                                   A( 1, p ), 1 )
                                              WORK( p ) = WORK( p ) / CS                                               WORK( p ) = WORK( p ) / CS
                                              WORK( q ) = WORK( q )*CS                                               WORK( q ) = WORK( q )*CS
                                              IF( RSVEC ) THEN                                               IF( RSVEC ) THEN
                                                 CALL DAXPY( MVL,                                                  CALL DAXPY( MVL,
      +                                               T*APOAQ, V( 1, p ),       $                                               T*APOAQ, V( 1, p ),
      +                                               1, V( 1, q ), 1 )       $                                               1, V( 1, q ), 1 )
                                                 CALL DAXPY( MVL,                                                  CALL DAXPY( MVL,
      +                                               -CS*SN*AQOAP,       $                                               -CS*SN*AQOAP,
      +                                               V( 1, q ), 1,       $                                               V( 1, q ), 1,
      +                                               V( 1, p ), 1 )       $                                               V( 1, p ), 1 )
                                              END IF                                               END IF
                                           END IF                                            END IF
                                        END IF                                         END IF
Line 1268 Line 1363
                               ELSE                                ELSE
                                  IF( AAPP.GT.AAQQ ) THEN                                   IF( AAPP.GT.AAQQ ) THEN
                                     CALL DCOPY( M, A( 1, p ), 1,                                      CALL DCOPY( M, A( 1, p ), 1,
      +                                          WORK( N+1 ), 1 )       $                                          WORK( N+1 ), 1 )
                                     CALL DLASCL( 'G', 0, 0, AAPP, ONE,                                      CALL DLASCL( 'G', 0, 0, AAPP, ONE,
      +                                           M, 1, WORK( N+1 ), LDA,       $                                           M, 1, WORK( N+1 ), LDA,
      +                                           IERR )       $                                           IERR )
                                     CALL DLASCL( 'G', 0, 0, AAQQ, ONE,                                      CALL DLASCL( 'G', 0, 0, AAQQ, ONE,
      +                                           M, 1, A( 1, q ), LDA,       $                                           M, 1, A( 1, q ), LDA,
      +                                           IERR )       $                                           IERR )
                                     TEMP1 = -AAPQ*WORK( p ) / WORK( q )                                      TEMP1 = -AAPQ*WORK( p ) / WORK( q )
                                     CALL DAXPY( M, TEMP1, WORK( N+1 ),                                      CALL DAXPY( M, TEMP1, WORK( N+1 ),
      +                                          1, A( 1, q ), 1 )       $                                          1, A( 1, q ), 1 )
                                     CALL DLASCL( 'G', 0, 0, ONE, AAQQ,                                      CALL DLASCL( 'G', 0, 0, ONE, AAQQ,
      +                                           M, 1, A( 1, q ), LDA,       $                                           M, 1, A( 1, q ), LDA,
      +                                           IERR )       $                                           IERR )
                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,                                      SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
      +                                         ONE-AAPQ*AAPQ ) )       $                                         ONE-AAPQ*AAPQ ) )
                                     MXSINJ = DMAX1( MXSINJ, SFMIN )                                      MXSINJ = DMAX1( MXSINJ, SFMIN )
                                  ELSE                                   ELSE
                                     CALL DCOPY( M, A( 1, q ), 1,                                      CALL DCOPY( M, A( 1, q ), 1,
      +                                          WORK( N+1 ), 1 )       $                                          WORK( N+1 ), 1 )
                                     CALL DLASCL( 'G', 0, 0, AAQQ, ONE,                                      CALL DLASCL( 'G', 0, 0, AAQQ, ONE,
      +                                           M, 1, WORK( N+1 ), LDA,       $                                           M, 1, WORK( N+1 ), LDA,
      +                                           IERR )       $                                           IERR )
                                     CALL DLASCL( 'G', 0, 0, AAPP, ONE,                                      CALL DLASCL( 'G', 0, 0, AAPP, ONE,
      +                                           M, 1, A( 1, p ), LDA,       $                                           M, 1, A( 1, p ), LDA,
      +                                           IERR )       $                                           IERR )
                                     TEMP1 = -AAPQ*WORK( q ) / WORK( p )                                      TEMP1 = -AAPQ*WORK( q ) / WORK( p )
                                     CALL DAXPY( M, TEMP1, WORK( N+1 ),                                      CALL DAXPY( M, TEMP1, WORK( N+1 ),
      +                                          1, A( 1, p ), 1 )       $                                          1, A( 1, p ), 1 )
                                     CALL DLASCL( 'G', 0, 0, ONE, AAPP,                                      CALL DLASCL( 'G', 0, 0, ONE, AAPP,
      +                                           M, 1, A( 1, p ), LDA,       $                                           M, 1, A( 1, p ), LDA,
      +                                           IERR )       $                                           IERR )
                                     SVA( p ) = AAPP*DSQRT( DMAX1( ZERO,                                      SVA( p ) = AAPP*DSQRT( DMAX1( ZERO,
      +                                         ONE-AAPQ*AAPQ ) )       $                                         ONE-AAPQ*AAPQ ) )
                                     MXSINJ = DMAX1( MXSINJ, SFMIN )                                      MXSINJ = DMAX1( MXSINJ, SFMIN )
                                  END IF                                   END IF
                               END IF                                END IF
Line 1309 Line 1404
 *           In the case of cancellation in updating SVA(q)  *           In the case of cancellation in updating SVA(q)
 *           .. recompute SVA(q)  *           .. recompute SVA(q)
                               IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )                                IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
      +                            THEN       $                            THEN
                                  IF( ( AAQQ.LT.ROOTBIG ) .AND.                                   IF( ( AAQQ.LT.ROOTBIG ) .AND.
      +                               ( AAQQ.GT.ROOTSFMIN ) ) THEN       $                               ( AAQQ.GT.ROOTSFMIN ) ) THEN
                                     SVA( q ) = DNRM2( M, A( 1, q ), 1 )*                                      SVA( q ) = DNRM2( M, A( 1, q ), 1 )*
      +                                         WORK( q )       $                                         WORK( q )
                                  ELSE                                   ELSE
                                     T = ZERO                                      T = ZERO
                                     AAQQ = ZERO                                      AAQQ = ONE
                                     CALL DLASSQ( M, A( 1, q ), 1, T,                                      CALL DLASSQ( M, A( 1, q ), 1, T,
      +                                           AAQQ )       $                                           AAQQ )
                                     SVA( q ) = T*DSQRT( AAQQ )*WORK( q )                                      SVA( q ) = T*DSQRT( AAQQ )*WORK( q )
                                  END IF                                   END IF
                               END IF                                END IF
                               IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN                                IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN
                                  IF( ( AAPP.LT.ROOTBIG ) .AND.                                   IF( ( AAPP.LT.ROOTBIG ) .AND.
      +                               ( AAPP.GT.ROOTSFMIN ) ) THEN       $                               ( AAPP.GT.ROOTSFMIN ) ) THEN
                                     AAPP = DNRM2( M, A( 1, p ), 1 )*                                      AAPP = DNRM2( M, A( 1, p ), 1 )*
      +                                     WORK( p )       $                                     WORK( p )
                                  ELSE                                   ELSE
                                     T = ZERO                                      T = ZERO
                                     AAPP = ZERO                                      AAPP = ONE
                                     CALL DLASSQ( M, A( 1, p ), 1, T,                                      CALL DLASSQ( M, A( 1, p ), 1, T,
      +                                           AAPP )       $                                           AAPP )
                                     AAPP = T*DSQRT( AAPP )*WORK( p )                                      AAPP = T*DSQRT( AAPP )*WORK( p )
                                  END IF                                   END IF
                                  SVA( p ) = AAPP                                   SVA( p ) = AAPP
Line 1350 Line 1445
                         END IF                          END IF
 *  *
                         IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) )                          IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) )
      +                      THEN       $                      THEN
                            SVA( p ) = AAPP                             SVA( p ) = AAPP
                            NOTROT = 0                             NOTROT = 0
                            GO TO 2011                             GO TO 2011
                         END IF                          END IF
                         IF( ( i.LE.SWBAND ) .AND.                          IF( ( i.LE.SWBAND ) .AND.
      +                      ( PSKIPPED.GT.ROWSKIP ) ) THEN       $                      ( PSKIPPED.GT.ROWSKIP ) ) THEN
                            AAPP = -AAPP                             AAPP = -AAPP
                            NOTROT = 0                             NOTROT = 0
                            GO TO 2203                             GO TO 2203
Line 1371 Line 1466
                   ELSE                    ELSE
 *  *
                      IF( AAPP.EQ.ZERO )NOTROT = NOTROT +                       IF( AAPP.EQ.ZERO )NOTROT = NOTROT +
      +                   MIN0( jgl+KBL-1, N ) - jgl + 1       $                   MIN0( jgl+KBL-1, N ) - jgl + 1
                      IF( AAPP.LT.ZERO )NOTROT = 0                       IF( AAPP.LT.ZERO )NOTROT = 0
 *  *
                   END IF                    END IF
Line 1391 Line 1486
 *  *
 *     .. update SVA(N)  *     .. update SVA(N)
          IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) )           IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) )
      +       THEN       $       THEN
             SVA( N ) = DNRM2( M, A( 1, N ), 1 )*WORK( N )              SVA( N ) = DNRM2( M, A( 1, N ), 1 )*WORK( N )
          ELSE           ELSE
             T = ZERO              T = ZERO
             AAPP = ZERO              AAPP = ONE
             CALL DLASSQ( M, A( 1, N ), 1, T, AAPP )              CALL DLASSQ( M, A( 1, N ), 1, T, AAPP )
             SVA( N ) = T*DSQRT( AAPP )*WORK( N )              SVA( N ) = T*DSQRT( AAPP )*WORK( N )
          END IF           END IF
Line 1403 Line 1498
 *     Additional steering devices  *     Additional steering devices
 *  *
          IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.           IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.
      +       ( ISWROT.LE.N ) ) )SWBAND = i       $       ( ISWROT.LE.N ) ) )SWBAND = i
 *  *
          IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.DSQRT( DBLE( N ) )*           IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.DSQRT( DBLE( N ) )*
      +       TOL ) .AND. ( DBLE( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN       $       TOL ) .AND. ( DBLE( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
             GO TO 1994              GO TO 1994
          END IF           END IF
 *  *
Line 1446 Line 1541
          END IF           END IF
          IF( SVA( p ).NE.ZERO ) THEN           IF( SVA( p ).NE.ZERO ) THEN
             N4 = N4 + 1              N4 = N4 + 1
             IF( SVA( p )*SCALE.GT.SFMIN )N2 = N2 + 1              IF( SVA( p )*SKL.GT.SFMIN )N2 = N2 + 1
          END IF           END IF
  5991 CONTINUE   5991 CONTINUE
       IF( SVA( N ).NE.ZERO ) THEN        IF( SVA( N ).NE.ZERO ) THEN
          N4 = N4 + 1           N4 = N4 + 1
          IF( SVA( N )*SCALE.GT.SFMIN )N2 = N2 + 1           IF( SVA( N )*SKL.GT.SFMIN )N2 = N2 + 1
       END IF        END IF
 *  *
 *     Normalize the left singular vectors.  *     Normalize the left singular vectors.
Line 1478 Line 1573
       END IF        END IF
 *  *
 *     Undo scaling, if necessary (and possible).  *     Undo scaling, if necessary (and possible).
       IF( ( ( SCALE.GT.ONE ) .AND. ( SVA( 1 ).LT.( BIG /        IF( ( ( SKL.GT.ONE ) .AND. ( SVA( 1 ).LT.( BIG / SKL) ) )
      +    SCALE ) ) ) .OR. ( ( SCALE.LT.ONE ) .AND. ( SVA( N2 ).GT.       $    .OR. ( ( SKL.LT.ONE ) .AND. ( SVA( MAX( N2, 1 ) ) .GT.
      +    ( SFMIN / SCALE ) ) ) ) THEN       $    ( SFMIN / SKL) ) ) ) THEN
          DO 2400 p = 1, N           DO 2400 p = 1, N
             SVA( p ) = SCALE*SVA( p )              SVA( P ) = SKL*SVA( P )
  2400    CONTINUE   2400    CONTINUE
          SCALE = ONE           SKL= ONE
       END IF        END IF
 *  *
       WORK( 1 ) = SCALE        WORK( 1 ) = SKL
 *     The singular values of A are SCALE*SVA(1:N). If SCALE.NE.ONE  *     The singular values of A are SKL*SVA(1:N). If SKL.NE.ONE
 *     then some of the singular values may overflow or underflow and  *     then some of the singular values may overflow or underflow and
 *     the spectrum is given in this factored representation.  *     the spectrum is given in this factored representation.
 *  *
Removed from v.1.1  
changed lines
  Added in v.1.11

CVSweb interface <joel.bertrand@systella.fr>