stgsja.f.html

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

</span>         DO 20 I = 1, L - 1
            DO 10 J = I + 1, L
<span class="comment">*</span><span class="comment">
</span>               A1 = ZERO
               A2 = ZERO
               A3 = ZERO
               IF( K+I.LE.M )
     $            A1 = A( K+I, N-L+I )
               IF( K+J.LE.M )
     $            A3 = A( K+J, N-L+J )
<span class="comment">*</span><span class="comment">
</span>               B1 = B( I, N-L+I )
               B3 = B( J, N-L+J )
<span class="comment">*</span><span class="comment">
</span>               IF( UPPER ) THEN
                  IF( K+I.LE.M )
     $               A2 = A( K+I, N-L+J )
                  B2 = B( I, N-L+J )
               ELSE
                  IF( K+J.LE.M )
     $               A2 = A( K+J, N-L+I )
                  B2 = B( J, N-L+I )
               END IF
<span class="comment">*</span><span class="comment">
</span>               CALL <a name="SLAGS2.366"></a><a href="slags2.f.html#SLAGS2.1">SLAGS2</a>( UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU,
     $                      CSV, SNV, CSQ, SNQ )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Update (K+I)-th and (K+J)-th rows of matrix A: U'*A
</span><span class="comment">*</span><span class="comment">
</span>               IF( K+J.LE.M )
     $            CALL SROT( L, A( K+J, N-L+1 ), LDA, A( K+I, N-L+1 ),
     $                       LDA, CSU, SNU )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Update I-th and J-th rows of matrix B: V'*B
</span><span class="comment">*</span><span class="comment">
</span>               CALL SROT( L, B( J, N-L+1 ), LDB, B( I, N-L+1 ), LDB,
     $                    CSV, SNV )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Update (N-L+I)-th and (N-L+J)-th columns of matrices
</span><span class="comment">*</span><span class="comment">              A and B: A*Q and B*Q
</span><span class="comment">*</span><span class="comment">
</span>               CALL SROT( MIN( K+L, M ), A( 1, N-L+J ), 1,
     $                    A( 1, N-L+I ), 1, CSQ, SNQ )
<span class="comment">*</span><span class="comment">
</span>               CALL SROT( L, B( 1, N-L+J ), 1, B( 1, N-L+I ), 1, CSQ,
     $                    SNQ )
<span class="comment">*</span><span class="comment">
</span>               IF( UPPER ) THEN
                  IF( K+I.LE.M )
     $               A( K+I, N-L+J ) = ZERO
                  B( I, N-L+J ) = ZERO
               ELSE
                  IF( K+J.LE.M )
     $               A( K+J, N-L+I ) = ZERO
                  B( J, N-L+I ) = ZERO
               END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Update orthogonal matrices U, V, Q, if desired.
</span><span class="comment">*</span><span class="comment">
</span>               IF( WANTU .AND. K+J.LE.M )
     $            CALL SROT( M, U( 1, K+J ), 1, U( 1, K+I ), 1, CSU,
     $                       SNU )
<span class="comment">*</span><span class="comment">
</span>               IF( WANTV )
     $            CALL SROT( P, V( 1, J ), 1, V( 1, I ), 1, CSV, SNV )
<span class="comment">*</span><span class="comment">
</span>               IF( WANTQ )
     $            CALL SROT( N, Q( 1, N-L+J ), 1, Q( 1, N-L+I ), 1, CSQ,
     $                       SNQ )
<span class="comment">*</span><span class="comment">
</span>   10       CONTINUE
   20    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         IF( .NOT.UPPER ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           The matrices A13 and B13 were lower triangular at the start
</span><span class="comment">*</span><span class="comment">           of the cycle, and are now upper triangular.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Convergence test: test the parallelism of the corresponding
</span><span class="comment">*</span><span class="comment">           rows of A and B.
</span><span class="comment">*</span><span class="comment">
</span>            ERROR = ZERO
            DO 30 I = 1, MIN( L, M-K )
               CALL SCOPY( L-I+1, A( K+I, N-L+I ), LDA, WORK, 1 )
               CALL SCOPY( L-I+1, B( I, N-L+I ), LDB, WORK( L+1 ), 1 )
               CALL <a name="SLAPLL.427"></a><a href="slapll.f.html#SLAPLL.1">SLAPLL</a>( L-I+1, WORK, 1, WORK( L+1 ), 1, SSMIN )
               ERROR = MAX( ERROR, SSMIN )
   30       CONTINUE
<span class="comment">*</span><span class="comment">
</span>            IF( ABS( ERROR ).LE.MIN( TOLA, TOLB ) )
     $         GO TO 50
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        End of cycle loop
</span><span class="comment">*</span><span class="comment">
</span>   40 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     The algorithm has not converged after MAXIT cycles.
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 1
      GO TO 100
<span class="comment">*</span><span class="comment">
</span>   50 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If ERROR &lt;= MIN(TOLA,TOLB), then the algorithm has converged.
</span><span class="comment">*</span><span class="comment">     Compute the generalized singular value pairs (ALPHA, BETA), and
</span><span class="comment">*</span><span class="comment">     set the triangular matrix R to array A.
</span><span class="comment">*</span><span class="comment">
</span>      DO 60 I = 1, K
         ALPHA( I ) = ONE
         BETA( I ) = ZERO
   60 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      DO 70 I = 1, MIN( L, M-K )
<span class="comment">*</span><span class="comment">
</span>         A1 = A( K+I, N-L+I )
         B1 = B( I, N-L+I )
<span class="comment">*</span><span class="comment">
</span>         IF( A1.NE.ZERO ) THEN
            GAMMA = B1 / A1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           change sign if necessary
</span><span class="comment">*</span><span class="comment">
</span>            IF( GAMMA.LT.ZERO ) THEN
               CALL SSCAL( L-I+1, -ONE, B( I, N-L+I ), LDB )
               IF( WANTV )
     $            CALL SSCAL( P, -ONE, V( 1, I ), 1 )
            END IF
<span class="comment">*</span><span class="comment">
</span>            CALL <a name="SLARTG.471"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>( ABS( GAMMA ), ONE, BETA( K+I ), ALPHA( K+I ),
     $                   RWK )
<span class="comment">*</span><span class="comment">
</span>            IF( ALPHA( K+I ).GE.BETA( K+I ) ) THEN
               CALL SSCAL( L-I+1, ONE / ALPHA( K+I ), A( K+I, N-L+I ),
     $                     LDA )
            ELSE
               CALL SSCAL( L-I+1, ONE / BETA( K+I ), B( I, N-L+I ),
     $                     LDB )
               CALL SCOPY( L-I+1, B( I, N-L+I ), LDB, A( K+I, N-L+I ),
     $                     LDA )
            END IF
<span class="comment">*</span><span class="comment">
</span>         ELSE
<span class="comment">*</span><span class="comment">
</span>            ALPHA( K+I ) = ZERO
            BETA( K+I ) = ONE
            CALL SCOPY( L-I+1, B( I, N-L+I ), LDB, A( K+I, N-L+I ),
     $                  LDA )
<span class="comment">*</span><span class="comment">
</span>         END IF
<span class="comment">*</span><span class="comment">
</span>   70 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Post-assignment
</span><span class="comment">*</span><span class="comment">
</span>      DO 80 I = M + 1, K + L
         ALPHA( I ) = ZERO
         BETA( I ) = ONE
   80 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      IF( K+L.LT.N ) THEN
         DO 90 I = K + L + 1, N
            ALPHA( I ) = ZERO
            BETA( I ) = ZERO
   90    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span>  100 CONTINUE
      NCYCLE = KCYCLE
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="STGSJA.513"></a><a href="stgsja.f.html#STGSJA.1">STGSJA</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>