dgglse.f.html

来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 291 行 · 第 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="DGEQRF.156"></a><a href="dgeqrf.f.html#DGEQRF.1">DGEQRF</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="DGERQF.157"></a><a href="dgerqf.f.html#DGERQF.1">DGERQF</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="DORMQR.158"></a><a href="dormqr.f.html#DORMQR.1">DORMQR</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="DORMRQ.159"></a><a href="dormrq.f.html#DORMRQ.1">DORMRQ</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="DGGLSE.172"></a><a href="dgglse.f.html#DGGLSE.1">DGGLSE</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">     orthogonal.
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="DGGRQF.192"></a><a href="dggrqf.f.html#DGGRQF.1">DGGRQF</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="DORMQR.199"></a><a href="dormqr.f.html#DORMQR.1">DORMQR</a>( <span class="string">'Left'</span>, <span class="string">'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="DTRTRS.206"></a><a href="dtrtrs.f.html#DTRTRS.1">DTRTRS</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 DCOPY( 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 DGEMV( <span class="string">'No transpose'</span>, N-P, P, -ONE, A( 1, N-P+1 ), LDA,
     $               D, 1, ONE, 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="DTRTRS.227"></a><a href="dtrtrs.f.html#DTRTRS.1">DTRTRS</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 DCOPY( 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 DGEMV( <span class="string">'No transpose'</span>, NR, N-M, -ONE, A( N-P+1, M+1 ),
     $                  LDA, D( NR+1 ), 1, ONE, C( N-P+1 ), 1 )
      ELSE
         NR = P
      END IF
      IF( NR.GT.0 ) THEN
         CALL DTRMV( <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 DAXPY( NR, -ONE, 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="DORMRQ.258"></a><a href="dormrq.f.html#DORMRQ.1">DORMRQ</a>( <span class="string">'Left'</span>, <span class="string">'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="DGGLSE.264"></a><a href="dgglse.f.html#DGGLSE.1">DGGLSE</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>