LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
stgsja.f
Go to the documentation of this file.
00001 *> \brief \b STGSJA
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download STGSJA + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/stgsja.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/stgsja.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/stgsja.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE STGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B,
00022 *                          LDB, TOLA, TOLB, ALPHA, BETA, U, LDU, V, LDV,
00023 *                          Q, LDQ, WORK, NCYCLE, INFO )
00024 * 
00025 *       .. Scalar Arguments ..
00026 *       CHARACTER          JOBQ, JOBU, JOBV
00027 *       INTEGER            INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N,
00028 *      $                   NCYCLE, P
00029 *       REAL               TOLA, TOLB
00030 *       ..
00031 *       .. Array Arguments ..
00032 *       REAL               A( LDA, * ), ALPHA( * ), B( LDB, * ),
00033 *      $                   BETA( * ), Q( LDQ, * ), U( LDU, * ),
00034 *      $                   V( LDV, * ), WORK( * )
00035 *       ..
00036 *  
00037 *
00038 *> \par Purpose:
00039 *  =============
00040 *>
00041 *> \verbatim
00042 *>
00043 *> STGSJA computes the generalized singular value decomposition (GSVD)
00044 *> of two real upper triangular (or trapezoidal) matrices A and B.
00045 *>
00046 *> On entry, it is assumed that matrices A and B have the following
00047 *> forms, which may be obtained by the preprocessing subroutine SGGSVP
00048 *> from a general M-by-N matrix A and P-by-N matrix B:
00049 *>
00050 *>              N-K-L  K    L
00051 *>    A =    K ( 0    A12  A13 ) if M-K-L >= 0;
00052 *>           L ( 0     0   A23 )
00053 *>       M-K-L ( 0     0    0  )
00054 *>
00055 *>            N-K-L  K    L
00056 *>    A =  K ( 0    A12  A13 ) if M-K-L < 0;
00057 *>       M-K ( 0     0   A23 )
00058 *>
00059 *>            N-K-L  K    L
00060 *>    B =  L ( 0     0   B13 )
00061 *>       P-L ( 0     0    0  )
00062 *>
00063 *> where the K-by-K matrix A12 and L-by-L matrix B13 are nonsingular
00064 *> upper triangular; A23 is L-by-L upper triangular if M-K-L >= 0,
00065 *> otherwise A23 is (M-K)-by-L upper trapezoidal.
00066 *>
00067 *> On exit,
00068 *>
00069 *>        U**T *A*Q = D1*( 0 R ),    V**T *B*Q = D2*( 0 R ),
00070 *>
00071 *> where U, V and Q are orthogonal matrices.
00072 *> R is a nonsingular upper triangular matrix, and D1 and D2 are
00073 *> ``diagonal'' matrices, which are of the following structures:
00074 *>
00075 *> If M-K-L >= 0,
00076 *>
00077 *>                     K  L
00078 *>        D1 =     K ( I  0 )
00079 *>                 L ( 0  C )
00080 *>             M-K-L ( 0  0 )
00081 *>
00082 *>                   K  L
00083 *>        D2 = L   ( 0  S )
00084 *>             P-L ( 0  0 )
00085 *>
00086 *>                N-K-L  K    L
00087 *>   ( 0 R ) = K (  0   R11  R12 ) K
00088 *>             L (  0    0   R22 ) L
00089 *>
00090 *> where
00091 *>
00092 *>   C = diag( ALPHA(K+1), ... , ALPHA(K+L) ),
00093 *>   S = diag( BETA(K+1),  ... , BETA(K+L) ),
00094 *>   C**2 + S**2 = I.
00095 *>
00096 *>   R is stored in A(1:K+L,N-K-L+1:N) on exit.
00097 *>
00098 *> If M-K-L < 0,
00099 *>
00100 *>                K M-K K+L-M
00101 *>     D1 =   K ( I  0    0   )
00102 *>          M-K ( 0  C    0   )
00103 *>
00104 *>                  K M-K K+L-M
00105 *>     D2 =   M-K ( 0  S    0   )
00106 *>          K+L-M ( 0  0    I   )
00107 *>            P-L ( 0  0    0   )
00108 *>
00109 *>                N-K-L  K   M-K  K+L-M
00110 *> ( 0 R ) =    K ( 0    R11  R12  R13  )
00111 *>           M-K ( 0     0   R22  R23  )
00112 *>         K+L-M ( 0     0    0   R33  )
00113 *>
00114 *> where
00115 *> C = diag( ALPHA(K+1), ... , ALPHA(M) ),
00116 *> S = diag( BETA(K+1),  ... , BETA(M) ),
00117 *> C**2 + S**2 = I.
00118 *>
00119 *> R = ( R11 R12 R13 ) is stored in A(1:M, N-K-L+1:N) and R33 is stored
00120 *>     (  0  R22 R23 )
00121 *> in B(M-K+1:L,N+M-K-L+1:N) on exit.
00122 *>
00123 *> The computation of the orthogonal transformation matrices U, V or Q
00124 *> is optional.  These matrices may either be formed explicitly, or they
00125 *> may be postmultiplied into input matrices U1, V1, or Q1.
00126 *> \endverbatim
00127 *
00128 *  Arguments:
00129 *  ==========
00130 *
00131 *> \param[in] JOBU
00132 *> \verbatim
00133 *>          JOBU is CHARACTER*1
00134 *>          = 'U':  U must contain an orthogonal matrix U1 on entry, and
00135 *>                  the product U1*U is returned;
00136 *>          = 'I':  U is initialized to the unit matrix, and the
00137 *>                  orthogonal matrix U is returned;
00138 *>          = 'N':  U is not computed.
00139 *> \endverbatim
00140 *>
00141 *> \param[in] JOBV
00142 *> \verbatim
00143 *>          JOBV is CHARACTER*1
00144 *>          = 'V':  V must contain an orthogonal matrix V1 on entry, and
00145 *>                  the product V1*V is returned;
00146 *>          = 'I':  V is initialized to the unit matrix, and the
00147 *>                  orthogonal matrix V is returned;
00148 *>          = 'N':  V is not computed.
00149 *> \endverbatim
00150 *>
00151 *> \param[in] JOBQ
00152 *> \verbatim
00153 *>          JOBQ is CHARACTER*1
00154 *>          = 'Q':  Q must contain an orthogonal matrix Q1 on entry, and
00155 *>                  the product Q1*Q is returned;
00156 *>          = 'I':  Q is initialized to the unit matrix, and the
00157 *>                  orthogonal matrix Q is returned;
00158 *>          = 'N':  Q is not computed.
00159 *> \endverbatim
00160 *>
00161 *> \param[in] M
00162 *> \verbatim
00163 *>          M is INTEGER
00164 *>          The number of rows of the matrix A.  M >= 0.
00165 *> \endverbatim
00166 *>
00167 *> \param[in] P
00168 *> \verbatim
00169 *>          P is INTEGER
00170 *>          The number of rows of the matrix B.  P >= 0.
00171 *> \endverbatim
00172 *>
00173 *> \param[in] N
00174 *> \verbatim
00175 *>          N is INTEGER
00176 *>          The number of columns of the matrices A and B.  N >= 0.
00177 *> \endverbatim
00178 *>
00179 *> \param[in] K
00180 *> \verbatim
00181 *>          K is INTEGER
00182 *> \endverbatim
00183 *>
00184 *> \param[in] L
00185 *> \verbatim
00186 *>          L is INTEGER
00187 *>
00188 *>          K and L specify the subblocks in the input matrices A and B:
00189 *>          A23 = A(K+1:MIN(K+L,M),N-L+1:N) and B13 = B(1:L,N-L+1:N)
00190 *>          of A and B, whose GSVD is going to be computed by STGSJA.
00191 *>          See Further Details.
00192 *> \endverbatim
00193 *>
00194 *> \param[in,out] A
00195 *> \verbatim
00196 *>          A is REAL array, dimension (LDA,N)
00197 *>          On entry, the M-by-N matrix A.
00198 *>          On exit, A(N-K+1:N,1:MIN(K+L,M) ) contains the triangular
00199 *>          matrix R or part of R.  See Purpose for details.
00200 *> \endverbatim
00201 *>
00202 *> \param[in] LDA
00203 *> \verbatim
00204 *>          LDA is INTEGER
00205 *>          The leading dimension of the array A. LDA >= max(1,M).
00206 *> \endverbatim
00207 *>
00208 *> \param[in,out] B
00209 *> \verbatim
00210 *>          B is REAL array, dimension (LDB,N)
00211 *>          On entry, the P-by-N matrix B.
00212 *>          On exit, if necessary, B(M-K+1:L,N+M-K-L+1:N) contains
00213 *>          a part of R.  See Purpose for details.
00214 *> \endverbatim
00215 *>
00216 *> \param[in] LDB
00217 *> \verbatim
00218 *>          LDB is INTEGER
00219 *>          The leading dimension of the array B. LDB >= max(1,P).
00220 *> \endverbatim
00221 *>
00222 *> \param[in] TOLA
00223 *> \verbatim
00224 *>          TOLA is REAL
00225 *> \endverbatim
00226 *>
00227 *> \param[in] TOLB
00228 *> \verbatim
00229 *>          TOLB is REAL
00230 *>
00231 *>          TOLA and TOLB are the convergence criteria for the Jacobi-
00232 *>          Kogbetliantz iteration procedure. Generally, they are the
00233 *>          same as used in the preprocessing step, say
00234 *>              TOLA = max(M,N)*norm(A)*MACHEPS,
00235 *>              TOLB = max(P,N)*norm(B)*MACHEPS.
00236 *> \endverbatim
00237 *>
00238 *> \param[out] ALPHA
00239 *> \verbatim
00240 *>          ALPHA is REAL array, dimension (N)
00241 *> \endverbatim
00242 *>
00243 *> \param[out] BETA
00244 *> \verbatim
00245 *>          BETA is REAL array, dimension (N)
00246 *>
00247 *>          On exit, ALPHA and BETA contain the generalized singular
00248 *>          value pairs of A and B;
00249 *>            ALPHA(1:K) = 1,
00250 *>            BETA(1:K)  = 0,
00251 *>          and if M-K-L >= 0,
00252 *>            ALPHA(K+1:K+L) = diag(C),
00253 *>            BETA(K+1:K+L)  = diag(S),
00254 *>          or if M-K-L < 0,
00255 *>            ALPHA(K+1:M)= C, ALPHA(M+1:K+L)= 0
00256 *>            BETA(K+1:M) = S, BETA(M+1:K+L) = 1.
00257 *>          Furthermore, if K+L < N,
00258 *>            ALPHA(K+L+1:N) = 0 and
00259 *>            BETA(K+L+1:N)  = 0.
00260 *> \endverbatim
00261 *>
00262 *> \param[in,out] U
00263 *> \verbatim
00264 *>          U is REAL array, dimension (LDU,M)
00265 *>          On entry, if JOBU = 'U', U must contain a matrix U1 (usually
00266 *>          the orthogonal matrix returned by SGGSVP).
00267 *>          On exit,
00268 *>          if JOBU = 'I', U contains the orthogonal matrix U;
00269 *>          if JOBU = 'U', U contains the product U1*U.
00270 *>          If JOBU = 'N', U is not referenced.
00271 *> \endverbatim
00272 *>
00273 *> \param[in] LDU
00274 *> \verbatim
00275 *>          LDU is INTEGER
00276 *>          The leading dimension of the array U. LDU >= max(1,M) if
00277 *>          JOBU = 'U'; LDU >= 1 otherwise.
00278 *> \endverbatim
00279 *>
00280 *> \param[in,out] V
00281 *> \verbatim
00282 *>          V is REAL array, dimension (LDV,P)
00283 *>          On entry, if JOBV = 'V', V must contain a matrix V1 (usually
00284 *>          the orthogonal matrix returned by SGGSVP).
00285 *>          On exit,
00286 *>          if JOBV = 'I', V contains the orthogonal matrix V;
00287 *>          if JOBV = 'V', V contains the product V1*V.
00288 *>          If JOBV = 'N', V is not referenced.
00289 *> \endverbatim
00290 *>
00291 *> \param[in] LDV
00292 *> \verbatim
00293 *>          LDV is INTEGER
00294 *>          The leading dimension of the array V. LDV >= max(1,P) if
00295 *>          JOBV = 'V'; LDV >= 1 otherwise.
00296 *> \endverbatim
00297 *>
00298 *> \param[in,out] Q
00299 *> \verbatim
00300 *>          Q is REAL array, dimension (LDQ,N)
00301 *>          On entry, if JOBQ = 'Q', Q must contain a matrix Q1 (usually
00302 *>          the orthogonal matrix returned by SGGSVP).
00303 *>          On exit,
00304 *>          if JOBQ = 'I', Q contains the orthogonal matrix Q;
00305 *>          if JOBQ = 'Q', Q contains the product Q1*Q.
00306 *>          If JOBQ = 'N', Q is not referenced.
00307 *> \endverbatim
00308 *>
00309 *> \param[in] LDQ
00310 *> \verbatim
00311 *>          LDQ is INTEGER
00312 *>          The leading dimension of the array Q. LDQ >= max(1,N) if
00313 *>          JOBQ = 'Q'; LDQ >= 1 otherwise.
00314 *> \endverbatim
00315 *>
00316 *> \param[out] WORK
00317 *> \verbatim
00318 *>          WORK is REAL array, dimension (2*N)
00319 *> \endverbatim
00320 *>
00321 *> \param[out] NCYCLE
00322 *> \verbatim
00323 *>          NCYCLE is INTEGER
00324 *>          The number of cycles required for convergence.
00325 *> \endverbatim
00326 *>
00327 *> \param[out] INFO
00328 *> \verbatim
00329 *>          INFO is INTEGER
00330 *>          = 0:  successful exit
00331 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
00332 *>          = 1:  the procedure does not converge after MAXIT cycles.
00333 *> \endverbatim
00334 *>
00335 *> \verbatim
00336 *>  Internal Parameters
00337 *>  ===================
00338 *>
00339 *>  MAXIT   INTEGER
00340 *>          MAXIT specifies the total loops that the iterative procedure
00341 *>          may take. If after MAXIT cycles, the routine fails to
00342 *>          converge, we return INFO = 1.
00343 *> \endverbatim
00344 *
00345 *  Authors:
00346 *  ========
00347 *
00348 *> \author Univ. of Tennessee 
00349 *> \author Univ. of California Berkeley 
00350 *> \author Univ. of Colorado Denver 
00351 *> \author NAG Ltd. 
00352 *
00353 *> \date November 2011
00354 *
00355 *> \ingroup realOTHERcomputational
00356 *
00357 *> \par Further Details:
00358 *  =====================
00359 *>
00360 *> \verbatim
00361 *>
00362 *>  STGSJA essentially uses a variant of Kogbetliantz algorithm to reduce
00363 *>  min(L,M-K)-by-L triangular (or trapezoidal) matrix A23 and L-by-L
00364 *>  matrix B13 to the form:
00365 *>
00366 *>           U1**T *A13*Q1 = C1*R1; V1**T *B13*Q1 = S1*R1,
00367 *>
00368 *>  where U1, V1 and Q1 are orthogonal matrix, and Z**T is the transpose
00369 *>  of Z.  C1 and S1 are diagonal matrices satisfying
00370 *>
00371 *>                C1**2 + S1**2 = I,
00372 *>
00373 *>  and R1 is an L-by-L nonsingular upper triangular matrix.
00374 *> \endverbatim
00375 *>
00376 *  =====================================================================
00377       SUBROUTINE STGSJA( JOBU, JOBV, JOBQ, M, P, N, K, L, A, LDA, B,
00378      $                   LDB, TOLA, TOLB, ALPHA, BETA, U, LDU, V, LDV,
00379      $                   Q, LDQ, WORK, NCYCLE, INFO )
00380 *
00381 *  -- LAPACK computational routine (version 3.4.0) --
00382 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00383 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00384 *     November 2011
00385 *
00386 *     .. Scalar Arguments ..
00387       CHARACTER          JOBQ, JOBU, JOBV
00388       INTEGER            INFO, K, L, LDA, LDB, LDQ, LDU, LDV, M, N,
00389      $                   NCYCLE, P
00390       REAL               TOLA, TOLB
00391 *     ..
00392 *     .. Array Arguments ..
00393       REAL               A( LDA, * ), ALPHA( * ), B( LDB, * ),
00394      $                   BETA( * ), Q( LDQ, * ), U( LDU, * ),
00395      $                   V( LDV, * ), WORK( * )
00396 *     ..
00397 *
00398 *  =====================================================================
00399 *
00400 *     .. Parameters ..
00401       INTEGER            MAXIT
00402       PARAMETER          ( MAXIT = 40 )
00403       REAL               ZERO, ONE
00404       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00405 *     ..
00406 *     .. Local Scalars ..
00407 *
00408       LOGICAL            INITQ, INITU, INITV, UPPER, WANTQ, WANTU, WANTV
00409       INTEGER            I, J, KCYCLE
00410       REAL               A1, A2, A3, B1, B2, B3, CSQ, CSU, CSV, ERROR,
00411      $                   GAMMA, RWK, SNQ, SNU, SNV, SSMIN
00412 *     ..
00413 *     .. External Functions ..
00414       LOGICAL            LSAME
00415       EXTERNAL           LSAME
00416 *     ..
00417 *     .. External Subroutines ..
00418       EXTERNAL           SCOPY, SLAGS2, SLAPLL, SLARTG, SLASET, SROT,
00419      $                   SSCAL, XERBLA
00420 *     ..
00421 *     .. Intrinsic Functions ..
00422       INTRINSIC          ABS, MAX, MIN
00423 *     ..
00424 *     .. Executable Statements ..
00425 *
00426 *     Decode and test the input parameters
00427 *
00428       INITU = LSAME( JOBU, 'I' )
00429       WANTU = INITU .OR. LSAME( JOBU, 'U' )
00430 *
00431       INITV = LSAME( JOBV, 'I' )
00432       WANTV = INITV .OR. LSAME( JOBV, 'V' )
00433 *
00434       INITQ = LSAME( JOBQ, 'I' )
00435       WANTQ = INITQ .OR. LSAME( JOBQ, 'Q' )
00436 *
00437       INFO = 0
00438       IF( .NOT.( INITU .OR. WANTU .OR. LSAME( JOBU, 'N' ) ) ) THEN
00439          INFO = -1
00440       ELSE IF( .NOT.( INITV .OR. WANTV .OR. LSAME( JOBV, 'N' ) ) ) THEN
00441          INFO = -2
00442       ELSE IF( .NOT.( INITQ .OR. WANTQ .OR. LSAME( JOBQ, 'N' ) ) ) THEN
00443          INFO = -3
00444       ELSE IF( M.LT.0 ) THEN
00445          INFO = -4
00446       ELSE IF( P.LT.0 ) THEN
00447          INFO = -5
00448       ELSE IF( N.LT.0 ) THEN
00449          INFO = -6
00450       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00451          INFO = -10
00452       ELSE IF( LDB.LT.MAX( 1, P ) ) THEN
00453          INFO = -12
00454       ELSE IF( LDU.LT.1 .OR. ( WANTU .AND. LDU.LT.M ) ) THEN
00455          INFO = -18
00456       ELSE IF( LDV.LT.1 .OR. ( WANTV .AND. LDV.LT.P ) ) THEN
00457          INFO = -20
00458       ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.N ) ) THEN
00459          INFO = -22
00460       END IF
00461       IF( INFO.NE.0 ) THEN
00462          CALL XERBLA( 'STGSJA', -INFO )
00463          RETURN
00464       END IF
00465 *
00466 *     Initialize U, V and Q, if necessary
00467 *
00468       IF( INITU )
00469      $   CALL SLASET( 'Full', M, M, ZERO, ONE, U, LDU )
00470       IF( INITV )
00471      $   CALL SLASET( 'Full', P, P, ZERO, ONE, V, LDV )
00472       IF( INITQ )
00473      $   CALL SLASET( 'Full', N, N, ZERO, ONE, Q, LDQ )
00474 *
00475 *     Loop until convergence
00476 *
00477       UPPER = .FALSE.
00478       DO 40 KCYCLE = 1, MAXIT
00479 *
00480          UPPER = .NOT.UPPER
00481 *
00482          DO 20 I = 1, L - 1
00483             DO 10 J = I + 1, L
00484 *
00485                A1 = ZERO
00486                A2 = ZERO
00487                A3 = ZERO
00488                IF( K+I.LE.M )
00489      $            A1 = A( K+I, N-L+I )
00490                IF( K+J.LE.M )
00491      $            A3 = A( K+J, N-L+J )
00492 *
00493                B1 = B( I, N-L+I )
00494                B3 = B( J, N-L+J )
00495 *
00496                IF( UPPER ) THEN
00497                   IF( K+I.LE.M )
00498      $               A2 = A( K+I, N-L+J )
00499                   B2 = B( I, N-L+J )
00500                ELSE
00501                   IF( K+J.LE.M )
00502      $               A2 = A( K+J, N-L+I )
00503                   B2 = B( J, N-L+I )
00504                END IF
00505 *
00506                CALL SLAGS2( UPPER, A1, A2, A3, B1, B2, B3, CSU, SNU,
00507      $                      CSV, SNV, CSQ, SNQ )
00508 *
00509 *              Update (K+I)-th and (K+J)-th rows of matrix A: U**T *A
00510 *
00511                IF( K+J.LE.M )
00512      $            CALL SROT( L, A( K+J, N-L+1 ), LDA, A( K+I, N-L+1 ),
00513      $                       LDA, CSU, SNU )
00514 *
00515 *              Update I-th and J-th rows of matrix B: V**T *B
00516 *
00517                CALL SROT( L, B( J, N-L+1 ), LDB, B( I, N-L+1 ), LDB,
00518      $                    CSV, SNV )
00519 *
00520 *              Update (N-L+I)-th and (N-L+J)-th columns of matrices
00521 *              A and B: A*Q and B*Q
00522 *
00523                CALL SROT( MIN( K+L, M ), A( 1, N-L+J ), 1,
00524      $                    A( 1, N-L+I ), 1, CSQ, SNQ )
00525 *
00526                CALL SROT( L, B( 1, N-L+J ), 1, B( 1, N-L+I ), 1, CSQ,
00527      $                    SNQ )
00528 *
00529                IF( UPPER ) THEN
00530                   IF( K+I.LE.M )
00531      $               A( K+I, N-L+J ) = ZERO
00532                   B( I, N-L+J ) = ZERO
00533                ELSE
00534                   IF( K+J.LE.M )
00535      $               A( K+J, N-L+I ) = ZERO
00536                   B( J, N-L+I ) = ZERO
00537                END IF
00538 *
00539 *              Update orthogonal matrices U, V, Q, if desired.
00540 *
00541                IF( WANTU .AND. K+J.LE.M )
00542      $            CALL SROT( M, U( 1, K+J ), 1, U( 1, K+I ), 1, CSU,
00543      $                       SNU )
00544 *
00545                IF( WANTV )
00546      $            CALL SROT( P, V( 1, J ), 1, V( 1, I ), 1, CSV, SNV )
00547 *
00548                IF( WANTQ )
00549      $            CALL SROT( N, Q( 1, N-L+J ), 1, Q( 1, N-L+I ), 1, CSQ,
00550      $                       SNQ )
00551 *
00552    10       CONTINUE
00553    20    CONTINUE
00554 *
00555          IF( .NOT.UPPER ) THEN
00556 *
00557 *           The matrices A13 and B13 were lower triangular at the start
00558 *           of the cycle, and are now upper triangular.
00559 *
00560 *           Convergence test: test the parallelism of the corresponding
00561 *           rows of A and B.
00562 *
00563             ERROR = ZERO
00564             DO 30 I = 1, MIN( L, M-K )
00565                CALL SCOPY( L-I+1, A( K+I, N-L+I ), LDA, WORK, 1 )
00566                CALL SCOPY( L-I+1, B( I, N-L+I ), LDB, WORK( L+1 ), 1 )
00567                CALL SLAPLL( L-I+1, WORK, 1, WORK( L+1 ), 1, SSMIN )
00568                ERROR = MAX( ERROR, SSMIN )
00569    30       CONTINUE
00570 *
00571             IF( ABS( ERROR ).LE.MIN( TOLA, TOLB ) )
00572      $         GO TO 50
00573          END IF
00574 *
00575 *        End of cycle loop
00576 *
00577    40 CONTINUE
00578 *
00579 *     The algorithm has not converged after MAXIT cycles.
00580 *
00581       INFO = 1
00582       GO TO 100
00583 *
00584    50 CONTINUE
00585 *
00586 *     If ERROR <= MIN(TOLA,TOLB), then the algorithm has converged.
00587 *     Compute the generalized singular value pairs (ALPHA, BETA), and
00588 *     set the triangular matrix R to array A.
00589 *
00590       DO 60 I = 1, K
00591          ALPHA( I ) = ONE
00592          BETA( I ) = ZERO
00593    60 CONTINUE
00594 *
00595       DO 70 I = 1, MIN( L, M-K )
00596 *
00597          A1 = A( K+I, N-L+I )
00598          B1 = B( I, N-L+I )
00599 *
00600          IF( A1.NE.ZERO ) THEN
00601             GAMMA = B1 / A1
00602 *
00603 *           change sign if necessary
00604 *
00605             IF( GAMMA.LT.ZERO ) THEN
00606                CALL SSCAL( L-I+1, -ONE, B( I, N-L+I ), LDB )
00607                IF( WANTV )
00608      $            CALL SSCAL( P, -ONE, V( 1, I ), 1 )
00609             END IF
00610 *
00611             CALL SLARTG( ABS( GAMMA ), ONE, BETA( K+I ), ALPHA( K+I ),
00612      $                   RWK )
00613 *
00614             IF( ALPHA( K+I ).GE.BETA( K+I ) ) THEN
00615                CALL SSCAL( L-I+1, ONE / ALPHA( K+I ), A( K+I, N-L+I ),
00616      $                     LDA )
00617             ELSE
00618                CALL SSCAL( L-I+1, ONE / BETA( K+I ), B( I, N-L+I ),
00619      $                     LDB )
00620                CALL SCOPY( L-I+1, B( I, N-L+I ), LDB, A( K+I, N-L+I ),
00621      $                     LDA )
00622             END IF
00623 *
00624          ELSE
00625 *
00626             ALPHA( K+I ) = ZERO
00627             BETA( K+I ) = ONE
00628             CALL SCOPY( L-I+1, B( I, N-L+I ), LDB, A( K+I, N-L+I ),
00629      $                  LDA )
00630 *
00631          END IF
00632 *
00633    70 CONTINUE
00634 *
00635 *     Post-assignment
00636 *
00637       DO 80 I = M + 1, K + L
00638          ALPHA( I ) = ZERO
00639          BETA( I ) = ONE
00640    80 CONTINUE
00641 *
00642       IF( K+L.LT.N ) THEN
00643          DO 90 I = K + L + 1, N
00644             ALPHA( I ) = ZERO
00645             BETA( I ) = ZERO
00646    90    CONTINUE
00647       END IF
00648 *
00649   100 CONTINUE
00650       NCYCLE = KCYCLE
00651       RETURN
00652 *
00653 *     End of STGSJA
00654 *
00655       END
 All Files Functions