LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
slalsa.f
Go to the documentation of this file.
00001 *> \brief \b SLALSA
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download SLALSA + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slalsa.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slalsa.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slalsa.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE SLALSA( ICOMPQ, SMLSIZ, N, NRHS, B, LDB, BX, LDBX, U,
00022 *                          LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR,
00023 *                          GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK,
00024 *                          IWORK, INFO )
00025 * 
00026 *       .. Scalar Arguments ..
00027 *       INTEGER            ICOMPQ, INFO, LDB, LDBX, LDGCOL, LDU, N, NRHS,
00028 *      $                   SMLSIZ
00029 *       ..
00030 *       .. Array Arguments ..
00031 *       INTEGER            GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ),
00032 *      $                   K( * ), PERM( LDGCOL, * )
00033 *       REAL               B( LDB, * ), BX( LDBX, * ), C( * ),
00034 *      $                   DIFL( LDU, * ), DIFR( LDU, * ),
00035 *      $                   GIVNUM( LDU, * ), POLES( LDU, * ), S( * ),
00036 *      $                   U( LDU, * ), VT( LDU, * ), WORK( * ),
00037 *      $                   Z( LDU, * )
00038 *       ..
00039 *  
00040 *
00041 *> \par Purpose:
00042 *  =============
00043 *>
00044 *> \verbatim
00045 *>
00046 *> SLALSA is an itermediate step in solving the least squares problem
00047 *> by computing the SVD of the coefficient matrix in compact form (The
00048 *> singular vectors are computed as products of simple orthorgonal
00049 *> matrices.).
00050 *>
00051 *> If ICOMPQ = 0, SLALSA applies the inverse of the left singular vector
00052 *> matrix of an upper bidiagonal matrix to the right hand side; and if
00053 *> ICOMPQ = 1, SLALSA applies the right singular vector matrix to the
00054 *> right hand side. The singular vector matrices were generated in
00055 *> compact form by SLALSA.
00056 *> \endverbatim
00057 *
00058 *  Arguments:
00059 *  ==========
00060 *
00061 *> \param[in] ICOMPQ
00062 *> \verbatim
00063 *>          ICOMPQ is INTEGER
00064 *>         Specifies whether the left or the right singular vector
00065 *>         matrix is involved.
00066 *>         = 0: Left singular vector matrix
00067 *>         = 1: Right singular vector matrix
00068 *> \endverbatim
00069 *>
00070 *> \param[in] SMLSIZ
00071 *> \verbatim
00072 *>          SMLSIZ is INTEGER
00073 *>         The maximum size of the subproblems at the bottom of the
00074 *>         computation tree.
00075 *> \endverbatim
00076 *>
00077 *> \param[in] N
00078 *> \verbatim
00079 *>          N is INTEGER
00080 *>         The row and column dimensions of the upper bidiagonal matrix.
00081 *> \endverbatim
00082 *>
00083 *> \param[in] NRHS
00084 *> \verbatim
00085 *>          NRHS is INTEGER
00086 *>         The number of columns of B and BX. NRHS must be at least 1.
00087 *> \endverbatim
00088 *>
00089 *> \param[in,out] B
00090 *> \verbatim
00091 *>          B is REAL array, dimension ( LDB, NRHS )
00092 *>         On input, B contains the right hand sides of the least
00093 *>         squares problem in rows 1 through M.
00094 *>         On output, B contains the solution X in rows 1 through N.
00095 *> \endverbatim
00096 *>
00097 *> \param[in] LDB
00098 *> \verbatim
00099 *>          LDB is INTEGER
00100 *>         The leading dimension of B in the calling subprogram.
00101 *>         LDB must be at least max(1,MAX( M, N ) ).
00102 *> \endverbatim
00103 *>
00104 *> \param[out] BX
00105 *> \verbatim
00106 *>          BX is REAL array, dimension ( LDBX, NRHS )
00107 *>         On exit, the result of applying the left or right singular
00108 *>         vector matrix to B.
00109 *> \endverbatim
00110 *>
00111 *> \param[in] LDBX
00112 *> \verbatim
00113 *>          LDBX is INTEGER
00114 *>         The leading dimension of BX.
00115 *> \endverbatim
00116 *>
00117 *> \param[in] U
00118 *> \verbatim
00119 *>          U is REAL array, dimension ( LDU, SMLSIZ ).
00120 *>         On entry, U contains the left singular vector matrices of all
00121 *>         subproblems at the bottom level.
00122 *> \endverbatim
00123 *>
00124 *> \param[in] LDU
00125 *> \verbatim
00126 *>          LDU is INTEGER, LDU = > N.
00127 *>         The leading dimension of arrays U, VT, DIFL, DIFR,
00128 *>         POLES, GIVNUM, and Z.
00129 *> \endverbatim
00130 *>
00131 *> \param[in] VT
00132 *> \verbatim
00133 *>          VT is REAL array, dimension ( LDU, SMLSIZ+1 ).
00134 *>         On entry, VT**T contains the right singular vector matrices of
00135 *>         all subproblems at the bottom level.
00136 *> \endverbatim
00137 *>
00138 *> \param[in] K
00139 *> \verbatim
00140 *>          K is INTEGER array, dimension ( N ).
00141 *> \endverbatim
00142 *>
00143 *> \param[in] DIFL
00144 *> \verbatim
00145 *>          DIFL is REAL array, dimension ( LDU, NLVL ).
00146 *>         where NLVL = INT(log_2 (N/(SMLSIZ+1))) + 1.
00147 *> \endverbatim
00148 *>
00149 *> \param[in] DIFR
00150 *> \verbatim
00151 *>          DIFR is REAL array, dimension ( LDU, 2 * NLVL ).
00152 *>         On entry, DIFL(*, I) and DIFR(*, 2 * I -1) record
00153 *>         distances between singular values on the I-th level and
00154 *>         singular values on the (I -1)-th level, and DIFR(*, 2 * I)
00155 *>         record the normalizing factors of the right singular vectors
00156 *>         matrices of subproblems on I-th level.
00157 *> \endverbatim
00158 *>
00159 *> \param[in] Z
00160 *> \verbatim
00161 *>          Z is REAL array, dimension ( LDU, NLVL ).
00162 *>         On entry, Z(1, I) contains the components of the deflation-
00163 *>         adjusted updating row vector for subproblems on the I-th
00164 *>         level.
00165 *> \endverbatim
00166 *>
00167 *> \param[in] POLES
00168 *> \verbatim
00169 *>          POLES is REAL array, dimension ( LDU, 2 * NLVL ).
00170 *>         On entry, POLES(*, 2 * I -1: 2 * I) contains the new and old
00171 *>         singular values involved in the secular equations on the I-th
00172 *>         level.
00173 *> \endverbatim
00174 *>
00175 *> \param[in] GIVPTR
00176 *> \verbatim
00177 *>          GIVPTR is INTEGER array, dimension ( N ).
00178 *>         On entry, GIVPTR( I ) records the number of Givens
00179 *>         rotations performed on the I-th problem on the computation
00180 *>         tree.
00181 *> \endverbatim
00182 *>
00183 *> \param[in] GIVCOL
00184 *> \verbatim
00185 *>          GIVCOL is INTEGER array, dimension ( LDGCOL, 2 * NLVL ).
00186 *>         On entry, for each I, GIVCOL(*, 2 * I - 1: 2 * I) records the
00187 *>         locations of Givens rotations performed on the I-th level on
00188 *>         the computation tree.
00189 *> \endverbatim
00190 *>
00191 *> \param[in] LDGCOL
00192 *> \verbatim
00193 *>          LDGCOL is INTEGER, LDGCOL = > N.
00194 *>         The leading dimension of arrays GIVCOL and PERM.
00195 *> \endverbatim
00196 *>
00197 *> \param[in] PERM
00198 *> \verbatim
00199 *>          PERM is INTEGER array, dimension ( LDGCOL, NLVL ).
00200 *>         On entry, PERM(*, I) records permutations done on the I-th
00201 *>         level of the computation tree.
00202 *> \endverbatim
00203 *>
00204 *> \param[in] GIVNUM
00205 *> \verbatim
00206 *>          GIVNUM is REAL array, dimension ( LDU, 2 * NLVL ).
00207 *>         On entry, GIVNUM(*, 2 *I -1 : 2 * I) records the C- and S-
00208 *>         values of Givens rotations performed on the I-th level on the
00209 *>         computation tree.
00210 *> \endverbatim
00211 *>
00212 *> \param[in] C
00213 *> \verbatim
00214 *>          C is REAL array, dimension ( N ).
00215 *>         On entry, if the I-th subproblem is not square,
00216 *>         C( I ) contains the C-value of a Givens rotation related to
00217 *>         the right null space of the I-th subproblem.
00218 *> \endverbatim
00219 *>
00220 *> \param[in] S
00221 *> \verbatim
00222 *>          S is REAL array, dimension ( N ).
00223 *>         On entry, if the I-th subproblem is not square,
00224 *>         S( I ) contains the S-value of a Givens rotation related to
00225 *>         the right null space of the I-th subproblem.
00226 *> \endverbatim
00227 *>
00228 *> \param[out] WORK
00229 *> \verbatim
00230 *>          WORK is REAL array.
00231 *>         The dimension must be at least N.
00232 *> \endverbatim
00233 *>
00234 *> \param[out] IWORK
00235 *> \verbatim
00236 *>          IWORK is INTEGER array.
00237 *>         The dimension must be at least 3 * N
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 realOTHERcomputational
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 SLALSA( ICOMPQ, SMLSIZ, N, NRHS, B, LDB, BX, LDBX, U,
00268      $                   LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR,
00269      $                   GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK,
00270      $                   IWORK, INFO )
00271 *
00272 *  -- LAPACK computational routine (version 3.4.0) --
00273 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00274 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00275 *     November 2011
00276 *
00277 *     .. Scalar Arguments ..
00278       INTEGER            ICOMPQ, INFO, LDB, LDBX, LDGCOL, LDU, N, NRHS,
00279      $                   SMLSIZ
00280 *     ..
00281 *     .. Array Arguments ..
00282       INTEGER            GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ),
00283      $                   K( * ), PERM( LDGCOL, * )
00284       REAL               B( LDB, * ), BX( LDBX, * ), C( * ),
00285      $                   DIFL( LDU, * ), DIFR( LDU, * ),
00286      $                   GIVNUM( LDU, * ), POLES( LDU, * ), S( * ),
00287      $                   U( LDU, * ), VT( LDU, * ), WORK( * ),
00288      $                   Z( LDU, * )
00289 *     ..
00290 *
00291 *  =====================================================================
00292 *
00293 *     .. Parameters ..
00294       REAL               ZERO, ONE
00295       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
00296 *     ..
00297 *     .. Local Scalars ..
00298       INTEGER            I, I1, IC, IM1, INODE, J, LF, LL, LVL, LVL2,
00299      $                   ND, NDB1, NDIML, NDIMR, NL, NLF, NLP1, NLVL,
00300      $                   NR, NRF, NRP1, SQRE
00301 *     ..
00302 *     .. External Subroutines ..
00303       EXTERNAL           SCOPY, SGEMM, SLALS0, SLASDT, XERBLA
00304 *     ..
00305 *     .. Executable Statements ..
00306 *
00307 *     Test the input parameters.
00308 *
00309       INFO = 0
00310 *
00311       IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
00312          INFO = -1
00313       ELSE IF( SMLSIZ.LT.3 ) THEN
00314          INFO = -2
00315       ELSE IF( N.LT.SMLSIZ ) THEN
00316          INFO = -3
00317       ELSE IF( NRHS.LT.1 ) THEN
00318          INFO = -4
00319       ELSE IF( LDB.LT.N ) THEN
00320          INFO = -6
00321       ELSE IF( LDBX.LT.N ) THEN
00322          INFO = -8
00323       ELSE IF( LDU.LT.N ) THEN
00324          INFO = -10
00325       ELSE IF( LDGCOL.LT.N ) THEN
00326          INFO = -19
00327       END IF
00328       IF( INFO.NE.0 ) THEN
00329          CALL XERBLA( 'SLALSA', -INFO )
00330          RETURN
00331       END IF
00332 *
00333 *     Book-keeping and  setting up the computation tree.
00334 *
00335       INODE = 1
00336       NDIML = INODE + N
00337       NDIMR = NDIML + N
00338 *
00339       CALL SLASDT( N, NLVL, ND, IWORK( INODE ), IWORK( NDIML ),
00340      $             IWORK( NDIMR ), SMLSIZ )
00341 *
00342 *     The following code applies back the left singular vector factors.
00343 *     For applying back the right singular vector factors, go to 50.
00344 *
00345       IF( ICOMPQ.EQ.1 ) THEN
00346          GO TO 50
00347       END IF
00348 *
00349 *     The nodes on the bottom level of the tree were solved
00350 *     by SLASDQ. The corresponding left and right singular vector
00351 *     matrices are in explicit form. First apply back the left
00352 *     singular vector matrices.
00353 *
00354       NDB1 = ( ND+1 ) / 2
00355       DO 10 I = NDB1, ND
00356 *
00357 *        IC : center row of each node
00358 *        NL : number of rows of left  subproblem
00359 *        NR : number of rows of right subproblem
00360 *        NLF: starting row of the left   subproblem
00361 *        NRF: starting row of the right  subproblem
00362 *
00363          I1 = I - 1
00364          IC = IWORK( INODE+I1 )
00365          NL = IWORK( NDIML+I1 )
00366          NR = IWORK( NDIMR+I1 )
00367          NLF = IC - NL
00368          NRF = IC + 1
00369          CALL SGEMM( 'T', 'N', NL, NRHS, NL, ONE, U( NLF, 1 ), LDU,
00370      $               B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX )
00371          CALL SGEMM( 'T', 'N', NR, NRHS, NR, ONE, U( NRF, 1 ), LDU,
00372      $               B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX )
00373    10 CONTINUE
00374 *
00375 *     Next copy the rows of B that correspond to unchanged rows
00376 *     in the bidiagonal matrix to BX.
00377 *
00378       DO 20 I = 1, ND
00379          IC = IWORK( INODE+I-1 )
00380          CALL SCOPY( NRHS, B( IC, 1 ), LDB, BX( IC, 1 ), LDBX )
00381    20 CONTINUE
00382 *
00383 *     Finally go through the left singular vector matrices of all
00384 *     the other subproblems bottom-up on the tree.
00385 *
00386       J = 2**NLVL
00387       SQRE = 0
00388 *
00389       DO 40 LVL = NLVL, 1, -1
00390          LVL2 = 2*LVL - 1
00391 *
00392 *        find the first node LF and last node LL on
00393 *        the current level LVL
00394 *
00395          IF( LVL.EQ.1 ) THEN
00396             LF = 1
00397             LL = 1
00398          ELSE
00399             LF = 2**( LVL-1 )
00400             LL = 2*LF - 1
00401          END IF
00402          DO 30 I = LF, LL
00403             IM1 = I - 1
00404             IC = IWORK( INODE+IM1 )
00405             NL = IWORK( NDIML+IM1 )
00406             NR = IWORK( NDIMR+IM1 )
00407             NLF = IC - NL
00408             NRF = IC + 1
00409             J = J - 1
00410             CALL SLALS0( ICOMPQ, NL, NR, SQRE, NRHS, BX( NLF, 1 ), LDBX,
00411      $                   B( NLF, 1 ), LDB, PERM( NLF, LVL ),
00412      $                   GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL,
00413      $                   GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ),
00414      $                   DIFL( NLF, LVL ), DIFR( NLF, LVL2 ),
00415      $                   Z( NLF, LVL ), K( J ), C( J ), S( J ), WORK,
00416      $                   INFO )
00417    30    CONTINUE
00418    40 CONTINUE
00419       GO TO 90
00420 *
00421 *     ICOMPQ = 1: applying back the right singular vector factors.
00422 *
00423    50 CONTINUE
00424 *
00425 *     First now go through the right singular vector matrices of all
00426 *     the tree nodes top-down.
00427 *
00428       J = 0
00429       DO 70 LVL = 1, NLVL
00430          LVL2 = 2*LVL - 1
00431 *
00432 *        Find the first node LF and last node LL on
00433 *        the current level LVL.
00434 *
00435          IF( LVL.EQ.1 ) THEN
00436             LF = 1
00437             LL = 1
00438          ELSE
00439             LF = 2**( LVL-1 )
00440             LL = 2*LF - 1
00441          END IF
00442          DO 60 I = LL, LF, -1
00443             IM1 = I - 1
00444             IC = IWORK( INODE+IM1 )
00445             NL = IWORK( NDIML+IM1 )
00446             NR = IWORK( NDIMR+IM1 )
00447             NLF = IC - NL
00448             NRF = IC + 1
00449             IF( I.EQ.LL ) THEN
00450                SQRE = 0
00451             ELSE
00452                SQRE = 1
00453             END IF
00454             J = J + 1
00455             CALL SLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B( NLF, 1 ), LDB,
00456      $                   BX( NLF, 1 ), LDBX, PERM( NLF, LVL ),
00457      $                   GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL,
00458      $                   GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ),
00459      $                   DIFL( NLF, LVL ), DIFR( NLF, LVL2 ),
00460      $                   Z( NLF, LVL ), K( J ), C( J ), S( J ), WORK,
00461      $                   INFO )
00462    60    CONTINUE
00463    70 CONTINUE
00464 *
00465 *     The nodes on the bottom level of the tree were solved
00466 *     by SLASDQ. The corresponding right singular vector
00467 *     matrices are in explicit form. Apply them back.
00468 *
00469       NDB1 = ( ND+1 ) / 2
00470       DO 80 I = NDB1, ND
00471          I1 = I - 1
00472          IC = IWORK( INODE+I1 )
00473          NL = IWORK( NDIML+I1 )
00474          NR = IWORK( NDIMR+I1 )
00475          NLP1 = NL + 1
00476          IF( I.EQ.ND ) THEN
00477             NRP1 = NR
00478          ELSE
00479             NRP1 = NR + 1
00480          END IF
00481          NLF = IC - NL
00482          NRF = IC + 1
00483          CALL SGEMM( 'T', 'N', NLP1, NRHS, NLP1, ONE, VT( NLF, 1 ), LDU,
00484      $               B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX )
00485          CALL SGEMM( 'T', 'N', NRP1, NRHS, NRP1, ONE, VT( NRF, 1 ), LDU,
00486      $               B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX )
00487    80 CONTINUE
00488 *
00489    90 CONTINUE
00490 *
00491       RETURN
00492 *
00493 *     End of SLALSA
00494 *
00495       END
 All Files Functions