SRC/dlalsd.f

   1 *> \brief \b DLALSD uses the singular value decomposition of A to solve the least squares problem.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DLALSD + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlalsd.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlalsd.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlalsd.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DLALSD( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
  22 *                          RANK, WORK, IWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          UPLO
  26 *       INTEGER            INFO, LDB, N, NRHS, RANK, SMLSIZ
  27 *       DOUBLE PRECISION   RCOND
  28 *       ..
  29 *       .. Array Arguments ..
  30 *       INTEGER            IWORK( * )
  31 *       DOUBLE PRECISION   B( LDB, * ), D( * ), E( * ), WORK( * )
  32 *       ..
  33 *
  34 *
  35 *> \par Purpose:
  36 *  =============
  37 *>
  38 *> \verbatim
  39 *>
  40 *> DLALSD uses the singular value decomposition of A to solve the least
  41 *> squares problem of finding X to minimize the Euclidean norm of each
  42 *> column of A*X-B, where A is N-by-N upper bidiagonal, and X and B
  43 *> are N-by-NRHS. The solution X overwrites B.
  44 *>
  45 *> The singular values of A smaller than RCOND times the largest
  46 *> singular value are treated as zero in solving the least squares
  47 *> problem; in this case a minimum norm solution is returned.
  48 *> The actual singular values are returned in D in ascending order.
  49 *>
  50 *> This code makes very mild assumptions about floating point
  51 *> arithmetic. It will work on machines with a guard digit in
  52 *> add/subtract, or on those binary machines without guard digits
  53 *> which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
  54 *> It could conceivably fail on hexadecimal or decimal machines
  55 *> without guard digits, but we know of none.
  56 *> \endverbatim
  57 *
  58 *  Arguments:
  59 *  ==========
  60 *
  61 *> \param[in] UPLO
  62 *> \verbatim
  63 *>          UPLO is CHARACTER*1
  64 *>         = 'U': D and E define an upper bidiagonal matrix.
  65 *>         = 'L': D and E define a  lower bidiagonal matrix.
  66 *> \endverbatim
  67 *>
  68 *> \param[in] SMLSIZ
  69 *> \verbatim
  70 *>          SMLSIZ is INTEGER
  71 *>         The maximum size of the subproblems at the bottom of the
  72 *>         computation tree.
  73 *> \endverbatim
  74 *>
  75 *> \param[in] N
  76 *> \verbatim
  77 *>          N is INTEGER
  78 *>         The dimension of the  bidiagonal matrix.  N >= 0.
  79 *> \endverbatim
  80 *>
  81 *> \param[in] NRHS
  82 *> \verbatim
  83 *>          NRHS is INTEGER
  84 *>         The number of columns of B. NRHS must be at least 1.
  85 *> \endverbatim
  86 *>
  87 *> \param[in,out] D
  88 *> \verbatim
  89 *>          D is DOUBLE PRECISION array, dimension (N)
  90 *>         On entry D contains the main diagonal of the bidiagonal
  91 *>         matrix. On exit, if INFO = 0, D contains its singular values.
  92 *> \endverbatim
  93 *>
  94 *> \param[in,out] E
  95 *> \verbatim
  96 *>          E is DOUBLE PRECISION array, dimension (N-1)
  97 *>         Contains the super-diagonal entries of the bidiagonal matrix.
  98 *>         On exit, E has been destroyed.
  99 *> \endverbatim
 100 *>
 101 *> \param[in,out] B
 102 *> \verbatim
 103 *>          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
 104 *>         On input, B contains the right hand sides of the least
 105 *>         squares problem. On output, B contains the solution X.
 106 *> \endverbatim
 107 *>
 108 *> \param[in] LDB
 109 *> \verbatim
 110 *>          LDB is INTEGER
 111 *>         The leading dimension of B in the calling subprogram.
 112 *>         LDB must be at least max(1,N).
 113 *> \endverbatim
 114 *>
 115 *> \param[in] RCOND
 116 *> \verbatim
 117 *>          RCOND is DOUBLE PRECISION
 118 *>         The singular values of A less than or equal to RCOND times
 119 *>         the largest singular value are treated as zero in solving
 120 *>         the least squares problem. If RCOND is negative,
 121 *>         machine precision is used instead.
 122 *>         For example, if diag(S)*X=B were the least squares problem,
 123 *>         where diag(S) is a diagonal matrix of singular values, the
 124 *>         solution would be X(i) = B(i) / S(i) if S(i) is greater than
 125 *>         RCOND*max(S), and X(i) = 0 if S(i) is less than or equal to
 126 *>         RCOND*max(S).
 127 *> \endverbatim
 128 *>
 129 *> \param[out] RANK
 130 *> \verbatim
 131 *>          RANK is INTEGER
 132 *>         The number of singular values of A greater than RCOND times
 133 *>         the largest singular value.
 134 *> \endverbatim
 135 *>
 136 *> \param[out] WORK
 137 *> \verbatim
 138 *>          WORK is DOUBLE PRECISION array, dimension at least
 139 *>         (9*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2),
 140 *>         where NLVL = max(0, INT(log_2 (N/(SMLSIZ+1))) + 1).
 141 *> \endverbatim
 142 *>
 143 *> \param[out] IWORK
 144 *> \verbatim
 145 *>          IWORK is INTEGER array, dimension at least
 146 *>         (3*N*NLVL + 11*N)
 147 *> \endverbatim
 148 *>
 149 *> \param[out] INFO
 150 *> \verbatim
 151 *>          INFO is INTEGER
 152 *>         = 0:  successful exit.
 153 *>         < 0:  if INFO = -i, the i-th argument had an illegal value.
 154 *>         > 0:  The algorithm failed to compute a singular value while
 155 *>               working on the submatrix lying in rows and columns
 156 *>               INFO/(N+1) through MOD(INFO,N+1).
 157 *> \endverbatim
 158 *
 159 *  Authors:
 160 *  ========
 161 *
 162 *> \author Univ. of Tennessee
 163 *> \author Univ. of California Berkeley
 164 *> \author Univ. of Colorado Denver
 165 *> \author NAG Ltd.
 166 *
 167 *> \date September 2012
 168 *
 169 *> \ingroup doubleOTHERcomputational
 170 *
 171 *> \par Contributors:
 172 *  ==================
 173 *>
 174 *>     Ming Gu and Ren-Cang Li, Computer Science Division, University of
 175 *>       California at Berkeley, USA \n
 176 *>     Osni Marques, LBNL/NERSC, USA \n
 177 *
 178 *  =====================================================================
 179       SUBROUTINE DLALSD( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
 180      $                   RANK, WORK, IWORK, INFO )
 181 *
 182 *  -- LAPACK computational routine (version 3.4.2) --
 183 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 184 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 185 *     September 2012
 186 *
 187 *     .. Scalar Arguments ..
 188       CHARACTER          UPLO
 189       INTEGER            INFO, LDB, N, NRHS, RANK, SMLSIZ
 190       DOUBLE PRECISION   RCOND
 191 *     ..
 192 *     .. Array Arguments ..
 193       INTEGER            IWORK( * )
 194       DOUBLE PRECISION   B( LDB, * ), D( * ), E( * ), WORK( * )
 195 *     ..
 196 *
 197 *  =====================================================================
 198 *
 199 *     .. Parameters ..
 200       DOUBLE PRECISION   ZERO, ONE, TWO
 201       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
 202 *     ..
 203 *     .. Local Scalars ..
 204       INTEGER            BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM,
 205      $                   GIVPTR, I, ICMPQ1, ICMPQ2, IWK, J, K, NLVL,
 206      $                   NM1, NSIZE, NSUB, NWORK, PERM, POLES, S, SIZEI,
 207      $                   SMLSZP, SQRE, ST, ST1, U, VT, Z
 208       DOUBLE PRECISION   CS, EPS, ORGNRM, R, RCND, SN, TOL
 209 *     ..
 210 *     .. External Functions ..
 211       INTEGER            IDAMAX
 212       DOUBLE PRECISION   DLAMCH, DLANST
 213       EXTERNAL           IDAMAX, DLAMCH, DLANST
 214 *     ..
 215 *     .. External Subroutines ..
 216       EXTERNAL           DCOPY, DGEMM, DLACPY, DLALSA, DLARTG, DLASCL,
 217      $                   DLASDA, DLASDQ, DLASET, DLASRT, DROT, XERBLA
 218 *     ..
 219 *     .. Intrinsic Functions ..
 220       INTRINSIC          ABS, DBLE, INT, LOG, SIGN
 221 *     ..
 222 *     .. Executable Statements ..
 223 *
 224 *     Test the input parameters.
 225 *
 226       INFO = 0
 227 *
 228       IF( N.LT.0 ) THEN
 229          INFO = -3
 230       ELSE IF( NRHS.LT.1 ) THEN
 231          INFO = -4
 232       ELSE IF( ( LDB.LT.1 ) .OR. ( LDB.LT.N ) ) THEN
 233          INFO = -8
 234       END IF
 235       IF( INFO.NE.0 ) THEN
 236          CALL XERBLA( 'DLALSD', -INFO )
 237          RETURN
 238       END IF
 239 *
 240       EPS = DLAMCH( 'Epsilon' )
 241 *
 242 *     Set up the tolerance.
 243 *
 244       IF( ( RCOND.LE.ZERO ) .OR. ( RCOND.GE.ONE ) ) THEN
 245          RCND = EPS
 246       ELSE
 247          RCND = RCOND
 248       END IF
 249 *
 250       RANK = 0
 251 *
 252 *     Quick return if possible.
 253 *
 254       IF( N.EQ.0 ) THEN
 255          RETURN
 256       ELSE IF( N.EQ.1 ) THEN
 257          IF( D( 1 ).EQ.ZERO ) THEN
 258             CALL DLASET( 'A', 1, NRHS, ZERO, ZERO, B, LDB )
 259          ELSE
 260             RANK = 1
 261             CALL DLASCL( 'G', 0, 0, D( 1 ), ONE, 1, NRHS, B, LDB, INFO )
 262             D( 1 ) = ABS( D( 1 ) )
 263          END IF
 264          RETURN
 265       END IF
 266 *
 267 *     Rotate the matrix if it is lower bidiagonal.
 268 *
 269       IF( UPLO.EQ.'L' ) THEN
 270          DO 10 I = 1, N - 1
 271             CALL DLARTG( D( I ), E( I ), CS, SN, R )
 272             D( I ) = R
 273             E( I ) = SN*D( I+1 )
 274             D( I+1 ) = CS*D( I+1 )
 275             IF( NRHS.EQ.1 ) THEN
 276                CALL DROT( 1, B( I, 1 ), 1, B( I+1, 1 ), 1, CS, SN )
 277             ELSE
 278                WORK( I*2-1 ) = CS
 279                WORK( I*2 ) = SN
 280             END IF
 281    10    CONTINUE
 282          IF( NRHS.GT.1 ) THEN
 283             DO 30 I = 1, NRHS
 284                DO 20 J = 1, N - 1
 285                   CS = WORK( J*2-1 )
 286                   SN = WORK( J*2 )
 287                   CALL DROT( 1, B( J, I ), 1, B( J+1, I ), 1, CS, SN )
 288    20          CONTINUE
 289    30       CONTINUE
 290          END IF
 291       END IF
 292 *
 293 *     Scale.
 294 *
 295       NM1 = N - 1
 296       ORGNRM = DLANST( 'M', N, D, E )
 297       IF( ORGNRM.EQ.ZERO ) THEN
 298          CALL DLASET( 'A', N, NRHS, ZERO, ZERO, B, LDB )
 299          RETURN
 300       END IF
 301 *
 302       CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
 303       CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, INFO )
 304 *
 305 *     If N is smaller than the minimum divide size SMLSIZ, then solve
 306 *     the problem with another solver.
 307 *
 308       IF( N.LE.SMLSIZ ) THEN
 309          NWORK = 1 + N*N
 310          CALL DLASET( 'A', N, N, ZERO, ONE, WORK, N )
 311          CALL DLASDQ( 'U', 0, N, N, 0, NRHS, D, E, WORK, N, WORK, N, B,
 312      $                LDB, WORK( NWORK ), INFO )
 313          IF( INFO.NE.0 ) THEN
 314             RETURN
 315          END IF
 316          TOL = RCND*ABS( D( IDAMAX( N, D, 1 ) ) )
 317          DO 40 I = 1, N
 318             IF( D( I ).LE.TOL ) THEN
 319                CALL DLASET( 'A', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB )
 320             ELSE
 321                CALL DLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS, B( I, 1 ),
 322      $                      LDB, INFO )
 323                RANK = RANK + 1
 324             END IF
 325    40    CONTINUE
 326          CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, WORK, N, B, LDB, ZERO,
 327      $               WORK( NWORK ), N )
 328          CALL DLACPY( 'A', N, NRHS, WORK( NWORK ), N, B, LDB )
 329 *
 330 *        Unscale.
 331 *
 332          CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
 333          CALL DLASRT( 'D', N, D, INFO )
 334          CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
 335 *
 336          RETURN
 337       END IF
 338 *
 339 *     Book-keeping and setting up some constants.
 340 *
 341       NLVL = INT( LOG( DBLE( N ) / DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1
 342 *
 343       SMLSZP = SMLSIZ + 1
 344 *
 345       U = 1
 346       VT = 1 + SMLSIZ*N
 347       DIFL = VT + SMLSZP*N
 348       DIFR = DIFL + NLVL*N
 349       Z = DIFR + NLVL*N*2
 350       C = Z + NLVL*N
 351       S = C + N
 352       POLES = S + N
 353       GIVNUM = POLES + 2*NLVL*N
 354       BX = GIVNUM + 2*NLVL*N
 355       NWORK = BX + N*NRHS
 356 *
 357       SIZEI = 1 + N
 358       K = SIZEI + N
 359       GIVPTR = K + N
 360       PERM = GIVPTR + N
 361       GIVCOL = PERM + NLVL*N
 362       IWK = GIVCOL + NLVL*N*2
 363 *
 364       ST = 1
 365       SQRE = 0
 366       ICMPQ1 = 1
 367       ICMPQ2 = 0
 368       NSUB = 0
 369 *
 370       DO 50 I = 1, N
 371          IF( ABS( D( I ) ).LT.EPS ) THEN
 372             D( I ) = SIGN( EPS, D( I ) )
 373          END IF
 374    50 CONTINUE
 375 *
 376       DO 60 I = 1, NM1
 377          IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN
 378             NSUB = NSUB + 1
 379             IWORK( NSUB ) = ST
 380 *
 381 *           Subproblem found. First determine its size and then
 382 *           apply divide and conquer on it.
 383 *
 384             IF( I.LT.NM1 ) THEN
 385 *
 386 *              A subproblem with E(I) small for I < NM1.
 387 *
 388                NSIZE = I - ST + 1
 389                IWORK( SIZEI+NSUB-1 ) = NSIZE
 390             ELSE IF( ABS( E( I ) ).GE.EPS ) THEN
 391 *
 392 *              A subproblem with E(NM1) not too small but I = NM1.
 393 *
 394                NSIZE = N - ST + 1
 395                IWORK( SIZEI+NSUB-1 ) = NSIZE
 396             ELSE
 397 *
 398 *              A subproblem with E(NM1) small. This implies an
 399 *              1-by-1 subproblem at D(N), which is not solved
 400 *              explicitly.
 401 *
 402                NSIZE = I - ST + 1
 403                IWORK( SIZEI+NSUB-1 ) = NSIZE
 404                NSUB = NSUB + 1
 405                IWORK( NSUB ) = N
 406                IWORK( SIZEI+NSUB-1 ) = 1
 407                CALL DCOPY( NRHS, B( N, 1 ), LDB, WORK( BX+NM1 ), N )
 408             END IF
 409             ST1 = ST - 1
 410             IF( NSIZE.EQ.1 ) THEN
 411 *
 412 *              This is a 1-by-1 subproblem and is not solved
 413 *              explicitly.
 414 *
 415                CALL DCOPY( NRHS, B( ST, 1 ), LDB, WORK( BX+ST1 ), N )
 416             ELSE IF( NSIZE.LE.SMLSIZ ) THEN
 417 *
 418 *              This is a small subproblem and is solved by DLASDQ.
 419 *
 420                CALL DLASET( 'A', NSIZE, NSIZE, ZERO, ONE,
 421      $                      WORK( VT+ST1 ), N )
 422                CALL DLASDQ( 'U', 0, NSIZE, NSIZE, 0, NRHS, D( ST ),
 423      $                      E( ST ), WORK( VT+ST1 ), N, WORK( NWORK ),
 424      $                      N, B( ST, 1 ), LDB, WORK( NWORK ), INFO )
 425                IF( INFO.NE.0 ) THEN
 426                   RETURN
 427                END IF
 428                CALL DLACPY( 'A', NSIZE, NRHS, B( ST, 1 ), LDB,
 429      $                      WORK( BX+ST1 ), N )
 430             ELSE
 431 *
 432 *              A large problem. Solve it using divide and conquer.
 433 *
 434                CALL DLASDA( ICMPQ1, SMLSIZ, NSIZE, SQRE, D( ST ),
 435      $                      E( ST ), WORK( U+ST1 ), N, WORK( VT+ST1 ),
 436      $                      IWORK( K+ST1 ), WORK( DIFL+ST1 ),
 437      $                      WORK( DIFR+ST1 ), WORK( Z+ST1 ),
 438      $                      WORK( POLES+ST1 ), IWORK( GIVPTR+ST1 ),
 439      $                      IWORK( GIVCOL+ST1 ), N, IWORK( PERM+ST1 ),
 440      $                      WORK( GIVNUM+ST1 ), WORK( C+ST1 ),
 441      $                      WORK( S+ST1 ), WORK( NWORK ), IWORK( IWK ),
 442      $                      INFO )
 443                IF( INFO.NE.0 ) THEN
 444                   RETURN
 445                END IF
 446                BXST = BX + ST1
 447                CALL DLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, B( ST, 1 ),
 448      $                      LDB, WORK( BXST ), N, WORK( U+ST1 ), N,
 449      $                      WORK( VT+ST1 ), IWORK( K+ST1 ),
 450      $                      WORK( DIFL+ST1 ), WORK( DIFR+ST1 ),
 451      $                      WORK( Z+ST1 ), WORK( POLES+ST1 ),
 452      $                      IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
 453      $                      IWORK( PERM+ST1 ), WORK( GIVNUM+ST1 ),
 454      $                      WORK( C+ST1 ), WORK( S+ST1 ), WORK( NWORK ),
 455      $                      IWORK( IWK ), INFO )
 456                IF( INFO.NE.0 ) THEN
 457                   RETURN
 458                END IF
 459             END IF
 460             ST = I + 1
 461          END IF
 462    60 CONTINUE
 463 *
 464 *     Apply the singular values and treat the tiny ones as zero.
 465 *
 466       TOL = RCND*ABS( D( IDAMAX( N, D, 1 ) ) )
 467 *
 468       DO 70 I = 1, N
 469 *
 470 *        Some of the elements in D can be negative because 1-by-1
 471 *        subproblems were not solved explicitly.
 472 *
 473          IF( ABS( D( I ) ).LE.TOL ) THEN
 474             CALL DLASET( 'A', 1, NRHS, ZERO, ZERO, WORK( BX+I-1 ), N )
 475          ELSE
 476             RANK = RANK + 1
 477             CALL DLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS,
 478      $                   WORK( BX+I-1 ), N, INFO )
 479          END IF
 480          D( I ) = ABS( D( I ) )
 481    70 CONTINUE
 482 *
 483 *     Now apply back the right singular vectors.
 484 *
 485       ICMPQ2 = 1
 486       DO 80 I = 1, NSUB
 487          ST = IWORK( I )
 488          ST1 = ST - 1
 489          NSIZE = IWORK( SIZEI+I-1 )
 490          BXST = BX + ST1
 491          IF( NSIZE.EQ.1 ) THEN
 492             CALL DCOPY( NRHS, WORK( BXST ), N, B( ST, 1 ), LDB )
 493          ELSE IF( NSIZE.LE.SMLSIZ ) THEN
 494             CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
 495      $                  WORK( VT+ST1 ), N, WORK( BXST ), N, ZERO,
 496      $                  B( ST, 1 ), LDB )
 497          ELSE
 498             CALL DLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, WORK( BXST ), N,
 499      $                   B( ST, 1 ), LDB, WORK( U+ST1 ), N,
 500      $                   WORK( VT+ST1 ), IWORK( K+ST1 ),
 501      $                   WORK( DIFL+ST1 ), WORK( DIFR+ST1 ),
 502      $                   WORK( Z+ST1 ), WORK( POLES+ST1 ),
 503      $                   IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
 504      $                   IWORK( PERM+ST1 ), WORK( GIVNUM+ST1 ),
 505      $                   WORK( C+ST1 ), WORK( S+ST1 ), WORK( NWORK ),
 506      $                   IWORK( IWK ), INFO )
 507             IF( INFO.NE.0 ) THEN
 508                RETURN
 509             END IF
 510          END IF
 511    80 CONTINUE
 512 *
 513 *     Unscale and sort the singular values.
 514 *
 515       CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
 516       CALL DLASRT( 'D', N, D, INFO )
 517       CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
 518 *
 519       RETURN
 520 *
 521 *     End of DLALSD
 522 *
 523       END