stgsja.f.html

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

</span><span class="comment">*</span><span class="comment">          The leading dimension of the array B. LDB &gt;= max(1,P).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  TOLA    (input) REAL
</span><span class="comment">*</span><span class="comment">  TOLB    (input) REAL
</span><span class="comment">*</span><span class="comment">          TOLA and TOLB are the convergence criteria for the Jacobi-
</span><span class="comment">*</span><span class="comment">          Kogbetliantz iteration procedure. Generally, they are the
</span><span class="comment">*</span><span class="comment">          same as used in the preprocessing step, say
</span><span class="comment">*</span><span class="comment">              TOLA = max(M,N)*norm(A)*MACHEPS,
</span><span class="comment">*</span><span class="comment">              TOLB = max(P,N)*norm(B)*MACHEPS.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  ALPHA   (output) REAL array, dimension (N)
</span><span class="comment">*</span><span class="comment">  BETA    (output) REAL array, dimension (N)
</span><span class="comment">*</span><span class="comment">          On exit, ALPHA and BETA contain the generalized singular
</span><span class="comment">*</span><span class="comment">          value pairs of A and B;
</span><span class="comment">*</span><span class="comment">            ALPHA(1:K) = 1,
</span><span class="comment">*</span><span class="comment">            BETA(1:K)  = 0,
</span><span class="comment">*</span><span class="comment">          and if M-K-L &gt;= 0,
</span><span class="comment">*</span><span class="comment">            ALPHA(K+1:K+L) = diag(C),
</span><span class="comment">*</span><span class="comment">            BETA(K+1:K+L)  = diag(S),
</span><span class="comment">*</span><span class="comment">          or if M-K-L &lt; 0,
</span><span class="comment">*</span><span class="comment">            ALPHA(K+1:M)= C, ALPHA(M+1:K+L)= 0
</span><span class="comment">*</span><span class="comment">            BETA(K+1:M) = S, BETA(M+1:K+L) = 1.
</span><span class="comment">*</span><span class="comment">          Furthermore, if K+L &lt; N,
</span><span class="comment">*</span><span class="comment">            ALPHA(K+L+1:N) = 0 and
</span><span class="comment">*</span><span class="comment">            BETA(K+L+1:N)  = 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  U       (input/output) REAL array, dimension (LDU,M)
</span><span class="comment">*</span><span class="comment">          On entry, if JOBU = 'U', U must contain a matrix U1 (usually
</span><span class="comment">*</span><span class="comment">          the orthogonal matrix returned by <a name="SGGSVP.190"></a><a href="sggsvp.f.html#SGGSVP.1">SGGSVP</a>).
</span><span class="comment">*</span><span class="comment">          On exit,
</span><span class="comment">*</span><span class="comment">          if JOBU = 'I', U contains the orthogonal matrix U;
</span><span class="comment">*</span><span class="comment">          if JOBU = 'U', U contains the product U1*U.
</span><span class="comment">*</span><span class="comment">          If JOBU = 'N', U is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDU     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array U. LDU &gt;= max(1,M) if
</span><span class="comment">*</span><span class="comment">          JOBU = 'U'; LDU &gt;= 1 otherwise.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  V       (input/output) REAL array, dimension (LDV,P)
</span><span class="comment">*</span><span class="comment">          On entry, if JOBV = 'V', V must contain a matrix V1 (usually
</span><span class="comment">*</span><span class="comment">          the orthogonal matrix returned by <a name="SGGSVP.202"></a><a href="sggsvp.f.html#SGGSVP.1">SGGSVP</a>).
</span><span class="comment">*</span><span class="comment">          On exit,
</span><span class="comment">*</span><span class="comment">          if JOBV = 'I', V contains the orthogonal matrix V;
</span><span class="comment">*</span><span class="comment">          if JOBV = 'V', V contains the product V1*V.
</span><span class="comment">*</span><span class="comment">          If JOBV = 'N', V is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDV     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array V. LDV &gt;= max(1,P) if
</span><span class="comment">*</span><span class="comment">          JOBV = 'V'; LDV &gt;= 1 otherwise.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Q       (input/output) REAL array, dimension (LDQ,N)
</span><span class="comment">*</span><span class="comment">          On entry, if JOBQ = 'Q', Q must contain a matrix Q1 (usually
</span><span class="comment">*</span><span class="comment">          the orthogonal matrix returned by <a name="SGGSVP.214"></a><a href="sggsvp.f.html#SGGSVP.1">SGGSVP</a>).
</span><span class="comment">*</span><span class="comment">          On exit,
</span><span class="comment">*</span><span class="comment">          if JOBQ = 'I', Q contains the orthogonal matrix Q;
</span><span class="comment">*</span><span class="comment">          if JOBQ = 'Q', Q contains the product Q1*Q.
</span><span class="comment">*</span><span class="comment">          If JOBQ = 'N', Q is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDQ     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array Q. LDQ &gt;= max(1,N) if
</span><span class="comment">*</span><span class="comment">          JOBQ = 'Q'; LDQ &gt;= 1 otherwise.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  WORK    (workspace) REAL array, dimension (2*N)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  NCYCLE  (output) INTEGER
</span><span class="comment">*</span><span class="comment">          The number of cycles required for convergence.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INFO    (output) INTEGER
</span><span class="comment">*</span><span class="comment">          = 0:  successful exit
</span><span class="comment">*</span><span class="comment">          &lt; 0:  if INFO = -i, the i-th argument had an illegal value.
</span><span class="comment">*</span><span class="comment">          = 1:  the procedure does not converge after MAXIT cycles.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Internal Parameters
</span><span class="comment">*</span><span class="comment">  ===================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  MAXIT   INTEGER
</span><span class="comment">*</span><span class="comment">          MAXIT specifies the total loops that the iterative procedure
</span><span class="comment">*</span><span class="comment">          may take. If after MAXIT cycles, the routine fails to
</span><span class="comment">*</span><span class="comment">          converge, we return INFO = 1.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Further Details
</span><span class="comment">*</span><span class="comment">  ===============
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  <a name="STGSJA.245"></a><a href="stgsja.f.html#STGSJA.1">STGSJA</a> essentially uses a variant of Kogbetliantz algorithm to reduce
</span><span class="comment">*</span><span class="comment">  min(L,M-K)-by-L triangular (or trapezoidal) matrix A23 and L-by-L
</span><span class="comment">*</span><span class="comment">  matrix B13 to the form:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           U1'*A13*Q1 = C1*R1; V1'*B13*Q1 = S1*R1,
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  where U1, V1 and Q1 are orthogonal matrix, and Z' is the transpose
</span><span class="comment">*</span><span class="comment">  of Z.  C1 and S1 are diagonal matrices satisfying
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                C1**2 + S1**2 = I,
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  and R1 is an L-by-L nonsingular upper triangular matrix.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      INTEGER            MAXIT
      PARAMETER          ( MAXIT = 40 )
      REAL               ZERO, ONE
      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span><span class="comment">*</span><span class="comment">
</span>      LOGICAL            INITQ, INITU, INITV, UPPER, WANTQ, WANTU, WANTV
      INTEGER            I, J, KCYCLE
      REAL               A1, A2, A3, B1, B2, B3, CSQ, CSU, CSV, ERROR,
     $                   GAMMA, RWK, SNQ, SNU, SNV, SSMIN
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.274"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      EXTERNAL           <a name="LSAME.275"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           SCOPY, <a name="SLAGS2.278"></a><a href="slags2.f.html#SLAGS2.1">SLAGS2</a>, <a name="SLAPLL.278"></a><a href="slapll.f.html#SLAPLL.1">SLAPLL</a>, <a name="SLARTG.278"></a><a href="slartg.f.html#SLARTG.1">SLARTG</a>, <a name="SLASET.278"></a><a href="slaset.f.html#SLASET.1">SLASET</a>, SROT,
     $                   SSCAL, <a name="XERBLA.279"></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">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, 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">     Decode and test the input parameters
</span><span class="comment">*</span><span class="comment">
</span>      INITU = <a name="LSAME.288"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOBU, <span class="string">'I'</span> )
      WANTU = INITU .OR. <a name="LSAME.289"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOBU, <span class="string">'U'</span> )
<span class="comment">*</span><span class="comment">
</span>      INITV = <a name="LSAME.291"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOBV, <span class="string">'I'</span> )
      WANTV = INITV .OR. <a name="LSAME.292"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOBV, <span class="string">'V'</span> )
<span class="comment">*</span><span class="comment">
</span>      INITQ = <a name="LSAME.294"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOBQ, <span class="string">'I'</span> )
      WANTQ = INITQ .OR. <a name="LSAME.295"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOBQ, <span class="string">'Q'</span> )
<span class="comment">*</span><span class="comment">
</span>      INFO = 0
      IF( .NOT.( INITU .OR. WANTU .OR. <a name="LSAME.298"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOBU, <span class="string">'N'</span> ) ) ) THEN
         INFO = -1
      ELSE IF( .NOT.( INITV .OR. WANTV .OR. <a name="LSAME.300"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOBV, <span class="string">'N'</span> ) ) ) THEN
         INFO = -2
      ELSE IF( .NOT.( INITQ .OR. WANTQ .OR. <a name="LSAME.302"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( JOBQ, <span class="string">'N'</span> ) ) ) THEN
         INFO = -3
      ELSE IF( M.LT.0 ) THEN
         INFO = -4
      ELSE IF( P.LT.0 ) THEN
         INFO = -5
      ELSE IF( N.LT.0 ) THEN
         INFO = -6
      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
         INFO = -10
      ELSE IF( LDB.LT.MAX( 1, P ) ) THEN
         INFO = -12
      ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN
         INFO = -18
      ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN
         INFO = -20
      ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
         INFO = -22
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.322"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="STGSJA.322"></a><a href="stgsja.f.html#STGSJA.1">STGSJA</a>'</span>, -INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Initialize U, V and Q, if necessary
</span><span class="comment">*</span><span class="comment">
</span>      IF( INITU )
     $   CALL <a name="SLASET.329"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'Full'</span>, M, M, ZERO, ONE, U, LDU )
      IF( INITV )
     $   CALL <a name="SLASET.331"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'Full'</span>, P, P, ZERO, ONE, V, LDV )
      IF( INITQ )
     $   CALL <a name="SLASET.333"></a><a href="slaset.f.html#SLASET.1">SLASET</a>( <span class="string">'Full'</span>, N, N, ZERO, ONE, Q, LDQ )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Loop until convergence
</span><span class="comment">*</span><span class="comment">
</span>      UPPER = .FALSE.
      DO 40 KCYCLE = 1, MAXIT
<span class="comment">*</span><span class="comment">
</span>         UPPER = .NOT.UPPER
<span class="comment">*</span><span class="comment">