dgelsd.f.html

来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 553 行 · 第 1/3 页
HTML
553 行
      LQUERY = ( LWORK.EQ.-1 )
      IF( M.LT.0 ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -2
      ELSE IF( NRHS.LT.0 ) THEN
         INFO = -3
      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
         INFO = -5
      ELSE IF( LDB.LT.MAX( 1, MAXMN ) ) THEN
         INFO = -7
      END IF
<span class="comment">*</span><span class="comment">
</span>      SMLSIZ = <a name="ILAENV.180"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 9, <span class="string">'<a name="DGELSD.180"></a><a href="dgelsd.f.html#DGELSD.1">DGELSD</a>'</span>, <span class="string">' '</span>, 0, 0, 0, 0 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute workspace.
</span><span class="comment">*</span><span class="comment">     (Note: Comments in the code beginning &quot;Workspace:&quot; describe the
</span><span class="comment">*</span><span class="comment">     minimal amount of workspace needed at that point in the code,
</span><span class="comment">*</span><span class="comment">     as well as the preferred amount for good performance.
</span><span class="comment">*</span><span class="comment">     NB refers to the optimal block size for the immediately
</span><span class="comment">*</span><span class="comment">     following subroutine, as returned by <a name="ILAENV.187"></a><a href="hfy-index.html#ILAENV">ILAENV</a>.)
</span><span class="comment">*</span><span class="comment">
</span>      MINWRK = 1
      MINMN = MAX( 1, MINMN )
      NLVL = MAX( INT( LOG( DBLE( MINMN ) / DBLE( SMLSIZ+1 ) ) /
     $       LOG( TWO ) ) + 1, 0 )
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.EQ.0 ) THEN
         MAXWRK = 0
         MM = M
         IF( M.GE.N .AND. 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
            MAXWRK = MAX( MAXWRK, N+N*<a name="ILAENV.202"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="DGEQRF.202"></a><a href="dgeqrf.f.html#DGEQRF.1">DGEQRF</a>'</span>, <span class="string">' '</span>, M, N,
     $               -1, -1 ) )
            MAXWRK = MAX( MAXWRK, N+NRHS*
     $               <a name="ILAENV.205"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="DORMQR.205"></a><a href="dormqr.f.html#DORMQR.1">DORMQR</a>'</span>, <span class="string">'LT'</span>, M, NRHS, N, -1 ) )
         END IF
         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>            MAXWRK = MAX( MAXWRK, 3*N+( MM+N )*
     $               <a name="ILAENV.212"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="DGEBRD.212"></a><a href="dgebrd.f.html#DGEBRD.1">DGEBRD</a>'</span>, <span class="string">' '</span>, MM, N, -1, -1 ) )
            MAXWRK = MAX( MAXWRK, 3*N+NRHS*
     $               <a name="ILAENV.214"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="DORMBR.214"></a><a href="dormbr.f.html#DORMBR.1">DORMBR</a>'</span>, <span class="string">'QLT'</span>, MM, NRHS, N, -1 ) )
            MAXWRK = MAX( MAXWRK, 3*N+( N-1 )*
     $               <a name="ILAENV.216"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="DORMBR.216"></a><a href="dormbr.f.html#DORMBR.1">DORMBR</a>'</span>, <span class="string">'PLN'</span>, N, NRHS, N, -1 ) )
            WLALSD = 9*N+2*N*SMLSIZ+8*N*NLVL+N*NRHS+(SMLSIZ+1)**2
            MAXWRK = MAX( MAXWRK, 3*N+WLALSD )
            MINWRK = MAX( 3*N+MM, 3*N+NRHS, 3*N+WLALSD )
         END IF
         IF( N.GT.M ) THEN
            WLALSD = 9*M+2*M*SMLSIZ+8*M*NLVL+M*NRHS+(SMLSIZ+1)**2
            IF( N.GE.MNTHR ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Path 2a - underdetermined, with many more columns
</span><span class="comment">*</span><span class="comment">              than rows.
</span><span class="comment">*</span><span class="comment">
</span>               MAXWRK = M + M*<a name="ILAENV.228"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="DGELQF.228"></a><a href="dgelqf.f.html#DGELQF.1">DGELQF</a>'</span>, <span class="string">' '</span>, M, N, -1, -1 )
               MAXWRK = MAX( MAXWRK, M*M+4*M+2*M*
     $                  <a name="ILAENV.230"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="DGEBRD.230"></a><a href="dgebrd.f.html#DGEBRD.1">DGEBRD</a>'</span>, <span class="string">' '</span>, M, M, -1, -1 ) )
               MAXWRK = MAX( MAXWRK, M*M+4*M+NRHS*
     $                  <a name="ILAENV.232"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="DORMBR.232"></a><a href="dormbr.f.html#DORMBR.1">DORMBR</a>'</span>, <span class="string">'QLT'</span>, M, NRHS, M, -1 ) )
               MAXWRK = MAX( MAXWRK, M*M+4*M+( M-1 )*
     $                  <a name="ILAENV.234"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="DORMBR.234"></a><a href="dormbr.f.html#DORMBR.1">DORMBR</a>'</span>, <span class="string">'PLN'</span>, M, NRHS, M, -1 ) )
               IF( NRHS.GT.1 ) THEN
                  MAXWRK = MAX( MAXWRK, M*M+M+M*NRHS )
               ELSE
                  MAXWRK = MAX( MAXWRK, M*M+2*M )
               END IF
               MAXWRK = MAX( MAXWRK, M+NRHS*
     $                  <a name="ILAENV.241"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="DORMLQ.241"></a><a href="dormlq.f.html#DORMLQ.1">DORMLQ</a>'</span>, <span class="string">'LT'</span>, N, NRHS, M, -1 ) )
               MAXWRK = MAX( MAXWRK, M*M+4*M+WLALSD )
            ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Path 2 - remaining underdetermined cases.
</span><span class="comment">*</span><span class="comment">
</span>               MAXWRK = 3*M + ( N+M )*<a name="ILAENV.247"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="DGEBRD.247"></a><a href="dgebrd.f.html#DGEBRD.1">DGEBRD</a>'</span>, <span class="string">' '</span>, M, N,
     $                  -1, -1 )
               MAXWRK = MAX( MAXWRK, 3*M+NRHS*
     $                  <a name="ILAENV.250"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="DORMBR.250"></a><a href="dormbr.f.html#DORMBR.1">DORMBR</a>'</span>, <span class="string">'QLT'</span>, M, NRHS, N, -1 ) )
               MAXWRK = MAX( MAXWRK, 3*M+M*
     $                  <a name="ILAENV.252"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="DORMBR.252"></a><a href="dormbr.f.html#DORMBR.1">DORMBR</a>'</span>, <span class="string">'PLN'</span>, N, NRHS, M, -1 ) )
               MAXWRK = MAX( MAXWRK, 3*M+WLALSD )
            END IF
            MINWRK = MAX( 3*M+NRHS, 3*M+M, 3*M+WLALSD )
         END IF
         MINWRK = MIN( MINWRK, MAXWRK )
         WORK( 1 ) = MAXWRK
         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.265"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="DGELSD.265"></a><a href="dgelsd.f.html#DGELSD.1">DGELSD</a>'</span>, -INFO )
         RETURN
      ELSE IF( LQUERY ) THEN
         GO TO 10
      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="DLAMCH.280"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'P'</span> )
      SFMIN = <a name="DLAMCH.281"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'S'</span> )
      SMLNUM = SFMIN / EPS
      BIGNUM = ONE / SMLNUM
      CALL <a name="DLABAD.284"></a><a href="dlabad.f.html#DLABAD.1">DLABAD</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="DLANGE.288"></a><a href="dlange.f.html#DLANGE.1">DLANGE</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="DLASCL.294"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</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="DLASCL.300"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</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="DLASET.306"></a><a href="dlaset.f.html#DLASET.1">DLASET</a>( <span class="string">'F'</span>, MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
         CALL <a name="DLASET.307"></a><a href="dlaset.f.html#DLASET.1">DLASET</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="DLANGE.314"></a><a href="dlange.f.html#DLANGE.1">DLANGE</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="DLASCL.320"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</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="DLASCL.326"></a><a href="dlascl.f.html#DLASCL.1">DLASCL</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="DLASET.333"></a><a href="dlaset.f.html#DLASET.1">DLASET</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)
dgelsd.f.html - 源码说明

本页面展示了「famous linear algebra library (LAPACK) ports to windows」中的 dgelsd.f.html 源码文件，采用 HTML 编程语言编写，共 553 行代码。您可以在线阅读完整代码内容，也可以返回资源详情页下载完整源码包进行本地学习和开发。
虫虫下载站收录了大量与algebra相关的技术资源，包括源代码、技术文档、电路图等，是电子工程师和嵌入式开发者的专业学习平台。
⌨️ 快捷键说明

复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?