LAPACK: SRC/dlals0.f Source File

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