LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
stgsyl.f
Go to the documentation of this file.
00001 *> \brief \b STGSYL
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download STGSYL + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/stgsyl.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/stgsyl.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/stgsyl.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE STGSYL( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
00022 *                          LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK,
00023 *                          IWORK, INFO )
00024 * 
00025 *       .. Scalar Arguments ..
00026 *       CHARACTER          TRANS
00027 *       INTEGER            IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF,
00028 *      $                   LWORK, M, N
00029 *       REAL               DIF, SCALE
00030 *       ..
00031 *       .. Array Arguments ..
00032 *       INTEGER            IWORK( * )
00033 *       REAL               A( LDA, * ), B( LDB, * ), C( LDC, * ),
00034 *      $                   D( LDD, * ), E( LDE, * ), F( LDF, * ),
00035 *      $                   WORK( * )
00036 *       ..
00037 *  
00038 *
00039 *> \par Purpose:
00040 *  =============
00041 *>
00042 *> \verbatim
00043 *>
00044 *> STGSYL solves the generalized Sylvester equation:
00045 *>
00046 *>             A * R - L * B = scale * C                 (1)
00047 *>             D * R - L * E = scale * F
00048 *>
00049 *> where R and L are unknown m-by-n matrices, (A, D), (B, E) and
00050 *> (C, F) are given matrix pairs of size m-by-m, n-by-n and m-by-n,
00051 *> respectively, with real entries. (A, D) and (B, E) must be in
00052 *> generalized (real) Schur canonical form, i.e. A, B are upper quasi
00053 *> triangular and D, E are upper triangular.
00054 *>
00055 *> The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1 is an output
00056 *> scaling factor chosen to avoid overflow.
00057 *>
00058 *> In matrix notation (1) is equivalent to solve  Zx = scale b, where
00059 *> Z is defined as
00060 *>
00061 *>            Z = [ kron(In, A)  -kron(B**T, Im) ]         (2)
00062 *>                [ kron(In, D)  -kron(E**T, Im) ].
00063 *>
00064 *> Here Ik is the identity matrix of size k and X**T is the transpose of
00065 *> X. kron(X, Y) is the Kronecker product between the matrices X and Y.
00066 *>
00067 *> If TRANS = 'T', STGSYL solves the transposed system Z**T*y = scale*b,
00068 *> which is equivalent to solve for R and L in
00069 *>
00070 *>             A**T * R + D**T * L = scale * C           (3)
00071 *>             R * B**T + L * E**T = scale * -F
00072 *>
00073 *> This case (TRANS = 'T') is used to compute an one-norm-based estimate
00074 *> of Dif[(A,D), (B,E)], the separation between the matrix pairs (A,D)
00075 *> and (B,E), using SLACON.
00076 *>
00077 *> If IJOB >= 1, STGSYL computes a Frobenius norm-based estimate
00078 *> of Dif[(A,D),(B,E)]. That is, the reciprocal of a lower bound on the
00079 *> reciprocal of the smallest singular value of Z. See [1-2] for more
00080 *> information.
00081 *>
00082 *> This is a level 3 BLAS algorithm.
00083 *> \endverbatim
00084 *
00085 *  Arguments:
00086 *  ==========
00087 *
00088 *> \param[in] TRANS
00089 *> \verbatim
00090 *>          TRANS is CHARACTER*1
00091 *>          = 'N', solve the generalized Sylvester equation (1).
00092 *>          = 'T', solve the 'transposed' system (3).
00093 *> \endverbatim
00094 *>
00095 *> \param[in] IJOB
00096 *> \verbatim
00097 *>          IJOB is INTEGER
00098 *>          Specifies what kind of functionality to be performed.
00099 *>           =0: solve (1) only.
00100 *>           =1: The functionality of 0 and 3.
00101 *>           =2: The functionality of 0 and 4.
00102 *>           =3: Only an estimate of Dif[(A,D), (B,E)] is computed.
00103 *>               (look ahead strategy IJOB  = 1 is used).
00104 *>           =4: Only an estimate of Dif[(A,D), (B,E)] is computed.
00105 *>               ( SGECON on sub-systems is used ).
00106 *>          Not referenced if TRANS = 'T'.
00107 *> \endverbatim
00108 *>
00109 *> \param[in] M
00110 *> \verbatim
00111 *>          M is INTEGER
00112 *>          The order of the matrices A and D, and the row dimension of
00113 *>          the matrices C, F, R and L.
00114 *> \endverbatim
00115 *>
00116 *> \param[in] N
00117 *> \verbatim
00118 *>          N is INTEGER
00119 *>          The order of the matrices B and E, and the column dimension
00120 *>          of the matrices C, F, R and L.
00121 *> \endverbatim
00122 *>
00123 *> \param[in] A
00124 *> \verbatim
00125 *>          A is REAL array, dimension (LDA, M)
00126 *>          The upper quasi triangular matrix A.
00127 *> \endverbatim
00128 *>
00129 *> \param[in] LDA
00130 *> \verbatim
00131 *>          LDA is INTEGER
00132 *>          The leading dimension of the array A. LDA >= max(1, M).
00133 *> \endverbatim
00134 *>
00135 *> \param[in] B
00136 *> \verbatim
00137 *>          B is REAL array, dimension (LDB, N)
00138 *>          The upper quasi triangular matrix B.
00139 *> \endverbatim
00140 *>
00141 *> \param[in] LDB
00142 *> \verbatim
00143 *>          LDB is INTEGER
00144 *>          The leading dimension of the array B. LDB >= max(1, N).
00145 *> \endverbatim
00146 *>
00147 *> \param[in,out] C
00148 *> \verbatim
00149 *>          C is REAL array, dimension (LDC, N)
00150 *>          On entry, C contains the right-hand-side of the first matrix
00151 *>          equation in (1) or (3).
00152 *>          On exit, if IJOB = 0, 1 or 2, C has been overwritten by
00153 *>          the solution R. If IJOB = 3 or 4 and TRANS = 'N', C holds R,
00154 *>          the solution achieved during the computation of the
00155 *>          Dif-estimate.
00156 *> \endverbatim
00157 *>
00158 *> \param[in] LDC
00159 *> \verbatim
00160 *>          LDC is INTEGER
00161 *>          The leading dimension of the array C. LDC >= max(1, M).
00162 *> \endverbatim
00163 *>
00164 *> \param[in] D
00165 *> \verbatim
00166 *>          D is REAL array, dimension (LDD, M)
00167 *>          The upper triangular matrix D.
00168 *> \endverbatim
00169 *>
00170 *> \param[in] LDD
00171 *> \verbatim
00172 *>          LDD is INTEGER
00173 *>          The leading dimension of the array D. LDD >= max(1, M).
00174 *> \endverbatim
00175 *>
00176 *> \param[in] E
00177 *> \verbatim
00178 *>          E is REAL array, dimension (LDE, N)
00179 *>          The upper triangular matrix E.
00180 *> \endverbatim
00181 *>
00182 *> \param[in] LDE
00183 *> \verbatim
00184 *>          LDE is INTEGER
00185 *>          The leading dimension of the array E. LDE >= max(1, N).
00186 *> \endverbatim
00187 *>
00188 *> \param[in,out] F
00189 *> \verbatim
00190 *>          F is REAL array, dimension (LDF, N)
00191 *>          On entry, F contains the right-hand-side of the second matrix
00192 *>          equation in (1) or (3).
00193 *>          On exit, if IJOB = 0, 1 or 2, F has been overwritten by
00194 *>          the solution L. If IJOB = 3 or 4 and TRANS = 'N', F holds L,
00195 *>          the solution achieved during the computation of the
00196 *>          Dif-estimate.
00197 *> \endverbatim
00198 *>
00199 *> \param[in] LDF
00200 *> \verbatim
00201 *>          LDF is INTEGER
00202 *>          The leading dimension of the array F. LDF >= max(1, M).
00203 *> \endverbatim
00204 *>
00205 *> \param[out] DIF
00206 *> \verbatim
00207 *>          DIF is REAL
00208 *>          On exit DIF is the reciprocal of a lower bound of the
00209 *>          reciprocal of the Dif-function, i.e. DIF is an upper bound of
00210 *>          Dif[(A,D), (B,E)] = sigma_min(Z), where Z as in (2).
00211 *>          IF IJOB = 0 or TRANS = 'T', DIF is not touched.
00212 *> \endverbatim
00213 *>
00214 *> \param[out] SCALE
00215 *> \verbatim
00216 *>          SCALE is REAL
00217 *>          On exit SCALE is the scaling factor in (1) or (3).
00218 *>          If 0 < SCALE < 1, C and F hold the solutions R and L, resp.,
00219 *>          to a slightly perturbed system but the input matrices A, B, D
00220 *>          and E have not been changed. If SCALE = 0, C and F hold the
00221 *>          solutions R and L, respectively, to the homogeneous system
00222 *>          with C = F = 0. Normally, SCALE = 1.
00223 *> \endverbatim
00224 *>
00225 *> \param[out] WORK
00226 *> \verbatim
00227 *>          WORK is REAL array, dimension (MAX(1,LWORK))
00228 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00229 *> \endverbatim
00230 *>
00231 *> \param[in] LWORK
00232 *> \verbatim
00233 *>          LWORK is INTEGER
00234 *>          The dimension of the array WORK. LWORK > = 1.
00235 *>          If IJOB = 1 or 2 and TRANS = 'N', LWORK >= max(1,2*M*N).
00236 *>
00237 *>          If LWORK = -1, then a workspace query is assumed; the routine
00238 *>          only calculates the optimal size of the WORK array, returns
00239 *>          this value as the first entry of the WORK array, and no error
00240 *>          message related to LWORK is issued by XERBLA.
00241 *> \endverbatim
00242 *>
00243 *> \param[out] IWORK
00244 *> \verbatim
00245 *>          IWORK is INTEGER array, dimension (M+N+6)
00246 *> \endverbatim
00247 *>
00248 *> \param[out] INFO
00249 *> \verbatim
00250 *>          INFO is INTEGER
00251 *>            =0: successful exit
00252 *>            <0: If INFO = -i, the i-th argument had an illegal value.
00253 *>            >0: (A, D) and (B, E) have common or close eigenvalues.
00254 *> \endverbatim
00255 *
00256 *  Authors:
00257 *  ========
00258 *
00259 *> \author Univ. of Tennessee 
00260 *> \author Univ. of California Berkeley 
00261 *> \author Univ. of Colorado Denver 
00262 *> \author NAG Ltd. 
00263 *
00264 *> \date November 2011
00265 *
00266 *> \ingroup realSYcomputational
00267 *
00268 *> \par Contributors:
00269 *  ==================
00270 *>
00271 *>     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
00272 *>     Umea University, S-901 87 Umea, Sweden.
00273 *
00274 *> \par References:
00275 *  ================
00276 *>
00277 *> \verbatim
00278 *>
00279 *>  [1] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
00280 *>      for Solving the Generalized Sylvester Equation and Estimating the
00281 *>      Separation between Regular Matrix Pairs, Report UMINF - 93.23,
00282 *>      Department of Computing Science, Umea University, S-901 87 Umea,
00283 *>      Sweden, December 1993, Revised April 1994, Also as LAPACK Working
00284 *>      Note 75.  To appear in ACM Trans. on Math. Software, Vol 22,
00285 *>      No 1, 1996.
00286 *>
00287 *>  [2] B. Kagstrom, A Perturbation Analysis of the Generalized Sylvester
00288 *>      Equation (AR - LB, DR - LE ) = (C, F), SIAM J. Matrix Anal.
00289 *>      Appl., 15(4):1045-1060, 1994
00290 *>
00291 *>  [3] B. Kagstrom and L. Westin, Generalized Schur Methods with
00292 *>      Condition Estimators for Solving the Generalized Sylvester
00293 *>      Equation, IEEE Transactions on Automatic Control, Vol. 34, No. 7,
00294 *>      July 1989, pp 745-751.
00295 *> \endverbatim
00296 *>
00297 *  =====================================================================
00298       SUBROUTINE STGSYL( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
00299      $                   LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK,
00300      $                   IWORK, INFO )
00301 *
00302 *  -- LAPACK computational routine (version 3.4.0) --
00303 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00304 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00305 *     November 2011
00306 *
00307 *     .. Scalar Arguments ..
00308       CHARACTER          TRANS
00309       INTEGER            IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF,
00310      $                   LWORK, M, N
00311       REAL               DIF, SCALE
00312 *     ..
00313 *     .. Array Arguments ..
00314       INTEGER            IWORK( * )
00315       REAL               A( LDA, * ), B( LDB, * ), C( LDC, * ),
00316      $                   D( LDD, * ), E( LDE, * ), F( LDF, * ),
00317      $                   WORK( * )
00318 *     ..
00319 *
00320 *  =====================================================================
00321 *  Replaced various illegal calls to SCOPY by calls to SLASET.
00322 *  Sven Hammarling, 1/5/02.
00323 *
00324 *     .. Parameters ..
00325       REAL               ZERO, ONE
00326       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00327 *     ..
00328 *     .. Local Scalars ..
00329       LOGICAL            LQUERY, NOTRAN
00330       INTEGER            I, IE, IFUNC, IROUND, IS, ISOLVE, J, JE, JS, K,
00331      $                   LINFO, LWMIN, MB, NB, P, PPQQ, PQ, Q
00332       REAL               DSCALE, DSUM, SCALE2, SCALOC
00333 *     ..
00334 *     .. External Functions ..
00335       LOGICAL            LSAME
00336       INTEGER            ILAENV
00337       EXTERNAL           LSAME, ILAENV
00338 *     ..
00339 *     .. External Subroutines ..
00340       EXTERNAL           SGEMM, SLACPY, SLASET, SSCAL, STGSY2, XERBLA
00341 *     ..
00342 *     .. Intrinsic Functions ..
00343       INTRINSIC          MAX, REAL, SQRT
00344 *     ..
00345 *     .. Executable Statements ..
00346 *
00347 *     Decode and test input parameters
00348 *
00349       INFO = 0
00350       NOTRAN = LSAME( TRANS, 'N' )
00351       LQUERY = ( LWORK.EQ.-1 )
00352 *
00353       IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN
00354          INFO = -1
00355       ELSE IF( NOTRAN ) THEN
00356          IF( ( IJOB.LT.0 ) .OR. ( IJOB.GT.4 ) ) THEN
00357             INFO = -2
00358          END IF
00359       END IF
00360       IF( INFO.EQ.0 ) THEN
00361          IF( M.LE.0 ) THEN
00362             INFO = -3
00363          ELSE IF( N.LE.0 ) THEN
00364             INFO = -4
00365          ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00366             INFO = -6
00367          ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00368             INFO = -8
00369          ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
00370             INFO = -10
00371          ELSE IF( LDD.LT.MAX( 1, M ) ) THEN
00372             INFO = -12
00373          ELSE IF( LDE.LT.MAX( 1, N ) ) THEN
00374             INFO = -14
00375          ELSE IF( LDF.LT.MAX( 1, M ) ) THEN
00376             INFO = -16
00377          END IF
00378       END IF
00379 *
00380       IF( INFO.EQ.0 ) THEN
00381          IF( NOTRAN ) THEN
00382             IF( IJOB.EQ.1 .OR. IJOB.EQ.2 ) THEN
00383                LWMIN = MAX( 1, 2*M*N )
00384             ELSE
00385                LWMIN = 1
00386             END IF
00387          ELSE
00388             LWMIN = 1
00389          END IF
00390          WORK( 1 ) = LWMIN
00391 *
00392          IF( LWORK.LT.LWMIN .AND. .NOT.LQUERY ) THEN
00393             INFO = -20
00394          END IF
00395       END IF
00396 *
00397       IF( INFO.NE.0 ) THEN
00398          CALL XERBLA( 'STGSYL', -INFO )
00399          RETURN
00400       ELSE IF( LQUERY ) THEN
00401          RETURN
00402       END IF
00403 *
00404 *     Quick return if possible
00405 *
00406       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
00407          SCALE = 1
00408          IF( NOTRAN ) THEN
00409             IF( IJOB.NE.0 ) THEN
00410                DIF = 0
00411             END IF
00412          END IF
00413          RETURN
00414       END IF
00415 *
00416 *     Determine optimal block sizes MB and NB
00417 *
00418       MB = ILAENV( 2, 'STGSYL', TRANS, M, N, -1, -1 )
00419       NB = ILAENV( 5, 'STGSYL', TRANS, M, N, -1, -1 )
00420 *
00421       ISOLVE = 1
00422       IFUNC = 0
00423       IF( NOTRAN ) THEN
00424          IF( IJOB.GE.3 ) THEN
00425             IFUNC = IJOB - 2
00426             CALL SLASET( 'F', M, N, ZERO, ZERO, C, LDC )
00427             CALL SLASET( 'F', M, N, ZERO, ZERO, F, LDF )
00428          ELSE IF( IJOB.GE.1 .AND. NOTRAN ) THEN
00429             ISOLVE = 2
00430          END IF
00431       END IF
00432 *
00433       IF( ( MB.LE.1 .AND. NB.LE.1 ) .OR. ( MB.GE.M .AND. NB.GE.N ) )
00434      $     THEN
00435 *
00436          DO 30 IROUND = 1, ISOLVE
00437 *
00438 *           Use unblocked Level 2 solver
00439 *
00440             DSCALE = ZERO
00441             DSUM = ONE
00442             PQ = 0
00443             CALL STGSY2( TRANS, IFUNC, M, N, A, LDA, B, LDB, C, LDC, D,
00444      $                   LDD, E, LDE, F, LDF, SCALE, DSUM, DSCALE,
00445      $                   IWORK, PQ, INFO )
00446             IF( DSCALE.NE.ZERO ) THEN
00447                IF( IJOB.EQ.1 .OR. IJOB.EQ.3 ) THEN
00448                   DIF = SQRT( REAL( 2*M*N ) ) / ( DSCALE*SQRT( DSUM ) )
00449                ELSE
00450                   DIF = SQRT( REAL( PQ ) ) / ( DSCALE*SQRT( DSUM ) )
00451                END IF
00452             END IF
00453 *
00454             IF( ISOLVE.EQ.2 .AND. IROUND.EQ.1 ) THEN
00455                IF( NOTRAN ) THEN
00456                   IFUNC = IJOB
00457                END IF
00458                SCALE2 = SCALE
00459                CALL SLACPY( 'F', M, N, C, LDC, WORK, M )
00460                CALL SLACPY( 'F', M, N, F, LDF, WORK( M*N+1 ), M )
00461                CALL SLASET( 'F', M, N, ZERO, ZERO, C, LDC )
00462                CALL SLASET( 'F', M, N, ZERO, ZERO, F, LDF )
00463             ELSE IF( ISOLVE.EQ.2 .AND. IROUND.EQ.2 ) THEN
00464                CALL SLACPY( 'F', M, N, WORK, M, C, LDC )
00465                CALL SLACPY( 'F', M, N, WORK( M*N+1 ), M, F, LDF )
00466                SCALE = SCALE2
00467             END IF
00468    30    CONTINUE
00469 *
00470          RETURN
00471       END IF
00472 *
00473 *     Determine block structure of A
00474 *
00475       P = 0
00476       I = 1
00477    40 CONTINUE
00478       IF( I.GT.M )
00479      $   GO TO 50
00480       P = P + 1
00481       IWORK( P ) = I
00482       I = I + MB
00483       IF( I.GE.M )
00484      $   GO TO 50
00485       IF( A( I, I-1 ).NE.ZERO )
00486      $   I = I + 1
00487       GO TO 40
00488    50 CONTINUE
00489 *
00490       IWORK( P+1 ) = M + 1
00491       IF( IWORK( P ).EQ.IWORK( P+1 ) )
00492      $   P = P - 1
00493 *
00494 *     Determine block structure of B
00495 *
00496       Q = P + 1
00497       J = 1
00498    60 CONTINUE
00499       IF( J.GT.N )
00500      $   GO TO 70
00501       Q = Q + 1
00502       IWORK( Q ) = J
00503       J = J + NB
00504       IF( J.GE.N )
00505      $   GO TO 70
00506       IF( B( J, J-1 ).NE.ZERO )
00507      $   J = J + 1
00508       GO TO 60
00509    70 CONTINUE
00510 *
00511       IWORK( Q+1 ) = N + 1
00512       IF( IWORK( Q ).EQ.IWORK( Q+1 ) )
00513      $   Q = Q - 1
00514 *
00515       IF( NOTRAN ) THEN
00516 *
00517          DO 150 IROUND = 1, ISOLVE
00518 *
00519 *           Solve (I, J)-subsystem
00520 *               A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
00521 *               D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
00522 *           for I = P, P - 1,..., 1; J = 1, 2,..., Q
00523 *
00524             DSCALE = ZERO
00525             DSUM = ONE
00526             PQ = 0
00527             SCALE = ONE
00528             DO 130 J = P + 2, Q
00529                JS = IWORK( J )
00530                JE = IWORK( J+1 ) - 1
00531                NB = JE - JS + 1
00532                DO 120 I = P, 1, -1
00533                   IS = IWORK( I )
00534                   IE = IWORK( I+1 ) - 1
00535                   MB = IE - IS + 1
00536                   PPQQ = 0
00537                   CALL STGSY2( TRANS, IFUNC, MB, NB, A( IS, IS ), LDA,
00538      $                         B( JS, JS ), LDB, C( IS, JS ), LDC,
00539      $                         D( IS, IS ), LDD, E( JS, JS ), LDE,
00540      $                         F( IS, JS ), LDF, SCALOC, DSUM, DSCALE,
00541      $                         IWORK( Q+2 ), PPQQ, LINFO )
00542                   IF( LINFO.GT.0 )
00543      $               INFO = LINFO
00544 *
00545                   PQ = PQ + PPQQ
00546                   IF( SCALOC.NE.ONE ) THEN
00547                      DO 80 K = 1, JS - 1
00548                         CALL SSCAL( M, SCALOC, C( 1, K ), 1 )
00549                         CALL SSCAL( M, SCALOC, F( 1, K ), 1 )
00550    80                CONTINUE
00551                      DO 90 K = JS, JE
00552                         CALL SSCAL( IS-1, SCALOC, C( 1, K ), 1 )
00553                         CALL SSCAL( IS-1, SCALOC, F( 1, K ), 1 )
00554    90                CONTINUE
00555                      DO 100 K = JS, JE
00556                         CALL SSCAL( M-IE, SCALOC, C( IE+1, K ), 1 )
00557                         CALL SSCAL( M-IE, SCALOC, F( IE+1, K ), 1 )
00558   100                CONTINUE
00559                      DO 110 K = JE + 1, N
00560                         CALL SSCAL( M, SCALOC, C( 1, K ), 1 )
00561                         CALL SSCAL( M, SCALOC, F( 1, K ), 1 )
00562   110                CONTINUE
00563                      SCALE = SCALE*SCALOC
00564                   END IF
00565 *
00566 *                 Substitute R(I, J) and L(I, J) into remaining
00567 *                 equation.
00568 *
00569                   IF( I.GT.1 ) THEN
00570                      CALL SGEMM( 'N', 'N', IS-1, NB, MB, -ONE,
00571      $                           A( 1, IS ), LDA, C( IS, JS ), LDC, ONE,
00572      $                           C( 1, JS ), LDC )
00573                      CALL SGEMM( 'N', 'N', IS-1, NB, MB, -ONE,
00574      $                           D( 1, IS ), LDD, C( IS, JS ), LDC, ONE,
00575      $                           F( 1, JS ), LDF )
00576                   END IF
00577                   IF( J.LT.Q ) THEN
00578                      CALL SGEMM( 'N', 'N', MB, N-JE, NB, ONE,
00579      $                           F( IS, JS ), LDF, B( JS, JE+1 ), LDB,
00580      $                           ONE, C( IS, JE+1 ), LDC )
00581                      CALL SGEMM( 'N', 'N', MB, N-JE, NB, ONE,
00582      $                           F( IS, JS ), LDF, E( JS, JE+1 ), LDE,
00583      $                           ONE, F( IS, JE+1 ), LDF )
00584                   END IF
00585   120          CONTINUE
00586   130       CONTINUE
00587             IF( DSCALE.NE.ZERO ) THEN
00588                IF( IJOB.EQ.1 .OR. IJOB.EQ.3 ) THEN
00589                   DIF = SQRT( REAL( 2*M*N ) ) / ( DSCALE*SQRT( DSUM ) )
00590                ELSE
00591                   DIF = SQRT( REAL( PQ ) ) / ( DSCALE*SQRT( DSUM ) )
00592                END IF
00593             END IF
00594             IF( ISOLVE.EQ.2 .AND. IROUND.EQ.1 ) THEN
00595                IF( NOTRAN ) THEN
00596                   IFUNC = IJOB
00597                END IF
00598                SCALE2 = SCALE
00599                CALL SLACPY( 'F', M, N, C, LDC, WORK, M )
00600                CALL SLACPY( 'F', M, N, F, LDF, WORK( M*N+1 ), M )
00601                CALL SLASET( 'F', M, N, ZERO, ZERO, C, LDC )
00602                CALL SLASET( 'F', M, N, ZERO, ZERO, F, LDF )
00603             ELSE IF( ISOLVE.EQ.2 .AND. IROUND.EQ.2 ) THEN
00604                CALL SLACPY( 'F', M, N, WORK, M, C, LDC )
00605                CALL SLACPY( 'F', M, N, WORK( M*N+1 ), M, F, LDF )
00606                SCALE = SCALE2
00607             END IF
00608   150    CONTINUE
00609 *
00610       ELSE
00611 *
00612 *        Solve transposed (I, J)-subsystem
00613 *             A(I, I)**T * R(I, J)  + D(I, I)**T * L(I, J)  =  C(I, J)
00614 *             R(I, J)  * B(J, J)**T + L(I, J)  * E(J, J)**T = -F(I, J)
00615 *        for I = 1,2,..., P; J = Q, Q-1,..., 1
00616 *
00617          SCALE = ONE
00618          DO 210 I = 1, P
00619             IS = IWORK( I )
00620             IE = IWORK( I+1 ) - 1
00621             MB = IE - IS + 1
00622             DO 200 J = Q, P + 2, -1
00623                JS = IWORK( J )
00624                JE = IWORK( J+1 ) - 1
00625                NB = JE - JS + 1
00626                CALL STGSY2( TRANS, IFUNC, MB, NB, A( IS, IS ), LDA,
00627      $                      B( JS, JS ), LDB, C( IS, JS ), LDC,
00628      $                      D( IS, IS ), LDD, E( JS, JS ), LDE,
00629      $                      F( IS, JS ), LDF, SCALOC, DSUM, DSCALE,
00630      $                      IWORK( Q+2 ), PPQQ, LINFO )
00631                IF( LINFO.GT.0 )
00632      $            INFO = LINFO
00633                IF( SCALOC.NE.ONE ) THEN
00634                   DO 160 K = 1, JS - 1
00635                      CALL SSCAL( M, SCALOC, C( 1, K ), 1 )
00636                      CALL SSCAL( M, SCALOC, F( 1, K ), 1 )
00637   160             CONTINUE
00638                   DO 170 K = JS, JE
00639                      CALL SSCAL( IS-1, SCALOC, C( 1, K ), 1 )
00640                      CALL SSCAL( IS-1, SCALOC, F( 1, K ), 1 )
00641   170             CONTINUE
00642                   DO 180 K = JS, JE
00643                      CALL SSCAL( M-IE, SCALOC, C( IE+1, K ), 1 )
00644                      CALL SSCAL( M-IE, SCALOC, F( IE+1, K ), 1 )
00645   180             CONTINUE
00646                   DO 190 K = JE + 1, N
00647                      CALL SSCAL( M, SCALOC, C( 1, K ), 1 )
00648                      CALL SSCAL( M, SCALOC, F( 1, K ), 1 )
00649   190             CONTINUE
00650                   SCALE = SCALE*SCALOC
00651                END IF
00652 *
00653 *              Substitute R(I, J) and L(I, J) into remaining equation.
00654 *
00655                IF( J.GT.P+2 ) THEN
00656                   CALL SGEMM( 'N', 'T', MB, JS-1, NB, ONE, C( IS, JS ),
00657      $                        LDC, B( 1, JS ), LDB, ONE, F( IS, 1 ),
00658      $                        LDF )
00659                   CALL SGEMM( 'N', 'T', MB, JS-1, NB, ONE, F( IS, JS ),
00660      $                        LDF, E( 1, JS ), LDE, ONE, F( IS, 1 ),
00661      $                        LDF )
00662                END IF
00663                IF( I.LT.P ) THEN
00664                   CALL SGEMM( 'T', 'N', M-IE, NB, MB, -ONE,
00665      $                        A( IS, IE+1 ), LDA, C( IS, JS ), LDC, ONE,
00666      $                        C( IE+1, JS ), LDC )
00667                   CALL SGEMM( 'T', 'N', M-IE, NB, MB, -ONE,
00668      $                        D( IS, IE+1 ), LDD, F( IS, JS ), LDF, ONE,
00669      $                        C( IE+1, JS ), LDC )
00670                END IF
00671   200       CONTINUE
00672   210    CONTINUE
00673 *
00674       END IF
00675 *
00676       WORK( 1 ) = LWMIN
00677 *
00678       RETURN
00679 *
00680 *     End of STGSYL
00681 *
00682       END
 All Files Functions