Annotation of rpl/lapack/lapack/zgesvj.f, revision 1.5

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

CVSweb interface <joel.bertrand@systella.fr>