LAPACK: SRC/clalsd.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CLALSD( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
00002      $                   RANK, WORK, RWORK, 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       REAL               RCOND
00013 *     ..
00014 *     .. Array Arguments ..
00015       INTEGER            IWORK( * )
00016       REAL               D( * ), E( * ), RWORK( * )
00017       COMPLEX            B( LDB, * ), WORK( * )
00018 *     ..
00019 *
00020 *  Purpose
00021 *  =======
00022 *
00023 *  CLALSD uses the singular value decomposition of A to solve the least
00024 *  squares problem of finding X to minimize the Euclidean norm of each
00025 *  column of A*X-B, where A is N-by-N upper bidiagonal, and X and B
00026 *  are N-by-NRHS. The solution X overwrites B.
00027 *
00028 *  The singular values of A smaller than RCOND times the largest
00029 *  singular value are treated as zero in solving the least squares
00030 *  problem; in this case a minimum norm solution is returned.
00031 *  The actual singular values are returned in D in ascending order.
00032 *
00033 *  This code makes very mild assumptions about floating point
00034 *  arithmetic. It will work on machines with a guard digit in
00035 *  add/subtract, or on those binary machines without guard digits
00036 *  which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
00037 *  It could conceivably fail on hexadecimal or decimal machines
00038 *  without guard digits, but we know of none.
00039 *
00040 *  Arguments
00041 *  =========
00042 *
00043 *  UPLO   (input) CHARACTER*1
00044 *         = 'U': D and E define an upper bidiagonal matrix.
00045 *         = 'L': D and E define a  lower bidiagonal matrix.
00046 *
00047 *  SMLSIZ (input) INTEGER
00048 *         The maximum size of the subproblems at the bottom of the
00049 *         computation tree.
00050 *
00051 *  N      (input) INTEGER
00052 *         The dimension of the  bidiagonal matrix.  N >= 0.
00053 *
00054 *  NRHS   (input) INTEGER
00055 *         The number of columns of B. NRHS must be at least 1.
00056 *
00057 *  D      (input/output) REAL array, dimension (N)
00058 *         On entry D contains the main diagonal of the bidiagonal
00059 *         matrix. On exit, if INFO = 0, D contains its singular values.
00060 *
00061 *  E      (input/output) REAL array, dimension (N-1)
00062 *         Contains the super-diagonal entries of the bidiagonal matrix.
00063 *         On exit, E has been destroyed.
00064 *
00065 *  B      (input/output) COMPLEX array, dimension (LDB,NRHS)
00066 *         On input, B contains the right hand sides of the least
00067 *         squares problem. On output, B contains the solution X.
00068 *
00069 *  LDB    (input) INTEGER
00070 *         The leading dimension of B in the calling subprogram.
00071 *         LDB must be at least max(1,N).
00072 *
00073 *  RCOND  (input) REAL
00074 *         The singular values of A less than or equal to RCOND times
00075 *         the largest singular value are treated as zero in solving
00076 *         the least squares problem. If RCOND is negative,
00077 *         machine precision is used instead.
00078 *         For example, if diag(S)*X=B were the least squares problem,
00079 *         where diag(S) is a diagonal matrix of singular values, the
00080 *         solution would be X(i) = B(i) / S(i) if S(i) is greater than
00081 *         RCOND*max(S), and X(i) = 0 if S(i) is less than or equal to
00082 *         RCOND*max(S).
00083 *
00084 *  RANK   (output) INTEGER
00085 *         The number of singular values of A greater than RCOND times
00086 *         the largest singular value.
00087 *
00088 *  WORK   (workspace) COMPLEX array, dimension (N * NRHS).
00089 *
00090 *  RWORK  (workspace) REAL array, dimension at least
00091 *         (9*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS +
00092 *         MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS ),
00093 *         where
00094 *         NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )
00095 *
00096 *  IWORK  (workspace) INTEGER array, dimension (3*N*NLVL + 11*N).
00097 *
00098 *  INFO   (output) INTEGER
00099 *         = 0:  successful exit.
00100 *         < 0:  if INFO = -i, the i-th argument had an illegal value.
00101 *         > 0:  The algorithm failed to compute a singular value while
00102 *               working on the submatrix lying in rows and columns
00103 *               INFO/(N+1) through MOD(INFO,N+1).
00104 *
00105 *  Further Details
00106 *  ===============
00107 *
00108 *  Based on contributions by
00109 *     Ming Gu and Ren-Cang Li, Computer Science Division, University of
00110 *       California at Berkeley, USA
00111 *     Osni Marques, LBNL/NERSC, USA
00112 *
00113 *  =====================================================================
00114 *
00115 *     .. Parameters ..
00116       REAL               ZERO, ONE, TWO
00117       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0 )
00118       COMPLEX            CZERO
00119       PARAMETER          ( CZERO = ( 0.0E0, 0.0E0 ) )
00120 *     ..
00121 *     .. Local Scalars ..
00122       INTEGER            BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM,
00123      $                   GIVPTR, I, ICMPQ1, ICMPQ2, IRWB, IRWIB, IRWRB,
00124      $                   IRWU, IRWVT, IRWWRK, IWK, J, JCOL, JIMAG,
00125      $                   JREAL, JROW, K, NLVL, NM1, NRWORK, NSIZE, NSUB,
00126      $                   PERM, POLES, S, SIZEI, SMLSZP, SQRE, ST, ST1,
00127      $                   U, VT, Z
00128       REAL               CS, EPS, ORGNRM, R, RCND, SN, TOL
00129 *     ..
00130 *     .. External Functions ..
00131       INTEGER            ISAMAX
00132       REAL               SLAMCH, SLANST
00133       EXTERNAL           ISAMAX, SLAMCH, SLANST
00134 *     ..
00135 *     .. External Subroutines ..
00136       EXTERNAL           CCOPY, CLACPY, CLALSA, CLASCL, CLASET, CSROT,
00137      $                   SGEMM, SLARTG, SLASCL, SLASDA, SLASDQ, SLASET,
00138      $                   SLASRT, XERBLA
00139 *     ..
00140 *     .. Intrinsic Functions ..
00141       INTRINSIC          ABS, AIMAG, CMPLX, INT, LOG, REAL, SIGN
00142 *     ..
00143 *     .. Executable Statements ..
00144 *
00145 *     Test the input parameters.
00146 *
00147       INFO = 0
00148 *
00149       IF( N.LT.0 ) THEN
00150          INFO = -3
00151       ELSE IF( NRHS.LT.1 ) THEN
00152          INFO = -4
00153       ELSE IF( ( LDB.LT.1 ) .OR. ( LDB.LT.N ) ) THEN
00154          INFO = -8
00155       END IF
00156       IF( INFO.NE.0 ) THEN
00157          CALL XERBLA( 'CLALSD', -INFO )
00158          RETURN
00159       END IF
00160 *
00161       EPS = SLAMCH( 'Epsilon' )
00162 *
00163 *     Set up the tolerance.
00164 *
00165       IF( ( RCOND.LE.ZERO ) .OR. ( RCOND.GE.ONE ) ) THEN
00166          RCND = EPS
00167       ELSE
00168          RCND = RCOND
00169       END IF
00170 *
00171       RANK = 0
00172 *
00173 *     Quick return if possible.
00174 *
00175       IF( N.EQ.0 ) THEN
00176          RETURN
00177       ELSE IF( N.EQ.1 ) THEN
00178          IF( D( 1 ).EQ.ZERO ) THEN
00179             CALL CLASET( 'A', 1, NRHS, CZERO, CZERO, B, LDB )
00180          ELSE
00181             RANK = 1
00182             CALL CLASCL( 'G', 0, 0, D( 1 ), ONE, 1, NRHS, B, LDB, INFO )
00183             D( 1 ) = ABS( D( 1 ) )
00184          END IF
00185          RETURN
00186       END IF
00187 *
00188 *     Rotate the matrix if it is lower bidiagonal.
00189 *
00190       IF( UPLO.EQ.'L' ) THEN
00191          DO 10 I = 1, N - 1
00192             CALL SLARTG( D( I ), E( I ), CS, SN, R )
00193             D( I ) = R
00194             E( I ) = SN*D( I+1 )
00195             D( I+1 ) = CS*D( I+1 )
00196             IF( NRHS.EQ.1 ) THEN
00197                CALL CSROT( 1, B( I, 1 ), 1, B( I+1, 1 ), 1, CS, SN )
00198             ELSE
00199                RWORK( I*2-1 ) = CS
00200                RWORK( I*2 ) = SN
00201             END IF
00202    10    CONTINUE
00203          IF( NRHS.GT.1 ) THEN
00204             DO 30 I = 1, NRHS
00205                DO 20 J = 1, N - 1
00206                   CS = RWORK( J*2-1 )
00207                   SN = RWORK( J*2 )
00208                   CALL CSROT( 1, B( J, I ), 1, B( J+1, I ), 1, CS, SN )
00209    20          CONTINUE
00210    30       CONTINUE
00211          END IF
00212       END IF
00213 *
00214 *     Scale.
00215 *
00216       NM1 = N - 1
00217       ORGNRM = SLANST( 'M', N, D, E )
00218       IF( ORGNRM.EQ.ZERO ) THEN
00219          CALL CLASET( 'A', N, NRHS, CZERO, CZERO, B, LDB )
00220          RETURN
00221       END IF
00222 *
00223       CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
00224       CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, INFO )
00225 *
00226 *     If N is smaller than the minimum divide size SMLSIZ, then solve
00227 *     the problem with another solver.
00228 *
00229       IF( N.LE.SMLSIZ ) THEN
00230          IRWU = 1
00231          IRWVT = IRWU + N*N
00232          IRWWRK = IRWVT + N*N
00233          IRWRB = IRWWRK
00234          IRWIB = IRWRB + N*NRHS
00235          IRWB = IRWIB + N*NRHS
00236          CALL SLASET( 'A', N, N, ZERO, ONE, RWORK( IRWU ), N )
00237          CALL SLASET( 'A', N, N, ZERO, ONE, RWORK( IRWVT ), N )
00238          CALL SLASDQ( 'U', 0, N, N, N, 0, D, E, RWORK( IRWVT ), N,
00239      $                RWORK( IRWU ), N, RWORK( IRWWRK ), 1,
00240      $                RWORK( IRWWRK ), INFO )
00241          IF( INFO.NE.0 ) THEN
00242             RETURN
00243          END IF
00244 *
00245 *        In the real version, B is passed to SLASDQ and multiplied
00246 *        internally by Q'. Here B is complex and that product is
00247 *        computed below in two steps (real and imaginary parts).
00248 *
00249          J = IRWB - 1
00250          DO 50 JCOL = 1, NRHS
00251             DO 40 JROW = 1, N
00252                J = J + 1
00253                RWORK( J ) = REAL( B( JROW, JCOL ) )
00254    40       CONTINUE
00255    50    CONTINUE
00256          CALL SGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWU ), N,
00257      $               RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N )
00258          J = IRWB - 1
00259          DO 70 JCOL = 1, NRHS
00260             DO 60 JROW = 1, N
00261                J = J + 1
00262                RWORK( J ) = AIMAG( B( JROW, JCOL ) )
00263    60       CONTINUE
00264    70    CONTINUE
00265          CALL SGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWU ), N,
00266      $               RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N )
00267          JREAL = IRWRB - 1
00268          JIMAG = IRWIB - 1
00269          DO 90 JCOL = 1, NRHS
00270             DO 80 JROW = 1, N
00271                JREAL = JREAL + 1
00272                JIMAG = JIMAG + 1
00273                B( JROW, JCOL ) = CMPLX( RWORK( JREAL ), RWORK( JIMAG ) )
00274    80       CONTINUE
00275    90    CONTINUE
00276 *
00277          TOL = RCND*ABS( D( ISAMAX( N, D, 1 ) ) )
00278          DO 100 I = 1, N
00279             IF( D( I ).LE.TOL ) THEN
00280                CALL CLASET( 'A', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
00281             ELSE
00282                CALL CLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS, B( I, 1 ),
00283      $                      LDB, INFO )
00284                RANK = RANK + 1
00285             END IF
00286   100    CONTINUE
00287 *
00288 *        Since B is complex, the following call to SGEMM is performed
00289 *        in two steps (real and imaginary parts). That is for V * B
00290 *        (in the real version of the code V' is stored in WORK).
00291 *
00292 *        CALL SGEMM( 'T', 'N', N, NRHS, N, ONE, WORK, N, B, LDB, ZERO,
00293 *    $               WORK( NWORK ), N )
00294 *
00295          J = IRWB - 1
00296          DO 120 JCOL = 1, NRHS
00297             DO 110 JROW = 1, N
00298                J = J + 1
00299                RWORK( J ) = REAL( B( JROW, JCOL ) )
00300   110       CONTINUE
00301   120    CONTINUE
00302          CALL SGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWVT ), N,
00303      $               RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N )
00304          J = IRWB - 1
00305          DO 140 JCOL = 1, NRHS
00306             DO 130 JROW = 1, N
00307                J = J + 1
00308                RWORK( J ) = AIMAG( B( JROW, JCOL ) )
00309   130       CONTINUE
00310   140    CONTINUE
00311          CALL SGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWVT ), N,
00312      $               RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N )
00313          JREAL = IRWRB - 1
00314          JIMAG = IRWIB - 1
00315          DO 160 JCOL = 1, NRHS
00316             DO 150 JROW = 1, N
00317                JREAL = JREAL + 1
00318                JIMAG = JIMAG + 1
00319                B( JROW, JCOL ) = CMPLX( RWORK( JREAL ), RWORK( JIMAG ) )
00320   150       CONTINUE
00321   160    CONTINUE
00322 *
00323 *        Unscale.
00324 *
00325          CALL SLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
00326          CALL SLASRT( 'D', N, D, INFO )
00327          CALL CLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
00328 *
00329          RETURN
00330       END IF
00331 *
00332 *     Book-keeping and setting up some constants.
00333 *
00334       NLVL = INT( LOG( REAL( N ) / REAL( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1
00335 *
00336       SMLSZP = SMLSIZ + 1
00337 *
00338       U = 1
00339       VT = 1 + SMLSIZ*N
00340       DIFL = VT + SMLSZP*N
00341       DIFR = DIFL + NLVL*N
00342       Z = DIFR + NLVL*N*2
00343       C = Z + NLVL*N
00344       S = C + N
00345       POLES = S + N
00346       GIVNUM = POLES + 2*NLVL*N
00347       NRWORK = GIVNUM + 2*NLVL*N
00348       BX = 1
00349 *
00350       IRWRB = NRWORK
00351       IRWIB = IRWRB + SMLSIZ*NRHS
00352       IRWB = IRWIB + SMLSIZ*NRHS
00353 *
00354       SIZEI = 1 + N
00355       K = SIZEI + N
00356       GIVPTR = K + N
00357       PERM = GIVPTR + N
00358       GIVCOL = PERM + NLVL*N
00359       IWK = GIVCOL + NLVL*N*2
00360 *
00361       ST = 1
00362       SQRE = 0
00363       ICMPQ1 = 1
00364       ICMPQ2 = 0
00365       NSUB = 0
00366 *
00367       DO 170 I = 1, N
00368          IF( ABS( D( I ) ).LT.EPS ) THEN
00369             D( I ) = SIGN( EPS, D( I ) )
00370          END IF
00371   170 CONTINUE
00372 *
00373       DO 240 I = 1, NM1
00374          IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN
00375             NSUB = NSUB + 1
00376             IWORK( NSUB ) = ST
00377 *
00378 *           Subproblem found. First determine its size and then
00379 *           apply divide and conquer on it.
00380 *
00381             IF( I.LT.NM1 ) THEN
00382 *
00383 *              A subproblem with E(I) small for I < NM1.
00384 *
00385                NSIZE = I - ST + 1
00386                IWORK( SIZEI+NSUB-1 ) = NSIZE
00387             ELSE IF( ABS( E( I ) ).GE.EPS ) THEN
00388 *
00389 *              A subproblem with E(NM1) not too small but I = NM1.
00390 *
00391                NSIZE = N - ST + 1
00392                IWORK( SIZEI+NSUB-1 ) = NSIZE
00393             ELSE
00394 *
00395 *              A subproblem with E(NM1) small. This implies an
00396 *              1-by-1 subproblem at D(N), which is not solved
00397 *              explicitly.
00398 *
00399                NSIZE = I - ST + 1
00400                IWORK( SIZEI+NSUB-1 ) = NSIZE
00401                NSUB = NSUB + 1
00402                IWORK( NSUB ) = N
00403                IWORK( SIZEI+NSUB-1 ) = 1
00404                CALL CCOPY( NRHS, B( N, 1 ), LDB, WORK( BX+NM1 ), N )
00405             END IF
00406             ST1 = ST - 1
00407             IF( NSIZE.EQ.1 ) THEN
00408 *
00409 *              This is a 1-by-1 subproblem and is not solved
00410 *              explicitly.
00411 *
00412                CALL CCOPY( NRHS, B( ST, 1 ), LDB, WORK( BX+ST1 ), N )
00413             ELSE IF( NSIZE.LE.SMLSIZ ) THEN
00414 *
00415 *              This is a small subproblem and is solved by SLASDQ.
00416 *
00417                CALL SLASET( 'A', NSIZE, NSIZE, ZERO, ONE,
00418      $                      RWORK( VT+ST1 ), N )
00419                CALL SLASET( 'A', NSIZE, NSIZE, ZERO, ONE,
00420      $                      RWORK( U+ST1 ), N )
00421                CALL SLASDQ( 'U', 0, NSIZE, NSIZE, NSIZE, 0, D( ST ),
00422      $                      E( ST ), RWORK( VT+ST1 ), N, RWORK( U+ST1 ),
00423      $                      N, RWORK( NRWORK ), 1, RWORK( NRWORK ),
00424      $                      INFO )
00425                IF( INFO.NE.0 ) THEN
00426                   RETURN
00427                END IF
00428 *
00429 *              In the real version, B is passed to SLASDQ and multiplied
00430 *              internally by Q'. Here B is complex and that product is
00431 *              computed below in two steps (real and imaginary parts).
00432 *
00433                J = IRWB - 1
00434                DO 190 JCOL = 1, NRHS
00435                   DO 180 JROW = ST, ST + NSIZE - 1
00436                      J = J + 1
00437                      RWORK( J ) = REAL( B( JROW, JCOL ) )
00438   180             CONTINUE
00439   190          CONTINUE
00440                CALL SGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
00441      $                     RWORK( U+ST1 ), N, RWORK( IRWB ), NSIZE,
00442      $                     ZERO, RWORK( IRWRB ), NSIZE )
00443                J = IRWB - 1
00444                DO 210 JCOL = 1, NRHS
00445                   DO 200 JROW = ST, ST + NSIZE - 1
00446                      J = J + 1
00447                      RWORK( J ) = AIMAG( B( JROW, JCOL ) )
00448   200             CONTINUE
00449   210          CONTINUE
00450                CALL SGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
00451      $                     RWORK( U+ST1 ), N, RWORK( IRWB ), NSIZE,
00452      $                     ZERO, RWORK( IRWIB ), NSIZE )
00453                JREAL = IRWRB - 1
00454                JIMAG = IRWIB - 1
00455                DO 230 JCOL = 1, NRHS
00456                   DO 220 JROW = ST, ST + NSIZE - 1
00457                      JREAL = JREAL + 1
00458                      JIMAG = JIMAG + 1
00459                      B( JROW, JCOL ) = CMPLX( RWORK( JREAL ),
00460      $                                 RWORK( JIMAG ) )
00461   220             CONTINUE
00462   230          CONTINUE
00463 *
00464                CALL CLACPY( 'A', NSIZE, NRHS, B( ST, 1 ), LDB,
00465      $                      WORK( BX+ST1 ), N )
00466             ELSE
00467 *
00468 *              A large problem. Solve it using divide and conquer.
00469 *
00470                CALL SLASDA( ICMPQ1, SMLSIZ, NSIZE, SQRE, D( ST ),
00471      $                      E( ST ), RWORK( U+ST1 ), N, RWORK( VT+ST1 ),
00472      $                      IWORK( K+ST1 ), RWORK( DIFL+ST1 ),
00473      $                      RWORK( DIFR+ST1 ), RWORK( Z+ST1 ),
00474      $                      RWORK( POLES+ST1 ), IWORK( GIVPTR+ST1 ),
00475      $                      IWORK( GIVCOL+ST1 ), N, IWORK( PERM+ST1 ),
00476      $                      RWORK( GIVNUM+ST1 ), RWORK( C+ST1 ),
00477      $                      RWORK( S+ST1 ), RWORK( NRWORK ),
00478      $                      IWORK( IWK ), INFO )
00479                IF( INFO.NE.0 ) THEN
00480                   RETURN
00481                END IF
00482                BXST = BX + ST1
00483                CALL CLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, B( ST, 1 ),
00484      $                      LDB, WORK( BXST ), N, RWORK( U+ST1 ), N,
00485      $                      RWORK( VT+ST1 ), IWORK( K+ST1 ),
00486      $                      RWORK( DIFL+ST1 ), RWORK( DIFR+ST1 ),
00487      $                      RWORK( Z+ST1 ), RWORK( POLES+ST1 ),
00488      $                      IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
00489      $                      IWORK( PERM+ST1 ), RWORK( GIVNUM+ST1 ),
00490      $                      RWORK( C+ST1 ), RWORK( S+ST1 ),
00491      $                      RWORK( NRWORK ), IWORK( IWK ), INFO )
00492                IF( INFO.NE.0 ) THEN
00493                   RETURN
00494                END IF
00495             END IF
00496             ST = I + 1
00497          END IF
00498   240 CONTINUE
00499 *
00500 *     Apply the singular values and treat the tiny ones as zero.
00501 *
00502       TOL = RCND*ABS( D( ISAMAX( N, D, 1 ) ) )
00503 *
00504       DO 250 I = 1, N
00505 *
00506 *        Some of the elements in D can be negative because 1-by-1
00507 *        subproblems were not solved explicitly.
00508 *
00509          IF( ABS( D( I ) ).LE.TOL ) THEN
00510             CALL CLASET( 'A', 1, NRHS, CZERO, CZERO, WORK( BX+I-1 ), N )
00511          ELSE
00512             RANK = RANK + 1
00513             CALL CLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS,
00514      $                   WORK( BX+I-1 ), N, INFO )
00515          END IF
00516          D( I ) = ABS( D( I ) )
00517   250 CONTINUE
00518 *
00519 *     Now apply back the right singular vectors.
00520 *
00521       ICMPQ2 = 1
00522       DO 320 I = 1, NSUB
00523          ST = IWORK( I )
00524          ST1 = ST - 1
00525          NSIZE = IWORK( SIZEI+I-1 )
00526          BXST = BX + ST1
00527          IF( NSIZE.EQ.1 ) THEN
00528             CALL CCOPY( NRHS, WORK( BXST ), N, B( ST, 1 ), LDB )
00529          ELSE IF( NSIZE.LE.SMLSIZ ) THEN
00530 *
00531 *           Since B and BX are complex, the following call to SGEMM
00532 *           is performed in two steps (real and imaginary parts).
00533 *
00534 *           CALL SGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
00535 *    $                  RWORK( VT+ST1 ), N, RWORK( BXST ), N, ZERO,
00536 *    $                  B( ST, 1 ), LDB )
00537 *
00538             J = BXST - N - 1
00539             JREAL = IRWB - 1
00540             DO 270 JCOL = 1, NRHS
00541                J = J + N
00542                DO 260 JROW = 1, NSIZE
00543                   JREAL = JREAL + 1
00544                   RWORK( JREAL ) = REAL( WORK( J+JROW ) )
00545   260          CONTINUE
00546   270       CONTINUE
00547             CALL SGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
00548      $                  RWORK( VT+ST1 ), N, RWORK( IRWB ), NSIZE, ZERO,
00549      $                  RWORK( IRWRB ), NSIZE )
00550             J = BXST - N - 1
00551             JIMAG = IRWB - 1
00552             DO 290 JCOL = 1, NRHS
00553                J = J + N
00554                DO 280 JROW = 1, NSIZE
00555                   JIMAG = JIMAG + 1
00556                   RWORK( JIMAG ) = AIMAG( WORK( J+JROW ) )
00557   280          CONTINUE
00558   290       CONTINUE
00559             CALL SGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
00560      $                  RWORK( VT+ST1 ), N, RWORK( IRWB ), NSIZE, ZERO,
00561      $                  RWORK( IRWIB ), NSIZE )
00562             JREAL = IRWRB - 1
00563             JIMAG = IRWIB - 1
00564             DO 310 JCOL = 1, NRHS
00565                DO 300 JROW = ST, ST + NSIZE - 1
00566                   JREAL = JREAL + 1
00567                   JIMAG = JIMAG + 1
00568                   B( JROW, JCOL ) = CMPLX( RWORK( JREAL ),
00569      $                              RWORK( JIMAG ) )
00570   300          CONTINUE
00571   310       CONTINUE
00572          ELSE
00573             CALL CLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, WORK( BXST ), N,
00574      $                   B( ST, 1 ), LDB, RWORK( U+ST1 ), N,
00575      $                   RWORK( VT+ST1 ), IWORK( K+ST1 ),
00576      $                   RWORK( DIFL+ST1 ), RWORK( DIFR+ST1 ),
00577      $                   RWORK( Z+ST1 ), RWORK( POLES+ST1 ),
00578      $                   IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
00579      $                   IWORK( PERM+ST1 ), RWORK( GIVNUM+ST1 ),
00580      $                   RWORK( C+ST1 ), RWORK( S+ST1 ),
00581      $                   RWORK( NRWORK ), IWORK( IWK ), INFO )
00582             IF( INFO.NE.0 ) THEN
00583                RETURN
00584             END IF
00585          END IF
00586   320 CONTINUE
00587 *
00588 *     Unscale and sort the singular values.
00589 *
00590       CALL SLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
00591       CALL SLASRT( 'D', N, D, INFO )
00592       CALL CLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
00593 *
00594       RETURN
00595 *
00596 *     End of CLALSD
00597 *
00598       END