cggglm.f.html

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

</span>      EXTERNAL           CCOPY, CGEMV, <a name="CGGQRF.124"></a><a href="cggqrf.f.html#CGGQRF.1">CGGQRF</a>, <a name="CTRTRS.124"></a><a href="ctrtrs.f.html#CTRTRS.1">CTRTRS</a>, <a name="CUNMQR.124"></a><a href="cunmqr.f.html#CUNMQR.1">CUNMQR</a>, <a name="CUNMRQ.124"></a><a href="cunmrq.f.html#CUNMRQ.1">CUNMRQ</a>,
     $                   <a name="XERBLA.125"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      INTEGER            <a name="ILAENV.128"></a><a href="hfy-index.html#ILAENV">ILAENV</a>
      EXTERNAL           <a name="ILAENV.129"></a><a href="hfy-index.html#ILAENV">ILAENV</a> 
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          INT, 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 parameters
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 0
      NP = MIN( N, P )
      LQUERY = ( LWORK.EQ.-1 )
      IF( N.LT.0 ) THEN
         INFO = -1
      ELSE IF( M.LT.0 .OR. M.GT.N ) THEN
         INFO = -2
      ELSE IF( P.LT.0 .OR. P.LT.N-M ) THEN
         INFO = -3
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -5
      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
         INFO = -7
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Calculate workspace
</span><span class="comment">*</span><span class="comment">
</span>      IF( INFO.EQ.0) THEN
         IF( N.EQ.0 ) THEN
            LWKMIN = 1
            LWKOPT = 1
         ELSE
            NB1 = <a name="ILAENV.160"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="CGEQRF.160"></a><a href="cgeqrf.f.html#CGEQRF.1">CGEQRF</a>'</span>, <span class="string">' '</span>, N, M, -1, -1 )
            NB2 = <a name="ILAENV.161"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="CGERQF.161"></a><a href="cgerqf.f.html#CGERQF.1">CGERQF</a>'</span>, <span class="string">' '</span>, N, M, -1, -1 )
            NB3 = <a name="ILAENV.162"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="CUNMQR.162"></a><a href="cunmqr.f.html#CUNMQR.1">CUNMQR</a>'</span>, <span class="string">' '</span>, N, M, P, -1 )
            NB4 = <a name="ILAENV.163"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="CUNMRQ.163"></a><a href="cunmrq.f.html#CUNMRQ.1">CUNMRQ</a>'</span>, <span class="string">' '</span>, N, M, P, -1 )
            NB = MAX( NB1, NB2, NB3, NB4 )
            LWKMIN = M + N + P
            LWKOPT = M + NP + MAX( N, P )*NB
         END IF
         WORK( 1 ) = LWKOPT
<span class="comment">*</span><span class="comment">
</span>         IF( LWORK.LT.LWKMIN .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.176"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="CGGGLM.176"></a><a href="cggglm.f.html#CGGGLM.1">CGGGLM</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( N.EQ.0 )
     $   RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute the GQR factorization of matrices A and B:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">            Q'*A = ( R11 ) M,    Q'*B*Z' = ( T11   T12 ) M
</span><span class="comment">*</span><span class="comment">                   (  0  ) N-M             (  0    T22 ) N-M
</span><span class="comment">*</span><span class="comment">                      M                     M+P-N  N-M
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     where R11 and T22 are upper triangular, and Q and Z are
</span><span class="comment">*</span><span class="comment">     unitary.
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="CGGQRF.196"></a><a href="cggqrf.f.html#CGGQRF.1">CGGQRF</a>( N, M, P, A, LDA, WORK, B, LDB, WORK( M+1 ),
     $             WORK( M+NP+1 ), LWORK-M-NP, INFO )
      LOPT = WORK( M+NP+1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Update left-hand-side vector d = Q'*d = ( d1 ) M
</span><span class="comment">*</span><span class="comment">                                             ( d2 ) N-M
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="CUNMQR.203"></a><a href="cunmqr.f.html#CUNMQR.1">CUNMQR</a>( <span class="string">'Left'</span>, <span class="string">'Conjugate transpose'</span>, N, 1, M, A, LDA, WORK,
     $             D, MAX( 1, N ), WORK( M+NP+1 ), LWORK-M-NP, INFO )
      LOPT = MAX( LOPT, INT( WORK( M+NP+1 ) ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Solve T22*y2 = d2 for y2
</span><span class="comment">*</span><span class="comment">
</span>      IF( N.GT.M ) THEN
         CALL <a name="CTRTRS.210"></a><a href="ctrtrs.f.html#CTRTRS.1">CTRTRS</a>( <span class="string">'Upper'</span>, <span class="string">'No transpose'</span>, <span class="string">'Non unit'</span>, N-M, 1,
     $                B( M+1, M+P-N+1 ), LDB, D( M+1 ), N-M, INFO )
<span class="comment">*</span><span class="comment">
</span>         IF( INFO.GT.0 ) THEN
            INFO = 1
            RETURN
         END IF
<span class="comment">*</span><span class="comment">
</span>         CALL CCOPY( N-M, D( M+1 ), 1, Y( M+P-N+1 ), 1 )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Set y1 = 0
</span><span class="comment">*</span><span class="comment">
</span>      DO 10 I = 1, M + P - N
         Y( I ) = CZERO
   10 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Update d1 = d1 - T12*y2
</span><span class="comment">*</span><span class="comment">
</span>      CALL CGEMV( <span class="string">'No transpose'</span>, M, N-M, -CONE, B( 1, M+P-N+1 ), LDB,
     $            Y( M+P-N+1 ), 1, CONE, D, 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Solve triangular system: R11*x = d1
</span><span class="comment">*</span><span class="comment">
</span>      IF( M.GT.0 ) THEN
         CALL <a name="CTRTRS.235"></a><a href="ctrtrs.f.html#CTRTRS.1">CTRTRS</a>( <span class="string">'Upper'</span>, <span class="string">'No Transpose'</span>, <span class="string">'Non unit'</span>, M, 1, A, LDA,
     $                D, M, INFO )
<span class="comment">*</span><span class="comment">
</span>         IF( INFO.GT.0 ) THEN
            INFO = 2
            RETURN
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Copy D to X
</span><span class="comment">*</span><span class="comment">
</span>         CALL CCOPY( M, D, 1, X, 1 )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Backward transformation y = Z'*y
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="CUNMRQ.250"></a><a href="cunmrq.f.html#CUNMRQ.1">CUNMRQ</a>( <span class="string">'Left'</span>, <span class="string">'Conjugate transpose'</span>, P, 1, NP,
     $             B( MAX( 1, N-P+1 ), 1 ), LDB, WORK( M+1 ), Y,
     $             MAX( 1, P ), WORK( M+NP+1 ), LWORK-M-NP, INFO )
      WORK( 1 ) = M + NP + MAX( LOPT, INT( WORK( M+NP+1 ) ) )
<span class="comment">*</span><span class="comment">
</span>      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="CGGGLM.257"></a><a href="cggglm.f.html#CGGGLM.1">CGGGLM</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>