LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
stgsyl.f
Go to the documentation of this file.
1 *> \brief \b STGSYL
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download STGSYL + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/stgsyl.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/stgsyl.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/stgsyl.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE STGSYL( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
22 * LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK,
23 * IWORK, INFO )
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER TRANS
27 * INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF,
28 * $ LWORK, M, N
29 * REAL DIF, SCALE
30 * ..
31 * .. Array Arguments ..
32 * INTEGER IWORK( * )
33 * REAL A( LDA, * ), B( LDB, * ), C( LDC, * ),
34 * $ D( LDD, * ), E( LDE, * ), F( LDF, * ),
35 * $ WORK( * )
36 * ..
37 *
38 *
39 *> \par Purpose:
40 * =============
41 *>
42 *> \verbatim
43 *>
44 *> STGSYL solves the generalized Sylvester equation:
45 *>
46 *> A * R - L * B = scale * C (1)
47 *> D * R - L * E = scale * F
48 *>
49 *> where R and L are unknown m-by-n matrices, (A, D), (B, E) and
50 *> (C, F) are given matrix pairs of size m-by-m, n-by-n and m-by-n,
51 *> respectively, with real entries. (A, D) and (B, E) must be in
52 *> generalized (real) Schur canonical form, i.e. A, B are upper quasi
53 *> triangular and D, E are upper triangular.
54 *>
55 *> The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1 is an output
56 *> scaling factor chosen to avoid overflow.
57 *>
58 *> In matrix notation (1) is equivalent to solve Zx = scale b, where
59 *> Z is defined as
60 *>
61 *> Z = [ kron(In, A) -kron(B**T, Im) ] (2)
62 *> [ kron(In, D) -kron(E**T, Im) ].
63 *>
64 *> Here Ik is the identity matrix of size k and X**T is the transpose of
65 *> X. kron(X, Y) is the Kronecker product between the matrices X and Y.
66 *>
67 *> If TRANS = 'T', STGSYL solves the transposed system Z**T*y = scale*b,
68 *> which is equivalent to solve for R and L in
69 *>
70 *> A**T * R + D**T * L = scale * C (3)
71 *> R * B**T + L * E**T = scale * -F
72 *>
73 *> This case (TRANS = 'T') is used to compute an one-norm-based estimate
74 *> of Dif[(A,D), (B,E)], the separation between the matrix pairs (A,D)
75 *> and (B,E), using SLACON.
76 *>
77 *> If IJOB >= 1, STGSYL computes a Frobenius norm-based estimate
78 *> of Dif[(A,D),(B,E)]. That is, the reciprocal of a lower bound on the
79 *> reciprocal of the smallest singular value of Z. See [1-2] for more
80 *> information.
81 *>
82 *> This is a level 3 BLAS algorithm.
83 *> \endverbatim
84 *
85 * Arguments:
86 * ==========
87 *
88 *> \param[in] TRANS
89 *> \verbatim
90 *> TRANS is CHARACTER*1
91 *> = 'N', solve the generalized Sylvester equation (1).
92 *> = 'T', solve the 'transposed' system (3).
93 *> \endverbatim
94 *>
95 *> \param[in] IJOB
96 *> \verbatim
97 *> IJOB is INTEGER
98 *> Specifies what kind of functionality to be performed.
99 *> =0: solve (1) only.
100 *> =1: The functionality of 0 and 3.
101 *> =2: The functionality of 0 and 4.
102 *> =3: Only an estimate of Dif[(A,D), (B,E)] is computed.
103 *> (look ahead strategy IJOB = 1 is used).
104 *> =4: Only an estimate of Dif[(A,D), (B,E)] is computed.
105 *> ( SGECON on sub-systems is used ).
106 *> Not referenced if TRANS = 'T'.
107 *> \endverbatim
108 *>
109 *> \param[in] M
110 *> \verbatim
111 *> M is INTEGER
112 *> The order of the matrices A and D, and the row dimension of
113 *> the matrices C, F, R and L.
114 *> \endverbatim
115 *>
116 *> \param[in] N
117 *> \verbatim
118 *> N is INTEGER
119 *> The order of the matrices B and E, and the column dimension
120 *> of the matrices C, F, R and L.
121 *> \endverbatim
122 *>
123 *> \param[in] A
124 *> \verbatim
125 *> A is REAL array, dimension (LDA, M)
126 *> The upper quasi triangular matrix A.
127 *> \endverbatim
128 *>
129 *> \param[in] LDA
130 *> \verbatim
131 *> LDA is INTEGER
132 *> The leading dimension of the array A. LDA >= max(1, M).
133 *> \endverbatim
134 *>
135 *> \param[in] B
136 *> \verbatim
137 *> B is REAL array, dimension (LDB, N)
138 *> The upper quasi triangular matrix B.
139 *> \endverbatim
140 *>
141 *> \param[in] LDB
142 *> \verbatim
143 *> LDB is INTEGER
144 *> The leading dimension of the array B. LDB >= max(1, N).
145 *> \endverbatim
146 *>
147 *> \param[in,out] C
148 *> \verbatim
149 *> C is REAL array, dimension (LDC, N)
150 *> On entry, C contains the right-hand-side of the first matrix
151 *> equation in (1) or (3).
152 *> On exit, if IJOB = 0, 1 or 2, C has been overwritten by
153 *> the solution R. If IJOB = 3 or 4 and TRANS = 'N', C holds R,
154 *> the solution achieved during the computation of the
155 *> Dif-estimate.
156 *> \endverbatim
157 *>
158 *> \param[in] LDC
159 *> \verbatim
160 *> LDC is INTEGER
161 *> The leading dimension of the array C. LDC >= max(1, M).
162 *> \endverbatim
163 *>
164 *> \param[in] D
165 *> \verbatim
166 *> D is REAL array, dimension (LDD, M)
167 *> The upper triangular matrix D.
168 *> \endverbatim
169 *>
170 *> \param[in] LDD
171 *> \verbatim
172 *> LDD is INTEGER
173 *> The leading dimension of the array D. LDD >= max(1, M).
174 *> \endverbatim
175 *>
176 *> \param[in] E
177 *> \verbatim
178 *> E is REAL array, dimension (LDE, N)
179 *> The upper triangular matrix E.
180 *> \endverbatim
181 *>
182 *> \param[in] LDE
183 *> \verbatim
184 *> LDE is INTEGER
185 *> The leading dimension of the array E. LDE >= max(1, N).
186 *> \endverbatim
187 *>
188 *> \param[in,out] F
189 *> \verbatim
190 *> F is REAL array, dimension (LDF, N)
191 *> On entry, F contains the right-hand-side of the second matrix
192 *> equation in (1) or (3).
193 *> On exit, if IJOB = 0, 1 or 2, F has been overwritten by
194 *> the solution L. If IJOB = 3 or 4 and TRANS = 'N', F holds L,
195 *> the solution achieved during the computation of the
196 *> Dif-estimate.
197 *> \endverbatim
198 *>
199 *> \param[in] LDF
200 *> \verbatim
201 *> LDF is INTEGER
202 *> The leading dimension of the array F. LDF >= max(1, M).
203 *> \endverbatim
204 *>
205 *> \param[out] DIF
206 *> \verbatim
207 *> DIF is REAL
208 *> On exit DIF is the reciprocal of a lower bound of the
209 *> reciprocal of the Dif-function, i.e. DIF is an upper bound of
210 *> Dif[(A,D), (B,E)] = sigma_min(Z), where Z as in (2).
211 *> IF IJOB = 0 or TRANS = 'T', DIF is not touched.
212 *> \endverbatim
213 *>
214 *> \param[out] SCALE
215 *> \verbatim
216 *> SCALE is REAL
217 *> On exit SCALE is the scaling factor in (1) or (3).
218 *> If 0 < SCALE < 1, C and F hold the solutions R and L, resp.,
219 *> to a slightly perturbed system but the input matrices A, B, D
220 *> and E have not been changed. If SCALE = 0, C and F hold the
221 *> solutions R and L, respectively, to the homogeneous system
222 *> with C = F = 0. Normally, SCALE = 1.
223 *> \endverbatim
224 *>
225 *> \param[out] WORK
226 *> \verbatim
227 *> WORK is REAL array, dimension (MAX(1,LWORK))
228 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
229 *> \endverbatim
230 *>
231 *> \param[in] LWORK
232 *> \verbatim
233 *> LWORK is INTEGER
234 *> The dimension of the array WORK. LWORK > = 1.
235 *> If IJOB = 1 or 2 and TRANS = 'N', LWORK >= max(1,2*M*N).
236 *>
237 *> If LWORK = -1, then a workspace query is assumed; the routine
238 *> only calculates the optimal size of the WORK array, returns
239 *> this value as the first entry of the WORK array, and no error
240 *> message related to LWORK is issued by XERBLA.
241 *> \endverbatim
242 *>
243 *> \param[out] IWORK
244 *> \verbatim
245 *> IWORK is INTEGER array, dimension (M+N+6)
246 *> \endverbatim
247 *>
248 *> \param[out] INFO
249 *> \verbatim
250 *> INFO is INTEGER
251 *> =0: successful exit
252 *> <0: If INFO = -i, the i-th argument had an illegal value.
253 *> >0: (A, D) and (B, E) have common or close eigenvalues.
254 *> \endverbatim
255 *
256 * Authors:
257 * ========
258 *
259 *> \author Univ. of Tennessee
260 *> \author Univ. of California Berkeley
261 *> \author Univ. of Colorado Denver
262 *> \author NAG Ltd.
263 *
264 *> \date November 2011
265 *
266 *> \ingroup realSYcomputational
267 *
268 *> \par Contributors:
269 * ==================
270 *>
271 *> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
272 *> Umea University, S-901 87 Umea, Sweden.
273 *
274 *> \par References:
275 * ================
276 *>
277 *> \verbatim
278 *>
279 *> [1] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software
280 *> for Solving the Generalized Sylvester Equation and Estimating the
281 *> Separation between Regular Matrix Pairs, Report UMINF - 93.23,
282 *> Department of Computing Science, Umea University, S-901 87 Umea,
283 *> Sweden, December 1993, Revised April 1994, Also as LAPACK Working
284 *> Note 75. To appear in ACM Trans. on Math. Software, Vol 22,
285 *> No 1, 1996.
286 *>
287 *> [2] B. Kagstrom, A Perturbation Analysis of the Generalized Sylvester
288 *> Equation (AR - LB, DR - LE ) = (C, F), SIAM J. Matrix Anal.
289 *> Appl., 15(4):1045-1060, 1994
290 *>
291 *> [3] B. Kagstrom and L. Westin, Generalized Schur Methods with
292 *> Condition Estimators for Solving the Generalized Sylvester
293 *> Equation, IEEE Transactions on Automatic Control, Vol. 34, No. 7,
294 *> July 1989, pp 745-751.
295 *> \endverbatim
296 *>
297 * =====================================================================
298  SUBROUTINE stgsyl( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
299  $ ldd, e, lde, f, ldf, scale, dif, work, lwork,
300  $ iwork, info )
301 *
302 * -- LAPACK computational routine (version 3.4.0) --
303 * -- LAPACK is a software package provided by Univ. of Tennessee, --
304 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
305 * November 2011
306 *
307 * .. Scalar Arguments ..
308  CHARACTER trans
309  INTEGER ijob, info, lda, ldb, ldc, ldd, lde, ldf,
310  $ lwork, m, n
311  REAL dif, scale
312 * ..
313 * .. Array Arguments ..
314  INTEGER iwork( * )
315  REAL a( lda, * ), b( ldb, * ), c( ldc, * ),
316  $ d( ldd, * ), e( lde, * ), f( ldf, * ),
317  $ work( * )
318 * ..
319 *
320 * =====================================================================
321 * Replaced various illegal calls to SCOPY by calls to SLASET.
322 * Sven Hammarling, 1/5/02.
323 *
324 * .. Parameters ..
325  REAL zero, one
326  parameter( zero = 0.0e+0, one = 1.0e+0 )
327 * ..
328 * .. Local Scalars ..
329  LOGICAL lquery, notran
330  INTEGER i, ie, ifunc, iround, is, isolve, j, je, js, k,
331  $ linfo, lwmin, mb, nb, p, ppqq, pq, q
332  REAL dscale, dsum, scale2, scaloc
333 * ..
334 * .. External Functions ..
335  LOGICAL lsame
336  INTEGER ilaenv
337  EXTERNAL lsame, ilaenv
338 * ..
339 * .. External Subroutines ..
340  EXTERNAL sgemm, slacpy, slaset, sscal, stgsy2, xerbla
341 * ..
342 * .. Intrinsic Functions ..
343  INTRINSIC max, REAL, sqrt
344 * ..
345 * .. Executable Statements ..
346 *
347 * Decode and test input parameters
348 *
349  info = 0
350  notran = lsame( trans, 'N' )
351  lquery = ( lwork.EQ.-1 )
352 *
353  IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) ) THEN
354  info = -1
355  ELSE IF( notran ) THEN
356  IF( ( ijob.LT.0 ) .OR. ( ijob.GT.4 ) ) THEN
357  info = -2
358  END IF
359  END IF
360  IF( info.EQ.0 ) THEN
361  IF( m.LE.0 ) THEN
362  info = -3
363  ELSE IF( n.LE.0 ) THEN
364  info = -4
365  ELSE IF( lda.LT.max( 1, m ) ) THEN
366  info = -6
367  ELSE IF( ldb.LT.max( 1, n ) ) THEN
368  info = -8
369  ELSE IF( ldc.LT.max( 1, m ) ) THEN
370  info = -10
371  ELSE IF( ldd.LT.max( 1, m ) ) THEN
372  info = -12
373  ELSE IF( lde.LT.max( 1, n ) ) THEN
374  info = -14
375  ELSE IF( ldf.LT.max( 1, m ) ) THEN
376  info = -16
377  END IF
378  END IF
379 *
380  IF( info.EQ.0 ) THEN
381  IF( notran ) THEN
382  IF( ijob.EQ.1 .OR. ijob.EQ.2 ) THEN
383  lwmin = max( 1, 2*m*n )
384  ELSE
385  lwmin = 1
386  END IF
387  ELSE
388  lwmin = 1
389  END IF
390  work( 1 ) = lwmin
391 *
392  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
393  info = -20
394  END IF
395  END IF
396 *
397  IF( info.NE.0 ) THEN
398  CALL xerbla( 'STGSYL', -info )
399  RETURN
400  ELSE IF( lquery ) THEN
401  RETURN
402  END IF
403 *
404 * Quick return if possible
405 *
406  IF( m.EQ.0 .OR. n.EQ.0 ) THEN
407  scale = 1
408  IF( notran ) THEN
409  IF( ijob.NE.0 ) THEN
410  dif = 0
411  END IF
412  END IF
413  RETURN
414  END IF
415 *
416 * Determine optimal block sizes MB and NB
417 *
418  mb = ilaenv( 2, 'STGSYL', trans, m, n, -1, -1 )
419  nb = ilaenv( 5, 'STGSYL', trans, m, n, -1, -1 )
420 *
421  isolve = 1
422  ifunc = 0
423  IF( notran ) THEN
424  IF( ijob.GE.3 ) THEN
425  ifunc = ijob - 2
426  CALL slaset( 'F', m, n, zero, zero, c, ldc )
427  CALL slaset( 'F', m, n, zero, zero, f, ldf )
428  ELSE IF( ijob.GE.1 .AND. notran ) THEN
429  isolve = 2
430  END IF
431  END IF
432 *
433  IF( ( mb.LE.1 .AND. nb.LE.1 ) .OR. ( mb.GE.m .AND. nb.GE.n ) )
434  $ THEN
435 *
436  DO 30 iround = 1, isolve
437 *
438 * Use unblocked Level 2 solver
439 *
440  dscale = zero
441  dsum = one
442  pq = 0
443  CALL stgsy2( trans, ifunc, m, n, a, lda, b, ldb, c, ldc, d,
444  $ ldd, e, lde, f, ldf, scale, dsum, dscale,
445  $ iwork, pq, info )
446  IF( dscale.NE.zero ) THEN
447  IF( ijob.EQ.1 .OR. ijob.EQ.3 ) THEN
448  dif = sqrt( REAL( 2*M*N ) ) / ( dscale*sqrt( dsum ) )
449  ELSE
450  dif = sqrt( REAL( PQ ) ) / ( dscale*sqrt( dsum ) )
451  END IF
452  END IF
453 *
454  IF( isolve.EQ.2 .AND. iround.EQ.1 ) THEN
455  IF( notran ) THEN
456  ifunc = ijob
457  END IF
458  scale2 = scale
459  CALL slacpy( 'F', m, n, c, ldc, work, m )
460  CALL slacpy( 'F', m, n, f, ldf, work( m*n+1 ), m )
461  CALL slaset( 'F', m, n, zero, zero, c, ldc )
462  CALL slaset( 'F', m, n, zero, zero, f, ldf )
463  ELSE IF( isolve.EQ.2 .AND. iround.EQ.2 ) THEN
464  CALL slacpy( 'F', m, n, work, m, c, ldc )
465  CALL slacpy( 'F', m, n, work( m*n+1 ), m, f, ldf )
466  scale = scale2
467  END IF
468  30 CONTINUE
469 *
470  RETURN
471  END IF
472 *
473 * Determine block structure of A
474 *
475  p = 0
476  i = 1
477  40 CONTINUE
478  IF( i.GT.m )
479  $ go to 50
480  p = p + 1
481  iwork( p ) = i
482  i = i + mb
483  IF( i.GE.m )
484  $ go to 50
485  IF( a( i, i-1 ).NE.zero )
486  $ i = i + 1
487  go to 40
488  50 CONTINUE
489 *
490  iwork( p+1 ) = m + 1
491  IF( iwork( p ).EQ.iwork( p+1 ) )
492  $ p = p - 1
493 *
494 * Determine block structure of B
495 *
496  q = p + 1
497  j = 1
498  60 CONTINUE
499  IF( j.GT.n )
500  $ go to 70
501  q = q + 1
502  iwork( q ) = j
503  j = j + nb
504  IF( j.GE.n )
505  $ go to 70
506  IF( b( j, j-1 ).NE.zero )
507  $ j = j + 1
508  go to 60
509  70 CONTINUE
510 *
511  iwork( q+1 ) = n + 1
512  IF( iwork( q ).EQ.iwork( q+1 ) )
513  $ q = q - 1
514 *
515  IF( notran ) THEN
516 *
517  DO 150 iround = 1, isolve
518 *
519 * Solve (I, J)-subsystem
520 * A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
521 * D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
522 * for I = P, P - 1,..., 1; J = 1, 2,..., Q
523 *
524  dscale = zero
525  dsum = one
526  pq = 0
527  scale = one
528  DO 130 j = p + 2, q
529  js = iwork( j )
530  je = iwork( j+1 ) - 1
531  nb = je - js + 1
532  DO 120 i = p, 1, -1
533  is = iwork( i )
534  ie = iwork( i+1 ) - 1
535  mb = ie - is + 1
536  ppqq = 0
537  CALL stgsy2( trans, ifunc, mb, nb, a( is, is ), lda,
538  $ b( js, js ), ldb, c( is, js ), ldc,
539  $ d( is, is ), ldd, e( js, js ), lde,
540  $ f( is, js ), ldf, scaloc, dsum, dscale,
541  $ iwork( q+2 ), ppqq, linfo )
542  IF( linfo.GT.0 )
543  $ info = linfo
544 *
545  pq = pq + ppqq
546  IF( scaloc.NE.one ) THEN
547  DO 80 k = 1, js - 1
548  CALL sscal( m, scaloc, c( 1, k ), 1 )
549  CALL sscal( m, scaloc, f( 1, k ), 1 )
550  80 CONTINUE
551  DO 90 k = js, je
552  CALL sscal( is-1, scaloc, c( 1, k ), 1 )
553  CALL sscal( is-1, scaloc, f( 1, k ), 1 )
554  90 CONTINUE
555  DO 100 k = js, je
556  CALL sscal( m-ie, scaloc, c( ie+1, k ), 1 )
557  CALL sscal( m-ie, scaloc, f( ie+1, k ), 1 )
558  100 CONTINUE
559  DO 110 k = je + 1, n
560  CALL sscal( m, scaloc, c( 1, k ), 1 )
561  CALL sscal( m, scaloc, f( 1, k ), 1 )
562  110 CONTINUE
563  scale = scale*scaloc
564  END IF
565 *
566 * Substitute R(I, J) and L(I, J) into remaining
567 * equation.
568 *
569  IF( i.GT.1 ) THEN
570  CALL sgemm( 'N', 'N', is-1, nb, mb, -one,
571  $ a( 1, is ), lda, c( is, js ), ldc, one,
572  $ c( 1, js ), ldc )
573  CALL sgemm( 'N', 'N', is-1, nb, mb, -one,
574  $ d( 1, is ), ldd, c( is, js ), ldc, one,
575  $ f( 1, js ), ldf )
576  END IF
577  IF( j.LT.q ) THEN
578  CALL sgemm( 'N', 'N', mb, n-je, nb, one,
579  $ f( is, js ), ldf, b( js, je+1 ), ldb,
580  $ one, c( is, je+1 ), ldc )
581  CALL sgemm( 'N', 'N', mb, n-je, nb, one,
582  $ f( is, js ), ldf, e( js, je+1 ), lde,
583  $ one, f( is, je+1 ), ldf )
584  END IF
585  120 CONTINUE
586  130 CONTINUE
587  IF( dscale.NE.zero ) THEN
588  IF( ijob.EQ.1 .OR. ijob.EQ.3 ) THEN
589  dif = sqrt( REAL( 2*M*N ) ) / ( dscale*sqrt( dsum ) )
590  ELSE
591  dif = sqrt( REAL( PQ ) ) / ( dscale*sqrt( dsum ) )
592  END IF
593  END IF
594  IF( isolve.EQ.2 .AND. iround.EQ.1 ) THEN
595  IF( notran ) THEN
596  ifunc = ijob
597  END IF
598  scale2 = scale
599  CALL slacpy( 'F', m, n, c, ldc, work, m )
600  CALL slacpy( 'F', m, n, f, ldf, work( m*n+1 ), m )
601  CALL slaset( 'F', m, n, zero, zero, c, ldc )
602  CALL slaset( 'F', m, n, zero, zero, f, ldf )
603  ELSE IF( isolve.EQ.2 .AND. iround.EQ.2 ) THEN
604  CALL slacpy( 'F', m, n, work, m, c, ldc )
605  CALL slacpy( 'F', m, n, work( m*n+1 ), m, f, ldf )
606  scale = scale2
607  END IF
608  150 CONTINUE
609 *
610  ELSE
611 *
612 * Solve transposed (I, J)-subsystem
613 * A(I, I)**T * R(I, J) + D(I, I)**T * L(I, J) = C(I, J)
614 * R(I, J) * B(J, J)**T + L(I, J) * E(J, J)**T = -F(I, J)
615 * for I = 1,2,..., P; J = Q, Q-1,..., 1
616 *
617  scale = one
618  DO 210 i = 1, p
619  is = iwork( i )
620  ie = iwork( i+1 ) - 1
621  mb = ie - is + 1
622  DO 200 j = q, p + 2, -1
623  js = iwork( j )
624  je = iwork( j+1 ) - 1
625  nb = je - js + 1
626  CALL stgsy2( trans, ifunc, mb, nb, a( is, is ), lda,
627  $ b( js, js ), ldb, c( is, js ), ldc,
628  $ d( is, is ), ldd, e( js, js ), lde,
629  $ f( is, js ), ldf, scaloc, dsum, dscale,
630  $ iwork( q+2 ), ppqq, linfo )
631  IF( linfo.GT.0 )
632  $ info = linfo
633  IF( scaloc.NE.one ) THEN
634  DO 160 k = 1, js - 1
635  CALL sscal( m, scaloc, c( 1, k ), 1 )
636  CALL sscal( m, scaloc, f( 1, k ), 1 )
637  160 CONTINUE
638  DO 170 k = js, je
639  CALL sscal( is-1, scaloc, c( 1, k ), 1 )
640  CALL sscal( is-1, scaloc, f( 1, k ), 1 )
641  170 CONTINUE
642  DO 180 k = js, je
643  CALL sscal( m-ie, scaloc, c( ie+1, k ), 1 )
644  CALL sscal( m-ie, scaloc, f( ie+1, k ), 1 )
645  180 CONTINUE
646  DO 190 k = je + 1, n
647  CALL sscal( m, scaloc, c( 1, k ), 1 )
648  CALL sscal( m, scaloc, f( 1, k ), 1 )
649  190 CONTINUE
650  scale = scale*scaloc
651  END IF
652 *
653 * Substitute R(I, J) and L(I, J) into remaining equation.
654 *
655  IF( j.GT.p+2 ) THEN
656  CALL sgemm( 'N', 'T', mb, js-1, nb, one, c( is, js ),
657  $ ldc, b( 1, js ), ldb, one, f( is, 1 ),
658  $ ldf )
659  CALL sgemm( 'N', 'T', mb, js-1, nb, one, f( is, js ),
660  $ ldf, e( 1, js ), lde, one, f( is, 1 ),
661  $ ldf )
662  END IF
663  IF( i.LT.p ) THEN
664  CALL sgemm( 'T', 'N', m-ie, nb, mb, -one,
665  $ a( is, ie+1 ), lda, c( is, js ), ldc, one,
666  $ c( ie+1, js ), ldc )
667  CALL sgemm( 'T', 'N', m-ie, nb, mb, -one,
668  $ d( is, ie+1 ), ldd, f( is, js ), ldf, one,
669  $ c( ie+1, js ), ldc )
670  END IF
671  200 CONTINUE
672  210 CONTINUE
673 *
674  END IF
675 *
676  work( 1 ) = lwmin
677 *
678  RETURN
679 *
680 * End of STGSYL
681 *
682  END