LAPACK: SRC/zlals0.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B, LDB, BX, LDBX,
00002      $                   PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
00003      $                   POLES, DIFL, DIFR, Z, K, C, S, RWORK, INFO )
00004 *
00005 *  -- LAPACK routine (version 3.2) --
00006 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00007 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00008 *     November 2006
00009 *
00010 *     .. Scalar Arguments ..
00011       INTEGER            GIVPTR, ICOMPQ, INFO, K, LDB, LDBX, LDGCOL,
00012      $                   LDGNUM, NL, NR, NRHS, SQRE
00013       DOUBLE PRECISION   C, S
00014 *     ..
00015 *     .. Array Arguments ..
00016       INTEGER            GIVCOL( LDGCOL, * ), PERM( * )
00017       DOUBLE PRECISION   DIFL( * ), DIFR( LDGNUM, * ),
00018      $                   GIVNUM( LDGNUM, * ), POLES( LDGNUM, * ),
00019      $                   RWORK( * ), Z( * )
00020       COMPLEX*16         B( LDB, * ), BX( LDBX, * )
00021 *     ..
00022 *
00023 *  Purpose
00024 *  =======
00025 *
00026 *  ZLALS0 applies back the multiplying factors of either the left or the
00027 *  right singular vector matrix of a diagonal matrix appended by a row
00028 *  to the right hand side matrix B in solving the least squares problem
00029 *  using the divide-and-conquer SVD approach.
00030 *
00031 *  For the left singular vector matrix, three types of orthogonal
00032 *  matrices are involved:
00033 *
00034 *  (1L) Givens rotations: the number of such rotations is GIVPTR; the
00035 *       pairs of columns/rows they were applied to are stored in GIVCOL;
00036 *       and the C- and S-values of these rotations are stored in GIVNUM.
00037 *
00038 *  (2L) Permutation. The (NL+1)-st row of B is to be moved to the first
00039 *       row, and for J=2:N, PERM(J)-th row of B is to be moved to the
00040 *       J-th row.
00041 *
00042 *  (3L) The left singular vector matrix of the remaining matrix.
00043 *
00044 *  For the right singular vector matrix, four types of orthogonal
00045 *  matrices are involved:
00046 *
00047 *  (1R) The right singular vector matrix of the remaining matrix.
00048 *
00049 *  (2R) If SQRE = 1, one extra Givens rotation to generate the right
00050 *       null space.
00051 *
00052 *  (3R) The inverse transformation of (2L).
00053 *
00054 *  (4R) The inverse transformation of (1L).
00055 *
00056 *  Arguments
00057 *  =========
00058 *
00059 *  ICOMPQ (input) INTEGER
00060 *         Specifies whether singular vectors are to be computed in
00061 *         factored form:
00062 *         = 0: Left singular vector matrix.
00063 *         = 1: Right singular vector matrix.
00064 *
00065 *  NL     (input) INTEGER
00066 *         The row dimension of the upper block. NL >= 1.
00067 *
00068 *  NR     (input) INTEGER
00069 *         The row dimension of the lower block. NR >= 1.
00070 *
00071 *  SQRE   (input) INTEGER
00072 *         = 0: the lower block is an NR-by-NR square matrix.
00073 *         = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
00074 *
00075 *         The bidiagonal matrix has row dimension N = NL + NR + 1,
00076 *         and column dimension M = N + SQRE.
00077 *
00078 *  NRHS   (input) INTEGER
00079 *         The number of columns of B and BX. NRHS must be at least 1.
00080 *
00081 *  B      (input/output) COMPLEX*16 array, dimension ( LDB, NRHS )
00082 *         On input, B contains the right hand sides of the least
00083 *         squares problem in rows 1 through M. On output, B contains
00084 *         the solution X in rows 1 through N.
00085 *
00086 *  LDB    (input) INTEGER
00087 *         The leading dimension of B. LDB must be at least
00088 *         max(1,MAX( M, N ) ).
00089 *
00090 *  BX     (workspace) COMPLEX*16 array, dimension ( LDBX, NRHS )
00091 *
00092 *  LDBX   (input) INTEGER
00093 *         The leading dimension of BX.
00094 *
00095 *  PERM   (input) INTEGER array, dimension ( N )
00096 *         The permutations (from deflation and sorting) applied
00097 *         to the two blocks.
00098 *
00099 *  GIVPTR (input) INTEGER
00100 *         The number of Givens rotations which took place in this
00101 *         subproblem.
00102 *
00103 *  GIVCOL (input) INTEGER array, dimension ( LDGCOL, 2 )
00104 *         Each pair of numbers indicates a pair of rows/columns
00105 *         involved in a Givens rotation.
00106 *
00107 *  LDGCOL (input) INTEGER
00108 *         The leading dimension of GIVCOL, must be at least N.
00109 *
00110 *  GIVNUM (input) DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
00111 *         Each number indicates the C or S value used in the
00112 *         corresponding Givens rotation.
00113 *
00114 *  LDGNUM (input) INTEGER
00115 *         The leading dimension of arrays DIFR, POLES and
00116 *         GIVNUM, must be at least K.
00117 *
00118 *  POLES  (input) DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
00119 *         On entry, POLES(1:K, 1) contains the new singular
00120 *         values obtained from solving the secular equation, and
00121 *         POLES(1:K, 2) is an array containing the poles in the secular
00122 *         equation.
00123 *
00124 *  DIFL   (input) DOUBLE PRECISION array, dimension ( K ).
00125 *         On entry, DIFL(I) is the distance between I-th updated
00126 *         (undeflated) singular value and the I-th (undeflated) old
00127 *         singular value.
00128 *
00129 *  DIFR   (input) DOUBLE PRECISION array, dimension ( LDGNUM, 2 ).
00130 *         On entry, DIFR(I, 1) contains the distances between I-th
00131 *         updated (undeflated) singular value and the I+1-th
00132 *         (undeflated) old singular value. And DIFR(I, 2) is the
00133 *         normalizing factor for the I-th right singular vector.
00134 *
00135 *  Z      (input) DOUBLE PRECISION array, dimension ( K )
00136 *         Contain the components of the deflation-adjusted updating row
00137 *         vector.
00138 *
00139 *  K      (input) INTEGER
00140 *         Contains the dimension of the non-deflated matrix,
00141 *         This is the order of the related secular equation. 1 <= K <=N.
00142 *
00143 *  C      (input) DOUBLE PRECISION
00144 *         C contains garbage if SQRE =0 and the C-value of a Givens
00145 *         rotation related to the right null space if SQRE = 1.
00146 *
00147 *  S      (input) DOUBLE PRECISION
00148 *         S contains garbage if SQRE =0 and the S-value of a Givens
00149 *         rotation related to the right null space if SQRE = 1.
00150 *
00151 *  RWORK  (workspace) DOUBLE PRECISION array, dimension
00152 *         ( K*(1+NRHS) + 2*NRHS )
00153 *
00154 *  INFO   (output) INTEGER
00155 *          = 0:  successful exit.
00156 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00157 *
00158 *  Further Details
00159 *  ===============
00160 *
00161 *  Based on contributions by
00162 *     Ming Gu and Ren-Cang Li, Computer Science Division, University of
00163 *       California at Berkeley, USA
00164 *     Osni Marques, LBNL/NERSC, USA
00165 *
00166 *  =====================================================================
00167 *
00168 *     .. Parameters ..
00169       DOUBLE PRECISION   ONE, ZERO, NEGONE
00170       PARAMETER          ( ONE = 1.0D0, ZERO = 0.0D0, NEGONE = -1.0D0 )
00171 *     ..
00172 *     .. Local Scalars ..
00173       INTEGER            I, J, JCOL, JROW, M, N, NLP1
00174       DOUBLE PRECISION   DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, TEMP
00175 *     ..
00176 *     .. External Subroutines ..
00177       EXTERNAL           DGEMV, XERBLA, ZCOPY, ZDROT, ZDSCAL, ZLACPY,
00178      $                   ZLASCL
00179 *     ..
00180 *     .. External Functions ..
00181       DOUBLE PRECISION   DLAMC3, DNRM2
00182       EXTERNAL           DLAMC3, DNRM2
00183 *     ..
00184 *     .. Intrinsic Functions ..
00185       INTRINSIC          DBLE, DCMPLX, DIMAG, MAX
00186 *     ..
00187 *     .. Executable Statements ..
00188 *
00189 *     Test the input parameters.
00190 *
00191       INFO = 0
00192 *
00193       IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
00194          INFO = -1
00195       ELSE IF( NL.LT.1 ) THEN
00196          INFO = -2
00197       ELSE IF( NR.LT.1 ) THEN
00198          INFO = -3
00199       ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
00200          INFO = -4
00201       END IF
00202 *
00203       N = NL + NR + 1
00204 *
00205       IF( NRHS.LT.1 ) THEN
00206          INFO = -5
00207       ELSE IF( LDB.LT.N ) THEN
00208          INFO = -7
00209       ELSE IF( LDBX.LT.N ) THEN
00210          INFO = -9
00211       ELSE IF( GIVPTR.LT.0 ) THEN
00212          INFO = -11
00213       ELSE IF( LDGCOL.LT.N ) THEN
00214          INFO = -13
00215       ELSE IF( LDGNUM.LT.N ) THEN
00216          INFO = -15
00217       ELSE IF( K.LT.1 ) THEN
00218          INFO = -20
00219       END IF
00220       IF( INFO.NE.0 ) THEN
00221          CALL XERBLA( 'ZLALS0', -INFO )
00222          RETURN
00223       END IF
00224 *
00225       M = N + SQRE
00226       NLP1 = NL + 1
00227 *
00228       IF( ICOMPQ.EQ.0 ) THEN
00229 *
00230 *        Apply back orthogonal transformations from the left.
00231 *
00232 *        Step (1L): apply back the Givens rotations performed.
00233 *
00234          DO 10 I = 1, GIVPTR
00235             CALL ZDROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB,
00236      $                  B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ),
00237      $                  GIVNUM( I, 1 ) )
00238    10    CONTINUE
00239 *
00240 *        Step (2L): permute rows of B.
00241 *
00242          CALL ZCOPY( NRHS, B( NLP1, 1 ), LDB, BX( 1, 1 ), LDBX )
00243          DO 20 I = 2, N
00244             CALL ZCOPY( NRHS, B( PERM( I ), 1 ), LDB, BX( I, 1 ), LDBX )
00245    20    CONTINUE
00246 *
00247 *        Step (3L): apply the inverse of the left singular vector
00248 *        matrix to BX.
00249 *
00250          IF( K.EQ.1 ) THEN
00251             CALL ZCOPY( NRHS, BX, LDBX, B, LDB )
00252             IF( Z( 1 ).LT.ZERO ) THEN
00253                CALL ZDSCAL( NRHS, NEGONE, B, LDB )
00254             END IF
00255          ELSE
00256             DO 100 J = 1, K
00257                DIFLJ = DIFL( J )
00258                DJ = POLES( J, 1 )
00259                DSIGJ = -POLES( J, 2 )
00260                IF( J.LT.K ) THEN
00261                   DIFRJ = -DIFR( J, 1 )
00262                   DSIGJP = -POLES( J+1, 2 )
00263                END IF
00264                IF( ( Z( J ).EQ.ZERO ) .OR. ( POLES( J, 2 ).EQ.ZERO ) )
00265      $              THEN
00266                   RWORK( J ) = ZERO
00267                ELSE
00268                   RWORK( J ) = -POLES( J, 2 )*Z( J ) / DIFLJ /
00269      $                         ( POLES( J, 2 )+DJ )
00270                END IF
00271                DO 30 I = 1, J - 1
00272                   IF( ( Z( I ).EQ.ZERO ) .OR.
00273      $                ( POLES( I, 2 ).EQ.ZERO ) ) THEN
00274                      RWORK( I ) = ZERO
00275                   ELSE
00276                      RWORK( I ) = POLES( I, 2 )*Z( I ) /
00277      $                            ( DLAMC3( POLES( I, 2 ), DSIGJ )-
00278      $                            DIFLJ ) / ( POLES( I, 2 )+DJ )
00279                   END IF
00280    30          CONTINUE
00281                DO 40 I = J + 1, K
00282                   IF( ( Z( I ).EQ.ZERO ) .OR.
00283      $                ( POLES( I, 2 ).EQ.ZERO ) ) THEN
00284                      RWORK( I ) = ZERO
00285                   ELSE
00286                      RWORK( I ) = POLES( I, 2 )*Z( I ) /
00287      $                            ( DLAMC3( POLES( I, 2 ), DSIGJP )+
00288      $                            DIFRJ ) / ( POLES( I, 2 )+DJ )
00289                   END IF
00290    40          CONTINUE
00291                RWORK( 1 ) = NEGONE
00292                TEMP = DNRM2( K, RWORK, 1 )
00293 *
00294 *              Since B and BX are complex, the following call to DGEMV
00295 *              is performed in two steps (real and imaginary parts).
00296 *
00297 *              CALL DGEMV( 'T', K, NRHS, ONE, BX, LDBX, WORK, 1, ZERO,
00298 *    $                     B( J, 1 ), LDB )
00299 *
00300                I = K + NRHS*2
00301                DO 60 JCOL = 1, NRHS
00302                   DO 50 JROW = 1, K
00303                      I = I + 1
00304                      RWORK( I ) = DBLE( BX( JROW, JCOL ) )
00305    50             CONTINUE
00306    60          CONTINUE
00307                CALL DGEMV( 'T', K, NRHS, ONE, RWORK( 1+K+NRHS*2 ), K,
00308      $                     RWORK( 1 ), 1, ZERO, RWORK( 1+K ), 1 )
00309                I = K + NRHS*2
00310                DO 80 JCOL = 1, NRHS
00311                   DO 70 JROW = 1, K
00312                      I = I + 1
00313                      RWORK( I ) = DIMAG( BX( JROW, JCOL ) )
00314    70             CONTINUE
00315    80          CONTINUE
00316                CALL DGEMV( 'T', K, NRHS, ONE, RWORK( 1+K+NRHS*2 ), K,
00317      $                     RWORK( 1 ), 1, ZERO, RWORK( 1+K+NRHS ), 1 )
00318                DO 90 JCOL = 1, NRHS
00319                   B( J, JCOL ) = DCMPLX( RWORK( JCOL+K ),
00320      $                           RWORK( JCOL+K+NRHS ) )
00321    90          CONTINUE
00322                CALL ZLASCL( 'G', 0, 0, TEMP, ONE, 1, NRHS, B( J, 1 ),
00323      $                      LDB, INFO )
00324   100       CONTINUE
00325          END IF
00326 *
00327 *        Move the deflated rows of BX to B also.
00328 *
00329          IF( K.LT.MAX( M, N ) )
00330      $      CALL ZLACPY( 'A', N-K, NRHS, BX( K+1, 1 ), LDBX,
00331      $                   B( K+1, 1 ), LDB )
00332       ELSE
00333 *
00334 *        Apply back the right orthogonal transformations.
00335 *
00336 *        Step (1R): apply back the new right singular vector matrix
00337 *        to B.
00338 *
00339          IF( K.EQ.1 ) THEN
00340             CALL ZCOPY( NRHS, B, LDB, BX, LDBX )
00341          ELSE
00342             DO 180 J = 1, K
00343                DSIGJ = POLES( J, 2 )
00344                IF( Z( J ).EQ.ZERO ) THEN
00345                   RWORK( J ) = ZERO
00346                ELSE
00347                   RWORK( J ) = -Z( J ) / DIFL( J ) /
00348      $                         ( DSIGJ+POLES( J, 1 ) ) / DIFR( J, 2 )
00349                END IF
00350                DO 110 I = 1, J - 1
00351                   IF( Z( J ).EQ.ZERO ) THEN
00352                      RWORK( I ) = ZERO
00353                   ELSE
00354                      RWORK( I ) = Z( J ) / ( DLAMC3( DSIGJ, -POLES( I+1,
00355      $                            2 ) )-DIFR( I, 1 ) ) /
00356      $                            ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 )
00357                   END IF
00358   110          CONTINUE
00359                DO 120 I = J + 1, K
00360                   IF( Z( J ).EQ.ZERO ) THEN
00361                      RWORK( I ) = ZERO
00362                   ELSE
00363                      RWORK( I ) = Z( J ) / ( DLAMC3( DSIGJ, -POLES( I,
00364      $                            2 ) )-DIFL( I ) ) /
00365      $                            ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 )
00366                   END IF
00367   120          CONTINUE
00368 *
00369 *              Since B and BX are complex, the following call to DGEMV
00370 *              is performed in two steps (real and imaginary parts).
00371 *
00372 *              CALL DGEMV( 'T', K, NRHS, ONE, B, LDB, WORK, 1, ZERO,
00373 *    $                     BX( J, 1 ), LDBX )
00374 *
00375                I = K + NRHS*2
00376                DO 140 JCOL = 1, NRHS
00377                   DO 130 JROW = 1, K
00378                      I = I + 1
00379                      RWORK( I ) = DBLE( B( JROW, JCOL ) )
00380   130             CONTINUE
00381   140          CONTINUE
00382                CALL DGEMV( 'T', K, NRHS, ONE, RWORK( 1+K+NRHS*2 ), K,
00383      $                     RWORK( 1 ), 1, ZERO, RWORK( 1+K ), 1 )
00384                I = K + NRHS*2
00385                DO 160 JCOL = 1, NRHS
00386                   DO 150 JROW = 1, K
00387                      I = I + 1
00388                      RWORK( I ) = DIMAG( B( JROW, JCOL ) )
00389   150             CONTINUE
00390   160          CONTINUE
00391                CALL DGEMV( 'T', K, NRHS, ONE, RWORK( 1+K+NRHS*2 ), K,
00392      $                     RWORK( 1 ), 1, ZERO, RWORK( 1+K+NRHS ), 1 )
00393                DO 170 JCOL = 1, NRHS
00394                   BX( J, JCOL ) = DCMPLX( RWORK( JCOL+K ),
00395      $                            RWORK( JCOL+K+NRHS ) )
00396   170          CONTINUE
00397   180       CONTINUE
00398          END IF
00399 *
00400 *        Step (2R): if SQRE = 1, apply back the rotation that is
00401 *        related to the right null space of the subproblem.
00402 *
00403          IF( SQRE.EQ.1 ) THEN
00404             CALL ZCOPY( NRHS, B( M, 1 ), LDB, BX( M, 1 ), LDBX )
00405             CALL ZDROT( NRHS, BX( 1, 1 ), LDBX, BX( M, 1 ), LDBX, C, S )
00406          END IF
00407          IF( K.LT.MAX( M, N ) )
00408      $      CALL ZLACPY( 'A', N-K, NRHS, B( K+1, 1 ), LDB, BX( K+1, 1 ),
00409      $                   LDBX )
00410 *
00411 *        Step (3R): permute rows of B.
00412 *
00413          CALL ZCOPY( NRHS, BX( 1, 1 ), LDBX, B( NLP1, 1 ), LDB )
00414          IF( SQRE.EQ.1 ) THEN
00415             CALL ZCOPY( NRHS, BX( M, 1 ), LDBX, B( M, 1 ), LDB )
00416          END IF
00417          DO 190 I = 2, N
00418             CALL ZCOPY( NRHS, BX( I, 1 ), LDBX, B( PERM( I ), 1 ), LDB )
00419   190    CONTINUE
00420 *
00421 *        Step (4R): apply back the Givens rotations performed.
00422 *
00423          DO 200 I = GIVPTR, 1, -1
00424             CALL ZDROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB,
00425      $                  B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ),
00426      $                  -GIVNUM( I, 1 ) )
00427   200    CONTINUE
00428       END IF
00429 *
00430       RETURN
00431 *
00432 *     End of ZLALS0
00433 *
00434       END