LAPACK  3.6.0
LAPACK: Linear Algebra PACKage
dtgsyl.f
Go to the documentation of this file.
1 *> \brief \b DTGSYL
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DTGSYL + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtgsyl.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtgsyl.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtgsyl.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DTGSYL( 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 * DOUBLE PRECISION DIF, SCALE
30 * ..
31 * .. Array Arguments ..
32 * INTEGER IWORK( * )
33 * DOUBLE PRECISION 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 *> DTGSYL 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', DTGSYL 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 DLACON.
76 *>
77 *> If IJOB >= 1, DTGSYL 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 *> ( DGECON 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
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 DOUBLE PRECISION
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 DOUBLE PRECISION 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 doubleSYcomputational
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 dtgsyl( 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  DOUBLE PRECISION DIF, SCALE
312 * ..
313 * .. Array Arguments ..
314  INTEGER IWORK( * )
315  DOUBLE PRECISION A( lda, * ), B( ldb, * ), C( ldc, * ),
316  $ d( ldd, * ), e( lde, * ), f( ldf, * ),
317  $ work( * )
318 * ..
319 *
320 * =====================================================================
321 * Replaced various illegal calls to DCOPY by calls to DLASET.
322 * Sven Hammarling, 1/5/02.
323 *
324 * .. Parameters ..
325  DOUBLE PRECISION ZERO, ONE
326  parameter( zero = 0.0d+0, one = 1.0d+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  DOUBLE PRECISION DSCALE, DSUM, SCALE2, SCALOC
333 * ..
334 * .. External Functions ..
335  LOGICAL LSAME
336  INTEGER ILAENV
337  EXTERNAL lsame, ilaenv
338 * ..
339 * .. External Subroutines ..
340  EXTERNAL dgemm, dlacpy, dlaset, dscal, dtgsy2, xerbla
341 * ..
342 * .. Intrinsic Functions ..
343  INTRINSIC dble, max, 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( 'DTGSYL', -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, 'DTGSYL', trans, m, n, -1, -1 )
419  nb = ilaenv( 5, 'DTGSYL', 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 dlaset( 'F', m, n, zero, zero, c, ldc )
427  CALL dlaset( 'F', m, n, zero, zero, f, ldf )
428  ELSE IF( ijob.GE.1 ) 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 dtgsy2( 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( dble( 2*m*n ) ) / ( dscale*sqrt( dsum ) )
449  ELSE
450  dif = sqrt( dble( 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 dlacpy( 'F', m, n, c, ldc, work, m )
460  CALL dlacpy( 'F', m, n, f, ldf, work( m*n+1 ), m )
461  CALL dlaset( 'F', m, n, zero, zero, c, ldc )
462  CALL dlaset( 'F', m, n, zero, zero, f, ldf )
463  ELSE IF( isolve.EQ.2 .AND. iround.EQ.2 ) THEN
464  CALL dlacpy( 'F', m, n, work, m, c, ldc )
465  CALL dlacpy( '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 dtgsy2( 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 dscal( m, scaloc, c( 1, k ), 1 )
549  CALL dscal( m, scaloc, f( 1, k ), 1 )
550  80 CONTINUE
551  DO 90 k = js, je
552  CALL dscal( is-1, scaloc, c( 1, k ), 1 )
553  CALL dscal( is-1, scaloc, f( 1, k ), 1 )
554  90 CONTINUE
555  DO 100 k = js, je
556  CALL dscal( m-ie, scaloc, c( ie+1, k ), 1 )
557  CALL dscal( m-ie, scaloc, f( ie+1, k ), 1 )
558  100 CONTINUE
559  DO 110 k = je + 1, n
560  CALL dscal( m, scaloc, c( 1, k ), 1 )
561  CALL dscal( 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 dgemm( 'N', 'N', is-1, nb, mb, -one,
571  $ a( 1, is ), lda, c( is, js ), ldc, one,
572  $ c( 1, js ), ldc )
573  CALL dgemm( '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 dgemm( '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 dgemm( '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( dble( 2*m*n ) ) / ( dscale*sqrt( dsum ) )
590  ELSE
591  dif = sqrt( dble( 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 dlacpy( 'F', m, n, c, ldc, work, m )
600  CALL dlacpy( 'F', m, n, f, ldf, work( m*n+1 ), m )
601  CALL dlaset( 'F', m, n, zero, zero, c, ldc )
602  CALL dlaset( 'F', m, n, zero, zero, f, ldf )
603  ELSE IF( isolve.EQ.2 .AND. iround.EQ.2 ) THEN
604  CALL dlacpy( 'F', m, n, work, m, c, ldc )
605  CALL dlacpy( '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 dtgsy2( 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 dscal( m, scaloc, c( 1, k ), 1 )
636  CALL dscal( m, scaloc, f( 1, k ), 1 )
637  160 CONTINUE
638  DO 170 k = js, je
639  CALL dscal( is-1, scaloc, c( 1, k ), 1 )
640  CALL dscal( is-1, scaloc, f( 1, k ), 1 )
641  170 CONTINUE
642  DO 180 k = js, je
643  CALL dscal( m-ie, scaloc, c( ie+1, k ), 1 )
644  CALL dscal( m-ie, scaloc, f( ie+1, k ), 1 )
645  180 CONTINUE
646  DO 190 k = je + 1, n
647  CALL dscal( m, scaloc, c( 1, k ), 1 )
648  CALL dscal( 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 dgemm( 'N', 'T', mb, js-1, nb, one, c( is, js ),
657  $ ldc, b( 1, js ), ldb, one, f( is, 1 ),
658  $ ldf )
659  CALL dgemm( '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 dgemm( '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 dgemm( '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 DTGSYL
681 *
682  END
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: dlaset.f:112
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:105
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:189
subroutine dtgsyl(TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK, IWORK, INFO)
DTGSYL
Definition: dtgsyl.f:301
subroutine dtgsy2(TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL, IWORK, PQ, INFO)
DTGSY2 solves the generalized Sylvester equation (unblocked algorithm).
Definition: dtgsy2.f:276
logical function lde(RI, RJ, LR)
Definition: dblat2.f:2945
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55