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

1.2       bertrand    1: *> \brief \b ZGSVJ1 pre-processor for the routine zgesvj, applies Jacobi rotations targeting only particular pivots.
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 ZGSVJ1 + dependencies
                     10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgsvj1.f">
                     11: *> [TGZ]</a>
                     12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgsvj1.f">
                     13: *> [ZIP]</a>
                     14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgsvj1.f">
1.1       bertrand   15: *> [TXT]</a>
1.4       bertrand   16: *> \endhtmlonly
1.1       bertrand   17: *
                     18: *  Definition:
                     19: *  ===========
                     20: *
                     21: *       SUBROUTINE ZGSVJ1( JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV,
                     22: *                          EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )
1.4       bertrand   23: *
1.1       bertrand   24: *       .. Scalar Arguments ..
                     25: *       DOUBLE PRECISION   EPS, SFMIN, TOL
                     26: *       INTEGER            INFO, LDA, LDV, LWORK, M, MV, N, N1, NSWEEP
                     27: *       CHARACTER*1        JOBV
                     28: *       ..
                     29: *       .. Array Arguments ..
1.4       bertrand   30: *       COMPLEX*16         A( LDA, * ), D( N ), V( LDV, * ), WORK( LWORK )
                     31: *       DOUBLE PRECISION   SVA( N )
1.1       bertrand   32: *       ..
1.4       bertrand   33: *
1.1       bertrand   34: *
                     35: *> \par Purpose:
                     36: *  =============
                     37: *>
                     38: *> \verbatim
                     39: *>
                     40: *> ZGSVJ1 is called from ZGESVJ as a pre-processor and that is its main
                     41: *> purpose. It applies Jacobi rotations in the same way as ZGESVJ does, but
                     42: *> it targets only particular pivots and it does not check convergence
                     43: *> (stopping criterion). Few tunning parameters (marked by [TP]) are
                     44: *> available for the implementer.
                     45: *>
                     46: *> Further Details
                     47: *> ~~~~~~~~~~~~~~~
                     48: *> ZGSVJ1 applies few sweeps of Jacobi rotations in the column space of
                     49: *> the input M-by-N matrix A. The pivot pairs are taken from the (1,2)
                     50: *> off-diagonal block in the corresponding N-by-N Gram matrix A^T * A. The
                     51: *> block-entries (tiles) of the (1,2) off-diagonal block are marked by the
                     52: *> [x]'s in the following scheme:
                     53: *>
                     54: *>    | *  *  * [x] [x] [x]|
                     55: *>    | *  *  * [x] [x] [x]|    Row-cycling in the nblr-by-nblc [x] blocks.
                     56: *>    | *  *  * [x] [x] [x]|    Row-cyclic pivoting inside each [x] block.
                     57: *>    |[x] [x] [x] *  *  * |
                     58: *>    |[x] [x] [x] *  *  * |
                     59: *>    |[x] [x] [x] *  *  * |
                     60: *>
                     61: *> In terms of the columns of A, the first N1 columns are rotated 'against'
                     62: *> the remaining N-N1 columns, trying to increase the angle between the
                     63: *> corresponding subspaces. The off-diagonal block is N1-by(N-N1) and it is
                     64: *> tiled using quadratic tiles of side KBL. Here, KBL is a tunning parmeter.
                     65: *> The number of sweeps is given in NSWEEP and the orthogonality threshold
                     66: *> is given in TOL.
                     67: *> \endverbatim
                     68: *
                     69: *  Arguments:
                     70: *  ==========
                     71: *
                     72: *> \param[in] JOBV
                     73: *> \verbatim
                     74: *>          JOBV is CHARACTER*1
                     75: *>          Specifies whether the output from this procedure is used
                     76: *>          to compute the matrix V:
                     77: *>          = 'V': the product of the Jacobi rotations is accumulated
                     78: *>                 by postmulyiplying the N-by-N array V.
                     79: *>                (See the description of V.)
                     80: *>          = 'A': the product of the Jacobi rotations is accumulated
                     81: *>                 by postmulyiplying the MV-by-N array V.
                     82: *>                (See the descriptions of MV and V.)
                     83: *>          = 'N': the Jacobi rotations are not accumulated.
                     84: *> \endverbatim
                     85: *>
                     86: *> \param[in] M
                     87: *> \verbatim
                     88: *>          M is INTEGER
                     89: *>          The number of rows of the input matrix A.  M >= 0.
                     90: *> \endverbatim
                     91: *>
                     92: *> \param[in] N
                     93: *> \verbatim
                     94: *>          N is INTEGER
                     95: *>          The number of columns of the input matrix A.
                     96: *>          M >= N >= 0.
                     97: *> \endverbatim
                     98: *>
                     99: *> \param[in] N1
                    100: *> \verbatim
                    101: *>          N1 is INTEGER
                    102: *>          N1 specifies the 2 x 2 block partition, the first N1 columns are
                    103: *>          rotated 'against' the remaining N-N1 columns of A.
                    104: *> \endverbatim
                    105: *>
                    106: *> \param[in,out] A
                    107: *> \verbatim
1.2       bertrand  108: *>          A is COMPLEX*16 array, dimension (LDA,N)
1.1       bertrand  109: *>          On entry, M-by-N matrix A, such that A*diag(D) represents
                    110: *>          the input matrix.
                    111: *>          On exit,
                    112: *>          A_onexit * D_onexit represents the input matrix A*diag(D)
                    113: *>          post-multiplied by a sequence of Jacobi rotations, where the
                    114: *>          rotation threshold and the total number of sweeps are given in
                    115: *>          TOL and NSWEEP, respectively.
                    116: *>          (See the descriptions of N1, D, TOL and NSWEEP.)
                    117: *> \endverbatim
                    118: *>
                    119: *> \param[in] LDA
                    120: *> \verbatim
                    121: *>          LDA is INTEGER
                    122: *>          The leading dimension of the array A.  LDA >= max(1,M).
                    123: *> \endverbatim
                    124: *>
                    125: *> \param[in,out] D
                    126: *> \verbatim
1.2       bertrand  127: *>          D is COMPLEX*16 array, dimension (N)
1.1       bertrand  128: *>          The array D accumulates the scaling factors from the fast scaled
                    129: *>          Jacobi rotations.
                    130: *>          On entry, A*diag(D) represents the input matrix.
                    131: *>          On exit, A_onexit*diag(D_onexit) represents the input matrix
                    132: *>          post-multiplied by a sequence of Jacobi rotations, where the
                    133: *>          rotation threshold and the total number of sweeps are given in
                    134: *>          TOL and NSWEEP, respectively.
                    135: *>          (See the descriptions of N1, A, TOL and NSWEEP.)
                    136: *> \endverbatim
                    137: *>
                    138: *> \param[in,out] SVA
                    139: *> \verbatim
                    140: *>          SVA is DOUBLE PRECISION array, dimension (N)
                    141: *>          On entry, SVA contains the Euclidean norms of the columns of
                    142: *>          the matrix A*diag(D).
                    143: *>          On exit, SVA contains the Euclidean norms of the columns of
                    144: *>          the matrix onexit*diag(D_onexit).
                    145: *> \endverbatim
                    146: *>
                    147: *> \param[in] MV
                    148: *> \verbatim
                    149: *>          MV is INTEGER
                    150: *>          If JOBV .EQ. 'A', then MV rows of V are post-multipled by a
                    151: *>                           sequence of Jacobi rotations.
                    152: *>          If JOBV = 'N',   then MV is not referenced.
                    153: *> \endverbatim
                    154: *>
                    155: *> \param[in,out] V
                    156: *> \verbatim
1.2       bertrand  157: *>          V is COMPLEX*16 array, dimension (LDV,N)
1.1       bertrand  158: *>          If JOBV .EQ. 'V' then N rows of V are post-multipled by a
                    159: *>                           sequence of Jacobi rotations.
                    160: *>          If JOBV .EQ. 'A' then MV rows of V are post-multipled by a
                    161: *>                           sequence of Jacobi rotations.
                    162: *>          If JOBV = 'N',   then V is not referenced.
                    163: *> \endverbatim
                    164: *>
                    165: *> \param[in] LDV
                    166: *> \verbatim
                    167: *>          LDV is INTEGER
                    168: *>          The leading dimension of the array V,  LDV >= 1.
                    169: *>          If JOBV = 'V', LDV .GE. N.
                    170: *>          If JOBV = 'A', LDV .GE. MV.
                    171: *> \endverbatim
                    172: *>
                    173: *> \param[in] EPS
                    174: *> \verbatim
                    175: *>          EPS is DOUBLE PRECISION
                    176: *>          EPS = DLAMCH('Epsilon')
                    177: *> \endverbatim
                    178: *>
                    179: *> \param[in] SFMIN
                    180: *> \verbatim
                    181: *>          SFMIN is DOUBLE PRECISION
                    182: *>          SFMIN = DLAMCH('Safe Minimum')
                    183: *> \endverbatim
                    184: *>
                    185: *> \param[in] TOL
                    186: *> \verbatim
                    187: *>          TOL is DOUBLE PRECISION
                    188: *>          TOL is the threshold for Jacobi rotations. For a pair
                    189: *>          A(:,p), A(:,q) of pivot columns, the Jacobi rotation is
                    190: *>          applied only if ABS(COS(angle(A(:,p),A(:,q)))) .GT. TOL.
                    191: *> \endverbatim
                    192: *>
                    193: *> \param[in] NSWEEP
                    194: *> \verbatim
                    195: *>          NSWEEP is INTEGER
                    196: *>          NSWEEP is the number of sweeps of Jacobi rotations to be
                    197: *>          performed.
                    198: *> \endverbatim
                    199: *>
                    200: *> \param[out] WORK
                    201: *> \verbatim
1.2       bertrand  202: *>          WORK is COMPLEX*16 array, dimension (LWORK)
1.1       bertrand  203: *> \endverbatim
                    204: *>
                    205: *> \param[in] LWORK
                    206: *> \verbatim
                    207: *>          LWORK is INTEGER
                    208: *>          LWORK is the dimension of WORK. LWORK .GE. M.
                    209: *> \endverbatim
                    210: *>
                    211: *> \param[out] INFO
                    212: *> \verbatim
                    213: *>          INFO is INTEGER
                    214: *>          = 0 : successful exit.
                    215: *>          < 0 : if INFO = -i, then the i-th argument had an illegal value
                    216: *> \endverbatim
                    217: *
                    218: *  Authors:
                    219: *  ========
                    220: *
1.4       bertrand  221: *> \author Univ. of Tennessee
                    222: *> \author Univ. of California Berkeley
                    223: *> \author Univ. of Colorado Denver
                    224: *> \author NAG Ltd.
1.1       bertrand  225: *
1.2       bertrand  226: *> \date June 2016
1.1       bertrand  227: *
                    228: *> \ingroup complex16OTHERcomputational
                    229: *
1.4       bertrand  230: *> \par Contributor:
1.1       bertrand  231: *  ==================
                    232: *>
1.4       bertrand  233: *> Zlatko Drmac (Zagreb, Croatia)
1.1       bertrand  234: *
                    235: *  =====================================================================
                    236:       SUBROUTINE ZGSVJ1( JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV,
                    237:      $                   EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )
                    238: *
1.4       bertrand  239: *  -- LAPACK computational routine (version 3.7.0) --
1.1       bertrand  240: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    241: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.2       bertrand  242: *     June 2016
1.1       bertrand  243: *
1.4       bertrand  244:       IMPLICIT NONE
1.1       bertrand  245: *     .. Scalar Arguments ..
                    246:       DOUBLE PRECISION   EPS, SFMIN, TOL
                    247:       INTEGER            INFO, LDA, LDV, LWORK, M, MV, N, N1, NSWEEP
                    248:       CHARACTER*1        JOBV
                    249: *     ..
                    250: *     .. Array Arguments ..
                    251:       COMPLEX*16         A( LDA, * ), D( N ), V( LDV, * ), WORK( LWORK )
1.4       bertrand  252:       DOUBLE PRECISION   SVA( N )
1.1       bertrand  253: *     ..
                    254: *
                    255: *  =====================================================================
                    256: *
                    257: *     .. Local Parameters ..
                    258:       DOUBLE PRECISION   ZERO, HALF, ONE
                    259:       PARAMETER          ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0)
                    260: *     ..
                    261: *     .. Local Scalars ..
                    262:       COMPLEX*16         AAPQ, OMPQ
                    263:       DOUBLE PRECISION   AAPP, AAPP0, AAPQ1, AAQQ, APOAQ, AQOAP, BIG,
1.4       bertrand  264:      $                   BIGTHETA, CS, MXAAPQ, MXSINJ, ROOTBIG,
1.1       bertrand  265:      $                   ROOTEPS, ROOTSFMIN, ROOTTOL, SMALL, SN, T,
                    266:      $                   TEMP1, THETA, THSIGN
                    267:       INTEGER            BLSKIP, EMPTSW, i, ibr, igl, IERR, IJBLSK,
                    268:      $                   ISWROT, jbc, jgl, KBL, MVL, NOTROT, nblc, nblr,
                    269:      $                   p, PSKIPPED, q, ROWSKIP, SWBAND
                    270:       LOGICAL            APPLV, ROTOK, RSVEC
                    271: *     ..
                    272: *     ..
                    273: *     .. Intrinsic Functions ..
1.4       bertrand  274:       INTRINSIC          ABS, CONJG, MAX, DBLE, MIN, SIGN, SQRT
1.1       bertrand  275: *     ..
                    276: *     .. External Functions ..
                    277:       DOUBLE PRECISION   DZNRM2
                    278:       COMPLEX*16         ZDOTC
                    279:       INTEGER            IDAMAX
                    280:       LOGICAL            LSAME
                    281:       EXTERNAL           IDAMAX, LSAME, ZDOTC, DZNRM2
                    282: *     ..
                    283: *     .. External Subroutines ..
1.4       bertrand  284: *     .. from BLAS
1.1       bertrand  285:       EXTERNAL           ZCOPY, ZROT, ZSWAP
                    286: *     .. from LAPACK
                    287:       EXTERNAL           ZLASCL, ZLASSQ, XERBLA
                    288: *     ..
                    289: *     .. Executable Statements ..
                    290: *
                    291: *     Test the input parameters.
                    292: *
                    293:       APPLV = LSAME( JOBV, 'A' )
                    294:       RSVEC = LSAME( JOBV, 'V' )
                    295:       IF( .NOT.( RSVEC .OR. APPLV .OR. LSAME( JOBV, 'N' ) ) ) THEN
                    296:          INFO = -1
                    297:       ELSE IF( M.LT.0 ) THEN
                    298:          INFO = -2
                    299:       ELSE IF( ( N.LT.0 ) .OR. ( N.GT.M ) ) THEN
                    300:          INFO = -3
                    301:       ELSE IF( N1.LT.0 ) THEN
                    302:          INFO = -4
                    303:       ELSE IF( LDA.LT.M ) THEN
                    304:          INFO = -6
                    305:       ELSE IF( ( RSVEC.OR.APPLV ) .AND. ( MV.LT.0 ) ) THEN
                    306:          INFO = -9
1.4       bertrand  307:       ELSE IF( ( RSVEC.AND.( LDV.LT.N ) ).OR.
1.1       bertrand  308:      $         ( APPLV.AND.( LDV.LT.MV ) )  ) THEN
                    309:          INFO = -11
                    310:       ELSE IF( TOL.LE.EPS ) THEN
                    311:          INFO = -14
                    312:       ELSE IF( NSWEEP.LT.0 ) THEN
                    313:          INFO = -15
                    314:       ELSE IF( LWORK.LT.M ) THEN
                    315:          INFO = -17
                    316:       ELSE
                    317:          INFO = 0
                    318:       END IF
                    319: *
                    320: *     #:(
                    321:       IF( INFO.NE.0 ) THEN
                    322:          CALL XERBLA( 'ZGSVJ1', -INFO )
                    323:          RETURN
                    324:       END IF
                    325: *
                    326:       IF( RSVEC ) THEN
                    327:          MVL = N
                    328:       ELSE IF( APPLV ) THEN
                    329:          MVL = MV
                    330:       END IF
                    331:       RSVEC = RSVEC .OR. APPLV
                    332: 
1.4       bertrand  333:       ROOTEPS = SQRT( EPS )
                    334:       ROOTSFMIN = SQRT( SFMIN )
1.1       bertrand  335:       SMALL = SFMIN / EPS
                    336:       BIG = ONE / SFMIN
                    337:       ROOTBIG = ONE / ROOTSFMIN
1.4       bertrand  338: *     LARGE = BIG / SQRT( DBLE( M*N ) )
1.1       bertrand  339:       BIGTHETA = ONE / ROOTEPS
1.4       bertrand  340:       ROOTTOL = SQRT( TOL )
1.1       bertrand  341: *
                    342: *     .. Initialize the right singular vector matrix ..
                    343: *
                    344: *     RSVEC = LSAME( JOBV, 'Y' )
                    345: *
                    346:       EMPTSW = N1*( N-N1 )
                    347:       NOTROT = 0
                    348: *
                    349: *     .. Row-cyclic pivot strategy with de Rijk's pivoting ..
                    350: *
1.4       bertrand  351:       KBL = MIN( 8, N )
1.1       bertrand  352:       NBLR = N1 / KBL
                    353:       IF( ( NBLR*KBL ).NE.N1 )NBLR = NBLR + 1
                    354: 
                    355: *     .. the tiling is nblr-by-nblc [tiles]
                    356: 
                    357:       NBLC = ( N-N1 ) / KBL
                    358:       IF( ( NBLC*KBL ).NE.( N-N1 ) )NBLC = NBLC + 1
                    359:       BLSKIP = ( KBL**2 ) + 1
                    360: *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
                    361: 
1.4       bertrand  362:       ROWSKIP = MIN( 5, KBL )
1.1       bertrand  363: *[TP] ROWSKIP is a tuning parameter.
                    364:       SWBAND = 0
                    365: *[TP] SWBAND is a tuning parameter. It is meaningful and effective
                    366: *     if ZGESVJ is used as a computational routine in the preconditioned
                    367: *     Jacobi SVD algorithm ZGEJSV.
                    368: *
                    369: *
                    370: *     | *   *   * [x] [x] [x]|
                    371: *     | *   *   * [x] [x] [x]|    Row-cycling in the nblr-by-nblc [x] blocks.
                    372: *     | *   *   * [x] [x] [x]|    Row-cyclic pivoting inside each [x] block.
                    373: *     |[x] [x] [x] *   *   * |
                    374: *     |[x] [x] [x] *   *   * |
                    375: *     |[x] [x] [x] *   *   * |
                    376: *
                    377: *
                    378:       DO 1993 i = 1, NSWEEP
                    379: *
                    380: *     .. go go go ...
                    381: *
                    382:          MXAAPQ = ZERO
                    383:          MXSINJ = ZERO
                    384:          ISWROT = 0
                    385: *
                    386:          NOTROT = 0
                    387:          PSKIPPED = 0
                    388: *
                    389: *     Each sweep is unrolled using KBL-by-KBL tiles over the pivot pairs
                    390: *     1 <= p < q <= N. This is the first step toward a blocked implementation
                    391: *     of the rotations. New implementation, based on block transformations,
                    392: *     is under development.
                    393: *
                    394:          DO 2000 ibr = 1, NBLR
                    395: *
                    396:             igl = ( ibr-1 )*KBL + 1
                    397: *
                    398: 
                    399: *
                    400: * ... go to the off diagonal blocks
                    401: *
                    402:             igl = ( ibr-1 )*KBL + 1
                    403: *
                    404: *            DO 2010 jbc = ibr + 1, NBL
1.4       bertrand  405:             DO 2010 jbc = 1, NBLC
1.1       bertrand  406: *
                    407:                jgl = ( jbc-1 )*KBL + N1 + 1
                    408: *
                    409: *        doing the block at ( ibr, jbc )
                    410: *
                    411:                IJBLSK = 0
1.4       bertrand  412:                DO 2100 p = igl, MIN( igl+KBL-1, N1 )
1.1       bertrand  413: *
                    414:                   AAPP = SVA( p )
                    415:                   IF( AAPP.GT.ZERO ) THEN
                    416: *
                    417:                      PSKIPPED = 0
                    418: *
1.4       bertrand  419:                      DO 2200 q = jgl, MIN( jgl+KBL-1, N )
1.1       bertrand  420: *
                    421:                         AAQQ = SVA( q )
                    422:                         IF( AAQQ.GT.ZERO ) THEN
                    423:                            AAPP0 = AAPP
                    424: *
                    425: *     .. M x 2 Jacobi SVD ..
                    426: *
                    427: *        Safe Gram matrix computation
                    428: *
                    429:                            IF( AAQQ.GE.ONE ) THEN
                    430:                               IF( AAPP.GE.AAQQ ) THEN
                    431:                                  ROTOK = ( SMALL*AAPP ).LE.AAQQ
                    432:                               ELSE
                    433:                                  ROTOK = ( SMALL*AAQQ ).LE.AAPP
                    434:                               END IF
                    435:                               IF( AAPP.LT.( BIG / AAQQ ) ) THEN
1.4       bertrand  436:                                  AAPQ = ( ZDOTC( M, A( 1, p ), 1,
1.1       bertrand  437:      $                                  A( 1, q ), 1 ) / AAQQ ) / AAPP
                    438:                               ELSE
                    439:                                  CALL ZCOPY( M, A( 1, p ), 1,
                    440:      $                                       WORK, 1 )
                    441:                                  CALL ZLASCL( 'G', 0, 0, AAPP,
                    442:      $                                        ONE, M, 1,
                    443:      $                                        WORK, LDA, IERR )
                    444:                                  AAPQ = ZDOTC( M, WORK, 1,
                    445:      $                                  A( 1, q ), 1 ) / AAQQ
                    446:                               END IF
                    447:                            ELSE
                    448:                               IF( AAPP.GE.AAQQ ) THEN
                    449:                                  ROTOK = AAPP.LE.( AAQQ / SMALL )
                    450:                               ELSE
                    451:                                  ROTOK = AAQQ.LE.( AAPP / SMALL )
                    452:                               END IF
                    453:                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
1.4       bertrand  454:                                  AAPQ = ( ZDOTC( M, A( 1, p ), 1,
                    455:      $                                 A( 1, q ), 1 ) / MAX(AAQQ,AAPP) )
                    456:      $                                               / MIN(AAQQ,AAPP)
1.1       bertrand  457:                               ELSE
                    458:                                  CALL ZCOPY( M, A( 1, q ), 1,
                    459:      $                                       WORK, 1 )
                    460:                                  CALL ZLASCL( 'G', 0, 0, AAQQ,
                    461:      $                                        ONE, M, 1,
                    462:      $                                        WORK, LDA, IERR )
                    463:                                  AAPQ = ZDOTC( M, A( 1, p ), 1,
                    464:      $                                  WORK, 1 ) / AAPP
                    465:                               END IF
                    466:                            END IF
                    467: *
1.4       bertrand  468: *                           AAPQ = AAPQ * CONJG(CWORK(p))*CWORK(q)
1.1       bertrand  469:                            AAPQ1  = -ABS(AAPQ)
1.4       bertrand  470:                            MXAAPQ = MAX( MXAAPQ, -AAPQ1 )
1.1       bertrand  471: *
                    472: *        TO rotate or NOT to rotate, THAT is the question ...
                    473: *
                    474:                            IF( ABS( AAPQ1 ).GT.TOL ) THEN
1.4       bertrand  475:                               OMPQ = AAPQ / ABS(AAPQ)
1.1       bertrand  476:                               NOTROT = 0
                    477: *[RTD]      ROTATED  = ROTATED + 1
                    478:                               PSKIPPED = 0
                    479:                               ISWROT = ISWROT + 1
                    480: *
                    481:                               IF( ROTOK ) THEN
                    482: *
                    483:                                  AQOAP = AAQQ / AAPP
                    484:                                  APOAQ = AAPP / AAQQ
                    485:                                  THETA = -HALF*ABS( AQOAP-APOAQ )/ AAPQ1
                    486:                                  IF( AAQQ.GT.AAPP0 )THETA = -THETA
                    487: *
                    488:                                  IF( ABS( THETA ).GT.BIGTHETA ) THEN
                    489:                                     T  = HALF / THETA
1.4       bertrand  490:                                     CS = ONE
1.1       bertrand  491:                                     CALL ZROT( M, A(1,p), 1, A(1,q), 1,
1.4       bertrand  492:      $                                          CS, CONJG(OMPQ)*T )
1.1       bertrand  493:                                     IF( RSVEC ) THEN
1.4       bertrand  494:                                         CALL ZROT( MVL, V(1,p), 1,
                    495:      $                                  V(1,q), 1, CS, CONJG(OMPQ)*T )
1.1       bertrand  496:                                     END IF
1.4       bertrand  497:                                     SVA( q ) = AAQQ*SQRT( MAX( ZERO,
1.1       bertrand  498:      $                                         ONE+T*APOAQ*AAPQ1 ) )
1.4       bertrand  499:                                     AAPP = AAPP*SQRT( MAX( ZERO,
1.1       bertrand  500:      $                                     ONE-T*AQOAP*AAPQ1 ) )
1.4       bertrand  501:                                     MXSINJ = MAX( MXSINJ, ABS( T ) )
1.1       bertrand  502:                                  ELSE
                    503: *
                    504: *                 .. choose correct signum for THETA and rotate
                    505: *
1.4       bertrand  506:                                     THSIGN = -SIGN( ONE, AAPQ1 )
1.1       bertrand  507:                                     IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN
                    508:                                     T = ONE / ( THETA+THSIGN*
1.4       bertrand  509:      $                                  SQRT( ONE+THETA*THETA ) )
                    510:                                     CS = SQRT( ONE / ( ONE+T*T ) )
1.1       bertrand  511:                                     SN = T*CS
1.4       bertrand  512:                                     MXSINJ = MAX( MXSINJ, ABS( SN ) )
                    513:                                     SVA( q ) = AAQQ*SQRT( MAX( ZERO,
1.1       bertrand  514:      $                                         ONE+T*APOAQ*AAPQ1 ) )
1.4       bertrand  515:                                     AAPP = AAPP*SQRT( MAX( ZERO,
1.1       bertrand  516:      $                                         ONE-T*AQOAP*AAPQ1 ) )
                    517: *
                    518:                                     CALL ZROT( M, A(1,p), 1, A(1,q), 1,
1.4       bertrand  519:      $                                          CS, CONJG(OMPQ)*SN )
1.1       bertrand  520:                                     IF( RSVEC ) THEN
1.4       bertrand  521:                                         CALL ZROT( MVL, V(1,p), 1,
                    522:      $                                  V(1,q), 1, CS, CONJG(OMPQ)*SN )
1.1       bertrand  523:                                     END IF
                    524:                                  END IF
                    525:                                  D(p) = -D(q) * OMPQ
                    526: *
                    527:                               ELSE
                    528: *              .. have to use modified Gram-Schmidt like transformation
                    529:                                IF( AAPP.GT.AAQQ ) THEN
                    530:                                     CALL ZCOPY( M, A( 1, p ), 1,
                    531:      $                                          WORK, 1 )
                    532:                                     CALL ZLASCL( 'G', 0, 0, AAPP, ONE,
                    533:      $                                           M, 1, WORK,LDA,
                    534:      $                                           IERR )
                    535:                                     CALL ZLASCL( 'G', 0, 0, AAQQ, ONE,
                    536:      $                                           M, 1, A( 1, q ), LDA,
                    537:      $                                           IERR )
                    538:                                     CALL ZAXPY( M, -AAPQ, WORK,
                    539:      $                                          1, A( 1, q ), 1 )
                    540:                                     CALL ZLASCL( 'G', 0, 0, ONE, AAQQ,
                    541:      $                                           M, 1, A( 1, q ), LDA,
                    542:      $                                           IERR )
1.4       bertrand  543:                                     SVA( q ) = AAQQ*SQRT( MAX( ZERO,
1.1       bertrand  544:      $                                         ONE-AAPQ1*AAPQ1 ) )
1.4       bertrand  545:                                     MXSINJ = MAX( MXSINJ, SFMIN )
1.1       bertrand  546:                                ELSE
                    547:                                    CALL ZCOPY( M, A( 1, q ), 1,
                    548:      $                                          WORK, 1 )
                    549:                                     CALL ZLASCL( 'G', 0, 0, AAQQ, ONE,
                    550:      $                                           M, 1, WORK,LDA,
                    551:      $                                           IERR )
                    552:                                     CALL ZLASCL( 'G', 0, 0, AAPP, ONE,
                    553:      $                                           M, 1, A( 1, p ), LDA,
                    554:      $                                           IERR )
1.4       bertrand  555:                                     CALL ZAXPY( M, -CONJG(AAPQ),
1.1       bertrand  556:      $                                   WORK, 1, A( 1, p ), 1 )
                    557:                                     CALL ZLASCL( 'G', 0, 0, ONE, AAPP,
                    558:      $                                           M, 1, A( 1, p ), LDA,
                    559:      $                                           IERR )
1.4       bertrand  560:                                     SVA( p ) = AAPP*SQRT( MAX( ZERO,
1.1       bertrand  561:      $                                         ONE-AAPQ1*AAPQ1 ) )
1.4       bertrand  562:                                     MXSINJ = MAX( MXSINJ, SFMIN )
1.1       bertrand  563:                                END IF
                    564:                               END IF
                    565: *           END IF ROTOK THEN ... ELSE
                    566: *
                    567: *           In the case of cancellation in updating SVA(q), SVA(p)
                    568: *           .. recompute SVA(q), SVA(p)
                    569:                               IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
                    570:      $                            THEN
                    571:                                  IF( ( AAQQ.LT.ROOTBIG ) .AND.
                    572:      $                               ( AAQQ.GT.ROOTSFMIN ) ) THEN
                    573:                                     SVA( q ) = DZNRM2( M, A( 1, q ), 1)
                    574:                                   ELSE
                    575:                                     T = ZERO
                    576:                                     AAQQ = ONE
                    577:                                     CALL ZLASSQ( M, A( 1, q ), 1, T,
                    578:      $                                           AAQQ )
1.4       bertrand  579:                                     SVA( q ) = T*SQRT( AAQQ )
1.1       bertrand  580:                                  END IF
                    581:                               END IF
                    582:                               IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN
                    583:                                  IF( ( AAPP.LT.ROOTBIG ) .AND.
                    584:      $                               ( AAPP.GT.ROOTSFMIN ) ) THEN
                    585:                                     AAPP = DZNRM2( M, A( 1, p ), 1 )
                    586:                                  ELSE
                    587:                                     T = ZERO
                    588:                                     AAPP = ONE
                    589:                                     CALL ZLASSQ( M, A( 1, p ), 1, T,
                    590:      $                                           AAPP )
1.4       bertrand  591:                                     AAPP = T*SQRT( AAPP )
1.1       bertrand  592:                                  END IF
                    593:                                  SVA( p ) = AAPP
                    594:                               END IF
                    595: *              end of OK rotation
                    596:                            ELSE
                    597:                               NOTROT = NOTROT + 1
                    598: *[RTD]      SKIPPED  = SKIPPED  + 1
                    599:                               PSKIPPED = PSKIPPED + 1
                    600:                               IJBLSK = IJBLSK + 1
                    601:                            END IF
                    602:                         ELSE
                    603:                            NOTROT = NOTROT + 1
                    604:                            PSKIPPED = PSKIPPED + 1
                    605:                            IJBLSK = IJBLSK + 1
                    606:                         END IF
                    607: *
                    608:                         IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) )
                    609:      $                      THEN
                    610:                            SVA( p ) = AAPP
                    611:                            NOTROT = 0
                    612:                            GO TO 2011
                    613:                         END IF
                    614:                         IF( ( i.LE.SWBAND ) .AND.
                    615:      $                      ( PSKIPPED.GT.ROWSKIP ) ) THEN
                    616:                            AAPP = -AAPP
                    617:                            NOTROT = 0
                    618:                            GO TO 2203
                    619:                         END IF
                    620: *
                    621:  2200                CONTINUE
                    622: *        end of the q-loop
                    623:  2203                CONTINUE
                    624: *
                    625:                      SVA( p ) = AAPP
                    626: *
                    627:                   ELSE
                    628: *
                    629:                      IF( AAPP.EQ.ZERO )NOTROT = NOTROT +
1.4       bertrand  630:      $                   MIN( jgl+KBL-1, N ) - jgl + 1
1.1       bertrand  631:                      IF( AAPP.LT.ZERO )NOTROT = 0
                    632: *
                    633:                   END IF
                    634: *
                    635:  2100          CONTINUE
                    636: *     end of the p-loop
                    637:  2010       CONTINUE
                    638: *     end of the jbc-loop
                    639:  2011       CONTINUE
                    640: *2011 bailed out of the jbc-loop
1.4       bertrand  641:             DO 2012 p = igl, MIN( igl+KBL-1, N )
1.1       bertrand  642:                SVA( p ) = ABS( SVA( p ) )
                    643:  2012       CONTINUE
                    644: ***
                    645:  2000    CONTINUE
                    646: *2000 :: end of the ibr-loop
                    647: *
                    648: *     .. update SVA(N)
                    649:          IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) )
                    650:      $       THEN
                    651:             SVA( N ) = DZNRM2( M, A( 1, N ), 1 )
                    652:          ELSE
                    653:             T = ZERO
                    654:             AAPP = ONE
                    655:             CALL ZLASSQ( M, A( 1, N ), 1, T, AAPP )
1.4       bertrand  656:             SVA( N ) = T*SQRT( AAPP )
1.1       bertrand  657:          END IF
                    658: *
                    659: *     Additional steering devices
                    660: *
                    661:          IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.
                    662:      $       ( ISWROT.LE.N ) ) )SWBAND = i
                    663: *
1.4       bertrand  664:          IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.SQRT( DBLE( N ) )*
1.2       bertrand  665:      $       TOL ) .AND. ( DBLE( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
1.1       bertrand  666:             GO TO 1994
                    667:          END IF
                    668: *
                    669:          IF( NOTROT.GE.EMPTSW )GO TO 1994
                    670: *
                    671:  1993 CONTINUE
                    672: *     end i=1:NSWEEP loop
                    673: *
                    674: * #:( Reaching this point means that the procedure has not converged.
                    675:       INFO = NSWEEP - 1
                    676:       GO TO 1995
                    677: *
                    678:  1994 CONTINUE
                    679: * #:) Reaching this point means numerical convergence after the i-th
                    680: *     sweep.
                    681: *
                    682:       INFO = 0
                    683: * #:) INFO = 0 confirms successful iterations.
                    684:  1995 CONTINUE
                    685: *
                    686: *     Sort the vector SVA() of column norms.
                    687:       DO 5991 p = 1, N - 1
                    688:          q = IDAMAX( N-p+1, SVA( p ), 1 ) + p - 1
                    689:          IF( p.NE.q ) THEN
                    690:             TEMP1 = SVA( p )
                    691:             SVA( p ) = SVA( q )
                    692:             SVA( q ) = TEMP1
                    693:             AAPQ = D( p )
                    694:             D( p ) = D( q )
                    695:             D( q ) = AAPQ
                    696:             CALL ZSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
                    697:             IF( RSVEC )CALL ZSWAP( MVL, V( 1, p ), 1, V( 1, q ), 1 )
                    698:          END IF
                    699:  5991 CONTINUE
                    700: *
                    701: *
                    702:       RETURN
                    703: *     ..
                    704: *     .. END OF ZGSVJ1
                    705: *     ..
                    706:       END

CVSweb interface <joel.bertrand@systella.fr>