LAPACK: SRC/dlalsd.f Source File

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