Annotation of rpl/lapack/lapack/dgesvj.f, revision 1.7

1.7     ! bertrand    1: *> \brief \b DGESVJ
        !             2: *
        !             3: *  =========== DOCUMENTATION ===========
        !             4: *
        !             5: * Online html documentation available at 
        !             6: *            http://www.netlib.org/lapack/explore-html/ 
        !             7: *
        !             8: *> \htmlonly
        !             9: *> Download DGESVJ + dependencies 
        !            10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgesvj.f"> 
        !            11: *> [TGZ]</a> 
        !            12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgesvj.f"> 
        !            13: *> [ZIP]</a> 
        !            14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgesvj.f"> 
        !            15: *> [TXT]</a>
        !            16: *> \endhtmlonly 
        !            17: *
        !            18: *  Definition:
        !            19: *  ===========
        !            20: *
        !            21: *       SUBROUTINE DGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
        !            22: *                          LDV, WORK, LWORK, INFO )
        !            23: * 
        !            24: *       .. Scalar Arguments ..
        !            25: *       INTEGER            INFO, LDA, LDV, LWORK, M, MV, N
        !            26: *       CHARACTER*1        JOBA, JOBU, JOBV
        !            27: *       ..
        !            28: *       .. Array Arguments ..
        !            29: *       DOUBLE PRECISION   A( LDA, * ), SVA( N ), V( LDV, * ),
        !            30: *      $                   WORK( LWORK )
        !            31: *       ..
        !            32: *  
        !            33: *
        !            34: *> \par Purpose:
        !            35: *  =============
        !            36: *>
        !            37: *> \verbatim
        !            38: *>
        !            39: *> DGESVJ computes the singular value decomposition (SVD) of a real
        !            40: *> M-by-N matrix A, where M >= N. The SVD of A is written as
        !            41: *>                                    [++]   [xx]   [x0]   [xx]
        !            42: *>              A = U * SIGMA * V^t,  [++] = [xx] * [ox] * [xx]
        !            43: *>                                    [++]   [xx]
        !            44: *> where SIGMA is an N-by-N diagonal matrix, U is an M-by-N orthonormal
        !            45: *> matrix, and V is an N-by-N orthogonal matrix. The diagonal elements
        !            46: *> of SIGMA are the singular values of A. The columns of U and V are the
        !            47: *> left and the right singular vectors of A, respectively.
        !            48: *> \endverbatim
        !            49: *
        !            50: *  Arguments:
        !            51: *  ==========
        !            52: *
        !            53: *> \param[in] JOBA
        !            54: *> \verbatim
        !            55: *>          JOBA is CHARACTER* 1
        !            56: *>          Specifies the structure of A.
        !            57: *>          = 'L': The input matrix A is lower triangular;
        !            58: *>          = 'U': The input matrix A is upper triangular;
        !            59: *>          = 'G': The input matrix A is general M-by-N matrix, M >= N.
        !            60: *> \endverbatim
        !            61: *>
        !            62: *> \param[in] JOBU
        !            63: *> \verbatim
        !            64: *>          JOBU is CHARACTER*1
        !            65: *>          Specifies whether to compute the left singular vectors
        !            66: *>          (columns of U):
        !            67: *>          = 'U': The left singular vectors corresponding to the nonzero
        !            68: *>                 singular values are computed and returned in the leading
        !            69: *>                 columns of A. See more details in the description of A.
        !            70: *>                 The default numerical orthogonality threshold is set to
        !            71: *>                 approximately TOL=CTOL*EPS, CTOL=DSQRT(M), EPS=DLAMCH('E').
        !            72: *>          = 'C': Analogous to JOBU='U', except that user can control the
        !            73: *>                 level of numerical orthogonality of the computed left
        !            74: *>                 singular vectors. TOL can be set to TOL = CTOL*EPS, where
        !            75: *>                 CTOL is given on input in the array WORK.
        !            76: *>                 No CTOL smaller than ONE is allowed. CTOL greater
        !            77: *>                 than 1 / EPS is meaningless. The option 'C'
        !            78: *>                 can be used if M*EPS is satisfactory orthogonality
        !            79: *>                 of the computed left singular vectors, so CTOL=M could
        !            80: *>                 save few sweeps of Jacobi rotations.
        !            81: *>                 See the descriptions of A and WORK(1).
        !            82: *>          = 'N': The matrix U is not computed. However, see the
        !            83: *>                 description of A.
        !            84: *> \endverbatim
        !            85: *>
        !            86: *> \param[in] JOBV
        !            87: *> \verbatim
        !            88: *>          JOBV is CHARACTER*1
        !            89: *>          Specifies whether to compute the right singular vectors, that
        !            90: *>          is, the matrix V:
        !            91: *>          = 'V' : the matrix V is computed and returned in the array V
        !            92: *>          = 'A' : the Jacobi rotations are applied to the MV-by-N
        !            93: *>                  array V. In other words, the right singular vector
        !            94: *>                  matrix V is not computed explicitly, instead it is
        !            95: *>                  applied to an MV-by-N matrix initially stored in the
        !            96: *>                  first MV rows of V.
        !            97: *>          = 'N' : the matrix V is not computed and the array V is not
        !            98: *>                  referenced
        !            99: *> \endverbatim
        !           100: *>
        !           101: *> \param[in] M
        !           102: *> \verbatim
        !           103: *>          M is INTEGER
        !           104: *>          The number of rows of the input matrix A. 1/DLAMCH('E') > M >= 0.  
        !           105: *> \endverbatim
        !           106: *>
        !           107: *> \param[in] N
        !           108: *> \verbatim
        !           109: *>          N is INTEGER
        !           110: *>          The number of columns of the input matrix A.
        !           111: *>          M >= N >= 0.
        !           112: *> \endverbatim
        !           113: *>
        !           114: *> \param[in,out] A
        !           115: *> \verbatim
        !           116: *>          A is DOUBLE PRECISION array, dimension (LDA,N)
        !           117: *>          On entry, the M-by-N matrix A.
        !           118: *>          On exit :
        !           119: *>          If JOBU .EQ. 'U' .OR. JOBU .EQ. 'C' :
        !           120: *>                 If INFO .EQ. 0 :
        !           121: *>                 RANKA orthonormal columns of U are returned in the
        !           122: *>                 leading RANKA columns of the array A. Here RANKA <= N
        !           123: *>                 is the number of computed singular values of A that are
        !           124: *>                 above the underflow threshold DLAMCH('S'). The singular
        !           125: *>                 vectors corresponding to underflowed or zero singular
        !           126: *>                 values are not computed. The value of RANKA is returned
        !           127: *>                 in the array WORK as RANKA=NINT(WORK(2)). Also see the
        !           128: *>                 descriptions of SVA and WORK. The computed columns of U
        !           129: *>                 are mutually numerically orthogonal up to approximately
        !           130: *>                 TOL=DSQRT(M)*EPS (default); or TOL=CTOL*EPS (JOBU.EQ.'C'),
        !           131: *>                 see the description of JOBU.
        !           132: *>                 If INFO .GT. 0 :
        !           133: *>                 the procedure DGESVJ did not converge in the given number
        !           134: *>                 of iterations (sweeps). In that case, the computed
        !           135: *>                 columns of U may not be orthogonal up to TOL. The output
        !           136: *>                 U (stored in A), SIGMA (given by the computed singular
        !           137: *>                 values in SVA(1:N)) and V is still a decomposition of the
        !           138: *>                 input matrix A in the sense that the residual
        !           139: *>                 ||A-SCALE*U*SIGMA*V^T||_2 / ||A||_2 is small.
        !           140: *>
        !           141: *>          If JOBU .EQ. 'N' :
        !           142: *>                 If INFO .EQ. 0 :
        !           143: *>                 Note that the left singular vectors are 'for free' in the
        !           144: *>                 one-sided Jacobi SVD algorithm. However, if only the
        !           145: *>                 singular values are needed, the level of numerical
        !           146: *>                 orthogonality of U is not an issue and iterations are
        !           147: *>                 stopped when the columns of the iterated matrix are
        !           148: *>                 numerically orthogonal up to approximately M*EPS. Thus,
        !           149: *>                 on exit, A contains the columns of U scaled with the
        !           150: *>                 corresponding singular values.
        !           151: *>                 If INFO .GT. 0 :
        !           152: *>                 the procedure DGESVJ did not converge in the given number
        !           153: *>                 of iterations (sweeps).
        !           154: *> \endverbatim
        !           155: *>
        !           156: *> \param[in] LDA
        !           157: *> \verbatim
        !           158: *>          LDA is INTEGER
        !           159: *>          The leading dimension of the array A.  LDA >= max(1,M).
        !           160: *> \endverbatim
        !           161: *>
        !           162: *> \param[out] SVA
        !           163: *> \verbatim
        !           164: *>          SVA is DOUBLE PRECISION array, dimension (N)
        !           165: *>          On exit :
        !           166: *>          If INFO .EQ. 0 :
        !           167: *>          depending on the value SCALE = WORK(1), we have:
        !           168: *>                 If SCALE .EQ. ONE :
        !           169: *>                 SVA(1:N) contains the computed singular values of A.
        !           170: *>                 During the computation SVA contains the Euclidean column
        !           171: *>                 norms of the iterated matrices in the array A.
        !           172: *>                 If SCALE .NE. ONE :
        !           173: *>                 The singular values of A are SCALE*SVA(1:N), and this
        !           174: *>                 factored representation is due to the fact that some of the
        !           175: *>                 singular values of A might underflow or overflow.
        !           176: *>          If INFO .GT. 0 :
        !           177: *>          the procedure DGESVJ did not converge in the given number of
        !           178: *>          iterations (sweeps) and SCALE*SVA(1:N) may not be accurate.
        !           179: *> \endverbatim
        !           180: *>
        !           181: *> \param[in] MV
        !           182: *> \verbatim
        !           183: *>          MV is INTEGER
        !           184: *>          If JOBV .EQ. 'A', then the product of Jacobi rotations in DGESVJ
        !           185: *>          is applied to the first MV rows of V. See the description of JOBV.
        !           186: *> \endverbatim
        !           187: *>
        !           188: *> \param[in,out] V
        !           189: *> \verbatim
        !           190: *>          V is DOUBLE PRECISION array, dimension (LDV,N)
        !           191: *>          If JOBV = 'V', then V contains on exit the N-by-N matrix of
        !           192: *>                         the right singular vectors;
        !           193: *>          If JOBV = 'A', then V contains the product of the computed right
        !           194: *>                         singular vector matrix and the initial matrix in
        !           195: *>                         the array V.
        !           196: *>          If JOBV = 'N', then V is not referenced.
        !           197: *> \endverbatim
        !           198: *>
        !           199: *> \param[in] LDV
        !           200: *> \verbatim
        !           201: *>          LDV is INTEGER
        !           202: *>          The leading dimension of the array V, LDV .GE. 1.
        !           203: *>          If JOBV .EQ. 'V', then LDV .GE. max(1,N).
        !           204: *>          If JOBV .EQ. 'A', then LDV .GE. max(1,MV) .
        !           205: *> \endverbatim
        !           206: *>
        !           207: *> \param[in,out] WORK
        !           208: *> \verbatim
        !           209: *>          WORK is DOUBLE PRECISION array, dimension max(4,M+N).
        !           210: *>          On entry :
        !           211: *>          If JOBU .EQ. 'C' :
        !           212: *>          WORK(1) = CTOL, where CTOL defines the threshold for convergence.
        !           213: *>                    The process stops if all columns of A are mutually
        !           214: *>                    orthogonal up to CTOL*EPS, EPS=DLAMCH('E').
        !           215: *>                    It is required that CTOL >= ONE, i.e. it is not
        !           216: *>                    allowed to force the routine to obtain orthogonality
        !           217: *>                    below EPS.
        !           218: *>          On exit :
        !           219: *>          WORK(1) = SCALE is the scaling factor such that SCALE*SVA(1:N)
        !           220: *>                    are the computed singular values of A.
        !           221: *>                    (See description of SVA().)
        !           222: *>          WORK(2) = NINT(WORK(2)) is the number of the computed nonzero
        !           223: *>                    singular values.
        !           224: *>          WORK(3) = NINT(WORK(3)) is the number of the computed singular
        !           225: *>                    values that are larger than the underflow threshold.
        !           226: *>          WORK(4) = NINT(WORK(4)) is the number of sweeps of Jacobi
        !           227: *>                    rotations needed for numerical convergence.
        !           228: *>          WORK(5) = max_{i.NE.j} |COS(A(:,i),A(:,j))| in the last sweep.
        !           229: *>                    This is useful information in cases when DGESVJ did
        !           230: *>                    not converge, as it can be used to estimate whether
        !           231: *>                    the output is stil useful and for post festum analysis.
        !           232: *>          WORK(6) = the largest absolute value over all sines of the
        !           233: *>                    Jacobi rotation angles in the last sweep. It can be
        !           234: *>                    useful for a post festum analysis.
        !           235: *> \endverbatim
        !           236: *>
        !           237: *> \param[in] LWORK
        !           238: *> \verbatim
        !           239: *>          LWORK is INTEGER
        !           240: *>          length of WORK, WORK >= MAX(6,M+N)
        !           241: *> \endverbatim
        !           242: *>
        !           243: *> \param[out] INFO
        !           244: *> \verbatim
        !           245: *>          INFO is INTEGER
        !           246: *>          = 0 : successful exit.
        !           247: *>          < 0 : if INFO = -i, then the i-th argument had an illegal value
        !           248: *>          > 0 : DGESVJ did not converge in the maximal allowed number (30)
        !           249: *>                of sweeps. The output may still be useful. See the
        !           250: *>                description of WORK.
        !           251: *> \endverbatim
        !           252: *
        !           253: *  Authors:
        !           254: *  ========
        !           255: *
        !           256: *> \author Univ. of Tennessee 
        !           257: *> \author Univ. of California Berkeley 
        !           258: *> \author Univ. of Colorado Denver 
        !           259: *> \author NAG Ltd. 
        !           260: *
        !           261: *> \date November 2011
        !           262: *
        !           263: *> \ingroup doubleGEcomputational
        !           264: *
        !           265: *> \par Further Details:
        !           266: *  =====================
        !           267: *>
        !           268: *> \verbatim
        !           269: *>
        !           270: *>  The orthogonal N-by-N matrix V is obtained as a product of Jacobi plane
        !           271: *>  rotations. The rotations are implemented as fast scaled rotations of
        !           272: *>  Anda and Park [1]. In the case of underflow of the Jacobi angle, a
        !           273: *>  modified Jacobi transformation of Drmac [4] is used. Pivot strategy uses
        !           274: *>  column interchanges of de Rijk [2]. The relative accuracy of the computed
        !           275: *>  singular values and the accuracy of the computed singular vectors (in
        !           276: *>  angle metric) is as guaranteed by the theory of Demmel and Veselic [3].
        !           277: *>  The condition number that determines the accuracy in the full rank case
        !           278: *>  is essentially min_{D=diag} kappa(A*D), where kappa(.) is the
        !           279: *>  spectral condition number. The best performance of this Jacobi SVD
        !           280: *>  procedure is achieved if used in an  accelerated version of Drmac and
        !           281: *>  Veselic [5,6], and it is the kernel routine in the SIGMA library [7].
        !           282: *>  Some tunning parameters (marked with [TP]) are available for the
        !           283: *>  implementer.
        !           284: *>  The computational range for the nonzero singular values is the  machine
        !           285: *>  number interval ( UNDERFLOW , OVERFLOW ). In extreme cases, even
        !           286: *>  denormalized singular values can be computed with the corresponding
        !           287: *>  gradual loss of accurate digits.
        !           288: *> \endverbatim
        !           289: *
        !           290: *> \par Contributors:
        !           291: *  ==================
        !           292: *>
        !           293: *> \verbatim
        !           294: *>
        !           295: *>  ============
        !           296: *>
        !           297: *>  Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
        !           298: *> \endverbatim
        !           299: *
        !           300: *> \par References:
        !           301: *  ================
        !           302: *>
        !           303: *> \verbatim
        !           304: *>
        !           305: *> [1] A. A. Anda and H. Park: Fast plane rotations with dynamic scaling.
        !           306: *>     SIAM J. matrix Anal. Appl., Vol. 15 (1994), pp. 162-174.
        !           307: *> [2] P. P. M. De Rijk: A one-sided Jacobi algorithm for computing the
        !           308: *>     singular value decomposition on a vector computer.
        !           309: *>     SIAM J. Sci. Stat. Comp., Vol. 10 (1998), pp. 359-371.
        !           310: *> [3] J. Demmel and K. Veselic: Jacobi method is more accurate than QR.
        !           311: *> [4] Z. Drmac: Implementation of Jacobi rotations for accurate singular
        !           312: *>     value computation in floating point arithmetic.
        !           313: *>     SIAM J. Sci. Comp., Vol. 18 (1997), pp. 1200-1222.
        !           314: *> [5] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm I.
        !           315: *>     SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1322-1342.
        !           316: *>     LAPACK Working note 169.
        !           317: *> [6] Z. Drmac and K. Veselic: New fast and accurate Jacobi SVD algorithm II.
        !           318: *>     SIAM J. Matrix Anal. Appl. Vol. 35, No. 2 (2008), pp. 1343-1362.
        !           319: *>     LAPACK Working note 170.
        !           320: *> [7] Z. Drmac: SIGMA - mathematical software library for accurate SVD, PSV,
        !           321: *>     QSVD, (H,K)-SVD computations.
        !           322: *>     Department of Mathematics, University of Zagreb, 2008.
        !           323: *> \endverbatim
        !           324: *
        !           325: *>  \par Bugs, examples and comments:
        !           326: *   =================================
        !           327: *>
        !           328: *> \verbatim
        !           329: *>  ===========================
        !           330: *>  Please report all bugs and send interesting test examples and comments to
        !           331: *>  drmac@math.hr. Thank you.
        !           332: *> \endverbatim
        !           333: *>
        !           334: *  =====================================================================
1.1       bertrand  335:       SUBROUTINE DGESVJ( JOBA, JOBU, JOBV, M, N, A, LDA, SVA, MV, V,
1.6       bertrand  336:      $                   LDV, WORK, LWORK, INFO )
1.1       bertrand  337: *
1.7     ! bertrand  338: *  -- LAPACK computational routine (version 3.4.0) --
1.1       bertrand  339: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    340: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.7     ! bertrand  341: *     November 2011
1.1       bertrand  342: *
                    343: *     .. Scalar Arguments ..
                    344:       INTEGER            INFO, LDA, LDV, LWORK, M, MV, N
                    345:       CHARACTER*1        JOBA, JOBU, JOBV
                    346: *     ..
                    347: *     .. Array Arguments ..
                    348:       DOUBLE PRECISION   A( LDA, * ), SVA( N ), V( LDV, * ),
1.6       bertrand  349:      $                   WORK( LWORK )
1.1       bertrand  350: *     ..
                    351: *
                    352: *  =====================================================================
                    353: *
                    354: *     .. Local Parameters ..
                    355:       DOUBLE PRECISION   ZERO, HALF, ONE, TWO
                    356:       PARAMETER          ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0,
1.6       bertrand  357:      $                   TWO = 2.0D0 )
1.1       bertrand  358:       INTEGER            NSWEEP
                    359:       PARAMETER          ( NSWEEP = 30 )
                    360: *     ..
                    361: *     .. Local Scalars ..
                    362:       DOUBLE PRECISION   AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG,
1.6       bertrand  363:      $                   BIGTHETA, CS, CTOL, EPSLN, LARGE, MXAAPQ,
                    364:      $                   MXSINJ, ROOTBIG, ROOTEPS, ROOTSFMIN, ROOTTOL,
                    365:      $                   SKL, SFMIN, SMALL, SN, T, TEMP1, THETA,
                    366:      $                   THSIGN, TOL
1.1       bertrand  367:       INTEGER            BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1,
1.6       bertrand  368:      $                   ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, N2, N34,
                    369:      $                   N4, NBL, NOTROT, p, PSKIPPED, q, ROWSKIP,
                    370:      $                   SWBAND
1.1       bertrand  371:       LOGICAL            APPLV, GOSCALE, LOWER, LSVEC, NOSCALE, ROTOK,
1.6       bertrand  372:      $                   RSVEC, UCTOL, UPPER
1.1       bertrand  373: *     ..
                    374: *     .. Local Arrays ..
                    375:       DOUBLE PRECISION   FASTR( 5 )
                    376: *     ..
                    377: *     .. Intrinsic Functions ..
                    378:       INTRINSIC          DABS, DMAX1, DMIN1, DBLE, MIN0, DSIGN, DSQRT
                    379: *     ..
                    380: *     .. External Functions ..
                    381: *     ..
                    382: *     from BLAS
                    383:       DOUBLE PRECISION   DDOT, DNRM2
                    384:       EXTERNAL           DDOT, DNRM2
                    385:       INTEGER            IDAMAX
                    386:       EXTERNAL           IDAMAX
                    387: *     from LAPACK
                    388:       DOUBLE PRECISION   DLAMCH
                    389:       EXTERNAL           DLAMCH
                    390:       LOGICAL            LSAME
                    391:       EXTERNAL           LSAME
                    392: *     ..
                    393: *     .. External Subroutines ..
                    394: *     ..
                    395: *     from BLAS
                    396:       EXTERNAL           DAXPY, DCOPY, DROTM, DSCAL, DSWAP
                    397: *     from LAPACK
                    398:       EXTERNAL           DLASCL, DLASET, DLASSQ, XERBLA
                    399: *
                    400:       EXTERNAL           DGSVJ0, DGSVJ1
                    401: *     ..
                    402: *     .. Executable Statements ..
                    403: *
                    404: *     Test the input arguments
                    405: *
                    406:       LSVEC = LSAME( JOBU, 'U' )
                    407:       UCTOL = LSAME( JOBU, 'C' )
                    408:       RSVEC = LSAME( JOBV, 'V' )
                    409:       APPLV = LSAME( JOBV, 'A' )
                    410:       UPPER = LSAME( JOBA, 'U' )
                    411:       LOWER = LSAME( JOBA, 'L' )
                    412: *
                    413:       IF( .NOT.( UPPER .OR. LOWER .OR. LSAME( JOBA, 'G' ) ) ) THEN
                    414:          INFO = -1
                    415:       ELSE IF( .NOT.( LSVEC .OR. UCTOL .OR. LSAME( JOBU, 'N' ) ) ) THEN
                    416:          INFO = -2
                    417:       ELSE IF( .NOT.( RSVEC .OR. APPLV .OR. LSAME( JOBV, 'N' ) ) ) THEN
                    418:          INFO = -3
                    419:       ELSE IF( M.LT.0 ) THEN
                    420:          INFO = -4
                    421:       ELSE IF( ( N.LT.0 ) .OR. ( N.GT.M ) ) THEN
                    422:          INFO = -5
                    423:       ELSE IF( LDA.LT.M ) THEN
                    424:          INFO = -7
                    425:       ELSE IF( MV.LT.0 ) THEN
                    426:          INFO = -9
                    427:       ELSE IF( ( RSVEC .AND. ( LDV.LT.N ) ) .OR.
1.6       bertrand  428:      $         ( APPLV .AND. ( LDV.LT.MV ) ) ) THEN
1.1       bertrand  429:          INFO = -11
                    430:       ELSE IF( UCTOL .AND. ( WORK( 1 ).LE.ONE ) ) THEN
                    431:          INFO = -12
                    432:       ELSE IF( LWORK.LT.MAX0( M+N, 6 ) ) THEN
                    433:          INFO = -13
                    434:       ELSE
                    435:          INFO = 0
                    436:       END IF
                    437: *
                    438: *     #:(
                    439:       IF( INFO.NE.0 ) THEN
                    440:          CALL XERBLA( 'DGESVJ', -INFO )
                    441:          RETURN
                    442:       END IF
                    443: *
                    444: * #:) Quick return for void matrix
                    445: *
                    446:       IF( ( M.EQ.0 ) .OR. ( N.EQ.0 ) )RETURN
                    447: *
                    448: *     Set numerical parameters
                    449: *     The stopping criterion for Jacobi rotations is
                    450: *
                    451: *     max_{i<>j}|A(:,i)^T * A(:,j)|/(||A(:,i)||*||A(:,j)||) < CTOL*EPS
                    452: *
                    453: *     where EPS is the round-off and CTOL is defined as follows:
                    454: *
                    455:       IF( UCTOL ) THEN
                    456: *        ... user controlled
                    457:          CTOL = WORK( 1 )
                    458:       ELSE
                    459: *        ... default
                    460:          IF( LSVEC .OR. RSVEC .OR. APPLV ) THEN
                    461:             CTOL = DSQRT( DBLE( M ) )
                    462:          ELSE
                    463:             CTOL = DBLE( M )
                    464:          END IF
                    465:       END IF
                    466: *     ... and the machine dependent parameters are
                    467: *[!]  (Make sure that DLAMCH() works properly on the target machine.)
                    468: *
1.4       bertrand  469:       EPSLN = DLAMCH( 'Epsilon' )
                    470:       ROOTEPS = DSQRT( EPSLN )
1.1       bertrand  471:       SFMIN = DLAMCH( 'SafeMinimum' )
                    472:       ROOTSFMIN = DSQRT( SFMIN )
1.4       bertrand  473:       SMALL = SFMIN / EPSLN
1.1       bertrand  474:       BIG = DLAMCH( 'Overflow' )
                    475: *     BIG         = ONE    / SFMIN
                    476:       ROOTBIG = ONE / ROOTSFMIN
                    477:       LARGE = BIG / DSQRT( DBLE( M*N ) )
                    478:       BIGTHETA = ONE / ROOTEPS
                    479: *
1.4       bertrand  480:       TOL = CTOL*EPSLN
1.1       bertrand  481:       ROOTTOL = DSQRT( TOL )
                    482: *
1.4       bertrand  483:       IF( DBLE( M )*EPSLN.GE.ONE ) THEN
1.6       bertrand  484:          INFO = -4
1.1       bertrand  485:          CALL XERBLA( 'DGESVJ', -INFO )
                    486:          RETURN
                    487:       END IF
                    488: *
                    489: *     Initialize the right singular vector matrix.
                    490: *
                    491:       IF( RSVEC ) THEN
                    492:          MVL = N
                    493:          CALL DLASET( 'A', MVL, N, ZERO, ONE, V, LDV )
                    494:       ELSE IF( APPLV ) THEN
                    495:          MVL = MV
                    496:       END IF
                    497:       RSVEC = RSVEC .OR. APPLV
                    498: *
                    499: *     Initialize SVA( 1:N ) = ( ||A e_i||_2, i = 1:N )
                    500: *(!)  If necessary, scale A to protect the largest singular value
                    501: *     from overflow. It is possible that saving the largest singular
                    502: *     value destroys the information about the small ones.
                    503: *     This initial scaling is almost minimal in the sense that the
                    504: *     goal is to make sure that no column norm overflows, and that
                    505: *     DSQRT(N)*max_i SVA(i) does not overflow. If INFinite entries
                    506: *     in A are detected, the procedure returns with INFO=-6.
                    507: *
1.4       bertrand  508:       SKL= ONE / DSQRT( DBLE( M )*DBLE( N ) )
1.1       bertrand  509:       NOSCALE = .TRUE.
                    510:       GOSCALE = .TRUE.
                    511: *
                    512:       IF( LOWER ) THEN
                    513: *        the input matrix is M-by-N lower triangular (trapezoidal)
                    514:          DO 1874 p = 1, N
                    515:             AAPP = ZERO
1.4       bertrand  516:             AAQQ = ONE
1.1       bertrand  517:             CALL DLASSQ( M-p+1, A( p, p ), 1, AAPP, AAQQ )
                    518:             IF( AAPP.GT.BIG ) THEN
                    519:                INFO = -6
                    520:                CALL XERBLA( 'DGESVJ', -INFO )
                    521:                RETURN
                    522:             END IF
                    523:             AAQQ = DSQRT( AAQQ )
                    524:             IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN
                    525:                SVA( p ) = AAPP*AAQQ
                    526:             ELSE
                    527:                NOSCALE = .FALSE.
1.4       bertrand  528:                SVA( p ) = AAPP*( AAQQ*SKL)
1.1       bertrand  529:                IF( GOSCALE ) THEN
                    530:                   GOSCALE = .FALSE.
                    531:                   DO 1873 q = 1, p - 1
1.4       bertrand  532:                      SVA( q ) = SVA( q )*SKL
1.1       bertrand  533:  1873             CONTINUE
                    534:                END IF
                    535:             END IF
                    536:  1874    CONTINUE
                    537:       ELSE IF( UPPER ) THEN
                    538: *        the input matrix is M-by-N upper triangular (trapezoidal)
                    539:          DO 2874 p = 1, N
                    540:             AAPP = ZERO
1.4       bertrand  541:             AAQQ = ONE
1.1       bertrand  542:             CALL DLASSQ( p, A( 1, p ), 1, AAPP, AAQQ )
                    543:             IF( AAPP.GT.BIG ) THEN
                    544:                INFO = -6
                    545:                CALL XERBLA( 'DGESVJ', -INFO )
                    546:                RETURN
                    547:             END IF
                    548:             AAQQ = DSQRT( AAQQ )
                    549:             IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN
                    550:                SVA( p ) = AAPP*AAQQ
                    551:             ELSE
                    552:                NOSCALE = .FALSE.
1.4       bertrand  553:                SVA( p ) = AAPP*( AAQQ*SKL)
1.1       bertrand  554:                IF( GOSCALE ) THEN
                    555:                   GOSCALE = .FALSE.
                    556:                   DO 2873 q = 1, p - 1
1.4       bertrand  557:                      SVA( q ) = SVA( q )*SKL
1.1       bertrand  558:  2873             CONTINUE
                    559:                END IF
                    560:             END IF
                    561:  2874    CONTINUE
                    562:       ELSE
                    563: *        the input matrix is M-by-N general dense
                    564:          DO 3874 p = 1, N
                    565:             AAPP = ZERO
1.4       bertrand  566:             AAQQ = ONE
1.1       bertrand  567:             CALL DLASSQ( M, A( 1, p ), 1, AAPP, AAQQ )
                    568:             IF( AAPP.GT.BIG ) THEN
                    569:                INFO = -6
                    570:                CALL XERBLA( 'DGESVJ', -INFO )
                    571:                RETURN
                    572:             END IF
                    573:             AAQQ = DSQRT( AAQQ )
                    574:             IF( ( AAPP.LT.( BIG / AAQQ ) ) .AND. NOSCALE ) THEN
                    575:                SVA( p ) = AAPP*AAQQ
                    576:             ELSE
                    577:                NOSCALE = .FALSE.
1.4       bertrand  578:                SVA( p ) = AAPP*( AAQQ*SKL)
1.1       bertrand  579:                IF( GOSCALE ) THEN
                    580:                   GOSCALE = .FALSE.
                    581:                   DO 3873 q = 1, p - 1
1.4       bertrand  582:                      SVA( q ) = SVA( q )*SKL
1.1       bertrand  583:  3873             CONTINUE
                    584:                END IF
                    585:             END IF
                    586:  3874    CONTINUE
                    587:       END IF
                    588: *
1.4       bertrand  589:       IF( NOSCALE )SKL= ONE
1.1       bertrand  590: *
                    591: *     Move the smaller part of the spectrum from the underflow threshold
                    592: *(!)  Start by determining the position of the nonzero entries of the
                    593: *     array SVA() relative to ( SFMIN, BIG ).
                    594: *
                    595:       AAPP = ZERO
                    596:       AAQQ = BIG
                    597:       DO 4781 p = 1, N
                    598:          IF( SVA( p ).NE.ZERO )AAQQ = DMIN1( AAQQ, SVA( p ) )
                    599:          AAPP = DMAX1( AAPP, SVA( p ) )
                    600:  4781 CONTINUE
                    601: *
                    602: * #:) Quick return for zero matrix
                    603: *
                    604:       IF( AAPP.EQ.ZERO ) THEN
                    605:          IF( LSVEC )CALL DLASET( 'G', M, N, ZERO, ONE, A, LDA )
                    606:          WORK( 1 ) = ONE
                    607:          WORK( 2 ) = ZERO
                    608:          WORK( 3 ) = ZERO
                    609:          WORK( 4 ) = ZERO
                    610:          WORK( 5 ) = ZERO
                    611:          WORK( 6 ) = ZERO
                    612:          RETURN
                    613:       END IF
                    614: *
                    615: * #:) Quick return for one-column matrix
                    616: *
                    617:       IF( N.EQ.1 ) THEN
1.4       bertrand  618:          IF( LSVEC )CALL DLASCL( 'G', 0, 0, SVA( 1 ), SKL, M, 1,
1.6       bertrand  619:      $                           A( 1, 1 ), LDA, IERR )
1.4       bertrand  620:          WORK( 1 ) = ONE / SKL
1.1       bertrand  621:          IF( SVA( 1 ).GE.SFMIN ) THEN
                    622:             WORK( 2 ) = ONE
                    623:          ELSE
                    624:             WORK( 2 ) = ZERO
                    625:          END IF
                    626:          WORK( 3 ) = ZERO
                    627:          WORK( 4 ) = ZERO
                    628:          WORK( 5 ) = ZERO
                    629:          WORK( 6 ) = ZERO
                    630:          RETURN
                    631:       END IF
                    632: *
                    633: *     Protect small singular values from underflow, and try to
                    634: *     avoid underflows/overflows in computing Jacobi rotations.
                    635: *
1.4       bertrand  636:       SN = DSQRT( SFMIN / EPSLN )
1.1       bertrand  637:       TEMP1 = DSQRT( BIG / DBLE( N ) )
                    638:       IF( ( AAPP.LE.SN ) .OR. ( AAQQ.GE.TEMP1 ) .OR.
1.6       bertrand  639:      $    ( ( SN.LE.AAQQ ) .AND. ( AAPP.LE.TEMP1 ) ) ) THEN
1.1       bertrand  640:          TEMP1 = DMIN1( BIG, TEMP1 / AAPP )
                    641: *         AAQQ  = AAQQ*TEMP1
                    642: *         AAPP  = AAPP*TEMP1
                    643:       ELSE IF( ( AAQQ.LE.SN ) .AND. ( AAPP.LE.TEMP1 ) ) THEN
                    644:          TEMP1 = DMIN1( SN / AAQQ, BIG / ( AAPP*DSQRT( DBLE( N ) ) ) )
                    645: *         AAQQ  = AAQQ*TEMP1
                    646: *         AAPP  = AAPP*TEMP1
                    647:       ELSE IF( ( AAQQ.GE.SN ) .AND. ( AAPP.GE.TEMP1 ) ) THEN
                    648:          TEMP1 = DMAX1( SN / AAQQ, TEMP1 / AAPP )
                    649: *         AAQQ  = AAQQ*TEMP1
                    650: *         AAPP  = AAPP*TEMP1
                    651:       ELSE IF( ( AAQQ.LE.SN ) .AND. ( AAPP.GE.TEMP1 ) ) THEN
                    652:          TEMP1 = DMIN1( SN / AAQQ, BIG / ( DSQRT( DBLE( N ) )*AAPP ) )
                    653: *         AAQQ  = AAQQ*TEMP1
                    654: *         AAPP  = AAPP*TEMP1
                    655:       ELSE
                    656:          TEMP1 = ONE
                    657:       END IF
                    658: *
                    659: *     Scale, if necessary
                    660: *
                    661:       IF( TEMP1.NE.ONE ) THEN
                    662:          CALL DLASCL( 'G', 0, 0, ONE, TEMP1, N, 1, SVA, N, IERR )
                    663:       END IF
1.4       bertrand  664:       SKL= TEMP1*SKL
                    665:       IF( SKL.NE.ONE ) THEN
                    666:          CALL DLASCL( JOBA, 0, 0, ONE, SKL, M, N, A, LDA, IERR )
                    667:          SKL= ONE / SKL
1.1       bertrand  668:       END IF
                    669: *
                    670: *     Row-cyclic Jacobi SVD algorithm with column pivoting
                    671: *
                    672:       EMPTSW = ( N*( N-1 ) ) / 2
                    673:       NOTROT = 0
                    674:       FASTR( 1 ) = ZERO
                    675: *
                    676: *     A is represented in factored form A = A * diag(WORK), where diag(WORK)
                    677: *     is initialized to identity. WORK is updated during fast scaled
                    678: *     rotations.
                    679: *
                    680:       DO 1868 q = 1, N
                    681:          WORK( q ) = ONE
                    682:  1868 CONTINUE
                    683: *
                    684: *
                    685:       SWBAND = 3
                    686: *[TP] SWBAND is a tuning parameter [TP]. It is meaningful and effective
                    687: *     if DGESVJ is used as a computational routine in the preconditioned
                    688: *     Jacobi SVD algorithm DGESVJ. For sweeps i=1:SWBAND the procedure
                    689: *     works on pivots inside a band-like region around the diagonal.
                    690: *     The boundaries are determined dynamically, based on the number of
                    691: *     pivots above a threshold.
                    692: *
                    693:       KBL = MIN0( 8, N )
                    694: *[TP] KBL is a tuning parameter that defines the tile size in the
                    695: *     tiling of the p-q loops of pivot pairs. In general, an optimal
                    696: *     value of KBL depends on the matrix dimensions and on the
                    697: *     parameters of the computer's memory.
                    698: *
                    699:       NBL = N / KBL
                    700:       IF( ( NBL*KBL ).NE.N )NBL = NBL + 1
                    701: *
                    702:       BLSKIP = KBL**2
                    703: *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
                    704: *
                    705:       ROWSKIP = MIN0( 5, KBL )
                    706: *[TP] ROWSKIP is a tuning parameter.
                    707: *
                    708:       LKAHEAD = 1
                    709: *[TP] LKAHEAD is a tuning parameter.
                    710: *
                    711: *     Quasi block transformations, using the lower (upper) triangular
                    712: *     structure of the input matrix. The quasi-block-cycling usually
                    713: *     invokes cubic convergence. Big part of this cycle is done inside
                    714: *     canonical subspaces of dimensions less than M.
                    715: *
                    716:       IF( ( LOWER .OR. UPPER ) .AND. ( N.GT.MAX0( 64, 4*KBL ) ) ) THEN
                    717: *[TP] The number of partition levels and the actual partition are
                    718: *     tuning parameters.
                    719:          N4 = N / 4
                    720:          N2 = N / 2
                    721:          N34 = 3*N4
                    722:          IF( APPLV ) THEN
                    723:             q = 0
                    724:          ELSE
                    725:             q = 1
                    726:          END IF
                    727: *
                    728:          IF( LOWER ) THEN
                    729: *
                    730: *     This works very well on lower triangular matrices, in particular
                    731: *     in the framework of the preconditioned Jacobi SVD (xGEJSV).
                    732: *     The idea is simple:
                    733: *     [+ 0 0 0]   Note that Jacobi transformations of [0 0]
                    734: *     [+ + 0 0]                                       [0 0]
                    735: *     [+ + x 0]   actually work on [x 0]              [x 0]
                    736: *     [+ + x x]                    [x x].             [x x]
                    737: *
                    738:             CALL DGSVJ0( JOBV, M-N34, N-N34, A( N34+1, N34+1 ), LDA,
1.6       bertrand  739:      $                   WORK( N34+1 ), SVA( N34+1 ), MVL,
                    740:      $                   V( N34*q+1, N34+1 ), LDV, EPSLN, SFMIN, TOL,
                    741:      $                   2, WORK( N+1 ), LWORK-N, IERR )
1.1       bertrand  742: *
                    743:             CALL DGSVJ0( JOBV, M-N2, N34-N2, A( N2+1, N2+1 ), LDA,
1.6       bertrand  744:      $                   WORK( N2+1 ), SVA( N2+1 ), MVL,
                    745:      $                   V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 2,
                    746:      $                   WORK( N+1 ), LWORK-N, IERR )
1.1       bertrand  747: *
                    748:             CALL DGSVJ1( JOBV, M-N2, N-N2, N4, A( N2+1, N2+1 ), LDA,
1.6       bertrand  749:      $                   WORK( N2+1 ), SVA( N2+1 ), MVL,
                    750:      $                   V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 1,
                    751:      $                   WORK( N+1 ), LWORK-N, IERR )
1.1       bertrand  752: *
                    753:             CALL DGSVJ0( JOBV, M-N4, N2-N4, A( N4+1, N4+1 ), LDA,
1.6       bertrand  754:      $                   WORK( N4+1 ), SVA( N4+1 ), MVL,
                    755:      $                   V( N4*q+1, N4+1 ), LDV, EPSLN, SFMIN, TOL, 1,
                    756:      $                   WORK( N+1 ), LWORK-N, IERR )
1.1       bertrand  757: *
                    758:             CALL DGSVJ0( JOBV, M, N4, A, LDA, WORK, SVA, MVL, V, LDV,
1.6       bertrand  759:      $                   EPSLN, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,
                    760:      $                   IERR )
1.1       bertrand  761: *
                    762:             CALL DGSVJ1( JOBV, M, N2, N4, A, LDA, WORK, SVA, MVL, V,
1.6       bertrand  763:      $                   LDV, EPSLN, SFMIN, TOL, 1, WORK( N+1 ),
                    764:      $                   LWORK-N, IERR )
1.1       bertrand  765: *
                    766: *
                    767:          ELSE IF( UPPER ) THEN
                    768: *
                    769: *
                    770:             CALL DGSVJ0( JOBV, N4, N4, A, LDA, WORK, SVA, MVL, V, LDV,
1.6       bertrand  771:      $                   EPSLN, SFMIN, TOL, 2, WORK( N+1 ), LWORK-N,
                    772:      $                   IERR )
1.1       bertrand  773: *
                    774:             CALL DGSVJ0( JOBV, N2, N4, A( 1, N4+1 ), LDA, WORK( N4+1 ),
1.6       bertrand  775:      $                   SVA( N4+1 ), MVL, V( N4*q+1, N4+1 ), LDV,
                    776:      $                   EPSLN, SFMIN, TOL, 1, WORK( N+1 ), LWORK-N,
                    777:      $                   IERR )
1.1       bertrand  778: *
                    779:             CALL DGSVJ1( JOBV, N2, N2, N4, A, LDA, WORK, SVA, MVL, V,
1.6       bertrand  780:      $                   LDV, EPSLN, SFMIN, TOL, 1, WORK( N+1 ),
                    781:      $                   LWORK-N, IERR )
1.1       bertrand  782: *
                    783:             CALL DGSVJ0( JOBV, N2+N4, N4, A( 1, N2+1 ), LDA,
1.6       bertrand  784:      $                   WORK( N2+1 ), SVA( N2+1 ), MVL,
                    785:      $                   V( N2*q+1, N2+1 ), LDV, EPSLN, SFMIN, TOL, 1,
                    786:      $                   WORK( N+1 ), LWORK-N, IERR )
1.1       bertrand  787: 
                    788:          END IF
                    789: *
                    790:       END IF
                    791: *
                    792: *     .. Row-cyclic pivot strategy with de Rijk's pivoting ..
                    793: *
                    794:       DO 1993 i = 1, NSWEEP
                    795: *
                    796: *     .. go go go ...
                    797: *
                    798:          MXAAPQ = ZERO
                    799:          MXSINJ = ZERO
                    800:          ISWROT = 0
                    801: *
                    802:          NOTROT = 0
                    803:          PSKIPPED = 0
                    804: *
                    805: *     Each sweep is unrolled using KBL-by-KBL tiles over the pivot pairs
                    806: *     1 <= p < q <= N. This is the first step toward a blocked implementation
                    807: *     of the rotations. New implementation, based on block transformations,
                    808: *     is under development.
                    809: *
                    810:          DO 2000 ibr = 1, NBL
                    811: *
                    812:             igl = ( ibr-1 )*KBL + 1
                    813: *
                    814:             DO 1002 ir1 = 0, MIN0( LKAHEAD, NBL-ibr )
                    815: *
                    816:                igl = igl + ir1*KBL
                    817: *
                    818:                DO 2001 p = igl, MIN0( igl+KBL-1, N-1 )
                    819: *
                    820: *     .. de Rijk's pivoting
                    821: *
                    822:                   q = IDAMAX( N-p+1, SVA( p ), 1 ) + p - 1
                    823:                   IF( p.NE.q ) THEN
                    824:                      CALL DSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
                    825:                      IF( RSVEC )CALL DSWAP( MVL, V( 1, p ), 1,
1.6       bertrand  826:      $                                      V( 1, q ), 1 )
1.1       bertrand  827:                      TEMP1 = SVA( p )
                    828:                      SVA( p ) = SVA( q )
                    829:                      SVA( q ) = TEMP1
                    830:                      TEMP1 = WORK( p )
                    831:                      WORK( p ) = WORK( q )
                    832:                      WORK( q ) = TEMP1
                    833:                   END IF
                    834: *
                    835:                   IF( ir1.EQ.0 ) THEN
                    836: *
                    837: *        Column norms are periodically updated by explicit
                    838: *        norm computation.
                    839: *        Caveat:
                    840: *        Unfortunately, some BLAS implementations compute DNRM2(M,A(1,p),1)
                    841: *        as DSQRT(DDOT(M,A(1,p),1,A(1,p),1)), which may cause the result to
                    842: *        overflow for ||A(:,p)||_2 > DSQRT(overflow_threshold), and to
                    843: *        underflow for ||A(:,p)||_2 < DSQRT(underflow_threshold).
                    844: *        Hence, DNRM2 cannot be trusted, not even in the case when
                    845: *        the true norm is far from the under(over)flow boundaries.
                    846: *        If properly implemented DNRM2 is available, the IF-THEN-ELSE
                    847: *        below should read "AAPP = DNRM2( M, A(1,p), 1 ) * WORK(p)".
                    848: *
                    849:                      IF( ( SVA( p ).LT.ROOTBIG ) .AND.
1.6       bertrand  850:      $                   ( SVA( p ).GT.ROOTSFMIN ) ) THEN
1.1       bertrand  851:                         SVA( p ) = DNRM2( M, A( 1, p ), 1 )*WORK( p )
                    852:                      ELSE
                    853:                         TEMP1 = ZERO
1.4       bertrand  854:                         AAPP = ONE
1.1       bertrand  855:                         CALL DLASSQ( M, A( 1, p ), 1, TEMP1, AAPP )
                    856:                         SVA( p ) = TEMP1*DSQRT( AAPP )*WORK( p )
                    857:                      END IF
                    858:                      AAPP = SVA( p )
                    859:                   ELSE
                    860:                      AAPP = SVA( p )
                    861:                   END IF
                    862: *
                    863:                   IF( AAPP.GT.ZERO ) THEN
                    864: *
                    865:                      PSKIPPED = 0
                    866: *
                    867:                      DO 2002 q = p + 1, MIN0( igl+KBL-1, N )
                    868: *
                    869:                         AAQQ = SVA( q )
                    870: *
                    871:                         IF( AAQQ.GT.ZERO ) THEN
                    872: *
                    873:                            AAPP0 = AAPP
                    874:                            IF( AAQQ.GE.ONE ) THEN
                    875:                               ROTOK = ( SMALL*AAPP ).LE.AAQQ
                    876:                               IF( AAPP.LT.( BIG / AAQQ ) ) THEN
                    877:                                  AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
1.6       bertrand  878:      $                                  q ), 1 )*WORK( p )*WORK( q ) /
                    879:      $                                  AAQQ ) / AAPP
1.1       bertrand  880:                               ELSE
                    881:                                  CALL DCOPY( M, A( 1, p ), 1,
1.6       bertrand  882:      $                                       WORK( N+1 ), 1 )
1.1       bertrand  883:                                  CALL DLASCL( 'G', 0, 0, AAPP,
1.6       bertrand  884:      $                                        WORK( p ), M, 1,
                    885:      $                                        WORK( N+1 ), LDA, IERR )
1.1       bertrand  886:                                  AAPQ = DDOT( M, WORK( N+1 ), 1,
1.6       bertrand  887:      $                                  A( 1, q ), 1 )*WORK( q ) / AAQQ
1.1       bertrand  888:                               END IF
                    889:                            ELSE
                    890:                               ROTOK = AAPP.LE.( AAQQ / SMALL )
                    891:                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
                    892:                                  AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
1.6       bertrand  893:      $                                  q ), 1 )*WORK( p )*WORK( q ) /
                    894:      $                                  AAQQ ) / AAPP
1.1       bertrand  895:                               ELSE
                    896:                                  CALL DCOPY( M, A( 1, q ), 1,
1.6       bertrand  897:      $                                       WORK( N+1 ), 1 )
1.1       bertrand  898:                                  CALL DLASCL( 'G', 0, 0, AAQQ,
1.6       bertrand  899:      $                                        WORK( q ), M, 1,
                    900:      $                                        WORK( N+1 ), LDA, IERR )
1.1       bertrand  901:                                  AAPQ = DDOT( M, WORK( N+1 ), 1,
1.6       bertrand  902:      $                                  A( 1, p ), 1 )*WORK( p ) / AAPP
1.1       bertrand  903:                               END IF
                    904:                            END IF
                    905: *
                    906:                            MXAAPQ = DMAX1( MXAAPQ, DABS( AAPQ ) )
                    907: *
                    908: *        TO rotate or NOT to rotate, THAT is the question ...
                    909: *
                    910:                            IF( DABS( AAPQ ).GT.TOL ) THEN
                    911: *
                    912: *           .. rotate
                    913: *[RTD]      ROTATED = ROTATED + ONE
                    914: *
                    915:                               IF( ir1.EQ.0 ) THEN
                    916:                                  NOTROT = 0
                    917:                                  PSKIPPED = 0
                    918:                                  ISWROT = ISWROT + 1
                    919:                               END IF
                    920: *
                    921:                               IF( ROTOK ) THEN
                    922: *
                    923:                                  AQOAP = AAQQ / AAPP
                    924:                                  APOAQ = AAPP / AAQQ
1.6       bertrand  925:                                  THETA = -HALF*DABS(AQOAP-APOAQ)/AAPQ
1.1       bertrand  926: *
                    927:                                  IF( DABS( THETA ).GT.BIGTHETA ) THEN
                    928: *
                    929:                                     T = HALF / THETA
                    930:                                     FASTR( 3 ) = T*WORK( p ) / WORK( q )
                    931:                                     FASTR( 4 ) = -T*WORK( q ) /
1.6       bertrand  932:      $                                           WORK( p )
1.1       bertrand  933:                                     CALL DROTM( M, A( 1, p ), 1,
1.6       bertrand  934:      $                                          A( 1, q ), 1, FASTR )
1.1       bertrand  935:                                     IF( RSVEC )CALL DROTM( MVL,
1.6       bertrand  936:      $                                              V( 1, p ), 1,
                    937:      $                                              V( 1, q ), 1,
                    938:      $                                              FASTR )
1.1       bertrand  939:                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
1.6       bertrand  940:      $                                         ONE+T*APOAQ*AAPQ ) )
1.4       bertrand  941:                                     AAPP = AAPP*DSQRT( DMAX1( ZERO,
1.6       bertrand  942:      $                                     ONE-T*AQOAP*AAPQ ) )
1.1       bertrand  943:                                     MXSINJ = DMAX1( MXSINJ, DABS( T ) )
                    944: *
                    945:                                  ELSE
                    946: *
                    947: *                 .. choose correct signum for THETA and rotate
                    948: *
                    949:                                     THSIGN = -DSIGN( ONE, AAPQ )
                    950:                                     T = ONE / ( THETA+THSIGN*
1.6       bertrand  951:      $                                  DSQRT( ONE+THETA*THETA ) )
1.1       bertrand  952:                                     CS = DSQRT( ONE / ( ONE+T*T ) )
                    953:                                     SN = T*CS
                    954: *
                    955:                                     MXSINJ = DMAX1( MXSINJ, DABS( SN ) )
                    956:                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
1.6       bertrand  957:      $                                         ONE+T*APOAQ*AAPQ ) )
1.1       bertrand  958:                                     AAPP = AAPP*DSQRT( DMAX1( ZERO,
1.6       bertrand  959:      $                                     ONE-T*AQOAP*AAPQ ) )
1.1       bertrand  960: *
                    961:                                     APOAQ = WORK( p ) / WORK( q )
                    962:                                     AQOAP = WORK( q ) / WORK( p )
                    963:                                     IF( WORK( p ).GE.ONE ) THEN
                    964:                                        IF( WORK( q ).GE.ONE ) THEN
                    965:                                           FASTR( 3 ) = T*APOAQ
                    966:                                           FASTR( 4 ) = -T*AQOAP
                    967:                                           WORK( p ) = WORK( p )*CS
                    968:                                           WORK( q ) = WORK( q )*CS
                    969:                                           CALL DROTM( M, A( 1, p ), 1,
1.6       bertrand  970:      $                                                A( 1, q ), 1,
                    971:      $                                                FASTR )
1.1       bertrand  972:                                           IF( RSVEC )CALL DROTM( MVL,
1.6       bertrand  973:      $                                        V( 1, p ), 1, V( 1, q ),
                    974:      $                                        1, FASTR )
1.1       bertrand  975:                                        ELSE
                    976:                                           CALL DAXPY( M, -T*AQOAP,
1.6       bertrand  977:      $                                                A( 1, q ), 1,
                    978:      $                                                A( 1, p ), 1 )
1.1       bertrand  979:                                           CALL DAXPY( M, CS*SN*APOAQ,
1.6       bertrand  980:      $                                                A( 1, p ), 1,
                    981:      $                                                A( 1, q ), 1 )
1.1       bertrand  982:                                           WORK( p ) = WORK( p )*CS
                    983:                                           WORK( q ) = WORK( q ) / CS
                    984:                                           IF( RSVEC ) THEN
                    985:                                              CALL DAXPY( MVL, -T*AQOAP,
1.6       bertrand  986:      $                                                   V( 1, q ), 1,
                    987:      $                                                   V( 1, p ), 1 )
1.1       bertrand  988:                                              CALL DAXPY( MVL,
1.6       bertrand  989:      $                                                   CS*SN*APOAQ,
                    990:      $                                                   V( 1, p ), 1,
                    991:      $                                                   V( 1, q ), 1 )
1.1       bertrand  992:                                           END IF
                    993:                                        END IF
                    994:                                     ELSE
                    995:                                        IF( WORK( q ).GE.ONE ) THEN
                    996:                                           CALL DAXPY( M, T*APOAQ,
1.6       bertrand  997:      $                                                A( 1, p ), 1,
                    998:      $                                                A( 1, q ), 1 )
1.1       bertrand  999:                                           CALL DAXPY( M, -CS*SN*AQOAP,
1.6       bertrand 1000:      $                                                A( 1, q ), 1,
                   1001:      $                                                A( 1, p ), 1 )
1.1       bertrand 1002:                                           WORK( p ) = WORK( p ) / CS
                   1003:                                           WORK( q ) = WORK( q )*CS
                   1004:                                           IF( RSVEC ) THEN
                   1005:                                              CALL DAXPY( MVL, T*APOAQ,
1.6       bertrand 1006:      $                                                   V( 1, p ), 1,
                   1007:      $                                                   V( 1, q ), 1 )
1.1       bertrand 1008:                                              CALL DAXPY( MVL,
1.6       bertrand 1009:      $                                                   -CS*SN*AQOAP,
                   1010:      $                                                   V( 1, q ), 1,
                   1011:      $                                                   V( 1, p ), 1 )
1.1       bertrand 1012:                                           END IF
                   1013:                                        ELSE
                   1014:                                           IF( WORK( p ).GE.WORK( q ) )
1.6       bertrand 1015:      $                                        THEN
1.1       bertrand 1016:                                              CALL DAXPY( M, -T*AQOAP,
1.6       bertrand 1017:      $                                                   A( 1, q ), 1,
                   1018:      $                                                   A( 1, p ), 1 )
1.1       bertrand 1019:                                              CALL DAXPY( M, CS*SN*APOAQ,
1.6       bertrand 1020:      $                                                   A( 1, p ), 1,
                   1021:      $                                                   A( 1, q ), 1 )
1.1       bertrand 1022:                                              WORK( p ) = WORK( p )*CS
                   1023:                                              WORK( q ) = WORK( q ) / CS
                   1024:                                              IF( RSVEC ) THEN
                   1025:                                                 CALL DAXPY( MVL,
1.6       bertrand 1026:      $                                               -T*AQOAP,
                   1027:      $                                               V( 1, q ), 1,
                   1028:      $                                               V( 1, p ), 1 )
1.1       bertrand 1029:                                                 CALL DAXPY( MVL,
1.6       bertrand 1030:      $                                               CS*SN*APOAQ,
                   1031:      $                                               V( 1, p ), 1,
                   1032:      $                                               V( 1, q ), 1 )
1.1       bertrand 1033:                                              END IF
                   1034:                                           ELSE
                   1035:                                              CALL DAXPY( M, T*APOAQ,
1.6       bertrand 1036:      $                                                   A( 1, p ), 1,
                   1037:      $                                                   A( 1, q ), 1 )
1.1       bertrand 1038:                                              CALL DAXPY( M,
1.6       bertrand 1039:      $                                                   -CS*SN*AQOAP,
                   1040:      $                                                   A( 1, q ), 1,
                   1041:      $                                                   A( 1, p ), 1 )
1.1       bertrand 1042:                                              WORK( p ) = WORK( p ) / CS
                   1043:                                              WORK( q ) = WORK( q )*CS
                   1044:                                              IF( RSVEC ) THEN
                   1045:                                                 CALL DAXPY( MVL,
1.6       bertrand 1046:      $                                               T*APOAQ, V( 1, p ),
                   1047:      $                                               1, V( 1, q ), 1 )
1.1       bertrand 1048:                                                 CALL DAXPY( MVL,
1.6       bertrand 1049:      $                                               -CS*SN*AQOAP,
                   1050:      $                                               V( 1, q ), 1,
                   1051:      $                                               V( 1, p ), 1 )
1.1       bertrand 1052:                                              END IF
                   1053:                                           END IF
                   1054:                                        END IF
                   1055:                                     END IF
                   1056:                                  END IF
                   1057: *
                   1058:                               ELSE
                   1059: *              .. have to use modified Gram-Schmidt like transformation
                   1060:                                  CALL DCOPY( M, A( 1, p ), 1,
1.6       bertrand 1061:      $                                       WORK( N+1 ), 1 )
1.1       bertrand 1062:                                  CALL DLASCL( 'G', 0, 0, AAPP, ONE, M,
1.6       bertrand 1063:      $                                        1, WORK( N+1 ), LDA,
                   1064:      $                                        IERR )
1.1       bertrand 1065:                                  CALL DLASCL( 'G', 0, 0, AAQQ, ONE, M,
1.6       bertrand 1066:      $                                        1, A( 1, q ), LDA, IERR )
1.1       bertrand 1067:                                  TEMP1 = -AAPQ*WORK( p ) / WORK( q )
                   1068:                                  CALL DAXPY( M, TEMP1, WORK( N+1 ), 1,
1.6       bertrand 1069:      $                                       A( 1, q ), 1 )
1.1       bertrand 1070:                                  CALL DLASCL( 'G', 0, 0, ONE, AAQQ, M,
1.6       bertrand 1071:      $                                        1, A( 1, q ), LDA, IERR )
1.1       bertrand 1072:                                  SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
1.6       bertrand 1073:      $                                      ONE-AAPQ*AAPQ ) )
1.1       bertrand 1074:                                  MXSINJ = DMAX1( MXSINJ, SFMIN )
                   1075:                               END IF
                   1076: *           END IF ROTOK THEN ... ELSE
                   1077: *
                   1078: *           In the case of cancellation in updating SVA(q), SVA(p)
                   1079: *           recompute SVA(q), SVA(p).
                   1080: *
                   1081:                               IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
1.6       bertrand 1082:      $                            THEN
1.1       bertrand 1083:                                  IF( ( AAQQ.LT.ROOTBIG ) .AND.
1.6       bertrand 1084:      $                               ( AAQQ.GT.ROOTSFMIN ) ) THEN
1.1       bertrand 1085:                                     SVA( q ) = DNRM2( M, A( 1, q ), 1 )*
1.6       bertrand 1086:      $                                         WORK( q )
1.1       bertrand 1087:                                  ELSE
                   1088:                                     T = ZERO
1.4       bertrand 1089:                                     AAQQ = ONE
1.1       bertrand 1090:                                     CALL DLASSQ( M, A( 1, q ), 1, T,
1.6       bertrand 1091:      $                                           AAQQ )
1.1       bertrand 1092:                                     SVA( q ) = T*DSQRT( AAQQ )*WORK( q )
                   1093:                                  END IF
                   1094:                               END IF
                   1095:                               IF( ( AAPP / AAPP0 ).LE.ROOTEPS ) THEN
                   1096:                                  IF( ( AAPP.LT.ROOTBIG ) .AND.
1.6       bertrand 1097:      $                               ( AAPP.GT.ROOTSFMIN ) ) THEN
1.1       bertrand 1098:                                     AAPP = DNRM2( M, A( 1, p ), 1 )*
1.6       bertrand 1099:      $                                     WORK( p )
1.1       bertrand 1100:                                  ELSE
                   1101:                                     T = ZERO
1.4       bertrand 1102:                                     AAPP = ONE
1.1       bertrand 1103:                                     CALL DLASSQ( M, A( 1, p ), 1, T,
1.6       bertrand 1104:      $                                           AAPP )
1.1       bertrand 1105:                                     AAPP = T*DSQRT( AAPP )*WORK( p )
                   1106:                                  END IF
                   1107:                                  SVA( p ) = AAPP
                   1108:                               END IF
                   1109: *
                   1110:                            ELSE
                   1111: *        A(:,p) and A(:,q) already numerically orthogonal
                   1112:                               IF( ir1.EQ.0 )NOTROT = NOTROT + 1
                   1113: *[RTD]      SKIPPED  = SKIPPED  + 1
                   1114:                               PSKIPPED = PSKIPPED + 1
                   1115:                            END IF
                   1116:                         ELSE
                   1117: *        A(:,q) is zero column
                   1118:                            IF( ir1.EQ.0 )NOTROT = NOTROT + 1
                   1119:                            PSKIPPED = PSKIPPED + 1
                   1120:                         END IF
                   1121: *
                   1122:                         IF( ( i.LE.SWBAND ) .AND.
1.6       bertrand 1123:      $                      ( PSKIPPED.GT.ROWSKIP ) ) THEN
1.1       bertrand 1124:                            IF( ir1.EQ.0 )AAPP = -AAPP
                   1125:                            NOTROT = 0
                   1126:                            GO TO 2103
                   1127:                         END IF
                   1128: *
                   1129:  2002                CONTINUE
                   1130: *     END q-LOOP
                   1131: *
                   1132:  2103                CONTINUE
                   1133: *     bailed out of q-loop
                   1134: *
                   1135:                      SVA( p ) = AAPP
                   1136: *
                   1137:                   ELSE
                   1138:                      SVA( p ) = AAPP
                   1139:                      IF( ( ir1.EQ.0 ) .AND. ( AAPP.EQ.ZERO ) )
1.6       bertrand 1140:      $                   NOTROT = NOTROT + MIN0( igl+KBL-1, N ) - p
1.1       bertrand 1141:                   END IF
                   1142: *
                   1143:  2001          CONTINUE
                   1144: *     end of the p-loop
                   1145: *     end of doing the block ( ibr, ibr )
                   1146:  1002       CONTINUE
                   1147: *     end of ir1-loop
                   1148: *
                   1149: * ... go to the off diagonal blocks
                   1150: *
                   1151:             igl = ( ibr-1 )*KBL + 1
                   1152: *
                   1153:             DO 2010 jbc = ibr + 1, NBL
                   1154: *
                   1155:                jgl = ( jbc-1 )*KBL + 1
                   1156: *
                   1157: *        doing the block at ( ibr, jbc )
                   1158: *
                   1159:                IJBLSK = 0
                   1160:                DO 2100 p = igl, MIN0( igl+KBL-1, N )
                   1161: *
                   1162:                   AAPP = SVA( p )
                   1163:                   IF( AAPP.GT.ZERO ) THEN
                   1164: *
                   1165:                      PSKIPPED = 0
                   1166: *
                   1167:                      DO 2200 q = jgl, MIN0( jgl+KBL-1, N )
                   1168: *
                   1169:                         AAQQ = SVA( q )
                   1170:                         IF( AAQQ.GT.ZERO ) THEN
                   1171:                            AAPP0 = AAPP
                   1172: *
                   1173: *     .. M x 2 Jacobi SVD ..
                   1174: *
                   1175: *        Safe Gram matrix computation
                   1176: *
                   1177:                            IF( AAQQ.GE.ONE ) THEN
                   1178:                               IF( AAPP.GE.AAQQ ) THEN
                   1179:                                  ROTOK = ( SMALL*AAPP ).LE.AAQQ
                   1180:                               ELSE
                   1181:                                  ROTOK = ( SMALL*AAQQ ).LE.AAPP
                   1182:                               END IF
                   1183:                               IF( AAPP.LT.( BIG / AAQQ ) ) THEN
                   1184:                                  AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
1.6       bertrand 1185:      $                                  q ), 1 )*WORK( p )*WORK( q ) /
                   1186:      $                                  AAQQ ) / AAPP
1.1       bertrand 1187:                               ELSE
                   1188:                                  CALL DCOPY( M, A( 1, p ), 1,
1.6       bertrand 1189:      $                                       WORK( N+1 ), 1 )
1.1       bertrand 1190:                                  CALL DLASCL( 'G', 0, 0, AAPP,
1.6       bertrand 1191:      $                                        WORK( p ), M, 1,
                   1192:      $                                        WORK( N+1 ), LDA, IERR )
1.1       bertrand 1193:                                  AAPQ = DDOT( M, WORK( N+1 ), 1,
1.6       bertrand 1194:      $                                  A( 1, q ), 1 )*WORK( q ) / AAQQ
1.1       bertrand 1195:                               END IF
                   1196:                            ELSE
                   1197:                               IF( AAPP.GE.AAQQ ) THEN
                   1198:                                  ROTOK = AAPP.LE.( AAQQ / SMALL )
                   1199:                               ELSE
                   1200:                                  ROTOK = AAQQ.LE.( AAPP / SMALL )
                   1201:                               END IF
                   1202:                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
                   1203:                                  AAPQ = ( DDOT( M, A( 1, p ), 1, A( 1,
1.6       bertrand 1204:      $                                  q ), 1 )*WORK( p )*WORK( q ) /
                   1205:      $                                  AAQQ ) / AAPP
1.1       bertrand 1206:                               ELSE
                   1207:                                  CALL DCOPY( M, A( 1, q ), 1,
1.6       bertrand 1208:      $                                       WORK( N+1 ), 1 )
1.1       bertrand 1209:                                  CALL DLASCL( 'G', 0, 0, AAQQ,
1.6       bertrand 1210:      $                                        WORK( q ), M, 1,
                   1211:      $                                        WORK( N+1 ), LDA, IERR )
1.1       bertrand 1212:                                  AAPQ = DDOT( M, WORK( N+1 ), 1,
1.6       bertrand 1213:      $                                  A( 1, p ), 1 )*WORK( p ) / AAPP
1.1       bertrand 1214:                               END IF
                   1215:                            END IF
                   1216: *
                   1217:                            MXAAPQ = DMAX1( MXAAPQ, DABS( AAPQ ) )
                   1218: *
                   1219: *        TO rotate or NOT to rotate, THAT is the question ...
                   1220: *
                   1221:                            IF( DABS( AAPQ ).GT.TOL ) THEN
                   1222:                               NOTROT = 0
                   1223: *[RTD]      ROTATED  = ROTATED + 1
                   1224:                               PSKIPPED = 0
                   1225:                               ISWROT = ISWROT + 1
                   1226: *
                   1227:                               IF( ROTOK ) THEN
                   1228: *
                   1229:                                  AQOAP = AAQQ / AAPP
                   1230:                                  APOAQ = AAPP / AAQQ
1.6       bertrand 1231:                                  THETA = -HALF*DABS(AQOAP-APOAQ)/AAPQ
1.1       bertrand 1232:                                  IF( AAQQ.GT.AAPP0 )THETA = -THETA
                   1233: *
                   1234:                                  IF( DABS( THETA ).GT.BIGTHETA ) THEN
                   1235:                                     T = HALF / THETA
                   1236:                                     FASTR( 3 ) = T*WORK( p ) / WORK( q )
                   1237:                                     FASTR( 4 ) = -T*WORK( q ) /
1.6       bertrand 1238:      $                                           WORK( p )
1.1       bertrand 1239:                                     CALL DROTM( M, A( 1, p ), 1,
1.6       bertrand 1240:      $                                          A( 1, q ), 1, FASTR )
1.1       bertrand 1241:                                     IF( RSVEC )CALL DROTM( MVL,
1.6       bertrand 1242:      $                                              V( 1, p ), 1,
                   1243:      $                                              V( 1, q ), 1,
                   1244:      $                                              FASTR )
1.1       bertrand 1245:                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
1.6       bertrand 1246:      $                                         ONE+T*APOAQ*AAPQ ) )
1.1       bertrand 1247:                                     AAPP = AAPP*DSQRT( DMAX1( ZERO,
1.6       bertrand 1248:      $                                     ONE-T*AQOAP*AAPQ ) )
1.1       bertrand 1249:                                     MXSINJ = DMAX1( MXSINJ, DABS( T ) )
                   1250:                                  ELSE
                   1251: *
                   1252: *                 .. choose correct signum for THETA and rotate
                   1253: *
                   1254:                                     THSIGN = -DSIGN( ONE, AAPQ )
                   1255:                                     IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN
                   1256:                                     T = ONE / ( THETA+THSIGN*
1.6       bertrand 1257:      $                                  DSQRT( ONE+THETA*THETA ) )
1.1       bertrand 1258:                                     CS = DSQRT( ONE / ( ONE+T*T ) )
                   1259:                                     SN = T*CS
                   1260:                                     MXSINJ = DMAX1( MXSINJ, DABS( SN ) )
                   1261:                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
1.6       bertrand 1262:      $                                         ONE+T*APOAQ*AAPQ ) )
1.4       bertrand 1263:                                     AAPP = AAPP*DSQRT( DMAX1( ZERO, 
1.6       bertrand 1264:      $                                     ONE-T*AQOAP*AAPQ ) )
1.1       bertrand 1265: *
                   1266:                                     APOAQ = WORK( p ) / WORK( q )
                   1267:                                     AQOAP = WORK( q ) / WORK( p )
                   1268:                                     IF( WORK( p ).GE.ONE ) THEN
                   1269: *
                   1270:                                        IF( WORK( q ).GE.ONE ) THEN
                   1271:                                           FASTR( 3 ) = T*APOAQ
                   1272:                                           FASTR( 4 ) = -T*AQOAP
                   1273:                                           WORK( p ) = WORK( p )*CS
                   1274:                                           WORK( q ) = WORK( q )*CS
                   1275:                                           CALL DROTM( M, A( 1, p ), 1,
1.6       bertrand 1276:      $                                                A( 1, q ), 1,
                   1277:      $                                                FASTR )
1.1       bertrand 1278:                                           IF( RSVEC )CALL DROTM( MVL,
1.6       bertrand 1279:      $                                        V( 1, p ), 1, V( 1, q ),
                   1280:      $                                        1, FASTR )
1.1       bertrand 1281:                                        ELSE
                   1282:                                           CALL DAXPY( M, -T*AQOAP,
1.6       bertrand 1283:      $                                                A( 1, q ), 1,
                   1284:      $                                                A( 1, p ), 1 )
1.1       bertrand 1285:                                           CALL DAXPY( M, CS*SN*APOAQ,
1.6       bertrand 1286:      $                                                A( 1, p ), 1,
                   1287:      $                                                A( 1, q ), 1 )
1.1       bertrand 1288:                                           IF( RSVEC ) THEN
                   1289:                                              CALL DAXPY( MVL, -T*AQOAP,
1.6       bertrand 1290:      $                                                   V( 1, q ), 1,
                   1291:      $                                                   V( 1, p ), 1 )
1.1       bertrand 1292:                                              CALL DAXPY( MVL,
1.6       bertrand 1293:      $                                                   CS*SN*APOAQ,
                   1294:      $                                                   V( 1, p ), 1,
                   1295:      $                                                   V( 1, q ), 1 )
1.1       bertrand 1296:                                           END IF
                   1297:                                           WORK( p ) = WORK( p )*CS
                   1298:                                           WORK( q ) = WORK( q ) / CS
                   1299:                                        END IF
                   1300:                                     ELSE
                   1301:                                        IF( WORK( q ).GE.ONE ) THEN
                   1302:                                           CALL DAXPY( M, T*APOAQ,
1.6       bertrand 1303:      $                                                A( 1, p ), 1,
                   1304:      $                                                A( 1, q ), 1 )
1.1       bertrand 1305:                                           CALL DAXPY( M, -CS*SN*AQOAP,
1.6       bertrand 1306:      $                                                A( 1, q ), 1,
                   1307:      $                                                A( 1, p ), 1 )
1.1       bertrand 1308:                                           IF( RSVEC ) THEN
                   1309:                                              CALL DAXPY( MVL, T*APOAQ,
1.6       bertrand 1310:      $                                                   V( 1, p ), 1,
                   1311:      $                                                   V( 1, q ), 1 )
1.1       bertrand 1312:                                              CALL DAXPY( MVL,
1.6       bertrand 1313:      $                                                   -CS*SN*AQOAP,
                   1314:      $                                                   V( 1, q ), 1,
                   1315:      $                                                   V( 1, p ), 1 )
1.1       bertrand 1316:                                           END IF
                   1317:                                           WORK( p ) = WORK( p ) / CS
                   1318:                                           WORK( q ) = WORK( q )*CS
                   1319:                                        ELSE
                   1320:                                           IF( WORK( p ).GE.WORK( q ) )
1.6       bertrand 1321:      $                                        THEN
1.1       bertrand 1322:                                              CALL DAXPY( M, -T*AQOAP,
1.6       bertrand 1323:      $                                                   A( 1, q ), 1,
                   1324:      $                                                   A( 1, p ), 1 )
1.1       bertrand 1325:                                              CALL DAXPY( M, CS*SN*APOAQ,
1.6       bertrand 1326:      $                                                   A( 1, p ), 1,
                   1327:      $                                                   A( 1, q ), 1 )
1.1       bertrand 1328:                                              WORK( p ) = WORK( p )*CS
                   1329:                                              WORK( q ) = WORK( q ) / CS
                   1330:                                              IF( RSVEC ) THEN
                   1331:                                                 CALL DAXPY( MVL,
1.6       bertrand 1332:      $                                               -T*AQOAP,
                   1333:      $                                               V( 1, q ), 1,
                   1334:      $                                               V( 1, p ), 1 )
1.1       bertrand 1335:                                                 CALL DAXPY( MVL,
1.6       bertrand 1336:      $                                               CS*SN*APOAQ,
                   1337:      $                                               V( 1, p ), 1,
                   1338:      $                                               V( 1, q ), 1 )
1.1       bertrand 1339:                                              END IF
                   1340:                                           ELSE
                   1341:                                              CALL DAXPY( M, T*APOAQ,
1.6       bertrand 1342:      $                                                   A( 1, p ), 1,
                   1343:      $                                                   A( 1, q ), 1 )
1.1       bertrand 1344:                                              CALL DAXPY( M,
1.6       bertrand 1345:      $                                                   -CS*SN*AQOAP,
                   1346:      $                                                   A( 1, q ), 1,
                   1347:      $                                                   A( 1, p ), 1 )
1.1       bertrand 1348:                                              WORK( p ) = WORK( p ) / CS
                   1349:                                              WORK( q ) = WORK( q )*CS
                   1350:                                              IF( RSVEC ) THEN
                   1351:                                                 CALL DAXPY( MVL,
1.6       bertrand 1352:      $                                               T*APOAQ, V( 1, p ),
                   1353:      $                                               1, V( 1, q ), 1 )
1.1       bertrand 1354:                                                 CALL DAXPY( MVL,
1.6       bertrand 1355:      $                                               -CS*SN*AQOAP,
                   1356:      $                                               V( 1, q ), 1,
                   1357:      $                                               V( 1, p ), 1 )
1.1       bertrand 1358:                                              END IF
                   1359:                                           END IF
                   1360:                                        END IF
                   1361:                                     END IF
                   1362:                                  END IF
                   1363: *
                   1364:                               ELSE
                   1365:                                  IF( AAPP.GT.AAQQ ) THEN
                   1366:                                     CALL DCOPY( M, A( 1, p ), 1,
1.6       bertrand 1367:      $                                          WORK( N+1 ), 1 )
1.1       bertrand 1368:                                     CALL DLASCL( 'G', 0, 0, AAPP, ONE,
1.6       bertrand 1369:      $                                           M, 1, WORK( N+1 ), LDA,
                   1370:      $                                           IERR )
1.1       bertrand 1371:                                     CALL DLASCL( 'G', 0, 0, AAQQ, ONE,
1.6       bertrand 1372:      $                                           M, 1, A( 1, q ), LDA,
                   1373:      $                                           IERR )
1.1       bertrand 1374:                                     TEMP1 = -AAPQ*WORK( p ) / WORK( q )
                   1375:                                     CALL DAXPY( M, TEMP1, WORK( N+1 ),
1.6       bertrand 1376:      $                                          1, A( 1, q ), 1 )
1.1       bertrand 1377:                                     CALL DLASCL( 'G', 0, 0, ONE, AAQQ,
1.6       bertrand 1378:      $                                           M, 1, A( 1, q ), LDA,
                   1379:      $                                           IERR )
1.1       bertrand 1380:                                     SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO,
1.6       bertrand 1381:      $                                         ONE-AAPQ*AAPQ ) )
1.1       bertrand 1382:                                     MXSINJ = DMAX1( MXSINJ, SFMIN )
                   1383:                                  ELSE
                   1384:                                     CALL DCOPY( M, A( 1, q ), 1,
1.6       bertrand 1385:      $                                          WORK( N+1 ), 1 )
1.1       bertrand 1386:                                     CALL DLASCL( 'G', 0, 0, AAQQ, ONE,
1.6       bertrand 1387:      $                                           M, 1, WORK( N+1 ), LDA,
                   1388:      $                                           IERR )
1.1       bertrand 1389:                                     CALL DLASCL( 'G', 0, 0, AAPP, ONE,
1.6       bertrand 1390:      $                                           M, 1, A( 1, p ), LDA,
                   1391:      $                                           IERR )
1.1       bertrand 1392:                                     TEMP1 = -AAPQ*WORK( q ) / WORK( p )
                   1393:                                     CALL DAXPY( M, TEMP1, WORK( N+1 ),
1.6       bertrand 1394:      $                                          1, A( 1, p ), 1 )
1.1       bertrand 1395:                                     CALL DLASCL( 'G', 0, 0, ONE, AAPP,
1.6       bertrand 1396:      $                                           M, 1, A( 1, p ), LDA,
                   1397:      $                                           IERR )
1.1       bertrand 1398:                                     SVA( p ) = AAPP*DSQRT( DMAX1( ZERO,
1.6       bertrand 1399:      $                                         ONE-AAPQ*AAPQ ) )
1.1       bertrand 1400:                                     MXSINJ = DMAX1( MXSINJ, SFMIN )
                   1401:                                  END IF
                   1402:                               END IF
                   1403: *           END IF ROTOK THEN ... ELSE
                   1404: *
                   1405: *           In the case of cancellation in updating SVA(q)
                   1406: *           .. recompute SVA(q)
                   1407:                               IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
1.6       bertrand 1408:      $                            THEN
1.1       bertrand 1409:                                  IF( ( AAQQ.LT.ROOTBIG ) .AND.
1.6       bertrand 1410:      $                               ( AAQQ.GT.ROOTSFMIN ) ) THEN
1.1       bertrand 1411:                                     SVA( q ) = DNRM2( M, A( 1, q ), 1 )*
1.6       bertrand 1412:      $                                         WORK( q )
1.1       bertrand 1413:                                  ELSE
                   1414:                                     T = ZERO
1.4       bertrand 1415:                                     AAQQ = ONE
1.1       bertrand 1416:                                     CALL DLASSQ( M, A( 1, q ), 1, T,
1.6       bertrand 1417:      $                                           AAQQ )
1.1       bertrand 1418:                                     SVA( q ) = T*DSQRT( AAQQ )*WORK( q )
                   1419:                                  END IF
                   1420:                               END IF
                   1421:                               IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN
                   1422:                                  IF( ( AAPP.LT.ROOTBIG ) .AND.
1.6       bertrand 1423:      $                               ( AAPP.GT.ROOTSFMIN ) ) THEN
1.1       bertrand 1424:                                     AAPP = DNRM2( M, A( 1, p ), 1 )*
1.6       bertrand 1425:      $                                     WORK( p )
1.1       bertrand 1426:                                  ELSE
                   1427:                                     T = ZERO
1.4       bertrand 1428:                                     AAPP = ONE
1.1       bertrand 1429:                                     CALL DLASSQ( M, A( 1, p ), 1, T,
1.6       bertrand 1430:      $                                           AAPP )
1.1       bertrand 1431:                                     AAPP = T*DSQRT( AAPP )*WORK( p )
                   1432:                                  END IF
                   1433:                                  SVA( p ) = AAPP
                   1434:                               END IF
                   1435: *              end of OK rotation
                   1436:                            ELSE
                   1437:                               NOTROT = NOTROT + 1
                   1438: *[RTD]      SKIPPED  = SKIPPED  + 1
                   1439:                               PSKIPPED = PSKIPPED + 1
                   1440:                               IJBLSK = IJBLSK + 1
                   1441:                            END IF
                   1442:                         ELSE
                   1443:                            NOTROT = NOTROT + 1
                   1444:                            PSKIPPED = PSKIPPED + 1
                   1445:                            IJBLSK = IJBLSK + 1
                   1446:                         END IF
                   1447: *
                   1448:                         IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) )
1.6       bertrand 1449:      $                      THEN
1.1       bertrand 1450:                            SVA( p ) = AAPP
                   1451:                            NOTROT = 0
                   1452:                            GO TO 2011
                   1453:                         END IF
                   1454:                         IF( ( i.LE.SWBAND ) .AND.
1.6       bertrand 1455:      $                      ( PSKIPPED.GT.ROWSKIP ) ) THEN
1.1       bertrand 1456:                            AAPP = -AAPP
                   1457:                            NOTROT = 0
                   1458:                            GO TO 2203
                   1459:                         END IF
                   1460: *
                   1461:  2200                CONTINUE
                   1462: *        end of the q-loop
                   1463:  2203                CONTINUE
                   1464: *
                   1465:                      SVA( p ) = AAPP
                   1466: *
                   1467:                   ELSE
                   1468: *
                   1469:                      IF( AAPP.EQ.ZERO )NOTROT = NOTROT +
1.6       bertrand 1470:      $                   MIN0( jgl+KBL-1, N ) - jgl + 1
1.1       bertrand 1471:                      IF( AAPP.LT.ZERO )NOTROT = 0
                   1472: *
                   1473:                   END IF
                   1474: *
                   1475:  2100          CONTINUE
                   1476: *     end of the p-loop
                   1477:  2010       CONTINUE
                   1478: *     end of the jbc-loop
                   1479:  2011       CONTINUE
                   1480: *2011 bailed out of the jbc-loop
                   1481:             DO 2012 p = igl, MIN0( igl+KBL-1, N )
                   1482:                SVA( p ) = DABS( SVA( p ) )
                   1483:  2012       CONTINUE
                   1484: ***
                   1485:  2000    CONTINUE
                   1486: *2000 :: end of the ibr-loop
                   1487: *
                   1488: *     .. update SVA(N)
                   1489:          IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) )
1.6       bertrand 1490:      $       THEN
1.1       bertrand 1491:             SVA( N ) = DNRM2( M, A( 1, N ), 1 )*WORK( N )
                   1492:          ELSE
                   1493:             T = ZERO
1.4       bertrand 1494:             AAPP = ONE
1.1       bertrand 1495:             CALL DLASSQ( M, A( 1, N ), 1, T, AAPP )
                   1496:             SVA( N ) = T*DSQRT( AAPP )*WORK( N )
                   1497:          END IF
                   1498: *
                   1499: *     Additional steering devices
                   1500: *
                   1501:          IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.
1.6       bertrand 1502:      $       ( ISWROT.LE.N ) ) )SWBAND = i
1.1       bertrand 1503: *
                   1504:          IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.DSQRT( DBLE( N ) )*
1.6       bertrand 1505:      $       TOL ) .AND. ( DBLE( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
1.1       bertrand 1506:             GO TO 1994
                   1507:          END IF
                   1508: *
                   1509:          IF( NOTROT.GE.EMPTSW )GO TO 1994
                   1510: *
                   1511:  1993 CONTINUE
                   1512: *     end i=1:NSWEEP loop
                   1513: *
                   1514: * #:( Reaching this point means that the procedure has not converged.
                   1515:       INFO = NSWEEP - 1
                   1516:       GO TO 1995
                   1517: *
                   1518:  1994 CONTINUE
                   1519: * #:) Reaching this point means numerical convergence after the i-th
                   1520: *     sweep.
                   1521: *
                   1522:       INFO = 0
                   1523: * #:) INFO = 0 confirms successful iterations.
                   1524:  1995 CONTINUE
                   1525: *
                   1526: *     Sort the singular values and find how many are above
                   1527: *     the underflow threshold.
                   1528: *
                   1529:       N2 = 0
                   1530:       N4 = 0
                   1531:       DO 5991 p = 1, N - 1
                   1532:          q = IDAMAX( N-p+1, SVA( p ), 1 ) + p - 1
                   1533:          IF( p.NE.q ) THEN
                   1534:             TEMP1 = SVA( p )
                   1535:             SVA( p ) = SVA( q )
                   1536:             SVA( q ) = TEMP1
                   1537:             TEMP1 = WORK( p )
                   1538:             WORK( p ) = WORK( q )
                   1539:             WORK( q ) = TEMP1
                   1540:             CALL DSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
                   1541:             IF( RSVEC )CALL DSWAP( MVL, V( 1, p ), 1, V( 1, q ), 1 )
                   1542:          END IF
                   1543:          IF( SVA( p ).NE.ZERO ) THEN
                   1544:             N4 = N4 + 1
1.4       bertrand 1545:             IF( SVA( p )*SKL.GT.SFMIN )N2 = N2 + 1
1.1       bertrand 1546:          END IF
                   1547:  5991 CONTINUE
                   1548:       IF( SVA( N ).NE.ZERO ) THEN
                   1549:          N4 = N4 + 1
1.4       bertrand 1550:          IF( SVA( N )*SKL.GT.SFMIN )N2 = N2 + 1
1.1       bertrand 1551:       END IF
                   1552: *
                   1553: *     Normalize the left singular vectors.
                   1554: *
                   1555:       IF( LSVEC .OR. UCTOL ) THEN
                   1556:          DO 1998 p = 1, N2
                   1557:             CALL DSCAL( M, WORK( p ) / SVA( p ), A( 1, p ), 1 )
                   1558:  1998    CONTINUE
                   1559:       END IF
                   1560: *
                   1561: *     Scale the product of Jacobi rotations (assemble the fast rotations).
                   1562: *
                   1563:       IF( RSVEC ) THEN
                   1564:          IF( APPLV ) THEN
                   1565:             DO 2398 p = 1, N
                   1566:                CALL DSCAL( MVL, WORK( p ), V( 1, p ), 1 )
                   1567:  2398       CONTINUE
                   1568:          ELSE
                   1569:             DO 2399 p = 1, N
                   1570:                TEMP1 = ONE / DNRM2( MVL, V( 1, p ), 1 )
                   1571:                CALL DSCAL( MVL, TEMP1, V( 1, p ), 1 )
                   1572:  2399       CONTINUE
                   1573:          END IF
                   1574:       END IF
                   1575: *
                   1576: *     Undo scaling, if necessary (and possible).
1.4       bertrand 1577:       IF( ( ( SKL.GT.ONE ) .AND. ( SVA( 1 ).LT.( BIG /
1.6       bertrand 1578:      $    SKL) ) ) .OR. ( ( SKL.LT.ONE ) .AND. ( SVA( N2 ).GT.
                   1579:      $    ( SFMIN / SKL) ) ) ) THEN
1.1       bertrand 1580:          DO 2400 p = 1, N
1.4       bertrand 1581:             SVA( p ) = SKL*SVA( p )
1.1       bertrand 1582:  2400    CONTINUE
1.4       bertrand 1583:          SKL= ONE
1.1       bertrand 1584:       END IF
                   1585: *
1.4       bertrand 1586:       WORK( 1 ) = SKL
                   1587: *     The singular values of A are SKL*SVA(1:N). If SKL.NE.ONE
1.1       bertrand 1588: *     then some of the singular values may overflow or underflow and
                   1589: *     the spectrum is given in this factored representation.
                   1590: *
                   1591:       WORK( 2 ) = DBLE( N4 )
                   1592: *     N4 is the number of computed nonzero singular values of A.
                   1593: *
                   1594:       WORK( 3 ) = DBLE( N2 )
                   1595: *     N2 is the number of singular values of A greater than SFMIN.
                   1596: *     If N2<N, SVA(N2:N) contains ZEROS and/or denormalized numbers
                   1597: *     that may carry some information.
                   1598: *
                   1599:       WORK( 4 ) = DBLE( i )
                   1600: *     i is the index of the last sweep before declaring convergence.
                   1601: *
                   1602:       WORK( 5 ) = MXAAPQ
                   1603: *     MXAAPQ is the largest absolute value of scaled pivots in the
                   1604: *     last sweep
                   1605: *
                   1606:       WORK( 6 ) = MXSINJ
                   1607: *     MXSINJ is the largest absolute value of the sines of Jacobi angles
                   1608: *     in the last sweep
                   1609: *
                   1610:       RETURN
                   1611: *     ..
                   1612: *     .. END OF DGESVJ
                   1613: *     ..
                   1614:       END

CVSweb interface <joel.bertrand@systella.fr>