LAPACK 3.3.1 Linear Algebra PACKage

# dlalsa.f

Go to the documentation of this file.
```00001       SUBROUTINE DLALSA( ICOMPQ, SMLSIZ, N, NRHS, B, LDB, BX, LDBX, U,
00002      \$                   LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR,
00003      \$                   GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK,
00004      \$                   IWORK, INFO )
00005 *
00006 *  -- LAPACK routine (version 3.2) --
00007 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00008 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00009 *     November 2006
00010 *
00011 *     .. Scalar Arguments ..
00012       INTEGER            ICOMPQ, INFO, LDB, LDBX, LDGCOL, LDU, N, NRHS,
00013      \$                   SMLSIZ
00014 *     ..
00015 *     .. Array Arguments ..
00016       INTEGER            GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ),
00017      \$                   K( * ), PERM( LDGCOL, * )
00018       DOUBLE PRECISION   B( LDB, * ), BX( LDBX, * ), C( * ),
00019      \$                   DIFL( LDU, * ), DIFR( LDU, * ),
00020      \$                   GIVNUM( LDU, * ), POLES( LDU, * ), S( * ),
00021      \$                   U( LDU, * ), VT( LDU, * ), WORK( * ),
00022      \$                   Z( LDU, * )
00023 *     ..
00024 *
00025 *  Purpose
00026 *  =======
00027 *
00028 *  DLALSA is an itermediate step in solving the least squares problem
00029 *  by computing the SVD of the coefficient matrix in compact form (The
00030 *  singular vectors are computed as products of simple orthorgonal
00031 *  matrices.).
00032 *
00033 *  If ICOMPQ = 0, DLALSA applies the inverse of the left singular vector
00034 *  matrix of an upper bidiagonal matrix to the right hand side; and if
00035 *  ICOMPQ = 1, DLALSA applies the right singular vector matrix to the
00036 *  right hand side. The singular vector matrices were generated in
00037 *  compact form by DLALSA.
00038 *
00039 *  Arguments
00040 *  =========
00041 *
00042 *
00043 *  ICOMPQ (input) INTEGER
00044 *         Specifies whether the left or the right singular vector
00045 *         matrix is involved.
00046 *         = 0: Left singular vector matrix
00047 *         = 1: Right singular vector matrix
00048 *
00049 *  SMLSIZ (input) INTEGER
00050 *         The maximum size of the subproblems at the bottom of the
00051 *         computation tree.
00052 *
00053 *  N      (input) INTEGER
00054 *         The row and column dimensions of the upper bidiagonal matrix.
00055 *
00056 *  NRHS   (input) INTEGER
00057 *         The number of columns of B and BX. NRHS must be at least 1.
00058 *
00059 *  B      (input/output) DOUBLE PRECISION array, dimension ( LDB, NRHS )
00060 *         On input, B contains the right hand sides of the least
00061 *         squares problem in rows 1 through M.
00062 *         On output, B contains the solution X in rows 1 through N.
00063 *
00064 *  LDB    (input) INTEGER
00065 *         The leading dimension of B in the calling subprogram.
00066 *         LDB must be at least max(1,MAX( M, N ) ).
00067 *
00068 *  BX     (output) DOUBLE PRECISION array, dimension ( LDBX, NRHS )
00069 *         On exit, the result of applying the left or right singular
00070 *         vector matrix to B.
00071 *
00072 *  LDBX   (input) INTEGER
00073 *         The leading dimension of BX.
00074 *
00075 *  U      (input) DOUBLE PRECISION array, dimension ( LDU, SMLSIZ ).
00076 *         On entry, U contains the left singular vector matrices of all
00077 *         subproblems at the bottom level.
00078 *
00079 *  LDU    (input) INTEGER, LDU = > N.
00080 *         The leading dimension of arrays U, VT, DIFL, DIFR,
00081 *         POLES, GIVNUM, and Z.
00082 *
00083 *  VT     (input) DOUBLE PRECISION array, dimension ( LDU, SMLSIZ+1 ).
00084 *         On entry, VT**T contains the right singular vector matrices of
00085 *         all subproblems at the bottom level.
00086 *
00087 *  K      (input) INTEGER array, dimension ( N ).
00088 *
00089 *  DIFL   (input) DOUBLE PRECISION array, dimension ( LDU, NLVL ).
00090 *         where NLVL = INT(log_2 (N/(SMLSIZ+1))) + 1.
00091 *
00092 *  DIFR   (input) DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ).
00093 *         On entry, DIFL(*, I) and DIFR(*, 2 * I -1) record
00094 *         distances between singular values on the I-th level and
00095 *         singular values on the (I -1)-th level, and DIFR(*, 2 * I)
00096 *         record the normalizing factors of the right singular vectors
00097 *         matrices of subproblems on I-th level.
00098 *
00099 *  Z      (input) DOUBLE PRECISION array, dimension ( LDU, NLVL ).
00100 *         On entry, Z(1, I) contains the components of the deflation-
00101 *         adjusted updating row vector for subproblems on the I-th
00102 *         level.
00103 *
00104 *  POLES  (input) DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ).
00105 *         On entry, POLES(*, 2 * I -1: 2 * I) contains the new and old
00106 *         singular values involved in the secular equations on the I-th
00107 *         level.
00108 *
00109 *  GIVPTR (input) INTEGER array, dimension ( N ).
00110 *         On entry, GIVPTR( I ) records the number of Givens
00111 *         rotations performed on the I-th problem on the computation
00112 *         tree.
00113 *
00114 *  GIVCOL (input) INTEGER array, dimension ( LDGCOL, 2 * NLVL ).
00115 *         On entry, for each I, GIVCOL(*, 2 * I - 1: 2 * I) records the
00116 *         locations of Givens rotations performed on the I-th level on
00117 *         the computation tree.
00118 *
00119 *  LDGCOL (input) INTEGER, LDGCOL = > N.
00120 *         The leading dimension of arrays GIVCOL and PERM.
00121 *
00122 *  PERM   (input) INTEGER array, dimension ( LDGCOL, NLVL ).
00123 *         On entry, PERM(*, I) records permutations done on the I-th
00124 *         level of the computation tree.
00125 *
00126 *  GIVNUM (input) DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ).
00127 *         On entry, GIVNUM(*, 2 *I -1 : 2 * I) records the C- and S-
00128 *         values of Givens rotations performed on the I-th level on the
00129 *         computation tree.
00130 *
00131 *  C      (input) DOUBLE PRECISION array, dimension ( N ).
00132 *         On entry, if the I-th subproblem is not square,
00133 *         C( I ) contains the C-value of a Givens rotation related to
00134 *         the right null space of the I-th subproblem.
00135 *
00136 *  S      (input) DOUBLE PRECISION array, dimension ( N ).
00137 *         On entry, if the I-th subproblem is not square,
00138 *         S( I ) contains the S-value of a Givens rotation related to
00139 *         the right null space of the I-th subproblem.
00140 *
00141 *  WORK   (workspace) DOUBLE PRECISION array.
00142 *         The dimension must be at least N.
00143 *
00144 *  IWORK  (workspace) INTEGER array.
00145 *         The dimension must be at least 3 * N
00146 *
00147 *  INFO   (output) INTEGER
00148 *          = 0:  successful exit.
00149 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00150 *
00151 *  Further Details
00152 *  ===============
00153 *
00154 *  Based on contributions by
00155 *     Ming Gu and Ren-Cang Li, Computer Science Division, University of
00156 *       California at Berkeley, USA
00157 *     Osni Marques, LBNL/NERSC, USA
00158 *
00159 *  =====================================================================
00160 *
00161 *     .. Parameters ..
00162       DOUBLE PRECISION   ZERO, ONE
00163       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
00164 *     ..
00165 *     .. Local Scalars ..
00166       INTEGER            I, I1, IC, IM1, INODE, J, LF, LL, LVL, LVL2,
00167      \$                   ND, NDB1, NDIML, NDIMR, NL, NLF, NLP1, NLVL,
00168      \$                   NR, NRF, NRP1, SQRE
00169 *     ..
00170 *     .. External Subroutines ..
00171       EXTERNAL           DCOPY, DGEMM, DLALS0, DLASDT, XERBLA
00172 *     ..
00173 *     .. Executable Statements ..
00174 *
00175 *     Test the input parameters.
00176 *
00177       INFO = 0
00178 *
00179       IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
00180          INFO = -1
00181       ELSE IF( SMLSIZ.LT.3 ) THEN
00182          INFO = -2
00183       ELSE IF( N.LT.SMLSIZ ) THEN
00184          INFO = -3
00185       ELSE IF( NRHS.LT.1 ) THEN
00186          INFO = -4
00187       ELSE IF( LDB.LT.N ) THEN
00188          INFO = -6
00189       ELSE IF( LDBX.LT.N ) THEN
00190          INFO = -8
00191       ELSE IF( LDU.LT.N ) THEN
00192          INFO = -10
00193       ELSE IF( LDGCOL.LT.N ) THEN
00194          INFO = -19
00195       END IF
00196       IF( INFO.NE.0 ) THEN
00197          CALL XERBLA( 'DLALSA', -INFO )
00198          RETURN
00199       END IF
00200 *
00201 *     Book-keeping and  setting up the computation tree.
00202 *
00203       INODE = 1
00204       NDIML = INODE + N
00205       NDIMR = NDIML + N
00206 *
00207       CALL DLASDT( N, NLVL, ND, IWORK( INODE ), IWORK( NDIML ),
00208      \$             IWORK( NDIMR ), SMLSIZ )
00209 *
00210 *     The following code applies back the left singular vector factors.
00211 *     For applying back the right singular vector factors, go to 50.
00212 *
00213       IF( ICOMPQ.EQ.1 ) THEN
00214          GO TO 50
00215       END IF
00216 *
00217 *     The nodes on the bottom level of the tree were solved
00218 *     by DLASDQ. The corresponding left and right singular vector
00219 *     matrices are in explicit form. First apply back the left
00220 *     singular vector matrices.
00221 *
00222       NDB1 = ( ND+1 ) / 2
00223       DO 10 I = NDB1, ND
00224 *
00225 *        IC : center row of each node
00226 *        NL : number of rows of left  subproblem
00227 *        NR : number of rows of right subproblem
00228 *        NLF: starting row of the left   subproblem
00229 *        NRF: starting row of the right  subproblem
00230 *
00231          I1 = I - 1
00232          IC = IWORK( INODE+I1 )
00233          NL = IWORK( NDIML+I1 )
00234          NR = IWORK( NDIMR+I1 )
00235          NLF = IC - NL
00236          NRF = IC + 1
00237          CALL DGEMM( 'T', 'N', NL, NRHS, NL, ONE, U( NLF, 1 ), LDU,
00238      \$               B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX )
00239          CALL DGEMM( 'T', 'N', NR, NRHS, NR, ONE, U( NRF, 1 ), LDU,
00240      \$               B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX )
00241    10 CONTINUE
00242 *
00243 *     Next copy the rows of B that correspond to unchanged rows
00244 *     in the bidiagonal matrix to BX.
00245 *
00246       DO 20 I = 1, ND
00247          IC = IWORK( INODE+I-1 )
00248          CALL DCOPY( NRHS, B( IC, 1 ), LDB, BX( IC, 1 ), LDBX )
00249    20 CONTINUE
00250 *
00251 *     Finally go through the left singular vector matrices of all
00252 *     the other subproblems bottom-up on the tree.
00253 *
00254       J = 2**NLVL
00255       SQRE = 0
00256 *
00257       DO 40 LVL = NLVL, 1, -1
00258          LVL2 = 2*LVL - 1
00259 *
00260 *        find the first node LF and last node LL on
00261 *        the current level LVL
00262 *
00263          IF( LVL.EQ.1 ) THEN
00264             LF = 1
00265             LL = 1
00266          ELSE
00267             LF = 2**( LVL-1 )
00268             LL = 2*LF - 1
00269          END IF
00270          DO 30 I = LF, LL
00271             IM1 = I - 1
00272             IC = IWORK( INODE+IM1 )
00273             NL = IWORK( NDIML+IM1 )
00274             NR = IWORK( NDIMR+IM1 )
00275             NLF = IC - NL
00276             NRF = IC + 1
00277             J = J - 1
00278             CALL DLALS0( ICOMPQ, NL, NR, SQRE, NRHS, BX( NLF, 1 ), LDBX,
00279      \$                   B( NLF, 1 ), LDB, PERM( NLF, LVL ),
00280      \$                   GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL,
00281      \$                   GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ),
00282      \$                   DIFL( NLF, LVL ), DIFR( NLF, LVL2 ),
00283      \$                   Z( NLF, LVL ), K( J ), C( J ), S( J ), WORK,
00284      \$                   INFO )
00285    30    CONTINUE
00286    40 CONTINUE
00287       GO TO 90
00288 *
00289 *     ICOMPQ = 1: applying back the right singular vector factors.
00290 *
00291    50 CONTINUE
00292 *
00293 *     First now go through the right singular vector matrices of all
00294 *     the tree nodes top-down.
00295 *
00296       J = 0
00297       DO 70 LVL = 1, NLVL
00298          LVL2 = 2*LVL - 1
00299 *
00300 *        Find the first node LF and last node LL on
00301 *        the current level LVL.
00302 *
00303          IF( LVL.EQ.1 ) THEN
00304             LF = 1
00305             LL = 1
00306          ELSE
00307             LF = 2**( LVL-1 )
00308             LL = 2*LF - 1
00309          END IF
00310          DO 60 I = LL, LF, -1
00311             IM1 = I - 1
00312             IC = IWORK( INODE+IM1 )
00313             NL = IWORK( NDIML+IM1 )
00314             NR = IWORK( NDIMR+IM1 )
00315             NLF = IC - NL
00316             NRF = IC + 1
00317             IF( I.EQ.LL ) THEN
00318                SQRE = 0
00319             ELSE
00320                SQRE = 1
00321             END IF
00322             J = J + 1
00323             CALL DLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B( NLF, 1 ), LDB,
00324      \$                   BX( NLF, 1 ), LDBX, PERM( NLF, LVL ),
00325      \$                   GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL,
00326      \$                   GIVNUM( NLF, LVL2 ), LDU, POLES( NLF, LVL2 ),
00327      \$                   DIFL( NLF, LVL ), DIFR( NLF, LVL2 ),
00328      \$                   Z( NLF, LVL ), K( J ), C( J ), S( J ), WORK,
00329      \$                   INFO )
00330    60    CONTINUE
00331    70 CONTINUE
00332 *
00333 *     The nodes on the bottom level of the tree were solved
00334 *     by DLASDQ. The corresponding right singular vector
00335 *     matrices are in explicit form. Apply them back.
00336 *
00337       NDB1 = ( ND+1 ) / 2
00338       DO 80 I = NDB1, ND
00339          I1 = I - 1
00340          IC = IWORK( INODE+I1 )
00341          NL = IWORK( NDIML+I1 )
00342          NR = IWORK( NDIMR+I1 )
00343          NLP1 = NL + 1
00344          IF( I.EQ.ND ) THEN
00345             NRP1 = NR
00346          ELSE
00347             NRP1 = NR + 1
00348          END IF
00349          NLF = IC - NL
00350          NRF = IC + 1
00351          CALL DGEMM( 'T', 'N', NLP1, NRHS, NLP1, ONE, VT( NLF, 1 ), LDU,
00352      \$               B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX )
00353          CALL DGEMM( 'T', 'N', NRP1, NRHS, NRP1, ONE, VT( NRF, 1 ), LDU,
00354      \$               B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX )
00355    80 CONTINUE
00356 *
00357    90 CONTINUE
00358 *
00359       RETURN
00360 *
00361 *     End of DLALSA
00362 *
00363       END
```