![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
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