cgelsd.f.html

来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 592 行 · 第 1/4 页

         END IF
         MINWRK = MIN( MINWRK, MAXWRK )
         WORK( 1 ) = MAXWRK
         IWORK( 1 ) = LIWORK
         RWORK( 1 ) = LRWORK
<span class="comment">*</span><span class="comment">
</span>         IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
            INFO = -12
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.290"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="CGELSD.290"></a><a href="cgelsd.f.html#CGELSD.1">CGELSD</a>'</span>, -INFO )
         RETURN
      ELSE IF( LQUERY ) THEN
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if possible.
</span><span class="comment">*</span><span class="comment">
</span>      IF( M.EQ.0 .OR. N.EQ.0 ) THEN
         RANK = 0
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Get machine parameters.
</span><span class="comment">*</span><span class="comment">
</span>      EPS = <a name="SLAMCH.305"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'P'</span> )
      SFMIN = <a name="SLAMCH.306"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'S'</span> )
      SMLNUM = SFMIN / EPS
      BIGNUM = ONE / SMLNUM
      CALL <a name="SLABAD.309"></a><a href="slabad.f.html#SLABAD.1">SLABAD</a>( SMLNUM, BIGNUM )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Scale A if max entry outside range [SMLNUM,BIGNUM].
</span><span class="comment">*</span><span class="comment">
</span>      ANRM = <a name="CLANGE.313"></a><a href="clange.f.html#CLANGE.1">CLANGE</a>( <span class="string">'M'</span>, M, N, A, LDA, RWORK )
      IASCL = 0
      IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Scale matrix norm up to SMLNUM
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="CLASCL.319"></a><a href="clascl.f.html#CLASCL.1">CLASCL</a>( <span class="string">'G'</span>, 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
         IASCL = 1
      ELSE IF( ANRM.GT.BIGNUM ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Scale matrix norm down to BIGNUM.
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="CLASCL.325"></a><a href="clascl.f.html#CLASCL.1">CLASCL</a>( <span class="string">'G'</span>, 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
         IASCL = 2
      ELSE IF( ANRM.EQ.ZERO ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Matrix all zero. Return zero solution.
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="CLASET.331"></a><a href="claset.f.html#CLASET.1">CLASET</a>( <span class="string">'F'</span>, MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
         CALL <a name="SLASET.332"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'F'</span>, MINMN, 1, ZERO, ZERO, S, 1 )
         RANK = 0
         GO TO 10
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Scale B if max entry outside range [SMLNUM,BIGNUM].
</span><span class="comment">*</span><span class="comment">
</span>      BNRM = <a name="CLANGE.339"></a><a href="clange.f.html#CLANGE.1">CLANGE</a>( <span class="string">'M'</span>, M, NRHS, B, LDB, RWORK )
      IBSCL = 0
      IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Scale matrix norm up to SMLNUM.
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="CLASCL.345"></a><a href="clascl.f.html#CLASCL.1">CLASCL</a>( <span class="string">'G'</span>, 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
         IBSCL = 1
      ELSE IF( BNRM.GT.BIGNUM ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Scale matrix norm down to BIGNUM.
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="CLASCL.351"></a><a href="clascl.f.html#CLASCL.1">CLASCL</a>( <span class="string">'G'</span>, 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
         IBSCL = 2
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If M &lt; N make sure B(M+1:N,:) = 0
</span><span class="comment">*</span><span class="comment">
</span>      IF( M.LT.N )
     $   CALL <a name="CLASET.358"></a><a href="claset.f.html#CLASET.1">CLASET</a>( <span class="string">'F'</span>, N-M, NRHS, CZERO, CZERO, B( M+1, 1 ), LDB )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Overdetermined case.
</span><span class="comment">*</span><span class="comment">
</span>      IF( M.GE.N ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Path 1 - overdetermined or exactly determined.
</span><span class="comment">*</span><span class="comment">
</span>         MM = M
         IF( M.GE.MNTHR ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Path 1a - overdetermined, with many more rows than columns
</span><span class="comment">*</span><span class="comment">
</span>            MM = N
            ITAU = 1
            NWORK = ITAU + N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute A=Q*R.
</span><span class="comment">*</span><span class="comment">           (RWorkspace: need N)
</span><span class="comment">*</span><span class="comment">           (CWorkspace: need N, prefer N*NB)
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="CGEQRF.379"></a><a href="cgeqrf.f.html#CGEQRF.1">CGEQRF</a>( M, N, A, LDA, WORK( ITAU ), WORK( NWORK ),
     $                   LWORK-NWORK+1, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Multiply B by transpose(Q).
</span><span class="comment">*</span><span class="comment">           (RWorkspace: need N)
</span><span class="comment">*</span><span class="comment">           (CWorkspace: need NRHS, prefer NRHS*NB)
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="CUNMQR.386"></a><a href="cunmqr.f.html#CUNMQR.1">CUNMQR</a>( <span class="string">'L'</span>, <span class="string">'C'</span>, M, NRHS, N, A, LDA, WORK( ITAU ), B,
     $                   LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Zero out below R.
</span><span class="comment">*</span><span class="comment">
</span>            IF( N.GT.1 ) THEN
               CALL <a name="CLASET.392"></a><a href="claset.f.html#CLASET.1">CLASET</a>( <span class="string">'L'</span>, N-1, N-1, CZERO, CZERO, A( 2, 1 ),
     $                      LDA )
            END IF
         END IF
<span class="comment">*</span><span class="comment">
</span>         ITAUQ = 1
         ITAUP = ITAUQ + N
         NWORK = ITAUP + N
         IE = 1
         NRWORK = IE + N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Bidiagonalize R in A.
</span><span class="comment">*</span><span class="comment">        (RWorkspace: need N)
</span><span class="comment">*</span><span class="comment">        (CWorkspace: need 2*N+MM, prefer 2*N+(MM+N)*NB)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="CGEBRD.407"></a><a href="cgebrd.f.html#CGEBRD.1">CGEBRD</a>( MM, N, A, LDA, S, RWORK( IE ), WORK( ITAUQ ),
     $                WORK( ITAUP ), WORK( NWORK ), LWORK-NWORK+1,
     $                INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Multiply B by transpose of left bidiagonalizing vectors of R.
</span><span class="comment">*</span><span class="comment">        (CWorkspace: need 2*N+NRHS, prefer 2*N+NRHS*NB)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="CUNMBR.414"></a><a href="cunmbr.f.html#CUNMBR.1">CUNMBR</a>( <span class="string">'Q'</span>, <span class="string">'L'</span>, <span class="string">'C'</span>, MM, NRHS, N, A, LDA, WORK( ITAUQ ),
     $                B, LDB, WORK( NWORK ), LWORK-NWORK+1, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Solve the bidiagonal least squares problem.
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="CLALSD.419"></a><a href="clalsd.f.html#CLALSD.1">CLALSD</a>( <span class="string">'U'</span>, SMLSIZ, N, NRHS, S, RWORK( IE ), B, LDB,
     $                RCOND, RANK, WORK( NWORK ), RWORK( NRWORK ),
     $                IWORK, INFO )
         IF( INFO.NE.0 ) THEN
            GO TO 10
         END IF
<span class="comment">*</span><span class="comment">