zgels.f.html

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

HTML
448
字号
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.133"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      INTEGER            <a name="ILAENV.134"></a><a href="hfy-index.html#ILAENV">ILAENV</a>
      DOUBLE PRECISION   <a name="DLAMCH.135"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="ZLANGE.135"></a><a href="zlange.f.html#ZLANGE.1">ZLANGE</a>
      EXTERNAL           <a name="LSAME.136"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, <a name="ILAENV.136"></a><a href="hfy-index.html#ILAENV">ILAENV</a>, <a name="DLAMCH.136"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="ZLANGE.136"></a><a href="zlange.f.html#ZLANGE.1">ZLANGE</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           <a name="DLABAD.139"></a><a href="dlabad.f.html#DLABAD.1">DLABAD</a>, <a name="XERBLA.139"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>, <a name="ZGELQF.139"></a><a href="zgelqf.f.html#ZGELQF.1">ZGELQF</a>, <a name="ZGEQRF.139"></a><a href="zgeqrf.f.html#ZGEQRF.1">ZGEQRF</a>, <a name="ZLASCL.139"></a><a href="zlascl.f.html#ZLASCL.1">ZLASCL</a>, <a name="ZLASET.139"></a><a href="zlaset.f.html#ZLASET.1">ZLASET</a>,
     $                   <a name="ZTRTRS.140"></a><a href="ztrtrs.f.html#ZTRTRS.1">ZTRTRS</a>, <a name="ZUNMLQ.140"></a><a href="zunmlq.f.html#ZUNMLQ.1">ZUNMLQ</a>, <a name="ZUNMQR.140"></a><a href="zunmqr.f.html#ZUNMQR.1">ZUNMQR</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          DBLE, MAX, MIN
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Test the input arguments.
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 0
      MN = MIN( M, N )
      LQUERY = ( LWORK.EQ.-1 )
      IF( .NOT.( <a name="LSAME.152"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( TRANS, <span class="string">'N'</span> ) .OR. <a name="LSAME.152"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( TRANS, <span class="string">'C'</span> ) ) ) THEN
         INFO = -1
      ELSE IF( M.LT.0 ) THEN
         INFO = -2
      ELSE IF( N.LT.0 ) THEN
         INFO = -3
      ELSE IF( NRHS.LT.0 ) THEN
         INFO = -4
      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
         INFO = -6
      ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
         INFO = -8
      ELSE IF( LWORK.LT.MAX( 1, MN+MAX( MN, NRHS ) ) .AND. .NOT.LQUERY )
     $          THEN
         INFO = -10
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Figure out optimal block size
</span><span class="comment">*</span><span class="comment">
</span>      IF( INFO.EQ.0 .OR. INFO.EQ.-10 ) THEN
<span class="comment">*</span><span class="comment">
</span>         TPSD = .TRUE.
         IF( <a name="LSAME.174"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( TRANS, <span class="string">'N'</span> ) )
     $      TPSD = .FALSE.
<span class="comment">*</span><span class="comment">
</span>         IF( M.GE.N ) THEN
            NB = <a name="ILAENV.178"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="ZGEQRF.178"></a><a href="zgeqrf.f.html#ZGEQRF.1">ZGEQRF</a>'</span>, <span class="string">' '</span>, M, N, -1, -1 )
            IF( TPSD ) THEN
               NB = MAX( NB, <a name="ILAENV.180"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="ZUNMQR.180"></a><a href="zunmqr.f.html#ZUNMQR.1">ZUNMQR</a>'</span>, <span class="string">'LN'</span>, M, NRHS, N,
     $              -1 ) )
            ELSE
               NB = MAX( NB, <a name="ILAENV.183"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="ZUNMQR.183"></a><a href="zunmqr.f.html#ZUNMQR.1">ZUNMQR</a>'</span>, <span class="string">'LC'</span>, M, NRHS, N,
     $              -1 ) )
            END IF
         ELSE
            NB = <a name="ILAENV.187"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="ZGELQF.187"></a><a href="zgelqf.f.html#ZGELQF.1">ZGELQF</a>'</span>, <span class="string">' '</span>, M, N, -1, -1 )
            IF( TPSD ) THEN
               NB = MAX( NB, <a name="ILAENV.189"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="ZUNMLQ.189"></a><a href="zunmlq.f.html#ZUNMLQ.1">ZUNMLQ</a>'</span>, <span class="string">'LC'</span>, N, NRHS, M,
     $              -1 ) )
            ELSE
               NB = MAX( NB, <a name="ILAENV.192"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="ZUNMLQ.192"></a><a href="zunmlq.f.html#ZUNMLQ.1">ZUNMLQ</a>'</span>, <span class="string">'LN'</span>, N, NRHS, M,
     $              -1 ) )
            END IF
         END IF
<span class="comment">*</span><span class="comment">
</span>         WSIZE = MAX( 1, MN+MAX( MN, NRHS )*NB )
         WORK( 1 ) = DBLE( WSIZE )
<span class="comment">*</span><span class="comment">
</span>      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.203"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="ZGELS.203"></a><a href="zgels.f.html#ZGELS.1">ZGELS</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( MIN( M, N, NRHS ).EQ.0 ) THEN
         CALL <a name="ZLASET.212"></a><a href="zlaset.f.html#ZLASET.1">ZLASET</a>( <span class="string">'Full'</span>, MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
         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>      SMLNUM = <a name="DLAMCH.218"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'S'</span> ) / <a name="DLAMCH.218"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'P'</span> )
      BIGNUM = ONE / SMLNUM
      CALL <a name="DLABAD.220"></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, B if max element outside range [SMLNUM,BIGNUM]
</span><span class="comment">*</span><span class="comment">
</span>      ANRM = <a name="ZLANGE.224"></a><a href="zlange.f.html#ZLANGE.1">ZLANGE</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="ZLASCL.230"></a><a href="zlascl.f.html#ZLASCL.1">ZLASCL</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="ZLASCL.236"></a><a href="zlascl.f.html#ZLASCL.1">ZLASCL</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="ZLASET.242"></a><a href="zlaset.f.html#ZLASET.1">ZLASET</a>( <span class="string">'F'</span>, MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
         GO TO 50
      END IF
<span class="comment">*</span><span class="comment">
</span>      BROW = M
      IF( TPSD )
     $   BROW = N
      BNRM = <a name="ZLANGE.249"></a><a href="zlange.f.html#ZLANGE.1">ZLANGE</a>( <span class="string">'M'</span>, BROW, 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="ZLASCL.255"></a><a href="zlascl.f.html#ZLASCL.1">ZLASCL</a>( <span class="string">'G'</span>, 0, 0, BNRM, SMLNUM, BROW, 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="ZLASCL.262"></a><a href="zlascl.f.html#ZLASCL.1">ZLASCL</a>( <span class="string">'G'</span>, 0, 0, BNRM, BIGNUM, BROW, NRHS, B, LDB,
     $                INFO )
         IBSCL = 2
      END IF
<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">        compute QR factorization of A
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="ZGEQRF.271"></a><a href="zgeqrf.f.html#ZGEQRF.1">ZGEQRF</a>( M, N, A, LDA, WORK( 1 ), WORK( MN+1 ), LWORK-MN,
     $                INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        workspace at least N, optimally N*NB
</span><span class="comment">*</span><span class="comment">
</span>         IF( .NOT.TPSD ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Least-Squares Problem min || A * X - B ||
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           B(1:M,1:NRHS) := Q' * B(1:M,1:NRHS)
</span><span class="comment">*</span><span class="comment">

⌨️ 快捷键说明

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