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