LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dlalsd.f
Go to the documentation of this file.
00001 *> \brief \b DLALSD
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DLALSD + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlalsd.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlalsd.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlalsd.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DLALSD( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
00022 *                          RANK, WORK, IWORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          UPLO
00026 *       INTEGER            INFO, LDB, N, NRHS, RANK, SMLSIZ
00027 *       DOUBLE PRECISION   RCOND
00028 *       ..
00029 *       .. Array Arguments ..
00030 *       INTEGER            IWORK( * )
00031 *       DOUBLE PRECISION   B( LDB, * ), D( * ), E( * ), WORK( * )
00032 *       ..
00033 *  
00034 *
00035 *> \par Purpose:
00036 *  =============
00037 *>
00038 *> \verbatim
00039 *>
00040 *> DLALSD uses the singular value decomposition of A to solve the least
00041 *> squares problem of finding X to minimize the Euclidean norm of each
00042 *> column of A*X-B, where A is N-by-N upper bidiagonal, and X and B
00043 *> are N-by-NRHS. The solution X overwrites B.
00044 *>
00045 *> The singular values of A smaller than RCOND times the largest
00046 *> singular value are treated as zero in solving the least squares
00047 *> problem; in this case a minimum norm solution is returned.
00048 *> The actual singular values are returned in D in ascending order.
00049 *>
00050 *> This code makes very mild assumptions about floating point
00051 *> arithmetic. It will work on machines with a guard digit in
00052 *> add/subtract, or on those binary machines without guard digits
00053 *> which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
00054 *> It could conceivably fail on hexadecimal or decimal machines
00055 *> without guard digits, but we know of none.
00056 *> \endverbatim
00057 *
00058 *  Arguments:
00059 *  ==========
00060 *
00061 *> \param[in] UPLO
00062 *> \verbatim
00063 *>          UPLO is CHARACTER*1
00064 *>         = 'U': D and E define an upper bidiagonal matrix.
00065 *>         = 'L': D and E define a  lower bidiagonal matrix.
00066 *> \endverbatim
00067 *>
00068 *> \param[in] SMLSIZ
00069 *> \verbatim
00070 *>          SMLSIZ is INTEGER
00071 *>         The maximum size of the subproblems at the bottom of the
00072 *>         computation tree.
00073 *> \endverbatim
00074 *>
00075 *> \param[in] N
00076 *> \verbatim
00077 *>          N is INTEGER
00078 *>         The dimension of the  bidiagonal matrix.  N >= 0.
00079 *> \endverbatim
00080 *>
00081 *> \param[in] NRHS
00082 *> \verbatim
00083 *>          NRHS is INTEGER
00084 *>         The number of columns of B. NRHS must be at least 1.
00085 *> \endverbatim
00086 *>
00087 *> \param[in,out] D
00088 *> \verbatim
00089 *>          D is DOUBLE PRECISION array, dimension (N)
00090 *>         On entry D contains the main diagonal of the bidiagonal
00091 *>         matrix. On exit, if INFO = 0, D contains its singular values.
00092 *> \endverbatim
00093 *>
00094 *> \param[in,out] E
00095 *> \verbatim
00096 *>          E is DOUBLE PRECISION array, dimension (N-1)
00097 *>         Contains the super-diagonal entries of the bidiagonal matrix.
00098 *>         On exit, E has been destroyed.
00099 *> \endverbatim
00100 *>
00101 *> \param[in,out] B
00102 *> \verbatim
00103 *>          B is DOUBLE PRECISION array, dimension (LDB,NRHS)
00104 *>         On input, B contains the right hand sides of the least
00105 *>         squares problem. On output, B contains the solution X.
00106 *> \endverbatim
00107 *>
00108 *> \param[in] LDB
00109 *> \verbatim
00110 *>          LDB is INTEGER
00111 *>         The leading dimension of B in the calling subprogram.
00112 *>         LDB must be at least max(1,N).
00113 *> \endverbatim
00114 *>
00115 *> \param[in] RCOND
00116 *> \verbatim
00117 *>          RCOND is DOUBLE PRECISION
00118 *>         The singular values of A less than or equal to RCOND times
00119 *>         the largest singular value are treated as zero in solving
00120 *>         the least squares problem. If RCOND is negative,
00121 *>         machine precision is used instead.
00122 *>         For example, if diag(S)*X=B were the least squares problem,
00123 *>         where diag(S) is a diagonal matrix of singular values, the
00124 *>         solution would be X(i) = B(i) / S(i) if S(i) is greater than
00125 *>         RCOND*max(S), and X(i) = 0 if S(i) is less than or equal to
00126 *>         RCOND*max(S).
00127 *> \endverbatim
00128 *>
00129 *> \param[out] RANK
00130 *> \verbatim
00131 *>          RANK is INTEGER
00132 *>         The number of singular values of A greater than RCOND times
00133 *>         the largest singular value.
00134 *> \endverbatim
00135 *>
00136 *> \param[out] WORK
00137 *> \verbatim
00138 *>          WORK is DOUBLE PRECISION array, dimension at least
00139 *>         (9*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2),
00140 *>         where NLVL = max(0, INT(log_2 (N/(SMLSIZ+1))) + 1).
00141 *> \endverbatim
00142 *>
00143 *> \param[out] IWORK
00144 *> \verbatim
00145 *>          IWORK is INTEGER array, dimension at least
00146 *>         (3*N*NLVL + 11*N)
00147 *> \endverbatim
00148 *>
00149 *> \param[out] INFO
00150 *> \verbatim
00151 *>          INFO is INTEGER
00152 *>         = 0:  successful exit.
00153 *>         < 0:  if INFO = -i, the i-th argument had an illegal value.
00154 *>         > 0:  The algorithm failed to compute a singular value while
00155 *>               working on the submatrix lying in rows and columns
00156 *>               INFO/(N+1) through MOD(INFO,N+1).
00157 *> \endverbatim
00158 *
00159 *  Authors:
00160 *  ========
00161 *
00162 *> \author Univ. of Tennessee 
00163 *> \author Univ. of California Berkeley 
00164 *> \author Univ. of Colorado Denver 
00165 *> \author NAG Ltd. 
00166 *
00167 *> \date November 2011
00168 *
00169 *> \ingroup doubleOTHERcomputational
00170 *
00171 *> \par Contributors:
00172 *  ==================
00173 *>
00174 *>     Ming Gu and Ren-Cang Li, Computer Science Division, University of
00175 *>       California at Berkeley, USA \n
00176 *>     Osni Marques, LBNL/NERSC, USA \n
00177 *
00178 *  =====================================================================
00179       SUBROUTINE DLALSD( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
00180      $                   RANK, WORK, IWORK, INFO )
00181 *
00182 *  -- LAPACK computational routine (version 3.4.0) --
00183 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00184 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00185 *     November 2011
00186 *
00187 *     .. Scalar Arguments ..
00188       CHARACTER          UPLO
00189       INTEGER            INFO, LDB, N, NRHS, RANK, SMLSIZ
00190       DOUBLE PRECISION   RCOND
00191 *     ..
00192 *     .. Array Arguments ..
00193       INTEGER            IWORK( * )
00194       DOUBLE PRECISION   B( LDB, * ), D( * ), E( * ), WORK( * )
00195 *     ..
00196 *
00197 *  =====================================================================
00198 *
00199 *     .. Parameters ..
00200       DOUBLE PRECISION   ZERO, ONE, TWO
00201       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
00202 *     ..
00203 *     .. Local Scalars ..
00204       INTEGER            BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM,
00205      $                   GIVPTR, I, ICMPQ1, ICMPQ2, IWK, J, K, NLVL,
00206      $                   NM1, NSIZE, NSUB, NWORK, PERM, POLES, S, SIZEI,
00207      $                   SMLSZP, SQRE, ST, ST1, U, VT, Z
00208       DOUBLE PRECISION   CS, EPS, ORGNRM, R, RCND, SN, TOL
00209 *     ..
00210 *     .. External Functions ..
00211       INTEGER            IDAMAX
00212       DOUBLE PRECISION   DLAMCH, DLANST
00213       EXTERNAL           IDAMAX, DLAMCH, DLANST
00214 *     ..
00215 *     .. External Subroutines ..
00216       EXTERNAL           DCOPY, DGEMM, DLACPY, DLALSA, DLARTG, DLASCL,
00217      $                   DLASDA, DLASDQ, DLASET, DLASRT, DROT, XERBLA
00218 *     ..
00219 *     .. Intrinsic Functions ..
00220       INTRINSIC          ABS, DBLE, INT, LOG, SIGN
00221 *     ..
00222 *     .. Executable Statements ..
00223 *
00224 *     Test the input parameters.
00225 *
00226       INFO = 0
00227 *
00228       IF( N.LT.0 ) THEN
00229          INFO = -3
00230       ELSE IF( NRHS.LT.1 ) THEN
00231          INFO = -4
00232       ELSE IF( ( LDB.LT.1 ) .OR. ( LDB.LT.N ) ) THEN
00233          INFO = -8
00234       END IF
00235       IF( INFO.NE.0 ) THEN
00236          CALL XERBLA( 'DLALSD', -INFO )
00237          RETURN
00238       END IF
00239 *
00240       EPS = DLAMCH( 'Epsilon' )
00241 *
00242 *     Set up the tolerance.
00243 *
00244       IF( ( RCOND.LE.ZERO ) .OR. ( RCOND.GE.ONE ) ) THEN
00245          RCND = EPS
00246       ELSE
00247          RCND = RCOND
00248       END IF
00249 *
00250       RANK = 0
00251 *
00252 *     Quick return if possible.
00253 *
00254       IF( N.EQ.0 ) THEN
00255          RETURN
00256       ELSE IF( N.EQ.1 ) THEN
00257          IF( D( 1 ).EQ.ZERO ) THEN
00258             CALL DLASET( 'A', 1, NRHS, ZERO, ZERO, B, LDB )
00259          ELSE
00260             RANK = 1
00261             CALL DLASCL( 'G', 0, 0, D( 1 ), ONE, 1, NRHS, B, LDB, INFO )
00262             D( 1 ) = ABS( D( 1 ) )
00263          END IF
00264          RETURN
00265       END IF
00266 *
00267 *     Rotate the matrix if it is lower bidiagonal.
00268 *
00269       IF( UPLO.EQ.'L' ) THEN
00270          DO 10 I = 1, N - 1
00271             CALL DLARTG( D( I ), E( I ), CS, SN, R )
00272             D( I ) = R
00273             E( I ) = SN*D( I+1 )
00274             D( I+1 ) = CS*D( I+1 )
00275             IF( NRHS.EQ.1 ) THEN
00276                CALL DROT( 1, B( I, 1 ), 1, B( I+1, 1 ), 1, CS, SN )
00277             ELSE
00278                WORK( I*2-1 ) = CS
00279                WORK( I*2 ) = SN
00280             END IF
00281    10    CONTINUE
00282          IF( NRHS.GT.1 ) THEN
00283             DO 30 I = 1, NRHS
00284                DO 20 J = 1, N - 1
00285                   CS = WORK( J*2-1 )
00286                   SN = WORK( J*2 )
00287                   CALL DROT( 1, B( J, I ), 1, B( J+1, I ), 1, CS, SN )
00288    20          CONTINUE
00289    30       CONTINUE
00290          END IF
00291       END IF
00292 *
00293 *     Scale.
00294 *
00295       NM1 = N - 1
00296       ORGNRM = DLANST( 'M', N, D, E )
00297       IF( ORGNRM.EQ.ZERO ) THEN
00298          CALL DLASET( 'A', N, NRHS, ZERO, ZERO, B, LDB )
00299          RETURN
00300       END IF
00301 *
00302       CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
00303       CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, INFO )
00304 *
00305 *     If N is smaller than the minimum divide size SMLSIZ, then solve
00306 *     the problem with another solver.
00307 *
00308       IF( N.LE.SMLSIZ ) THEN
00309          NWORK = 1 + N*N
00310          CALL DLASET( 'A', N, N, ZERO, ONE, WORK, N )
00311          CALL DLASDQ( 'U', 0, N, N, 0, NRHS, D, E, WORK, N, WORK, N, B,
00312      $                LDB, WORK( NWORK ), INFO )
00313          IF( INFO.NE.0 ) THEN
00314             RETURN
00315          END IF
00316          TOL = RCND*ABS( D( IDAMAX( N, D, 1 ) ) )
00317          DO 40 I = 1, N
00318             IF( D( I ).LE.TOL ) THEN
00319                CALL DLASET( 'A', 1, NRHS, ZERO, ZERO, B( I, 1 ), LDB )
00320             ELSE
00321                CALL DLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS, B( I, 1 ),
00322      $                      LDB, INFO )
00323                RANK = RANK + 1
00324             END IF
00325    40    CONTINUE
00326          CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, WORK, N, B, LDB, ZERO,
00327      $               WORK( NWORK ), N )
00328          CALL DLACPY( 'A', N, NRHS, WORK( NWORK ), N, B, LDB )
00329 *
00330 *        Unscale.
00331 *
00332          CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
00333          CALL DLASRT( 'D', N, D, INFO )
00334          CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
00335 *
00336          RETURN
00337       END IF
00338 *
00339 *     Book-keeping and setting up some constants.
00340 *
00341       NLVL = INT( LOG( DBLE( N ) / DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1
00342 *
00343       SMLSZP = SMLSIZ + 1
00344 *
00345       U = 1
00346       VT = 1 + SMLSIZ*N
00347       DIFL = VT + SMLSZP*N
00348       DIFR = DIFL + NLVL*N
00349       Z = DIFR + NLVL*N*2
00350       C = Z + NLVL*N
00351       S = C + N
00352       POLES = S + N
00353       GIVNUM = POLES + 2*NLVL*N
00354       BX = GIVNUM + 2*NLVL*N
00355       NWORK = BX + N*NRHS
00356 *
00357       SIZEI = 1 + N
00358       K = SIZEI + N
00359       GIVPTR = K + N
00360       PERM = GIVPTR + N
00361       GIVCOL = PERM + NLVL*N
00362       IWK = GIVCOL + NLVL*N*2
00363 *
00364       ST = 1
00365       SQRE = 0
00366       ICMPQ1 = 1
00367       ICMPQ2 = 0
00368       NSUB = 0
00369 *
00370       DO 50 I = 1, N
00371          IF( ABS( D( I ) ).LT.EPS ) THEN
00372             D( I ) = SIGN( EPS, D( I ) )
00373          END IF
00374    50 CONTINUE
00375 *
00376       DO 60 I = 1, NM1
00377          IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN
00378             NSUB = NSUB + 1
00379             IWORK( NSUB ) = ST
00380 *
00381 *           Subproblem found. First determine its size and then
00382 *           apply divide and conquer on it.
00383 *
00384             IF( I.LT.NM1 ) THEN
00385 *
00386 *              A subproblem with E(I) small for I < NM1.
00387 *
00388                NSIZE = I - ST + 1
00389                IWORK( SIZEI+NSUB-1 ) = NSIZE
00390             ELSE IF( ABS( E( I ) ).GE.EPS ) THEN
00391 *
00392 *              A subproblem with E(NM1) not too small but I = NM1.
00393 *
00394                NSIZE = N - ST + 1
00395                IWORK( SIZEI+NSUB-1 ) = NSIZE
00396             ELSE
00397 *
00398 *              A subproblem with E(NM1) small. This implies an
00399 *              1-by-1 subproblem at D(N), which is not solved
00400 *              explicitly.
00401 *
00402                NSIZE = I - ST + 1
00403                IWORK( SIZEI+NSUB-1 ) = NSIZE
00404                NSUB = NSUB + 1
00405                IWORK( NSUB ) = N
00406                IWORK( SIZEI+NSUB-1 ) = 1
00407                CALL DCOPY( NRHS, B( N, 1 ), LDB, WORK( BX+NM1 ), N )
00408             END IF
00409             ST1 = ST - 1
00410             IF( NSIZE.EQ.1 ) THEN
00411 *
00412 *              This is a 1-by-1 subproblem and is not solved
00413 *              explicitly.
00414 *
00415                CALL DCOPY( NRHS, B( ST, 1 ), LDB, WORK( BX+ST1 ), N )
00416             ELSE IF( NSIZE.LE.SMLSIZ ) THEN
00417 *
00418 *              This is a small subproblem and is solved by DLASDQ.
00419 *
00420                CALL DLASET( 'A', NSIZE, NSIZE, ZERO, ONE,
00421      $                      WORK( VT+ST1 ), N )
00422                CALL DLASDQ( 'U', 0, NSIZE, NSIZE, 0, NRHS, D( ST ),
00423      $                      E( ST ), WORK( VT+ST1 ), N, WORK( NWORK ),
00424      $                      N, B( ST, 1 ), LDB, WORK( NWORK ), INFO )
00425                IF( INFO.NE.0 ) THEN
00426                   RETURN
00427                END IF
00428                CALL DLACPY( 'A', NSIZE, NRHS, B( ST, 1 ), LDB,
00429      $                      WORK( BX+ST1 ), N )
00430             ELSE
00431 *
00432 *              A large problem. Solve it using divide and conquer.
00433 *
00434                CALL DLASDA( ICMPQ1, SMLSIZ, NSIZE, SQRE, D( ST ),
00435      $                      E( ST ), WORK( U+ST1 ), N, WORK( VT+ST1 ),
00436      $                      IWORK( K+ST1 ), WORK( DIFL+ST1 ),
00437      $                      WORK( DIFR+ST1 ), WORK( Z+ST1 ),
00438      $                      WORK( POLES+ST1 ), IWORK( GIVPTR+ST1 ),
00439      $                      IWORK( GIVCOL+ST1 ), N, IWORK( PERM+ST1 ),
00440      $                      WORK( GIVNUM+ST1 ), WORK( C+ST1 ),
00441      $                      WORK( S+ST1 ), WORK( NWORK ), IWORK( IWK ),
00442      $                      INFO )
00443                IF( INFO.NE.0 ) THEN
00444                   RETURN
00445                END IF
00446                BXST = BX + ST1
00447                CALL DLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, B( ST, 1 ),
00448      $                      LDB, WORK( BXST ), N, WORK( U+ST1 ), N,
00449      $                      WORK( VT+ST1 ), IWORK( K+ST1 ),
00450      $                      WORK( DIFL+ST1 ), WORK( DIFR+ST1 ),
00451      $                      WORK( Z+ST1 ), WORK( POLES+ST1 ),
00452      $                      IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
00453      $                      IWORK( PERM+ST1 ), WORK( GIVNUM+ST1 ),
00454      $                      WORK( C+ST1 ), WORK( S+ST1 ), WORK( NWORK ),
00455      $                      IWORK( IWK ), INFO )
00456                IF( INFO.NE.0 ) THEN
00457                   RETURN
00458                END IF
00459             END IF
00460             ST = I + 1
00461          END IF
00462    60 CONTINUE
00463 *
00464 *     Apply the singular values and treat the tiny ones as zero.
00465 *
00466       TOL = RCND*ABS( D( IDAMAX( N, D, 1 ) ) )
00467 *
00468       DO 70 I = 1, N
00469 *
00470 *        Some of the elements in D can be negative because 1-by-1
00471 *        subproblems were not solved explicitly.
00472 *
00473          IF( ABS( D( I ) ).LE.TOL ) THEN
00474             CALL DLASET( 'A', 1, NRHS, ZERO, ZERO, WORK( BX+I-1 ), N )
00475          ELSE
00476             RANK = RANK + 1
00477             CALL DLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS,
00478      $                   WORK( BX+I-1 ), N, INFO )
00479          END IF
00480          D( I ) = ABS( D( I ) )
00481    70 CONTINUE
00482 *
00483 *     Now apply back the right singular vectors.
00484 *
00485       ICMPQ2 = 1
00486       DO 80 I = 1, NSUB
00487          ST = IWORK( I )
00488          ST1 = ST - 1
00489          NSIZE = IWORK( SIZEI+I-1 )
00490          BXST = BX + ST1
00491          IF( NSIZE.EQ.1 ) THEN
00492             CALL DCOPY( NRHS, WORK( BXST ), N, B( ST, 1 ), LDB )
00493          ELSE IF( NSIZE.LE.SMLSIZ ) THEN
00494             CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
00495      $                  WORK( VT+ST1 ), N, WORK( BXST ), N, ZERO,
00496      $                  B( ST, 1 ), LDB )
00497          ELSE
00498             CALL DLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, WORK( BXST ), N,
00499      $                   B( ST, 1 ), LDB, WORK( U+ST1 ), N,
00500      $                   WORK( VT+ST1 ), IWORK( K+ST1 ),
00501      $                   WORK( DIFL+ST1 ), WORK( DIFR+ST1 ),
00502      $                   WORK( Z+ST1 ), WORK( POLES+ST1 ),
00503      $                   IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
00504      $                   IWORK( PERM+ST1 ), WORK( GIVNUM+ST1 ),
00505      $                   WORK( C+ST1 ), WORK( S+ST1 ), WORK( NWORK ),
00506      $                   IWORK( IWK ), INFO )
00507             IF( INFO.NE.0 ) THEN
00508                RETURN
00509             END IF
00510          END IF
00511    80 CONTINUE
00512 *
00513 *     Unscale and sort the singular values.
00514 *
00515       CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
00516       CALL DLASRT( 'D', N, D, INFO )
00517       CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
00518 *
00519       RETURN
00520 *
00521 *     End of DLALSD
00522 *
00523       END
 All Files Functions