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