cgglse.f.html

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

</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
      MN = MIN( M, N )
      LQUERY = ( LWORK.EQ.-1 )
      IF( M.LT.0 ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -2
      ELSE IF( P.LT.0 .OR. P.GT.N .OR. P.LT.N-M ) THEN
         INFO = -3
      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
         INFO = -5
      ELSE IF( LDB.LT.MAX( 1, P ) ) 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.156"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="CGEQRF.156"></a><a href="cgeqrf.f.html#CGEQRF.1">CGEQRF</a>'</span>, <span class="string">' '</span>, M, N, -1, -1 )
            NB2 = <a name="ILAENV.157"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="CGERQF.157"></a><a href="cgerqf.f.html#CGERQF.1">CGERQF</a>'</span>, <span class="string">' '</span>, M, N, -1, -1 )
            NB3 = <a name="ILAENV.158"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="CUNMQR.158"></a><a href="cunmqr.f.html#CUNMQR.1">CUNMQR</a>'</span>, <span class="string">' '</span>, M, N, P, -1 )
            NB4 = <a name="ILAENV.159"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 1, <span class="string">'<a name="CUNMRQ.159"></a><a href="cunmrq.f.html#CUNMRQ.1">CUNMRQ</a>'</span>, <span class="string">' '</span>, M, N, P, -1 )
            NB = MAX( NB1, NB2, NB3, NB4 )
            LWKMIN = M + N + P
            LWKOPT = P + MN + MAX( M, N )*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.172"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="CGGLSE.172"></a><a href="cgglse.f.html#CGGLSE.1">CGGLSE</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 GRQ factorization of matrices B and A:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">            B*Q' = (  0  T12 ) P   Z'*A*Q' = ( R11 R12 ) N-P
</span><span class="comment">*</span><span class="comment">                     N-P  P                  (  0  R22 ) M+P-N
</span><span class="comment">*</span><span class="comment">                                               N-P  P
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     where T12 and R11 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="CGGRQF.192"></a><a href="cggrqf.f.html#CGGRQF.1">CGGRQF</a>( P, M, N, B, LDB, WORK, A, LDA, WORK( P+1 ),
     $             WORK( P+MN+1 ), LWORK-P-MN, INFO )
      LOPT = WORK( P+MN+1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Update c = Z'*c = ( c1 ) N-P
</span><span class="comment">*</span><span class="comment">                       ( c2 ) M+P-N
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="CUNMQR.199"></a><a href="cunmqr.f.html#CUNMQR.1">CUNMQR</a>( <span class="string">'Left'</span>, <span class="string">'Conjugate Transpose'</span>, M, 1, MN, A, LDA,
     $             WORK( P+1 ), C, MAX( 1, M ), WORK( P+MN+1 ),
     $             LWORK-P-MN, INFO )
      LOPT = MAX( LOPT, INT( WORK( P+MN+1 ) ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Solve T12*x2 = d for x2
</span><span class="comment">*</span><span class="comment">
</span>      IF( P.GT.0 ) THEN
         CALL <a name="CTRTRS.207"></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>, P, 1,
     $                B( 1, N-P+1 ), LDB, D, P, 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><span class="comment">*</span><span class="comment">        Put the solution in X
</span><span class="comment">*</span><span class="comment">
</span>      CALL CCOPY( P, D, 1, X( N-P+1 ), 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Update c1
</span><span class="comment">*</span><span class="comment">
</span>         CALL CGEMV( <span class="string">'No transpose'</span>, N-P, P, -CONE, A( 1, N-P+1 ), LDA,
     $               D, 1, CONE, C, 1 )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Solve R11*x1 = c1 for x1
</span><span class="comment">*</span><span class="comment">
</span>      IF( N.GT.P ) THEN
         CALL <a name="CTRTRS.228"></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-P, 1,
     $                A, LDA, C, N-P, 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">        Put the solutions in X
</span><span class="comment">*</span><span class="comment">
</span>         CALL CCOPY( N-P, C, 1, X, 1 )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute the residual vector:
</span><span class="comment">*</span><span class="comment">
</span>      IF( M.LT.N ) THEN
         NR = M + P - N
         IF( NR.GT.0 )
     $      CALL CGEMV( <span class="string">'No transpose'</span>, NR, N-M, -CONE, A( N-P+1, M+1 ),
     $                  LDA, D( NR+1 ), 1, CONE, C( N-P+1 ), 1 )
      ELSE
         NR = P
      END IF
      IF( NR.GT.0 ) THEN
         CALL CTRMV( <span class="string">'Upper'</span>, <span class="string">'No transpose'</span>, <span class="string">'Non unit'</span>, NR,
     $               A( N-P+1, N-P+1 ), LDA, D, 1 )
         CALL CAXPY( NR, -CONE, D, 1, C( N-P+1 ), 1 )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Backward transformation x = Q'*x
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="CUNMRQ.259"></a><a href="cunmrq.f.html#CUNMRQ.1">CUNMRQ</a>( <span class="string">'Left'</span>, <span class="string">'Conjugate Transpose'</span>, N, 1, P, B, LDB,
     $             WORK( 1 ), X, N, WORK( P+MN+1 ), LWORK-P-MN, INFO )
      WORK( 1 ) = P + MN + MAX( LOPT, INT( WORK( P+MN+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="CGGLSE.265"></a><a href="cgglse.f.html#CGGLSE.1">CGGLSE</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>