Annotation of rpl/lapack/lapack/zlals0.f, revision 1.17

1.11      bertrand    1: *> \brief \b ZLALS0 applies back multiplying factors in solving the least squares problem using divide and conquer SVD approach. Used by sgelsd.
1.8       bertrand    2: *
                      3: *  =========== DOCUMENTATION ===========
                      4: *
1.16      bertrand    5: * Online html documentation available at
                      6: *            http://www.netlib.org/lapack/explore-html/
1.8       bertrand    7: *
                      8: *> \htmlonly
1.16      bertrand    9: *> Download ZLALS0 + dependencies
                     10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlals0.f">
                     11: *> [TGZ]</a>
                     12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlals0.f">
                     13: *> [ZIP]</a>
                     14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlals0.f">
1.8       bertrand   15: *> [TXT]</a>
1.16      bertrand   16: *> \endhtmlonly
1.8       bertrand   17: *
                     18: *  Definition:
                     19: *  ===========
                     20: *
                     21: *       SUBROUTINE ZLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B, LDB, BX, LDBX,
                     22: *                          PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
                     23: *                          POLES, DIFL, DIFR, Z, K, C, S, RWORK, INFO )
1.16      bertrand   24: *
1.8       bertrand   25: *       .. Scalar Arguments ..
                     26: *       INTEGER            GIVPTR, ICOMPQ, INFO, K, LDB, LDBX, LDGCOL,
                     27: *      $                   LDGNUM, NL, NR, NRHS, SQRE
                     28: *       DOUBLE PRECISION   C, S
                     29: *       ..
                     30: *       .. Array Arguments ..
                     31: *       INTEGER            GIVCOL( LDGCOL, * ), PERM( * )
                     32: *       DOUBLE PRECISION   DIFL( * ), DIFR( LDGNUM, * ),
                     33: *      $                   GIVNUM( LDGNUM, * ), POLES( LDGNUM, * ),
                     34: *      $                   RWORK( * ), Z( * )
                     35: *       COMPLEX*16         B( LDB, * ), BX( LDBX, * )
                     36: *       ..
1.16      bertrand   37: *
1.8       bertrand   38: *
                     39: *> \par Purpose:
                     40: *  =============
                     41: *>
                     42: *> \verbatim
                     43: *>
                     44: *> ZLALS0 applies back the multiplying factors of either the left or the
                     45: *> right singular vector matrix of a diagonal matrix appended by a row
                     46: *> to the right hand side matrix B in solving the least squares problem
                     47: *> using the divide-and-conquer SVD approach.
                     48: *>
                     49: *> For the left singular vector matrix, three types of orthogonal
                     50: *> matrices are involved:
                     51: *>
                     52: *> (1L) Givens rotations: the number of such rotations is GIVPTR; the
                     53: *>      pairs of columns/rows they were applied to are stored in GIVCOL;
                     54: *>      and the C- and S-values of these rotations are stored in GIVNUM.
                     55: *>
                     56: *> (2L) Permutation. The (NL+1)-st row of B is to be moved to the first
                     57: *>      row, and for J=2:N, PERM(J)-th row of B is to be moved to the
                     58: *>      J-th row.
                     59: *>
                     60: *> (3L) The left singular vector matrix of the remaining matrix.
                     61: *>
                     62: *> For the right singular vector matrix, four types of orthogonal
                     63: *> matrices are involved:
                     64: *>
                     65: *> (1R) The right singular vector matrix of the remaining matrix.
                     66: *>
                     67: *> (2R) If SQRE = 1, one extra Givens rotation to generate the right
                     68: *>      null space.
                     69: *>
                     70: *> (3R) The inverse transformation of (2L).
                     71: *>
                     72: *> (4R) The inverse transformation of (1L).
                     73: *> \endverbatim
                     74: *
                     75: *  Arguments:
                     76: *  ==========
                     77: *
                     78: *> \param[in] ICOMPQ
                     79: *> \verbatim
                     80: *>          ICOMPQ is INTEGER
                     81: *>         Specifies whether singular vectors are to be computed in
                     82: *>         factored form:
                     83: *>         = 0: Left singular vector matrix.
                     84: *>         = 1: Right singular vector matrix.
                     85: *> \endverbatim
                     86: *>
                     87: *> \param[in] NL
                     88: *> \verbatim
                     89: *>          NL is INTEGER
                     90: *>         The row dimension of the upper block. NL >= 1.
                     91: *> \endverbatim
                     92: *>
                     93: *> \param[in] NR
                     94: *> \verbatim
                     95: *>          NR is INTEGER
                     96: *>         The row dimension of the lower block. NR >= 1.
                     97: *> \endverbatim
                     98: *>
                     99: *> \param[in] SQRE
                    100: *> \verbatim
                    101: *>          SQRE is INTEGER
                    102: *>         = 0: the lower block is an NR-by-NR square matrix.
                    103: *>         = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
                    104: *>
                    105: *>         The bidiagonal matrix has row dimension N = NL + NR + 1,
                    106: *>         and column dimension M = N + SQRE.
                    107: *> \endverbatim
                    108: *>
                    109: *> \param[in] NRHS
                    110: *> \verbatim
                    111: *>          NRHS is INTEGER
                    112: *>         The number of columns of B and BX. NRHS must be at least 1.
                    113: *> \endverbatim
                    114: *>
                    115: *> \param[in,out] B
                    116: *> \verbatim
                    117: *>          B is COMPLEX*16 array, dimension ( LDB, NRHS )
                    118: *>         On input, B contains the right hand sides of the least
                    119: *>         squares problem in rows 1 through M. On output, B contains
                    120: *>         the solution X in rows 1 through N.
                    121: *> \endverbatim
                    122: *>
                    123: *> \param[in] LDB
                    124: *> \verbatim
                    125: *>          LDB is INTEGER
                    126: *>         The leading dimension of B. LDB must be at least
                    127: *>         max(1,MAX( M, N ) ).
                    128: *> \endverbatim
                    129: *>
                    130: *> \param[out] BX
                    131: *> \verbatim
                    132: *>          BX is COMPLEX*16 array, dimension ( LDBX, NRHS )
                    133: *> \endverbatim
                    134: *>
                    135: *> \param[in] LDBX
                    136: *> \verbatim
                    137: *>          LDBX is INTEGER
                    138: *>         The leading dimension of BX.
                    139: *> \endverbatim
                    140: *>
                    141: *> \param[in] PERM
                    142: *> \verbatim
                    143: *>          PERM is INTEGER array, dimension ( N )
                    144: *>         The permutations (from deflation and sorting) applied
                    145: *>         to the two blocks.
                    146: *> \endverbatim
                    147: *>
                    148: *> \param[in] GIVPTR
                    149: *> \verbatim
                    150: *>          GIVPTR is INTEGER
                    151: *>         The number of Givens rotations which took place in this
                    152: *>         subproblem.
                    153: *> \endverbatim
                    154: *>
                    155: *> \param[in] GIVCOL
                    156: *> \verbatim
                    157: *>          GIVCOL is INTEGER array, dimension ( LDGCOL, 2 )
                    158: *>         Each pair of numbers indicates a pair of rows/columns
                    159: *>         involved in a Givens rotation.
                    160: *> \endverbatim
                    161: *>
                    162: *> \param[in] LDGCOL
                    163: *> \verbatim
                    164: *>          LDGCOL is INTEGER
                    165: *>         The leading dimension of GIVCOL, must be at least N.
                    166: *> \endverbatim
                    167: *>
                    168: *> \param[in] GIVNUM
                    169: *> \verbatim
                    170: *>          GIVNUM is DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
                    171: *>         Each number indicates the C or S value used in the
                    172: *>         corresponding Givens rotation.
                    173: *> \endverbatim
                    174: *>
                    175: *> \param[in] LDGNUM
                    176: *> \verbatim
                    177: *>          LDGNUM is INTEGER
                    178: *>         The leading dimension of arrays DIFR, POLES and
                    179: *>         GIVNUM, must be at least K.
                    180: *> \endverbatim
                    181: *>
                    182: *> \param[in] POLES
                    183: *> \verbatim
                    184: *>          POLES is DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
                    185: *>         On entry, POLES(1:K, 1) contains the new singular
                    186: *>         values obtained from solving the secular equation, and
                    187: *>         POLES(1:K, 2) is an array containing the poles in the secular
                    188: *>         equation.
                    189: *> \endverbatim
                    190: *>
                    191: *> \param[in] DIFL
                    192: *> \verbatim
                    193: *>          DIFL is DOUBLE PRECISION array, dimension ( K ).
                    194: *>         On entry, DIFL(I) is the distance between I-th updated
                    195: *>         (undeflated) singular value and the I-th (undeflated) old
                    196: *>         singular value.
                    197: *> \endverbatim
                    198: *>
                    199: *> \param[in] DIFR
                    200: *> \verbatim
                    201: *>          DIFR is DOUBLE PRECISION array, dimension ( LDGNUM, 2 ).
                    202: *>         On entry, DIFR(I, 1) contains the distances between I-th
                    203: *>         updated (undeflated) singular value and the I+1-th
                    204: *>         (undeflated) old singular value. And DIFR(I, 2) is the
                    205: *>         normalizing factor for the I-th right singular vector.
                    206: *> \endverbatim
                    207: *>
                    208: *> \param[in] Z
                    209: *> \verbatim
                    210: *>          Z is DOUBLE PRECISION array, dimension ( K )
                    211: *>         Contain the components of the deflation-adjusted updating row
                    212: *>         vector.
                    213: *> \endverbatim
                    214: *>
                    215: *> \param[in] K
                    216: *> \verbatim
                    217: *>          K is INTEGER
                    218: *>         Contains the dimension of the non-deflated matrix,
                    219: *>         This is the order of the related secular equation. 1 <= K <=N.
                    220: *> \endverbatim
                    221: *>
                    222: *> \param[in] C
                    223: *> \verbatim
                    224: *>          C is DOUBLE PRECISION
                    225: *>         C contains garbage if SQRE =0 and the C-value of a Givens
                    226: *>         rotation related to the right null space if SQRE = 1.
                    227: *> \endverbatim
                    228: *>
                    229: *> \param[in] S
                    230: *> \verbatim
                    231: *>          S is DOUBLE PRECISION
                    232: *>         S contains garbage if SQRE =0 and the S-value of a Givens
                    233: *>         rotation related to the right null space if SQRE = 1.
                    234: *> \endverbatim
                    235: *>
                    236: *> \param[out] RWORK
                    237: *> \verbatim
                    238: *>          RWORK is DOUBLE PRECISION array, dimension
                    239: *>         ( K*(1+NRHS) + 2*NRHS )
                    240: *> \endverbatim
                    241: *>
                    242: *> \param[out] INFO
                    243: *> \verbatim
                    244: *>          INFO is INTEGER
                    245: *>          = 0:  successful exit.
                    246: *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
                    247: *> \endverbatim
                    248: *
                    249: *  Authors:
                    250: *  ========
                    251: *
1.16      bertrand  252: *> \author Univ. of Tennessee
                    253: *> \author Univ. of California Berkeley
                    254: *> \author Univ. of Colorado Denver
                    255: *> \author NAG Ltd.
1.8       bertrand  256: *
1.16      bertrand  257: *> \date December 2016
1.8       bertrand  258: *
                    259: *> \ingroup complex16OTHERcomputational
                    260: *
                    261: *> \par Contributors:
                    262: *  ==================
                    263: *>
                    264: *>     Ming Gu and Ren-Cang Li, Computer Science Division, University of
                    265: *>       California at Berkeley, USA \n
                    266: *>     Osni Marques, LBNL/NERSC, USA \n
                    267: *
                    268: *  =====================================================================
1.1       bertrand  269:       SUBROUTINE ZLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B, LDB, BX, LDBX,
                    270:      $                   PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
                    271:      $                   POLES, DIFL, DIFR, Z, K, C, S, RWORK, INFO )
                    272: *
1.16      bertrand  273: *  -- LAPACK computational routine (version 3.7.0) --
1.1       bertrand  274: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    275: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.16      bertrand  276: *     December 2016
1.1       bertrand  277: *
                    278: *     .. Scalar Arguments ..
                    279:       INTEGER            GIVPTR, ICOMPQ, INFO, K, LDB, LDBX, LDGCOL,
                    280:      $                   LDGNUM, NL, NR, NRHS, SQRE
                    281:       DOUBLE PRECISION   C, S
                    282: *     ..
                    283: *     .. Array Arguments ..
                    284:       INTEGER            GIVCOL( LDGCOL, * ), PERM( * )
                    285:       DOUBLE PRECISION   DIFL( * ), DIFR( LDGNUM, * ),
                    286:      $                   GIVNUM( LDGNUM, * ), POLES( LDGNUM, * ),
                    287:      $                   RWORK( * ), Z( * )
                    288:       COMPLEX*16         B( LDB, * ), BX( LDBX, * )
                    289: *     ..
                    290: *
                    291: *  =====================================================================
                    292: *
                    293: *     .. Parameters ..
                    294:       DOUBLE PRECISION   ONE, ZERO, NEGONE
                    295:       PARAMETER          ( ONE = 1.0D0, ZERO = 0.0D0, NEGONE = -1.0D0 )
                    296: *     ..
                    297: *     .. Local Scalars ..
                    298:       INTEGER            I, J, JCOL, JROW, M, N, NLP1
                    299:       DOUBLE PRECISION   DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, TEMP
                    300: *     ..
                    301: *     .. External Subroutines ..
                    302:       EXTERNAL           DGEMV, XERBLA, ZCOPY, ZDROT, ZDSCAL, ZLACPY,
                    303:      $                   ZLASCL
                    304: *     ..
                    305: *     .. External Functions ..
                    306:       DOUBLE PRECISION   DLAMC3, DNRM2
                    307:       EXTERNAL           DLAMC3, DNRM2
                    308: *     ..
                    309: *     .. Intrinsic Functions ..
                    310:       INTRINSIC          DBLE, DCMPLX, DIMAG, MAX
                    311: *     ..
                    312: *     .. Executable Statements ..
                    313: *
                    314: *     Test the input parameters.
                    315: *
                    316:       INFO = 0
1.14      bertrand  317:       N = NL + NR + 1
1.1       bertrand  318: *
                    319:       IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
                    320:          INFO = -1
                    321:       ELSE IF( NL.LT.1 ) THEN
                    322:          INFO = -2
                    323:       ELSE IF( NR.LT.1 ) THEN
                    324:          INFO = -3
                    325:       ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
                    326:          INFO = -4
1.14      bertrand  327:       ELSE IF( NRHS.LT.1 ) THEN
1.1       bertrand  328:          INFO = -5
                    329:       ELSE IF( LDB.LT.N ) THEN
                    330:          INFO = -7
                    331:       ELSE IF( LDBX.LT.N ) THEN
                    332:          INFO = -9
                    333:       ELSE IF( GIVPTR.LT.0 ) THEN
                    334:          INFO = -11
                    335:       ELSE IF( LDGCOL.LT.N ) THEN
                    336:          INFO = -13
                    337:       ELSE IF( LDGNUM.LT.N ) THEN
                    338:          INFO = -15
                    339:       ELSE IF( K.LT.1 ) THEN
                    340:          INFO = -20
                    341:       END IF
                    342:       IF( INFO.NE.0 ) THEN
                    343:          CALL XERBLA( 'ZLALS0', -INFO )
                    344:          RETURN
                    345:       END IF
                    346: *
                    347:       M = N + SQRE
                    348:       NLP1 = NL + 1
                    349: *
                    350:       IF( ICOMPQ.EQ.0 ) THEN
                    351: *
                    352: *        Apply back orthogonal transformations from the left.
                    353: *
                    354: *        Step (1L): apply back the Givens rotations performed.
                    355: *
                    356:          DO 10 I = 1, GIVPTR
                    357:             CALL ZDROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB,
                    358:      $                  B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ),
                    359:      $                  GIVNUM( I, 1 ) )
                    360:    10    CONTINUE
                    361: *
                    362: *        Step (2L): permute rows of B.
                    363: *
                    364:          CALL ZCOPY( NRHS, B( NLP1, 1 ), LDB, BX( 1, 1 ), LDBX )
                    365:          DO 20 I = 2, N
                    366:             CALL ZCOPY( NRHS, B( PERM( I ), 1 ), LDB, BX( I, 1 ), LDBX )
                    367:    20    CONTINUE
                    368: *
                    369: *        Step (3L): apply the inverse of the left singular vector
                    370: *        matrix to BX.
                    371: *
                    372:          IF( K.EQ.1 ) THEN
                    373:             CALL ZCOPY( NRHS, BX, LDBX, B, LDB )
                    374:             IF( Z( 1 ).LT.ZERO ) THEN
                    375:                CALL ZDSCAL( NRHS, NEGONE, B, LDB )
                    376:             END IF
                    377:          ELSE
                    378:             DO 100 J = 1, K
                    379:                DIFLJ = DIFL( J )
                    380:                DJ = POLES( J, 1 )
                    381:                DSIGJ = -POLES( J, 2 )
                    382:                IF( J.LT.K ) THEN
                    383:                   DIFRJ = -DIFR( J, 1 )
                    384:                   DSIGJP = -POLES( J+1, 2 )
                    385:                END IF
                    386:                IF( ( Z( J ).EQ.ZERO ) .OR. ( POLES( J, 2 ).EQ.ZERO ) )
                    387:      $              THEN
                    388:                   RWORK( J ) = ZERO
                    389:                ELSE
                    390:                   RWORK( J ) = -POLES( J, 2 )*Z( J ) / DIFLJ /
                    391:      $                         ( POLES( J, 2 )+DJ )
                    392:                END IF
                    393:                DO 30 I = 1, J - 1
                    394:                   IF( ( Z( I ).EQ.ZERO ) .OR.
                    395:      $                ( POLES( I, 2 ).EQ.ZERO ) ) THEN
                    396:                      RWORK( I ) = ZERO
                    397:                   ELSE
                    398:                      RWORK( I ) = POLES( I, 2 )*Z( I ) /
                    399:      $                            ( DLAMC3( POLES( I, 2 ), DSIGJ )-
                    400:      $                            DIFLJ ) / ( POLES( I, 2 )+DJ )
                    401:                   END IF
                    402:    30          CONTINUE
                    403:                DO 40 I = J + 1, K
                    404:                   IF( ( Z( I ).EQ.ZERO ) .OR.
                    405:      $                ( POLES( I, 2 ).EQ.ZERO ) ) THEN
                    406:                      RWORK( I ) = ZERO
                    407:                   ELSE
                    408:                      RWORK( I ) = POLES( I, 2 )*Z( I ) /
                    409:      $                            ( DLAMC3( POLES( I, 2 ), DSIGJP )+
                    410:      $                            DIFRJ ) / ( POLES( I, 2 )+DJ )
                    411:                   END IF
                    412:    40          CONTINUE
                    413:                RWORK( 1 ) = NEGONE
                    414:                TEMP = DNRM2( K, RWORK, 1 )
                    415: *
                    416: *              Since B and BX are complex, the following call to DGEMV
                    417: *              is performed in two steps (real and imaginary parts).
                    418: *
                    419: *              CALL DGEMV( 'T', K, NRHS, ONE, BX, LDBX, WORK, 1, ZERO,
                    420: *    $                     B( J, 1 ), LDB )
                    421: *
                    422:                I = K + NRHS*2
                    423:                DO 60 JCOL = 1, NRHS
                    424:                   DO 50 JROW = 1, K
                    425:                      I = I + 1
                    426:                      RWORK( I ) = DBLE( BX( JROW, JCOL ) )
                    427:    50             CONTINUE
                    428:    60          CONTINUE
                    429:                CALL DGEMV( 'T', K, NRHS, ONE, RWORK( 1+K+NRHS*2 ), K,
                    430:      $                     RWORK( 1 ), 1, ZERO, RWORK( 1+K ), 1 )
                    431:                I = K + NRHS*2
                    432:                DO 80 JCOL = 1, NRHS
                    433:                   DO 70 JROW = 1, K
                    434:                      I = I + 1
                    435:                      RWORK( I ) = DIMAG( BX( JROW, JCOL ) )
                    436:    70             CONTINUE
                    437:    80          CONTINUE
                    438:                CALL DGEMV( 'T', K, NRHS, ONE, RWORK( 1+K+NRHS*2 ), K,
                    439:      $                     RWORK( 1 ), 1, ZERO, RWORK( 1+K+NRHS ), 1 )
                    440:                DO 90 JCOL = 1, NRHS
                    441:                   B( J, JCOL ) = DCMPLX( RWORK( JCOL+K ),
                    442:      $                           RWORK( JCOL+K+NRHS ) )
                    443:    90          CONTINUE
                    444:                CALL ZLASCL( 'G', 0, 0, TEMP, ONE, 1, NRHS, B( J, 1 ),
                    445:      $                      LDB, INFO )
                    446:   100       CONTINUE
                    447:          END IF
                    448: *
                    449: *        Move the deflated rows of BX to B also.
                    450: *
                    451:          IF( K.LT.MAX( M, N ) )
                    452:      $      CALL ZLACPY( 'A', N-K, NRHS, BX( K+1, 1 ), LDBX,
                    453:      $                   B( K+1, 1 ), LDB )
                    454:       ELSE
                    455: *
                    456: *        Apply back the right orthogonal transformations.
                    457: *
                    458: *        Step (1R): apply back the new right singular vector matrix
                    459: *        to B.
                    460: *
                    461:          IF( K.EQ.1 ) THEN
                    462:             CALL ZCOPY( NRHS, B, LDB, BX, LDBX )
                    463:          ELSE
                    464:             DO 180 J = 1, K
                    465:                DSIGJ = POLES( J, 2 )
                    466:                IF( Z( J ).EQ.ZERO ) THEN
                    467:                   RWORK( J ) = ZERO
                    468:                ELSE
                    469:                   RWORK( J ) = -Z( J ) / DIFL( J ) /
                    470:      $                         ( DSIGJ+POLES( J, 1 ) ) / DIFR( J, 2 )
                    471:                END IF
                    472:                DO 110 I = 1, J - 1
                    473:                   IF( Z( J ).EQ.ZERO ) THEN
                    474:                      RWORK( I ) = ZERO
                    475:                   ELSE
                    476:                      RWORK( I ) = Z( J ) / ( DLAMC3( DSIGJ, -POLES( I+1,
                    477:      $                            2 ) )-DIFR( I, 1 ) ) /
                    478:      $                            ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 )
                    479:                   END IF
                    480:   110          CONTINUE
                    481:                DO 120 I = J + 1, K
                    482:                   IF( Z( J ).EQ.ZERO ) THEN
                    483:                      RWORK( I ) = ZERO
                    484:                   ELSE
                    485:                      RWORK( I ) = Z( J ) / ( DLAMC3( DSIGJ, -POLES( I,
                    486:      $                            2 ) )-DIFL( I ) ) /
                    487:      $                            ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 )
                    488:                   END IF
                    489:   120          CONTINUE
                    490: *
                    491: *              Since B and BX are complex, the following call to DGEMV
                    492: *              is performed in two steps (real and imaginary parts).
                    493: *
                    494: *              CALL DGEMV( 'T', K, NRHS, ONE, B, LDB, WORK, 1, ZERO,
                    495: *    $                     BX( J, 1 ), LDBX )
                    496: *
                    497:                I = K + NRHS*2
                    498:                DO 140 JCOL = 1, NRHS
                    499:                   DO 130 JROW = 1, K
                    500:                      I = I + 1
                    501:                      RWORK( I ) = DBLE( B( JROW, JCOL ) )
                    502:   130             CONTINUE
                    503:   140          CONTINUE
                    504:                CALL DGEMV( 'T', K, NRHS, ONE, RWORK( 1+K+NRHS*2 ), K,
                    505:      $                     RWORK( 1 ), 1, ZERO, RWORK( 1+K ), 1 )
                    506:                I = K + NRHS*2
                    507:                DO 160 JCOL = 1, NRHS
                    508:                   DO 150 JROW = 1, K
                    509:                      I = I + 1
                    510:                      RWORK( I ) = DIMAG( B( JROW, JCOL ) )
                    511:   150             CONTINUE
                    512:   160          CONTINUE
                    513:                CALL DGEMV( 'T', K, NRHS, ONE, RWORK( 1+K+NRHS*2 ), K,
                    514:      $                     RWORK( 1 ), 1, ZERO, RWORK( 1+K+NRHS ), 1 )
                    515:                DO 170 JCOL = 1, NRHS
                    516:                   BX( J, JCOL ) = DCMPLX( RWORK( JCOL+K ),
                    517:      $                            RWORK( JCOL+K+NRHS ) )
                    518:   170          CONTINUE
                    519:   180       CONTINUE
                    520:          END IF
                    521: *
                    522: *        Step (2R): if SQRE = 1, apply back the rotation that is
                    523: *        related to the right null space of the subproblem.
                    524: *
                    525:          IF( SQRE.EQ.1 ) THEN
                    526:             CALL ZCOPY( NRHS, B( M, 1 ), LDB, BX( M, 1 ), LDBX )
                    527:             CALL ZDROT( NRHS, BX( 1, 1 ), LDBX, BX( M, 1 ), LDBX, C, S )
                    528:          END IF
                    529:          IF( K.LT.MAX( M, N ) )
                    530:      $      CALL ZLACPY( 'A', N-K, NRHS, B( K+1, 1 ), LDB, BX( K+1, 1 ),
                    531:      $                   LDBX )
                    532: *
                    533: *        Step (3R): permute rows of B.
                    534: *
                    535:          CALL ZCOPY( NRHS, BX( 1, 1 ), LDBX, B( NLP1, 1 ), LDB )
                    536:          IF( SQRE.EQ.1 ) THEN
                    537:             CALL ZCOPY( NRHS, BX( M, 1 ), LDBX, B( M, 1 ), LDB )
                    538:          END IF
                    539:          DO 190 I = 2, N
                    540:             CALL ZCOPY( NRHS, BX( I, 1 ), LDBX, B( PERM( I ), 1 ), LDB )
                    541:   190    CONTINUE
                    542: *
                    543: *        Step (4R): apply back the Givens rotations performed.
                    544: *
                    545:          DO 200 I = GIVPTR, 1, -1
                    546:             CALL ZDROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB,
                    547:      $                  B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ),
                    548:      $                  -GIVNUM( I, 1 ) )
                    549:   200    CONTINUE
                    550:       END IF
                    551: *
                    552:       RETURN
                    553: *
                    554: *     End of ZLALS0
                    555: *
                    556:       END

CVSweb interface <joel.bertrand@systella.fr>