sgglse.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="SGEQRF.156"></a><a href="sgeqrf.f.html#SGEQRF.1">SGEQRF</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="SGERQF.157"></a><a href="sgerqf.f.html#SGERQF.1">SGERQF</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="SORMQR.158"></a><a href="sormqr.f.html#SORMQR.1">SORMQR</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="SORMRQ.159"></a><a href="sormrq.f.html#SORMRQ.1">SORMRQ</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="SGGLSE.172"></a><a href="sgglse.f.html#SGGLSE.1">SGGLSE</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="SGGRQF.192"></a><a href="sggrqf.f.html#SGGRQF.1">SGGRQF</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="SORMQR.199"></a><a href="sormqr.f.html#SORMQR.1">SORMQR</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="STRTRS.206"></a><a href="strtrs.f.html#STRTRS.1">STRTRS</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 SCOPY( 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 SGEMV( <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="STRTRS.227"></a><a href="strtrs.f.html#STRTRS.1">STRTRS</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 solution in X
</span><span class="comment">*</span><span class="comment">
</span>         CALL SCOPY( 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 SGEMV( <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 STRMV( <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 SAXPY( 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="SORMRQ.258"></a><a href="sormrq.f.html#SORMRQ.1">SORMRQ</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="SGGLSE.264"></a><a href="sgglse.f.html#SGGLSE.1">SGGLSE</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>