sgelsd.f.html

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

         MINWRK = MIN( MINWRK, MAXWRK )
         WORK( 1 ) = MAXWRK
         IWORK( 1 ) = LIWORK
<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.274"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="SGELSD.274"></a><a href="sgelsd.f.html#SGELSD.1">SGELSD</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.289"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'P'</span> )
      SFMIN = <a name="SLAMCH.290"></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.293"></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="SLANGE.297"></a><a href="slange.f.html#SLANGE.1">SLANGE</a>( <span class="string">'M'</span>, M, N, A, LDA, WORK )
      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="SLASCL.303"></a><a href="slascl.f.html#SLASCL.1">SLASCL</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="SLASCL.309"></a><a href="slascl.f.html#SLASCL.1">SLASCL</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="SLASET.315"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'F'</span>, MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
         CALL <a name="SLASET.316"></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="SLANGE.323"></a><a href="slange.f.html#SLANGE.1">SLANGE</a>( <span class="string">'M'</span>, M, NRHS, B, LDB, WORK )
      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="SLASCL.329"></a><a href="slascl.f.html#SLASCL.1">SLASCL</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="SLASCL.335"></a><a href="slascl.f.html#SLASCL.1">SLASCL</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 certain entries of B are zero.
</span><span class="comment">*</span><span class="comment">
</span>      IF( M.LT.N )
     $   CALL <a name="SLASET.342"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'F'</span>, N-M, NRHS, ZERO, ZERO, 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">           (Workspace: need 2*N, prefer N+N*NB)
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="SGEQRF.362"></a><a href="sgeqrf.f.html#SGEQRF.1">SGEQRF</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">           (Workspace: need N+NRHS, prefer N+NRHS*NB)
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="SORMQR.368"></a><a href="sormqr.f.html#SORMQR.1">SORMQR</a>( <span class="string">'L'</span>, <span class="string">'T'</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="SLASET.374"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'L'</span>, N-1, N-1, ZERO, ZERO, A( 2, 1 ), LDA )
            END IF
         END IF
<span class="comment">*</span><span class="comment">
</span>         IE = 1
         ITAUQ = IE + N
         ITAUP = ITAUQ + N
         NWORK = ITAUP + 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">        (Workspace: need 3*N+MM, prefer 3*N+(MM+N)*NB)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SGEBRD.386"></a><a href="sgebrd.f.html#SGEBRD.1">SGEBRD</a>( MM, N, A, LDA, S, WORK( 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">        (Workspace: need 3*N+NRHS, prefer 3*N+NRHS*NB)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SORMBR.393"></a><a href="sormbr.f.html#SORMBR.1">SORMBR</a>( <span class="string">'Q'</span>, <span class="string">'L'</span>, <span class="string">'T'</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="SLALSD.398"></a><a href="slalsd.f.html#SLALSD.1">SLALSD</a>( <span class="string">'U'</span>, SMLSIZ, N, NRHS, S, WORK( IE ), B, LDB,
     $                RCOND, RANK, WORK( NWORK ), IWORK, INFO )
         IF( INFO.NE.0 ) THEN
            GO TO 10
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Multiply B by right bidiagonalizing vectors of R.