SRC/dlasd1.f

   1 *> \brief \b DLASD1 computes the SVD of an upper bidiagonal matrix B of the specified size. Used by sbdsdc.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DLASD1 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasd1.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasd1.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasd1.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DLASD1( NL, NR, SQRE, D, ALPHA, BETA, U, LDU, VT, LDVT,
  22 *                          IDXQ, IWORK, WORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       INTEGER            INFO, LDU, LDVT, NL, NR, SQRE
  26 *       DOUBLE PRECISION   ALPHA, BETA
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       INTEGER            IDXQ( * ), IWORK( * )
  30 *       DOUBLE PRECISION   D( * ), U( LDU, * ), VT( LDVT, * ), WORK( * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> DLASD1 computes the SVD of an upper bidiagonal N-by-M matrix B,
  40 *> where N = NL + NR + 1 and M = N + SQRE. DLASD1 is called from DLASD0.
  41 *>
  42 *> A related subroutine DLASD7 handles the case in which the singular
  43 *> values (and the singular vectors in factored form) are desired.
  44 *>
  45 *> DLASD1 computes the SVD as follows:
  46 *>
  47 *>               ( D1(in)    0    0       0 )
  48 *>   B = U(in) * (   Z1**T   a   Z2**T    b ) * VT(in)
  49 *>               (   0       0   D2(in)   0 )
  50 *>
  51 *>     = U(out) * ( D(out) 0) * VT(out)
  52 *>
  53 *> where Z**T = (Z1**T a Z2**T b) = u**T VT**T, and u is a vector of dimension M
  54 *> with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros
  55 *> elsewhere; and the entry b is empty if SQRE = 0.
  56 *>
  57 *> The left singular vectors of the original matrix are stored in U, and
  58 *> the transpose of the right singular vectors are stored in VT, and the
  59 *> singular values are in D.  The algorithm consists of three stages:
  60 *>
  61 *>    The first stage consists of deflating the size of the problem
  62 *>    when there are multiple singular values or when there are zeros in
  63 *>    the Z vector.  For each such occurrence the dimension of the
  64 *>    secular equation problem is reduced by one.  This stage is
  65 *>    performed by the routine DLASD2.
  66 *>
  67 *>    The second stage consists of calculating the updated
  68 *>    singular values. This is done by finding the square roots of the
  69 *>    roots of the secular equation via the routine DLASD4 (as called
  70 *>    by DLASD3). This routine also calculates the singular vectors of
  71 *>    the current problem.
  72 *>
  73 *>    The final stage consists of computing the updated singular vectors
  74 *>    directly using the updated singular values.  The singular vectors
  75 *>    for the current problem are multiplied with the singular vectors
  76 *>    from the overall problem.
  77 *> \endverbatim
  78 *
  79 *  Arguments:
  80 *  ==========
  81 *
  82 *> \param[in] NL
  83 *> \verbatim
  84 *>          NL is INTEGER
  85 *>         The row dimension of the upper block.  NL >= 1.
  86 *> \endverbatim
  87 *>
  88 *> \param[in] NR
  89 *> \verbatim
  90 *>          NR is INTEGER
  91 *>         The row dimension of the lower block.  NR >= 1.
  92 *> \endverbatim
  93 *>
  94 *> \param[in] SQRE
  95 *> \verbatim
  96 *>          SQRE is INTEGER
  97 *>         = 0: the lower block is an NR-by-NR square matrix.
  98 *>         = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
  99 *>
 100 *>         The bidiagonal matrix has row dimension N = NL + NR + 1,
 101 *>         and column dimension M = N + SQRE.
 102 *> \endverbatim
 103 *>
 104 *> \param[in,out] D
 105 *> \verbatim
 106 *>          D is DOUBLE PRECISION array,
 107 *>                        dimension (N = NL+NR+1).
 108 *>         On entry D(1:NL,1:NL) contains the singular values of the
 109 *>         upper block; and D(NL+2:N) contains the singular values of
 110 *>         the lower block. On exit D(1:N) contains the singular values
 111 *>         of the modified matrix.
 112 *> \endverbatim
 113 *>
 114 *> \param[in,out] ALPHA
 115 *> \verbatim
 116 *>          ALPHA is DOUBLE PRECISION
 117 *>         Contains the diagonal element associated with the added row.
 118 *> \endverbatim
 119 *>
 120 *> \param[in,out] BETA
 121 *> \verbatim
 122 *>          BETA is DOUBLE PRECISION
 123 *>         Contains the off-diagonal element associated with the added
 124 *>         row.
 125 *> \endverbatim
 126 *>
 127 *> \param[in,out] U
 128 *> \verbatim
 129 *>          U is DOUBLE PRECISION array, dimension(LDU,N)
 130 *>         On entry U(1:NL, 1:NL) contains the left singular vectors of
 131 *>         the upper block; U(NL+2:N, NL+2:N) contains the left singular
 132 *>         vectors of the lower block. On exit U contains the left
 133 *>         singular vectors of the bidiagonal matrix.
 134 *> \endverbatim
 135 *>
 136 *> \param[in] LDU
 137 *> \verbatim
 138 *>          LDU is INTEGER
 139 *>         The leading dimension of the array U.  LDU >= max( 1, N ).
 140 *> \endverbatim
 141 *>
 142 *> \param[in,out] VT
 143 *> \verbatim
 144 *>          VT is DOUBLE PRECISION array, dimension(LDVT,M)
 145 *>         where M = N + SQRE.
 146 *>         On entry VT(1:NL+1, 1:NL+1)**T contains the right singular
 147 *>         vectors of the upper block; VT(NL+2:M, NL+2:M)**T contains
 148 *>         the right singular vectors of the lower block. On exit
 149 *>         VT**T contains the right singular vectors of the
 150 *>         bidiagonal matrix.
 151 *> \endverbatim
 152 *>
 153 *> \param[in] LDVT
 154 *> \verbatim
 155 *>          LDVT is INTEGER
 156 *>         The leading dimension of the array VT.  LDVT >= max( 1, M ).
 157 *> \endverbatim
 158 *>
 159 *> \param[in,out] IDXQ
 160 *> \verbatim
 161 *>          IDXQ is INTEGER array, dimension(N)
 162 *>         This contains the permutation which will reintegrate the
 163 *>         subproblem just solved back into sorted order, i.e.
 164 *>         D( IDXQ( I = 1, N ) ) will be in ascending order.
 165 *> \endverbatim
 166 *>
 167 *> \param[out] IWORK
 168 *> \verbatim
 169 *>          IWORK is INTEGER array, dimension( 4 * N )
 170 *> \endverbatim
 171 *>
 172 *> \param[out] WORK
 173 *> \verbatim
 174 *>          WORK is DOUBLE PRECISION array, dimension( 3*M**2 + 2*M )
 175 *> \endverbatim
 176 *>
 177 *> \param[out] INFO
 178 *> \verbatim
 179 *>          INFO is INTEGER
 180 *>          = 0:  successful exit.
 181 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 182 *>          > 0:  if INFO = 1, a singular value did not converge
 183 *> \endverbatim
 184 *
 185 *  Authors:
 186 *  ========
 187 *
 188 *> \author Univ. of Tennessee
 189 *> \author Univ. of California Berkeley
 190 *> \author Univ. of Colorado Denver
 191 *> \author NAG Ltd.
 192 *
 193 *> \date June 2016
 194 *
 195 *> \ingroup auxOTHERauxiliary
 196 *
 197 *> \par Contributors:
 198 *  ==================
 199 *>
 200 *>     Ming Gu and Huan Ren, Computer Science Division, University of
 201 *>     California at Berkeley, USA
 202 *>
 203 *  =====================================================================
 204       SUBROUTINE DLASD1( NL, NR, SQRE, D, ALPHA, BETA, U, LDU, VT, LDVT,
 205      $                   IDXQ, IWORK, WORK, INFO )
 206 *
 207 *  -- LAPACK auxiliary routine (version 3.6.1) --
 208 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 209 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 210 *     June 2016
 211 *
 212 *     .. Scalar Arguments ..
 213       INTEGER            INFO, LDU, LDVT, NL, NR, SQRE
 214       DOUBLE PRECISION   ALPHA, BETA
 215 *     ..
 216 *     .. Array Arguments ..
 217       INTEGER            IDXQ( * ), IWORK( * )
 218       DOUBLE PRECISION   D( * ), U( LDU, * ), VT( LDVT, * ), WORK( * )
 219 *     ..
 220 *
 221 *  =====================================================================
 222 *
 223 *     .. Parameters ..
 224 *
 225       DOUBLE PRECISION   ONE, ZERO
 226       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 227 *     ..
 228 *     .. Local Scalars ..
 229       INTEGER            COLTYP, I, IDX, IDXC, IDXP, IQ, ISIGMA, IU2,
 230      $                   IVT2, IZ, K, LDQ, LDU2, LDVT2, M, N, N1, N2
 231       DOUBLE PRECISION   ORGNRM
 232 *     ..
 233 *     .. External Subroutines ..
 234       EXTERNAL           DLAMRG, DLASCL, DLASD2, DLASD3, XERBLA
 235 *     ..
 236 *     .. Intrinsic Functions ..
 237       INTRINSIC          ABS, MAX
 238 *     ..
 239 *     .. Executable Statements ..
 240 *
 241 *     Test the input parameters.
 242 *
 243       INFO = 0
 244 *
 245       IF( NL.LT.1 ) THEN
 246          INFO = -1
 247       ELSE IF( NR.LT.1 ) THEN
 248          INFO = -2
 249       ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
 250          INFO = -3
 251       END IF
 252       IF( INFO.NE.0 ) THEN
 253          CALL XERBLA( 'DLASD1', -INFO )
 254          RETURN
 255       END IF
 256 *
 257       N = NL + NR + 1
 258       M = N + SQRE
 259 *
 260 *     The following values are for bookkeeping purposes only.  They are
 261 *     integer pointers which indicate the portion of the workspace
 262 *     used by a particular array in DLASD2 and DLASD3.
 263 *
 264       LDU2 = N
 265       LDVT2 = M
 266 *
 267       IZ = 1
 268       ISIGMA = IZ + M
 269       IU2 = ISIGMA + N
 270       IVT2 = IU2 + LDU2*N
 271       IQ = IVT2 + LDVT2*M
 272 *
 273       IDX = 1
 274       IDXC = IDX + N
 275       COLTYP = IDXC + N
 276       IDXP = COLTYP + N
 277 *
 278 *     Scale.
 279 *
 280       ORGNRM = MAX( ABS( ALPHA ), ABS( BETA ) )
 281       D( NL+1 ) = ZERO
 282       DO 10 I = 1, N
 283          IF( ABS( D( I ) ).GT.ORGNRM ) THEN
 284             ORGNRM = ABS( D( I ) )
 285          END IF
 286    10 CONTINUE
 287       CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
 288       ALPHA = ALPHA / ORGNRM
 289       BETA = BETA / ORGNRM
 290 *
 291 *     Deflate singular values.
 292 *
 293       CALL DLASD2( NL, NR, SQRE, K, D, WORK( IZ ), ALPHA, BETA, U, LDU,
 294      $             VT, LDVT, WORK( ISIGMA ), WORK( IU2 ), LDU2,
 295      $             WORK( IVT2 ), LDVT2, IWORK( IDXP ), IWORK( IDX ),
 296      $             IWORK( IDXC ), IDXQ, IWORK( COLTYP ), INFO )
 297 *
 298 *     Solve Secular Equation and update singular vectors.
 299 *
 300       LDQ = K
 301       CALL DLASD3( NL, NR, SQRE, K, D, WORK( IQ ), LDQ, WORK( ISIGMA ),
 302      $             U, LDU, WORK( IU2 ), LDU2, VT, LDVT, WORK( IVT2 ),
 303      $             LDVT2, IWORK( IDXC ), IWORK( COLTYP ), WORK( IZ ),
 304      $             INFO )
 305 *
 306 *     Report the convergence failure.
 307 *
 308       IF( INFO.NE.0 ) THEN
 309          RETURN
 310       END IF
 311 *
 312 *     Unscale.
 313 *
 314       CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
 315 *
 316 *     Prepare the IDXQ sorting permutation.
 317 *
 318       N1 = K
 319       N2 = N - K
 320       CALL DLAMRG( N1, N2, D, 1, -1, IDXQ )
 321 *
 322       RETURN
 323 *
 324 *     End of DLASD1
 325 *
 326       END