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