LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zlalsd.f
Go to the documentation of this file.
00001 *> \brief \b ZLALSD
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download ZLALSD + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlalsd.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlalsd.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlalsd.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE ZLALSD( 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 *       DOUBLE PRECISION   RCOND
00028 *       ..
00029 *       .. Array Arguments ..
00030 *       INTEGER            IWORK( * )
00031 *       DOUBLE PRECISION   D( * ), E( * ), RWORK( * )
00032 *       COMPLEX*16         B( LDB, * ), WORK( * )
00033 *       ..
00034 *  
00035 *
00036 *> \par Purpose:
00037 *  =============
00038 *>
00039 *> \verbatim
00040 *>
00041 *> ZLALSD 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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*16 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 DOUBLE PRECISION
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*16 array, dimension at least
00140 *>         (N * NRHS).
00141 *> \endverbatim
00142 *>
00143 *> \param[out] RWORK
00144 *> \verbatim
00145 *>          RWORK is DOUBLE PRECISION array, dimension at least
00146 *>         (9*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS +
00147 *>         MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS ),
00148 *>         where
00149 *>         NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )
00150 *> \endverbatim
00151 *>
00152 *> \param[out] IWORK
00153 *> \verbatim
00154 *>          IWORK is INTEGER array, dimension at least
00155 *>         (3*N*NLVL + 11*N).
00156 *> \endverbatim
00157 *>
00158 *> \param[out] INFO
00159 *> \verbatim
00160 *>          INFO is INTEGER
00161 *>         = 0:  successful exit.
00162 *>         < 0:  if INFO = -i, the i-th argument had an illegal value.
00163 *>         > 0:  The algorithm failed to compute a singular value while
00164 *>               working on the submatrix lying in rows and columns
00165 *>               INFO/(N+1) through MOD(INFO,N+1).
00166 *> \endverbatim
00167 *
00168 *  Authors:
00169 *  ========
00170 *
00171 *> \author Univ. of Tennessee 
00172 *> \author Univ. of California Berkeley 
00173 *> \author Univ. of Colorado Denver 
00174 *> \author NAG Ltd. 
00175 *
00176 *> \date November 2011
00177 *
00178 *> \ingroup complex16OTHERcomputational
00179 *
00180 *> \par Contributors:
00181 *  ==================
00182 *>
00183 *>     Ming Gu and Ren-Cang Li, Computer Science Division, University of
00184 *>       California at Berkeley, USA \n
00185 *>     Osni Marques, LBNL/NERSC, USA \n
00186 *
00187 *  =====================================================================
00188       SUBROUTINE ZLALSD( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
00189      $                   RANK, WORK, RWORK, IWORK, INFO )
00190 *
00191 *  -- LAPACK computational routine (version 3.4.0) --
00192 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00193 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00194 *     November 2011
00195 *
00196 *     .. Scalar Arguments ..
00197       CHARACTER          UPLO
00198       INTEGER            INFO, LDB, N, NRHS, RANK, SMLSIZ
00199       DOUBLE PRECISION   RCOND
00200 *     ..
00201 *     .. Array Arguments ..
00202       INTEGER            IWORK( * )
00203       DOUBLE PRECISION   D( * ), E( * ), RWORK( * )
00204       COMPLEX*16         B( LDB, * ), WORK( * )
00205 *     ..
00206 *
00207 *  =====================================================================
00208 *
00209 *     .. Parameters ..
00210       DOUBLE PRECISION   ZERO, ONE, TWO
00211       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
00212       COMPLEX*16         CZERO
00213       PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ) )
00214 *     ..
00215 *     .. Local Scalars ..
00216       INTEGER            BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM,
00217      $                   GIVPTR, I, ICMPQ1, ICMPQ2, IRWB, IRWIB, IRWRB,
00218      $                   IRWU, IRWVT, IRWWRK, IWK, J, JCOL, JIMAG,
00219      $                   JREAL, JROW, K, NLVL, NM1, NRWORK, NSIZE, NSUB,
00220      $                   PERM, POLES, S, SIZEI, SMLSZP, SQRE, ST, ST1,
00221      $                   U, VT, Z
00222       DOUBLE PRECISION   CS, EPS, ORGNRM, RCND, R, SN, TOL
00223 *     ..
00224 *     .. External Functions ..
00225       INTEGER            IDAMAX
00226       DOUBLE PRECISION   DLAMCH, DLANST
00227       EXTERNAL           IDAMAX, DLAMCH, DLANST
00228 *     ..
00229 *     .. External Subroutines ..
00230       EXTERNAL           DGEMM, DLARTG, DLASCL, DLASDA, DLASDQ, DLASET,
00231      $                   DLASRT, XERBLA, ZCOPY, ZDROT, ZLACPY, ZLALSA,
00232      $                   ZLASCL, ZLASET
00233 *     ..
00234 *     .. Intrinsic Functions ..
00235       INTRINSIC          ABS, DBLE, DCMPLX, DIMAG, INT, LOG, SIGN
00236 *     ..
00237 *     .. Executable Statements ..
00238 *
00239 *     Test the input parameters.
00240 *
00241       INFO = 0
00242 *
00243       IF( N.LT.0 ) THEN
00244          INFO = -3
00245       ELSE IF( NRHS.LT.1 ) THEN
00246          INFO = -4
00247       ELSE IF( ( LDB.LT.1 ) .OR. ( LDB.LT.N ) ) THEN
00248          INFO = -8
00249       END IF
00250       IF( INFO.NE.0 ) THEN
00251          CALL XERBLA( 'ZLALSD', -INFO )
00252          RETURN
00253       END IF
00254 *
00255       EPS = DLAMCH( 'Epsilon' )
00256 *
00257 *     Set up the tolerance.
00258 *
00259       IF( ( RCOND.LE.ZERO ) .OR. ( RCOND.GE.ONE ) ) THEN
00260          RCND = EPS
00261       ELSE
00262          RCND = RCOND
00263       END IF
00264 *
00265       RANK = 0
00266 *
00267 *     Quick return if possible.
00268 *
00269       IF( N.EQ.0 ) THEN
00270          RETURN
00271       ELSE IF( N.EQ.1 ) THEN
00272          IF( D( 1 ).EQ.ZERO ) THEN
00273             CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, B, LDB )
00274          ELSE
00275             RANK = 1
00276             CALL ZLASCL( 'G', 0, 0, D( 1 ), ONE, 1, NRHS, B, LDB, INFO )
00277             D( 1 ) = ABS( D( 1 ) )
00278          END IF
00279          RETURN
00280       END IF
00281 *
00282 *     Rotate the matrix if it is lower bidiagonal.
00283 *
00284       IF( UPLO.EQ.'L' ) THEN
00285          DO 10 I = 1, N - 1
00286             CALL DLARTG( D( I ), E( I ), CS, SN, R )
00287             D( I ) = R
00288             E( I ) = SN*D( I+1 )
00289             D( I+1 ) = CS*D( I+1 )
00290             IF( NRHS.EQ.1 ) THEN
00291                CALL ZDROT( 1, B( I, 1 ), 1, B( I+1, 1 ), 1, CS, SN )
00292             ELSE
00293                RWORK( I*2-1 ) = CS
00294                RWORK( I*2 ) = SN
00295             END IF
00296    10    CONTINUE
00297          IF( NRHS.GT.1 ) THEN
00298             DO 30 I = 1, NRHS
00299                DO 20 J = 1, N - 1
00300                   CS = RWORK( J*2-1 )
00301                   SN = RWORK( J*2 )
00302                   CALL ZDROT( 1, B( J, I ), 1, B( J+1, I ), 1, CS, SN )
00303    20          CONTINUE
00304    30       CONTINUE
00305          END IF
00306       END IF
00307 *
00308 *     Scale.
00309 *
00310       NM1 = N - 1
00311       ORGNRM = DLANST( 'M', N, D, E )
00312       IF( ORGNRM.EQ.ZERO ) THEN
00313          CALL ZLASET( 'A', N, NRHS, CZERO, CZERO, B, LDB )
00314          RETURN
00315       END IF
00316 *
00317       CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
00318       CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, INFO )
00319 *
00320 *     If N is smaller than the minimum divide size SMLSIZ, then solve
00321 *     the problem with another solver.
00322 *
00323       IF( N.LE.SMLSIZ ) THEN
00324          IRWU = 1
00325          IRWVT = IRWU + N*N
00326          IRWWRK = IRWVT + N*N
00327          IRWRB = IRWWRK
00328          IRWIB = IRWRB + N*NRHS
00329          IRWB = IRWIB + N*NRHS
00330          CALL DLASET( 'A', N, N, ZERO, ONE, RWORK( IRWU ), N )
00331          CALL DLASET( 'A', N, N, ZERO, ONE, RWORK( IRWVT ), N )
00332          CALL DLASDQ( 'U', 0, N, N, N, 0, D, E, RWORK( IRWVT ), N,
00333      $                RWORK( IRWU ), N, RWORK( IRWWRK ), 1,
00334      $                RWORK( IRWWRK ), INFO )
00335          IF( INFO.NE.0 ) THEN
00336             RETURN
00337          END IF
00338 *
00339 *        In the real version, B is passed to DLASDQ and multiplied
00340 *        internally by Q**H. Here B is complex and that product is
00341 *        computed below in two steps (real and imaginary parts).
00342 *
00343          J = IRWB - 1
00344          DO 50 JCOL = 1, NRHS
00345             DO 40 JROW = 1, N
00346                J = J + 1
00347                RWORK( J ) = DBLE( B( JROW, JCOL ) )
00348    40       CONTINUE
00349    50    CONTINUE
00350          CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWU ), N,
00351      $               RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N )
00352          J = IRWB - 1
00353          DO 70 JCOL = 1, NRHS
00354             DO 60 JROW = 1, N
00355                J = J + 1
00356                RWORK( J ) = DIMAG( B( JROW, JCOL ) )
00357    60       CONTINUE
00358    70    CONTINUE
00359          CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWU ), N,
00360      $               RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N )
00361          JREAL = IRWRB - 1
00362          JIMAG = IRWIB - 1
00363          DO 90 JCOL = 1, NRHS
00364             DO 80 JROW = 1, N
00365                JREAL = JREAL + 1
00366                JIMAG = JIMAG + 1
00367                B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
00368      $                           RWORK( JIMAG ) )
00369    80       CONTINUE
00370    90    CONTINUE
00371 *
00372          TOL = RCND*ABS( D( IDAMAX( N, D, 1 ) ) )
00373          DO 100 I = 1, N
00374             IF( D( I ).LE.TOL ) THEN
00375                CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
00376             ELSE
00377                CALL ZLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS, B( I, 1 ),
00378      $                      LDB, INFO )
00379                RANK = RANK + 1
00380             END IF
00381   100    CONTINUE
00382 *
00383 *        Since B is complex, the following call to DGEMM is performed
00384 *        in two steps (real and imaginary parts). That is for V * B
00385 *        (in the real version of the code V**H is stored in WORK).
00386 *
00387 *        CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, WORK, N, B, LDB, ZERO,
00388 *    $               WORK( NWORK ), N )
00389 *
00390          J = IRWB - 1
00391          DO 120 JCOL = 1, NRHS
00392             DO 110 JROW = 1, N
00393                J = J + 1
00394                RWORK( J ) = DBLE( B( JROW, JCOL ) )
00395   110       CONTINUE
00396   120    CONTINUE
00397          CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWVT ), N,
00398      $               RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N )
00399          J = IRWB - 1
00400          DO 140 JCOL = 1, NRHS
00401             DO 130 JROW = 1, N
00402                J = J + 1
00403                RWORK( J ) = DIMAG( B( JROW, JCOL ) )
00404   130       CONTINUE
00405   140    CONTINUE
00406          CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWVT ), N,
00407      $               RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N )
00408          JREAL = IRWRB - 1
00409          JIMAG = IRWIB - 1
00410          DO 160 JCOL = 1, NRHS
00411             DO 150 JROW = 1, N
00412                JREAL = JREAL + 1
00413                JIMAG = JIMAG + 1
00414                B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
00415      $                           RWORK( JIMAG ) )
00416   150       CONTINUE
00417   160    CONTINUE
00418 *
00419 *        Unscale.
00420 *
00421          CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
00422          CALL DLASRT( 'D', N, D, INFO )
00423          CALL ZLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
00424 *
00425          RETURN
00426       END IF
00427 *
00428 *     Book-keeping and setting up some constants.
00429 *
00430       NLVL = INT( LOG( DBLE( N ) / DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1
00431 *
00432       SMLSZP = SMLSIZ + 1
00433 *
00434       U = 1
00435       VT = 1 + SMLSIZ*N
00436       DIFL = VT + SMLSZP*N
00437       DIFR = DIFL + NLVL*N
00438       Z = DIFR + NLVL*N*2
00439       C = Z + NLVL*N
00440       S = C + N
00441       POLES = S + N
00442       GIVNUM = POLES + 2*NLVL*N
00443       NRWORK = GIVNUM + 2*NLVL*N
00444       BX = 1
00445 *
00446       IRWRB = NRWORK
00447       IRWIB = IRWRB + SMLSIZ*NRHS
00448       IRWB = IRWIB + SMLSIZ*NRHS
00449 *
00450       SIZEI = 1 + N
00451       K = SIZEI + N
00452       GIVPTR = K + N
00453       PERM = GIVPTR + N
00454       GIVCOL = PERM + NLVL*N
00455       IWK = GIVCOL + NLVL*N*2
00456 *
00457       ST = 1
00458       SQRE = 0
00459       ICMPQ1 = 1
00460       ICMPQ2 = 0
00461       NSUB = 0
00462 *
00463       DO 170 I = 1, N
00464          IF( ABS( D( I ) ).LT.EPS ) THEN
00465             D( I ) = SIGN( EPS, D( I ) )
00466          END IF
00467   170 CONTINUE
00468 *
00469       DO 240 I = 1, NM1
00470          IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN
00471             NSUB = NSUB + 1
00472             IWORK( NSUB ) = ST
00473 *
00474 *           Subproblem found. First determine its size and then
00475 *           apply divide and conquer on it.
00476 *
00477             IF( I.LT.NM1 ) THEN
00478 *
00479 *              A subproblem with E(I) small for I < NM1.
00480 *
00481                NSIZE = I - ST + 1
00482                IWORK( SIZEI+NSUB-1 ) = NSIZE
00483             ELSE IF( ABS( E( I ) ).GE.EPS ) THEN
00484 *
00485 *              A subproblem with E(NM1) not too small but I = NM1.
00486 *
00487                NSIZE = N - ST + 1
00488                IWORK( SIZEI+NSUB-1 ) = NSIZE
00489             ELSE
00490 *
00491 *              A subproblem with E(NM1) small. This implies an
00492 *              1-by-1 subproblem at D(N), which is not solved
00493 *              explicitly.
00494 *
00495                NSIZE = I - ST + 1
00496                IWORK( SIZEI+NSUB-1 ) = NSIZE
00497                NSUB = NSUB + 1
00498                IWORK( NSUB ) = N
00499                IWORK( SIZEI+NSUB-1 ) = 1
00500                CALL ZCOPY( NRHS, B( N, 1 ), LDB, WORK( BX+NM1 ), N )
00501             END IF
00502             ST1 = ST - 1
00503             IF( NSIZE.EQ.1 ) THEN
00504 *
00505 *              This is a 1-by-1 subproblem and is not solved
00506 *              explicitly.
00507 *
00508                CALL ZCOPY( NRHS, B( ST, 1 ), LDB, WORK( BX+ST1 ), N )
00509             ELSE IF( NSIZE.LE.SMLSIZ ) THEN
00510 *
00511 *              This is a small subproblem and is solved by DLASDQ.
00512 *
00513                CALL DLASET( 'A', NSIZE, NSIZE, ZERO, ONE,
00514      $                      RWORK( VT+ST1 ), N )
00515                CALL DLASET( 'A', NSIZE, NSIZE, ZERO, ONE,
00516      $                      RWORK( U+ST1 ), N )
00517                CALL DLASDQ( 'U', 0, NSIZE, NSIZE, NSIZE, 0, D( ST ),
00518      $                      E( ST ), RWORK( VT+ST1 ), N, RWORK( U+ST1 ),
00519      $                      N, RWORK( NRWORK ), 1, RWORK( NRWORK ),
00520      $                      INFO )
00521                IF( INFO.NE.0 ) THEN
00522                   RETURN
00523                END IF
00524 *
00525 *              In the real version, B is passed to DLASDQ and multiplied
00526 *              internally by Q**H. Here B is complex and that product is
00527 *              computed below in two steps (real and imaginary parts).
00528 *
00529                J = IRWB - 1
00530                DO 190 JCOL = 1, NRHS
00531                   DO 180 JROW = ST, ST + NSIZE - 1
00532                      J = J + 1
00533                      RWORK( J ) = DBLE( B( JROW, JCOL ) )
00534   180             CONTINUE
00535   190          CONTINUE
00536                CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
00537      $                     RWORK( U+ST1 ), N, RWORK( IRWB ), NSIZE,
00538      $                     ZERO, RWORK( IRWRB ), NSIZE )
00539                J = IRWB - 1
00540                DO 210 JCOL = 1, NRHS
00541                   DO 200 JROW = ST, ST + NSIZE - 1
00542                      J = J + 1
00543                      RWORK( J ) = DIMAG( B( JROW, JCOL ) )
00544   200             CONTINUE
00545   210          CONTINUE
00546                CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
00547      $                     RWORK( U+ST1 ), N, RWORK( IRWB ), NSIZE,
00548      $                     ZERO, RWORK( IRWIB ), NSIZE )
00549                JREAL = IRWRB - 1
00550                JIMAG = IRWIB - 1
00551                DO 230 JCOL = 1, NRHS
00552                   DO 220 JROW = ST, ST + NSIZE - 1
00553                      JREAL = JREAL + 1
00554                      JIMAG = JIMAG + 1
00555                      B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
00556      $                                 RWORK( JIMAG ) )
00557   220             CONTINUE
00558   230          CONTINUE
00559 *
00560                CALL ZLACPY( 'A', NSIZE, NRHS, B( ST, 1 ), LDB,
00561      $                      WORK( BX+ST1 ), N )
00562             ELSE
00563 *
00564 *              A large problem. Solve it using divide and conquer.
00565 *
00566                CALL DLASDA( ICMPQ1, SMLSIZ, NSIZE, SQRE, D( ST ),
00567      $                      E( ST ), RWORK( U+ST1 ), N, RWORK( VT+ST1 ),
00568      $                      IWORK( K+ST1 ), RWORK( DIFL+ST1 ),
00569      $                      RWORK( DIFR+ST1 ), RWORK( Z+ST1 ),
00570      $                      RWORK( POLES+ST1 ), IWORK( GIVPTR+ST1 ),
00571      $                      IWORK( GIVCOL+ST1 ), N, IWORK( PERM+ST1 ),
00572      $                      RWORK( GIVNUM+ST1 ), RWORK( C+ST1 ),
00573      $                      RWORK( S+ST1 ), RWORK( NRWORK ),
00574      $                      IWORK( IWK ), INFO )
00575                IF( INFO.NE.0 ) THEN
00576                   RETURN
00577                END IF
00578                BXST = BX + ST1
00579                CALL ZLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, B( ST, 1 ),
00580      $                      LDB, WORK( BXST ), N, RWORK( U+ST1 ), N,
00581      $                      RWORK( VT+ST1 ), IWORK( K+ST1 ),
00582      $                      RWORK( DIFL+ST1 ), RWORK( DIFR+ST1 ),
00583      $                      RWORK( Z+ST1 ), RWORK( POLES+ST1 ),
00584      $                      IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
00585      $                      IWORK( PERM+ST1 ), RWORK( GIVNUM+ST1 ),
00586      $                      RWORK( C+ST1 ), RWORK( S+ST1 ),
00587      $                      RWORK( NRWORK ), IWORK( IWK ), INFO )
00588                IF( INFO.NE.0 ) THEN
00589                   RETURN
00590                END IF
00591             END IF
00592             ST = I + 1
00593          END IF
00594   240 CONTINUE
00595 *
00596 *     Apply the singular values and treat the tiny ones as zero.
00597 *
00598       TOL = RCND*ABS( D( IDAMAX( N, D, 1 ) ) )
00599 *
00600       DO 250 I = 1, N
00601 *
00602 *        Some of the elements in D can be negative because 1-by-1
00603 *        subproblems were not solved explicitly.
00604 *
00605          IF( ABS( D( I ) ).LE.TOL ) THEN
00606             CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, WORK( BX+I-1 ), N )
00607          ELSE
00608             RANK = RANK + 1
00609             CALL ZLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS,
00610      $                   WORK( BX+I-1 ), N, INFO )
00611          END IF
00612          D( I ) = ABS( D( I ) )
00613   250 CONTINUE
00614 *
00615 *     Now apply back the right singular vectors.
00616 *
00617       ICMPQ2 = 1
00618       DO 320 I = 1, NSUB
00619          ST = IWORK( I )
00620          ST1 = ST - 1
00621          NSIZE = IWORK( SIZEI+I-1 )
00622          BXST = BX + ST1
00623          IF( NSIZE.EQ.1 ) THEN
00624             CALL ZCOPY( NRHS, WORK( BXST ), N, B( ST, 1 ), LDB )
00625          ELSE IF( NSIZE.LE.SMLSIZ ) THEN
00626 *
00627 *           Since B and BX are complex, the following call to DGEMM
00628 *           is performed in two steps (real and imaginary parts).
00629 *
00630 *           CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
00631 *    $                  RWORK( VT+ST1 ), N, RWORK( BXST ), N, ZERO,
00632 *    $                  B( ST, 1 ), LDB )
00633 *
00634             J = BXST - N - 1
00635             JREAL = IRWB - 1
00636             DO 270 JCOL = 1, NRHS
00637                J = J + N
00638                DO 260 JROW = 1, NSIZE
00639                   JREAL = JREAL + 1
00640                   RWORK( JREAL ) = DBLE( WORK( J+JROW ) )
00641   260          CONTINUE
00642   270       CONTINUE
00643             CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
00644      $                  RWORK( VT+ST1 ), N, RWORK( IRWB ), NSIZE, ZERO,
00645      $                  RWORK( IRWRB ), NSIZE )
00646             J = BXST - N - 1
00647             JIMAG = IRWB - 1
00648             DO 290 JCOL = 1, NRHS
00649                J = J + N
00650                DO 280 JROW = 1, NSIZE
00651                   JIMAG = JIMAG + 1
00652                   RWORK( JIMAG ) = DIMAG( WORK( J+JROW ) )
00653   280          CONTINUE
00654   290       CONTINUE
00655             CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
00656      $                  RWORK( VT+ST1 ), N, RWORK( IRWB ), NSIZE, ZERO,
00657      $                  RWORK( IRWIB ), NSIZE )
00658             JREAL = IRWRB - 1
00659             JIMAG = IRWIB - 1
00660             DO 310 JCOL = 1, NRHS
00661                DO 300 JROW = ST, ST + NSIZE - 1
00662                   JREAL = JREAL + 1
00663                   JIMAG = JIMAG + 1
00664                   B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
00665      $                              RWORK( JIMAG ) )
00666   300          CONTINUE
00667   310       CONTINUE
00668          ELSE
00669             CALL ZLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, WORK( BXST ), N,
00670      $                   B( ST, 1 ), LDB, RWORK( U+ST1 ), N,
00671      $                   RWORK( VT+ST1 ), IWORK( K+ST1 ),
00672      $                   RWORK( DIFL+ST1 ), RWORK( DIFR+ST1 ),
00673      $                   RWORK( Z+ST1 ), RWORK( POLES+ST1 ),
00674      $                   IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
00675      $                   IWORK( PERM+ST1 ), RWORK( GIVNUM+ST1 ),
00676      $                   RWORK( C+ST1 ), RWORK( S+ST1 ),
00677      $                   RWORK( NRWORK ), IWORK( IWK ), INFO )
00678             IF( INFO.NE.0 ) THEN
00679                RETURN
00680             END IF
00681          END IF
00682   320 CONTINUE
00683 *
00684 *     Unscale and sort the singular values.
00685 *
00686       CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
00687       CALL DLASRT( 'D', N, D, INFO )
00688       CALL ZLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
00689 *
00690       RETURN
00691 *
00692 *     End of ZLALSD
00693 *
00694       END
 All Files Functions