LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zgelss.f
Go to the documentation of this file.
00001 *> \brief <b> ZGELSS solves overdetermined or underdetermined systems for GE matrices</b>
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download ZGELSS + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgelss.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgelss.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgelss.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE ZGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
00022 *                          WORK, LWORK, RWORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
00026 *       DOUBLE PRECISION   RCOND
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       DOUBLE PRECISION   RWORK( * ), S( * )
00030 *       COMPLEX*16         A( LDA, * ), B( LDB, * ), WORK( * )
00031 *       ..
00032 *  
00033 *
00034 *> \par Purpose:
00035 *  =============
00036 *>
00037 *> \verbatim
00038 *>
00039 *> ZGELSS computes the minimum norm solution to a complex linear
00040 *> least squares problem:
00041 *>
00042 *> Minimize 2-norm(| b - A*x |).
00043 *>
00044 *> using the singular value decomposition (SVD) of A. A is an M-by-N
00045 *> matrix which may be rank-deficient.
00046 *>
00047 *> Several right hand side vectors b and solution vectors x can be
00048 *> handled in a single call; they are stored as the columns of the
00049 *> M-by-NRHS right hand side matrix B and the N-by-NRHS solution matrix
00050 *> X.
00051 *>
00052 *> The effective rank of A is determined by treating as zero those
00053 *> singular values which are less than RCOND times the largest singular
00054 *> value.
00055 *> \endverbatim
00056 *
00057 *  Arguments:
00058 *  ==========
00059 *
00060 *> \param[in] M
00061 *> \verbatim
00062 *>          M is INTEGER
00063 *>          The number of rows of the matrix A. M >= 0.
00064 *> \endverbatim
00065 *>
00066 *> \param[in] N
00067 *> \verbatim
00068 *>          N is INTEGER
00069 *>          The number of columns of the matrix A. N >= 0.
00070 *> \endverbatim
00071 *>
00072 *> \param[in] NRHS
00073 *> \verbatim
00074 *>          NRHS is INTEGER
00075 *>          The number of right hand sides, i.e., the number of columns
00076 *>          of the matrices B and X. NRHS >= 0.
00077 *> \endverbatim
00078 *>
00079 *> \param[in,out] A
00080 *> \verbatim
00081 *>          A is COMPLEX*16 array, dimension (LDA,N)
00082 *>          On entry, the M-by-N matrix A.
00083 *>          On exit, the first min(m,n) rows of A are overwritten with
00084 *>          its right singular vectors, stored rowwise.
00085 *> \endverbatim
00086 *>
00087 *> \param[in] LDA
00088 *> \verbatim
00089 *>          LDA is INTEGER
00090 *>          The leading dimension of the array A. LDA >= max(1,M).
00091 *> \endverbatim
00092 *>
00093 *> \param[in,out] B
00094 *> \verbatim
00095 *>          B is COMPLEX*16 array, dimension (LDB,NRHS)
00096 *>          On entry, the M-by-NRHS right hand side matrix B.
00097 *>          On exit, B is overwritten by the N-by-NRHS solution matrix X.
00098 *>          If m >= n and RANK = n, the residual sum-of-squares for
00099 *>          the solution in the i-th column is given by the sum of
00100 *>          squares of the modulus of elements n+1:m in that column.
00101 *> \endverbatim
00102 *>
00103 *> \param[in] LDB
00104 *> \verbatim
00105 *>          LDB is INTEGER
00106 *>          The leading dimension of the array B.  LDB >= max(1,M,N).
00107 *> \endverbatim
00108 *>
00109 *> \param[out] S
00110 *> \verbatim
00111 *>          S is DOUBLE PRECISION array, dimension (min(M,N))
00112 *>          The singular values of A in decreasing order.
00113 *>          The condition number of A in the 2-norm = S(1)/S(min(m,n)).
00114 *> \endverbatim
00115 *>
00116 *> \param[in] RCOND
00117 *> \verbatim
00118 *>          RCOND is DOUBLE PRECISION
00119 *>          RCOND is used to determine the effective rank of A.
00120 *>          Singular values S(i) <= RCOND*S(1) are treated as zero.
00121 *>          If RCOND < 0, machine precision is used instead.
00122 *> \endverbatim
00123 *>
00124 *> \param[out] RANK
00125 *> \verbatim
00126 *>          RANK is INTEGER
00127 *>          The effective rank of A, i.e., the number of singular values
00128 *>          which are greater than RCOND*S(1).
00129 *> \endverbatim
00130 *>
00131 *> \param[out] WORK
00132 *> \verbatim
00133 *>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
00134 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00135 *> \endverbatim
00136 *>
00137 *> \param[in] LWORK
00138 *> \verbatim
00139 *>          LWORK is INTEGER
00140 *>          The dimension of the array WORK. LWORK >= 1, and also:
00141 *>          LWORK >=  2*min(M,N) + max(M,N,NRHS)
00142 *>          For good performance, LWORK should generally be larger.
00143 *>
00144 *>          If LWORK = -1, then a workspace query is assumed; the routine
00145 *>          only calculates the optimal size of the WORK array, returns
00146 *>          this value as the first entry of the WORK array, and no error
00147 *>          message related to LWORK is issued by XERBLA.
00148 *> \endverbatim
00149 *>
00150 *> \param[out] RWORK
00151 *> \verbatim
00152 *>          RWORK is DOUBLE PRECISION array, dimension (5*min(M,N))
00153 *> \endverbatim
00154 *>
00155 *> \param[out] INFO
00156 *> \verbatim
00157 *>          INFO is INTEGER
00158 *>          = 0:  successful exit
00159 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
00160 *>          > 0:  the algorithm for computing the SVD failed to converge;
00161 *>                if INFO = i, i off-diagonal elements of an intermediate
00162 *>                bidiagonal form did not converge to zero.
00163 *> \endverbatim
00164 *
00165 *  Authors:
00166 *  ========
00167 *
00168 *> \author Univ. of Tennessee 
00169 *> \author Univ. of California Berkeley 
00170 *> \author Univ. of Colorado Denver 
00171 *> \author NAG Ltd. 
00172 *
00173 *> \date November 2011
00174 *
00175 *> \ingroup complex16GEsolve
00176 *
00177 *  =====================================================================
00178       SUBROUTINE ZGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
00179      $                   WORK, LWORK, RWORK, INFO )
00180 *
00181 *  -- LAPACK driver routine (version 3.4.0) --
00182 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00183 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00184 *     November 2011
00185 *
00186 *     .. Scalar Arguments ..
00187       INTEGER            INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
00188       DOUBLE PRECISION   RCOND
00189 *     ..
00190 *     .. Array Arguments ..
00191       DOUBLE PRECISION   RWORK( * ), S( * )
00192       COMPLEX*16         A( LDA, * ), B( LDB, * ), WORK( * )
00193 *     ..
00194 *
00195 *  =====================================================================
00196 *
00197 *     .. Parameters ..
00198       DOUBLE PRECISION   ZERO, ONE
00199       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00200       COMPLEX*16         CZERO, CONE
00201       PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
00202      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
00203 *     ..
00204 *     .. Local Scalars ..
00205       LOGICAL            LQUERY
00206       INTEGER            BL, CHUNK, I, IASCL, IBSCL, IE, IL, IRWORK,
00207      $                   ITAU, ITAUP, ITAUQ, IWORK, LDWORK, MAXMN,
00208      $                   MAXWRK, MINMN, MINWRK, MM, MNTHR
00209       INTEGER            LWORK_ZGEQRF, LWORK_ZUNMQR, LWORK_ZGEBRD,
00210      $                   LWORK_ZUNMBR, LWORK_ZUNGBR, LWORK_ZUNMLQ,
00211      $                   LWORK_ZGELQF
00212       DOUBLE PRECISION   ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM, THR
00213 *     ..
00214 *     .. Local Arrays ..
00215       COMPLEX*16         DUM( 1 )
00216 *     ..
00217 *     .. External Subroutines ..
00218       EXTERNAL           DLABAD, DLASCL, DLASET, XERBLA, ZBDSQR, ZCOPY,
00219      $                   ZDRSCL, ZGEBRD, ZGELQF, ZGEMM, ZGEMV, ZGEQRF,
00220      $                   ZLACPY, ZLASCL, ZLASET, ZUNGBR, ZUNMBR, ZUNMLQ,
00221      $                   ZUNMQR
00222 *     ..
00223 *     .. External Functions ..
00224       INTEGER            ILAENV
00225       DOUBLE PRECISION   DLAMCH, ZLANGE
00226       EXTERNAL           ILAENV, DLAMCH, ZLANGE
00227 *     ..
00228 *     .. Intrinsic Functions ..
00229       INTRINSIC          MAX, MIN
00230 *     ..
00231 *     .. Executable Statements ..
00232 *
00233 *     Test the input arguments
00234 *
00235       INFO = 0
00236       MINMN = MIN( M, N )
00237       MAXMN = MAX( M, N )
00238       LQUERY = ( LWORK.EQ.-1 )
00239       IF( M.LT.0 ) THEN
00240          INFO = -1
00241       ELSE IF( N.LT.0 ) THEN
00242          INFO = -2
00243       ELSE IF( NRHS.LT.0 ) THEN
00244          INFO = -3
00245       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00246          INFO = -5
00247       ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN
00248          INFO = -7
00249       END IF
00250 *
00251 *     Compute workspace
00252 *      (Note: Comments in the code beginning "Workspace:" describe the
00253 *       minimal amount of workspace needed at that point in the code,
00254 *       as well as the preferred amount for good performance.
00255 *       CWorkspace refers to complex workspace, and RWorkspace refers
00256 *       to real workspace. NB refers to the optimal block size for the
00257 *       immediately following subroutine, as returned by ILAENV.)
00258 *
00259       IF( INFO.EQ.0 ) THEN
00260          MINWRK = 1
00261          MAXWRK = 1
00262          IF( MINMN.GT.0 ) THEN
00263             MM = M
00264             MNTHR = ILAENV( 6, 'ZGELSS', ' ', M, N, NRHS, -1 )
00265             IF( M.GE.N .AND. M.GE.MNTHR ) THEN
00266 *
00267 *              Path 1a - overdetermined, with many more rows than
00268 *                        columns
00269 *
00270 *              Compute space needed for ZGEQRF
00271                CALL ZGEQRF( M, N, A, LDA, DUM(1), DUM(1), -1, INFO )
00272                LWORK_ZGEQRF=DUM(1)
00273 *              Compute space needed for ZUNMQR
00274                CALL ZUNMQR( 'L', 'C', M, NRHS, N, A, LDA, DUM(1), B,
00275      $                   LDB, DUM(1), -1, INFO )
00276                LWORK_ZUNMQR=DUM(1)
00277                MM = N
00278                MAXWRK = MAX( MAXWRK, N + N*ILAENV( 1, 'ZGEQRF', ' ', M,
00279      $                       N, -1, -1 ) )
00280                MAXWRK = MAX( MAXWRK, N + NRHS*ILAENV( 1, 'ZUNMQR', 'LC',
00281      $                       M, NRHS, N, -1 ) )
00282             END IF
00283             IF( M.GE.N ) THEN
00284 *
00285 *              Path 1 - overdetermined or exactly determined
00286 *
00287 *              Compute space needed for ZGEBRD
00288                CALL ZGEBRD( MM, N, A, LDA, S, DUM(1), DUM(1),
00289      $                      DUM(1), DUM(1), -1, INFO )
00290                LWORK_ZGEBRD=DUM(1)
00291 *              Compute space needed for ZUNMBR
00292                CALL ZUNMBR( 'Q', 'L', 'C', MM, NRHS, N, A, LDA, DUM(1),
00293      $                B, LDB, DUM(1), -1, INFO )
00294                LWORK_ZUNMBR=DUM(1)
00295 *              Compute space needed for ZUNGBR
00296                CALL ZUNGBR( 'P', N, N, N, A, LDA, DUM(1),
00297      $                   DUM(1), -1, INFO )
00298                LWORK_ZUNGBR=DUM(1)
00299 *              Compute total workspace needed 
00300                MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZGEBRD )
00301                MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNMBR )
00302                MAXWRK = MAX( MAXWRK, 2*N + LWORK_ZUNGBR )
00303                MAXWRK = MAX( MAXWRK, N*NRHS )
00304                MINWRK = 2*N + MAX( NRHS, M )
00305             END IF
00306             IF( N.GT.M ) THEN
00307                MINWRK = 2*M + MAX( NRHS, N )
00308                IF( N.GE.MNTHR ) THEN
00309 *
00310 *                 Path 2a - underdetermined, with many more columns
00311 *                 than rows
00312 *
00313 *                 Compute space needed for ZGELQF
00314                   CALL ZGELQF( M, N, A, LDA, DUM(1), DUM(1),
00315      $                -1, INFO )
00316                   LWORK_ZGELQF=DUM(1)
00317 *                 Compute space needed for ZGEBRD
00318                   CALL ZGEBRD( M, M, A, LDA, S, DUM(1), DUM(1),
00319      $                      DUM(1), DUM(1), -1, INFO )
00320                   LWORK_ZGEBRD=DUM(1)
00321 *                 Compute space needed for ZUNMBR
00322                   CALL ZUNMBR( 'Q', 'L', 'C', M, NRHS, N, A, LDA, 
00323      $                DUM(1), B, LDB, DUM(1), -1, INFO )
00324                   LWORK_ZUNMBR=DUM(1)
00325 *                 Compute space needed for ZUNGBR
00326                   CALL ZUNGBR( 'P', M, M, M, A, LDA, DUM(1),
00327      $                   DUM(1), -1, INFO )
00328                   LWORK_ZUNGBR=DUM(1)
00329 *                 Compute space needed for ZUNMLQ
00330                   CALL ZUNMLQ( 'L', 'C', N, NRHS, M, A, LDA, DUM(1),
00331      $                 B, LDB, DUM(1), -1, INFO )
00332                   LWORK_ZUNMLQ=DUM(1)
00333 *                 Compute total workspace needed 
00334                   MAXWRK = M + LWORK_ZGELQF
00335                   MAXWRK = MAX( MAXWRK, 3*M + M*M + LWORK_ZGEBRD )
00336                   MAXWRK = MAX( MAXWRK, 3*M + M*M + LWORK_ZUNMBR )
00337                   MAXWRK = MAX( MAXWRK, 3*M + M*M + LWORK_ZUNGBR )
00338                   IF( NRHS.GT.1 ) THEN
00339                      MAXWRK = MAX( MAXWRK, M*M + M + M*NRHS )
00340                   ELSE
00341                      MAXWRK = MAX( MAXWRK, M*M + 2*M )
00342                   END IF
00343                   MAXWRK = MAX( MAXWRK, M + LWORK_ZUNMLQ )
00344                ELSE
00345 *
00346 *                 Path 2 - underdetermined
00347 *
00348 *                 Compute space needed for ZGEBRD
00349                   CALL ZGEBRD( M, N, A, LDA, S, DUM(1), DUM(1),
00350      $                      DUM(1), DUM(1), -1, INFO )
00351                   LWORK_ZGEBRD=DUM(1)
00352 *                 Compute space needed for ZUNMBR
00353                   CALL ZUNMBR( 'Q', 'L', 'C', M, NRHS, M, A, LDA, 
00354      $                DUM(1), B, LDB, DUM(1), -1, INFO )
00355                   LWORK_ZUNMBR=DUM(1)
00356 *                 Compute space needed for ZUNGBR
00357                   CALL ZUNGBR( 'P', M, N, M, A, LDA, DUM(1),
00358      $                   DUM(1), -1, INFO )
00359                   LWORK_ZUNGBR=DUM(1)
00360                   MAXWRK = 2*M + LWORK_ZGEBRD
00361                   MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNMBR )
00362                   MAXWRK = MAX( MAXWRK, 2*M + LWORK_ZUNGBR )
00363                   MAXWRK = MAX( MAXWRK, N*NRHS )
00364                END IF
00365             END IF
00366             MAXWRK = MAX( MINWRK, MAXWRK )
00367          END IF
00368          WORK( 1 ) = MAXWRK
00369 *
00370          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY )
00371      $      INFO = -12
00372       END IF
00373 *
00374       IF( INFO.NE.0 ) THEN
00375          CALL XERBLA( 'ZGELSS', -INFO )
00376          RETURN
00377       ELSE IF( LQUERY ) THEN
00378          RETURN
00379       END IF
00380 *
00381 *     Quick return if possible
00382 *
00383       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
00384          RANK = 0
00385          RETURN
00386       END IF
00387 *
00388 *     Get machine parameters
00389 *
00390       EPS = DLAMCH( 'P' )
00391       SFMIN = DLAMCH( 'S' )
00392       SMLNUM = SFMIN / EPS
00393       BIGNUM = ONE / SMLNUM
00394       CALL DLABAD( SMLNUM, BIGNUM )
00395 *
00396 *     Scale A if max element outside range [SMLNUM,BIGNUM]
00397 *
00398       ANRM = ZLANGE( 'M', M, N, A, LDA, RWORK )
00399       IASCL = 0
00400       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
00401 *
00402 *        Scale matrix norm up to SMLNUM
00403 *
00404          CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
00405          IASCL = 1
00406       ELSE IF( ANRM.GT.BIGNUM ) THEN
00407 *
00408 *        Scale matrix norm down to BIGNUM
00409 *
00410          CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
00411          IASCL = 2
00412       ELSE IF( ANRM.EQ.ZERO ) THEN
00413 *
00414 *        Matrix all zero. Return zero solution.
00415 *
00416          CALL ZLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
00417          CALL DLASET( 'F', MINMN, 1, ZERO, ZERO, S, MINMN )
00418          RANK = 0
00419          GO TO 70
00420       END IF
00421 *
00422 *     Scale B if max element outside range [SMLNUM,BIGNUM]
00423 *
00424       BNRM = ZLANGE( 'M', M, NRHS, B, LDB, RWORK )
00425       IBSCL = 0
00426       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
00427 *
00428 *        Scale matrix norm up to SMLNUM
00429 *
00430          CALL ZLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
00431          IBSCL = 1
00432       ELSE IF( BNRM.GT.BIGNUM ) THEN
00433 *
00434 *        Scale matrix norm down to BIGNUM
00435 *
00436          CALL ZLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
00437          IBSCL = 2
00438       END IF
00439 *
00440 *     Overdetermined case
00441 *
00442       IF( M.GE.N ) THEN
00443 *
00444 *        Path 1 - overdetermined or exactly determined
00445 *
00446          MM = M
00447          IF( M.GE.MNTHR ) THEN
00448 *
00449 *           Path 1a - overdetermined, with many more rows than columns
00450 *
00451             MM = N
00452             ITAU = 1
00453             IWORK = ITAU + N
00454 *
00455 *           Compute A=Q*R
00456 *           (CWorkspace: need 2*N, prefer N+N*NB)
00457 *           (RWorkspace: none)
00458 *
00459             CALL ZGEQRF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
00460      $                   LWORK-IWORK+1, INFO )
00461 *
00462 *           Multiply B by transpose(Q)
00463 *           (CWorkspace: need N+NRHS, prefer N+NRHS*NB)
00464 *           (RWorkspace: none)
00465 *
00466             CALL ZUNMQR( 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAU ), B,
00467      $                   LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
00468 *
00469 *           Zero out below R
00470 *
00471             IF( N.GT.1 )
00472      $         CALL ZLASET( 'L', N-1, N-1, CZERO, CZERO, A( 2, 1 ),
00473      $                      LDA )
00474          END IF
00475 *
00476          IE = 1
00477          ITAUQ = 1
00478          ITAUP = ITAUQ + N
00479          IWORK = ITAUP + N
00480 *
00481 *        Bidiagonalize R in A
00482 *        (CWorkspace: need 2*N+MM, prefer 2*N+(MM+N)*NB)
00483 *        (RWorkspace: need N)
00484 *
00485          CALL ZGEBRD( MM, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
00486      $                WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
00487      $                INFO )
00488 *
00489 *        Multiply B by transpose of left bidiagonalizing vectors of R
00490 *        (CWorkspace: need 2*N+NRHS, prefer 2*N+NRHS*NB)
00491 *        (RWorkspace: none)
00492 *
00493          CALL ZUNMBR( 'Q', 'L', 'C', MM, NRHS, N, A, LDA, WORK( ITAUQ ),
00494      $                B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
00495 *
00496 *        Generate right bidiagonalizing vectors of R in A
00497 *        (CWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
00498 *        (RWorkspace: none)
00499 *
00500          CALL ZUNGBR( 'P', N, N, N, A, LDA, WORK( ITAUP ),
00501      $                WORK( IWORK ), LWORK-IWORK+1, INFO )
00502          IRWORK = IE + N
00503 *
00504 *        Perform bidiagonal QR iteration
00505 *          multiply B by transpose of left singular vectors
00506 *          compute right singular vectors in A
00507 *        (CWorkspace: none)
00508 *        (RWorkspace: need BDSPAC)
00509 *
00510          CALL ZBDSQR( 'U', N, N, 0, NRHS, S, RWORK( IE ), A, LDA, DUM,
00511      $                1, B, LDB, RWORK( IRWORK ), INFO )
00512          IF( INFO.NE.0 )
00513      $      GO TO 70
00514 *
00515 *        Multiply B by reciprocals of singular values
00516 *
00517          THR = MAX( RCOND*S( 1 ), SFMIN )
00518          IF( RCOND.LT.ZERO )
00519      $      THR = MAX( EPS*S( 1 ), SFMIN )
00520          RANK = 0
00521          DO 10 I = 1, N
00522             IF( S( I ).GT.THR ) THEN
00523                CALL ZDRSCL( NRHS, S( I ), B( I, 1 ), LDB )
00524                RANK = RANK + 1
00525             ELSE
00526                CALL ZLASET( 'F', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
00527             END IF
00528    10    CONTINUE
00529 *
00530 *        Multiply B by right singular vectors
00531 *        (CWorkspace: need N, prefer N*NRHS)
00532 *        (RWorkspace: none)
00533 *
00534          IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN
00535             CALL ZGEMM( 'C', 'N', N, NRHS, N, CONE, A, LDA, B, LDB,
00536      $                  CZERO, WORK, LDB )
00537             CALL ZLACPY( 'G', N, NRHS, WORK, LDB, B, LDB )
00538          ELSE IF( NRHS.GT.1 ) THEN
00539             CHUNK = LWORK / N
00540             DO 20 I = 1, NRHS, CHUNK
00541                BL = MIN( NRHS-I+1, CHUNK )
00542                CALL ZGEMM( 'C', 'N', N, BL, N, CONE, A, LDA, B( 1, I ),
00543      $                     LDB, CZERO, WORK, N )
00544                CALL ZLACPY( 'G', N, BL, WORK, N, B( 1, I ), LDB )
00545    20       CONTINUE
00546          ELSE
00547             CALL ZGEMV( 'C', N, N, CONE, A, LDA, B, 1, CZERO, WORK, 1 )
00548             CALL ZCOPY( N, WORK, 1, B, 1 )
00549          END IF
00550 *
00551       ELSE IF( N.GE.MNTHR .AND. LWORK.GE.3*M+M*M+MAX( M, NRHS, N-2*M ) )
00552      $          THEN
00553 *
00554 *        Underdetermined case, M much less than N
00555 *
00556 *        Path 2a - underdetermined, with many more columns than rows
00557 *        and sufficient workspace for an efficient algorithm
00558 *
00559          LDWORK = M
00560          IF( LWORK.GE.3*M+M*LDA+MAX( M, NRHS, N-2*M ) )
00561      $      LDWORK = LDA
00562          ITAU = 1
00563          IWORK = M + 1
00564 *
00565 *        Compute A=L*Q
00566 *        (CWorkspace: need 2*M, prefer M+M*NB)
00567 *        (RWorkspace: none)
00568 *
00569          CALL ZGELQF( M, N, A, LDA, WORK( ITAU ), WORK( IWORK ),
00570      $                LWORK-IWORK+1, INFO )
00571          IL = IWORK
00572 *
00573 *        Copy L to WORK(IL), zeroing out above it
00574 *
00575          CALL ZLACPY( 'L', M, M, A, LDA, WORK( IL ), LDWORK )
00576          CALL ZLASET( 'U', M-1, M-1, CZERO, CZERO, WORK( IL+LDWORK ),
00577      $                LDWORK )
00578          IE = 1
00579          ITAUQ = IL + LDWORK*M
00580          ITAUP = ITAUQ + M
00581          IWORK = ITAUP + M
00582 *
00583 *        Bidiagonalize L in WORK(IL)
00584 *        (CWorkspace: need M*M+4*M, prefer M*M+3*M+2*M*NB)
00585 *        (RWorkspace: need M)
00586 *
00587          CALL ZGEBRD( M, M, WORK( IL ), LDWORK, S, RWORK( IE ),
00588      $                WORK( ITAUQ ), WORK( ITAUP ), WORK( IWORK ),
00589      $                LWORK-IWORK+1, INFO )
00590 *
00591 *        Multiply B by transpose of left bidiagonalizing vectors of L
00592 *        (CWorkspace: need M*M+3*M+NRHS, prefer M*M+3*M+NRHS*NB)
00593 *        (RWorkspace: none)
00594 *
00595          CALL ZUNMBR( 'Q', 'L', 'C', M, NRHS, M, WORK( IL ), LDWORK,
00596      $                WORK( ITAUQ ), B, LDB, WORK( IWORK ),
00597      $                LWORK-IWORK+1, INFO )
00598 *
00599 *        Generate right bidiagonalizing vectors of R in WORK(IL)
00600 *        (CWorkspace: need M*M+4*M-1, prefer M*M+3*M+(M-1)*NB)
00601 *        (RWorkspace: none)
00602 *
00603          CALL ZUNGBR( 'P', M, M, M, WORK( IL ), LDWORK, WORK( ITAUP ),
00604      $                WORK( IWORK ), LWORK-IWORK+1, INFO )
00605          IRWORK = IE + M
00606 *
00607 *        Perform bidiagonal QR iteration, computing right singular
00608 *        vectors of L in WORK(IL) and multiplying B by transpose of
00609 *        left singular vectors
00610 *        (CWorkspace: need M*M)
00611 *        (RWorkspace: need BDSPAC)
00612 *
00613          CALL ZBDSQR( 'U', M, M, 0, NRHS, S, RWORK( IE ), WORK( IL ),
00614      $                LDWORK, A, LDA, B, LDB, RWORK( IRWORK ), INFO )
00615          IF( INFO.NE.0 )
00616      $      GO TO 70
00617 *
00618 *        Multiply B by reciprocals of singular values
00619 *
00620          THR = MAX( RCOND*S( 1 ), SFMIN )
00621          IF( RCOND.LT.ZERO )
00622      $      THR = MAX( EPS*S( 1 ), SFMIN )
00623          RANK = 0
00624          DO 30 I = 1, M
00625             IF( S( I ).GT.THR ) THEN
00626                CALL ZDRSCL( NRHS, S( I ), B( I, 1 ), LDB )
00627                RANK = RANK + 1
00628             ELSE
00629                CALL ZLASET( 'F', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
00630             END IF
00631    30    CONTINUE
00632          IWORK = IL + M*LDWORK
00633 *
00634 *        Multiply B by right singular vectors of L in WORK(IL)
00635 *        (CWorkspace: need M*M+2*M, prefer M*M+M+M*NRHS)
00636 *        (RWorkspace: none)
00637 *
00638          IF( LWORK.GE.LDB*NRHS+IWORK-1 .AND. NRHS.GT.1 ) THEN
00639             CALL ZGEMM( 'C', 'N', M, NRHS, M, CONE, WORK( IL ), LDWORK,
00640      $                  B, LDB, CZERO, WORK( IWORK ), LDB )
00641             CALL ZLACPY( 'G', M, NRHS, WORK( IWORK ), LDB, B, LDB )
00642          ELSE IF( NRHS.GT.1 ) THEN
00643             CHUNK = ( LWORK-IWORK+1 ) / M
00644             DO 40 I = 1, NRHS, CHUNK
00645                BL = MIN( NRHS-I+1, CHUNK )
00646                CALL ZGEMM( 'C', 'N', M, BL, M, CONE, WORK( IL ), LDWORK,
00647      $                     B( 1, I ), LDB, CZERO, WORK( IWORK ), M )
00648                CALL ZLACPY( 'G', M, BL, WORK( IWORK ), M, B( 1, I ),
00649      $                      LDB )
00650    40       CONTINUE
00651          ELSE
00652             CALL ZGEMV( 'C', M, M, CONE, WORK( IL ), LDWORK, B( 1, 1 ),
00653      $                  1, CZERO, WORK( IWORK ), 1 )
00654             CALL ZCOPY( M, WORK( IWORK ), 1, B( 1, 1 ), 1 )
00655          END IF
00656 *
00657 *        Zero out below first M rows of B
00658 *
00659          CALL ZLASET( 'F', N-M, NRHS, CZERO, CZERO, B( M+1, 1 ), LDB )
00660          IWORK = ITAU + M
00661 *
00662 *        Multiply transpose(Q) by B
00663 *        (CWorkspace: need M+NRHS, prefer M+NHRS*NB)
00664 *        (RWorkspace: none)
00665 *
00666          CALL ZUNMLQ( 'L', 'C', N, NRHS, M, A, LDA, WORK( ITAU ), B,
00667      $                LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
00668 *
00669       ELSE
00670 *
00671 *        Path 2 - remaining underdetermined cases
00672 *
00673          IE = 1
00674          ITAUQ = 1
00675          ITAUP = ITAUQ + M
00676          IWORK = ITAUP + M
00677 *
00678 *        Bidiagonalize A
00679 *        (CWorkspace: need 3*M, prefer 2*M+(M+N)*NB)
00680 *        (RWorkspace: need N)
00681 *
00682          CALL ZGEBRD( M, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
00683      $                WORK( ITAUP ), WORK( IWORK ), LWORK-IWORK+1,
00684      $                INFO )
00685 *
00686 *        Multiply B by transpose of left bidiagonalizing vectors
00687 *        (CWorkspace: need 2*M+NRHS, prefer 2*M+NRHS*NB)
00688 *        (RWorkspace: none)
00689 *
00690          CALL ZUNMBR( 'Q', 'L', 'C', M, NRHS, N, A, LDA, WORK( ITAUQ ),
00691      $                B, LDB, WORK( IWORK ), LWORK-IWORK+1, INFO )
00692 *
00693 *        Generate right bidiagonalizing vectors in A
00694 *        (CWorkspace: need 3*M, prefer 2*M+M*NB)
00695 *        (RWorkspace: none)
00696 *
00697          CALL ZUNGBR( 'P', M, N, M, A, LDA, WORK( ITAUP ),
00698      $                WORK( IWORK ), LWORK-IWORK+1, INFO )
00699          IRWORK = IE + M
00700 *
00701 *        Perform bidiagonal QR iteration,
00702 *           computing right singular vectors of A in A and
00703 *           multiplying B by transpose of left singular vectors
00704 *        (CWorkspace: none)
00705 *        (RWorkspace: need BDSPAC)
00706 *
00707          CALL ZBDSQR( 'L', M, N, 0, NRHS, S, RWORK( IE ), A, LDA, DUM,
00708      $                1, B, LDB, RWORK( IRWORK ), INFO )
00709          IF( INFO.NE.0 )
00710      $      GO TO 70
00711 *
00712 *        Multiply B by reciprocals of singular values
00713 *
00714          THR = MAX( RCOND*S( 1 ), SFMIN )
00715          IF( RCOND.LT.ZERO )
00716      $      THR = MAX( EPS*S( 1 ), SFMIN )
00717          RANK = 0
00718          DO 50 I = 1, M
00719             IF( S( I ).GT.THR ) THEN
00720                CALL ZDRSCL( NRHS, S( I ), B( I, 1 ), LDB )
00721                RANK = RANK + 1
00722             ELSE
00723                CALL ZLASET( 'F', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
00724             END IF
00725    50    CONTINUE
00726 *
00727 *        Multiply B by right singular vectors of A
00728 *        (CWorkspace: need N, prefer N*NRHS)
00729 *        (RWorkspace: none)
00730 *
00731          IF( LWORK.GE.LDB*NRHS .AND. NRHS.GT.1 ) THEN
00732             CALL ZGEMM( 'C', 'N', N, NRHS, M, CONE, A, LDA, B, LDB,
00733      $                  CZERO, WORK, LDB )
00734             CALL ZLACPY( 'G', N, NRHS, WORK, LDB, B, LDB )
00735          ELSE IF( NRHS.GT.1 ) THEN
00736             CHUNK = LWORK / N
00737             DO 60 I = 1, NRHS, CHUNK
00738                BL = MIN( NRHS-I+1, CHUNK )
00739                CALL ZGEMM( 'C', 'N', N, BL, M, CONE, A, LDA, B( 1, I ),
00740      $                     LDB, CZERO, WORK, N )
00741                CALL ZLACPY( 'F', N, BL, WORK, N, B( 1, I ), LDB )
00742    60       CONTINUE
00743          ELSE
00744             CALL ZGEMV( 'C', M, N, CONE, A, LDA, B, 1, CZERO, WORK, 1 )
00745             CALL ZCOPY( N, WORK, 1, B, 1 )
00746          END IF
00747       END IF
00748 *
00749 *     Undo scaling
00750 *
00751       IF( IASCL.EQ.1 ) THEN
00752          CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
00753          CALL DLASCL( 'G', 0, 0, SMLNUM, ANRM, MINMN, 1, S, MINMN,
00754      $                INFO )
00755       ELSE IF( IASCL.EQ.2 ) THEN
00756          CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
00757          CALL DLASCL( 'G', 0, 0, BIGNUM, ANRM, MINMN, 1, S, MINMN,
00758      $                INFO )
00759       END IF
00760       IF( IBSCL.EQ.1 ) THEN
00761          CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
00762       ELSE IF( IBSCL.EQ.2 ) THEN
00763          CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
00764       END IF
00765    70 CONTINUE
00766       WORK( 1 ) = MAXWRK
00767       RETURN
00768 *
00769 *     End of ZGELSS
00770 *
00771       END
 All Files Functions