LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
ctgsy2.f
Go to the documentation of this file.
1*> \brief \b CTGSY2 solves the generalized Sylvester equation (unblocked algorithm).
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ctgsy2.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ctgsy2.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ctgsy2.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE CTGSY2( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
22* LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL,
23* INFO )
24*
25* .. Scalar Arguments ..
26* CHARACTER TRANS
27* INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N
28* REAL RDSCAL, RDSUM, SCALE
29* ..
30* .. Array Arguments ..
31* COMPLEX A( LDA, * ), B( LDB, * ), C( LDC, * ),
32* \$ D( LDD, * ), E( LDE, * ), F( LDF, * )
33* ..
34*
35*
36*> \par Purpose:
37* =============
38*>
39*> \verbatim
40*>
41*> CTGSY2 solves the generalized Sylvester equation
42*>
43*> A * R - L * B = scale * C (1)
44*> D * R - L * E = scale * F
45*>
46*> using Level 1 and 2 BLAS, where R and L are unknown M-by-N matrices,
47*> (A, D), (B, E) and (C, F) are given matrix pairs of size M-by-M,
48*> N-by-N and M-by-N, respectively. A, B, D and E are upper triangular
49*> (i.e., (A,D) and (B,E) in generalized Schur form).
50*>
51*> The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1 is an output
52*> scaling factor chosen to avoid overflow.
53*>
54*> In matrix notation solving equation (1) corresponds to solve
55*> Zx = scale * b, where Z is defined as
56*>
57*> Z = [ kron(In, A) -kron(B**H, Im) ] (2)
58*> [ kron(In, D) -kron(E**H, Im) ],
59*>
60*> Ik is the identity matrix of size k and X**H is the transpose of X.
61*> kron(X, Y) is the Kronecker product between the matrices X and Y.
62*>
63*> If TRANS = 'C', y in the conjugate transposed system Z**H*y = scale*b
64*> is solved for, which is equivalent to solve for R and L in
65*>
66*> A**H * R + D**H * L = scale * C (3)
67*> R * B**H + L * E**H = scale * -F
68*>
69*> This case is used to compute an estimate of Dif[(A, D), (B, E)] =
70*> = sigma_min(Z) using reverse communication with CLACON.
71*>
72*> CTGSY2 also (IJOB >= 1) contributes to the computation in CTGSYL
73*> of an upper bound on the separation between to matrix pairs. Then
74*> the input (A, D), (B, E) are sub-pencils of two matrix pairs in
75*> CTGSYL.
76*> \endverbatim
77*
78* Arguments:
79* ==========
80*
81*> \param[in] TRANS
82*> \verbatim
83*> TRANS is CHARACTER*1
84*> = 'N': solve the generalized Sylvester equation (1).
85*> = 'T': solve the 'transposed' system (3).
86*> \endverbatim
87*>
88*> \param[in] IJOB
89*> \verbatim
90*> IJOB is INTEGER
91*> Specifies what kind of functionality to be performed.
92*> = 0: solve (1) only.
93*> = 1: A contribution from this subsystem to a Frobenius
94*> norm-based estimate of the separation between two matrix
95*> pairs is computed. (look ahead strategy is used).
96*> = 2: A contribution from this subsystem to a Frobenius
97*> norm-based estimate of the separation between two matrix
98*> pairs is computed. (SGECON on sub-systems is used.)
99*> Not referenced if TRANS = 'T'.
100*> \endverbatim
101*>
102*> \param[in] M
103*> \verbatim
104*> M is INTEGER
105*> On entry, M specifies the order of A and D, and the row
106*> dimension of C, F, R and L.
107*> \endverbatim
108*>
109*> \param[in] N
110*> \verbatim
111*> N is INTEGER
112*> On entry, N specifies the order of B and E, and the column
113*> dimension of C, F, R and L.
114*> \endverbatim
115*>
116*> \param[in] A
117*> \verbatim
118*> A is COMPLEX array, dimension (LDA, M)
119*> On entry, A contains an upper triangular matrix.
120*> \endverbatim
121*>
122*> \param[in] LDA
123*> \verbatim
124*> LDA is INTEGER
125*> The leading dimension of the matrix A. LDA >= max(1, M).
126*> \endverbatim
127*>
128*> \param[in] B
129*> \verbatim
130*> B is COMPLEX array, dimension (LDB, N)
131*> On entry, B contains an upper triangular matrix.
132*> \endverbatim
133*>
134*> \param[in] LDB
135*> \verbatim
136*> LDB is INTEGER
137*> The leading dimension of the matrix B. LDB >= max(1, N).
138*> \endverbatim
139*>
140*> \param[in,out] C
141*> \verbatim
142*> C is COMPLEX array, dimension (LDC, N)
143*> On entry, C contains the right-hand-side of the first matrix
144*> equation in (1).
145*> On exit, if IJOB = 0, C has been overwritten by the solution
146*> R.
147*> \endverbatim
148*>
149*> \param[in] LDC
150*> \verbatim
151*> LDC is INTEGER
152*> The leading dimension of the matrix C. LDC >= max(1, M).
153*> \endverbatim
154*>
155*> \param[in] D
156*> \verbatim
157*> D is COMPLEX array, dimension (LDD, M)
158*> On entry, D contains an upper triangular matrix.
159*> \endverbatim
160*>
161*> \param[in] LDD
162*> \verbatim
163*> LDD is INTEGER
164*> The leading dimension of the matrix D. LDD >= max(1, M).
165*> \endverbatim
166*>
167*> \param[in] E
168*> \verbatim
169*> E is COMPLEX array, dimension (LDE, N)
170*> On entry, E contains an upper triangular matrix.
171*> \endverbatim
172*>
173*> \param[in] LDE
174*> \verbatim
175*> LDE is INTEGER
176*> The leading dimension of the matrix E. LDE >= max(1, N).
177*> \endverbatim
178*>
179*> \param[in,out] F
180*> \verbatim
181*> F is COMPLEX array, dimension (LDF, N)
182*> On entry, F contains the right-hand-side of the second matrix
183*> equation in (1).
184*> On exit, if IJOB = 0, F has been overwritten by the solution
185*> L.
186*> \endverbatim
187*>
188*> \param[in] LDF
189*> \verbatim
190*> LDF is INTEGER
191*> The leading dimension of the matrix F. LDF >= max(1, M).
192*> \endverbatim
193*>
194*> \param[out] SCALE
195*> \verbatim
196*> SCALE is REAL
197*> On exit, 0 <= SCALE <= 1. If 0 < SCALE < 1, the solutions
198*> R and L (C and F on entry) will hold the solutions to a
199*> slightly perturbed system but the input matrices A, B, D and
200*> E have not been changed. If SCALE = 0, R and L will hold the
201*> solutions to the homogeneous system with C = F = 0.
202*> Normally, SCALE = 1.
203*> \endverbatim
204*>
205*> \param[in,out] RDSUM
206*> \verbatim
207*> RDSUM is REAL
208*> On entry, the sum of squares of computed contributions to
209*> the Dif-estimate under computation by CTGSYL, where the
210*> scaling factor RDSCAL (see below) has been factored out.
211*> On exit, the corresponding sum of squares updated with the
212*> contributions from the current sub-system.
213*> If TRANS = 'T' RDSUM is not touched.
214*> NOTE: RDSUM only makes sense when CTGSY2 is called by
215*> CTGSYL.
216*> \endverbatim
217*>
218*> \param[in,out] RDSCAL
219*> \verbatim
220*> RDSCAL is REAL
221*> On entry, scaling factor used to prevent overflow in RDSUM.
222*> On exit, RDSCAL is updated w.r.t. the current contributions
223*> in RDSUM.
224*> If TRANS = 'T', RDSCAL is not touched.
225*> NOTE: RDSCAL only makes sense when CTGSY2 is called by
226*> CTGSYL.
227*> \endverbatim
228*>
229*> \param[out] INFO
230*> \verbatim
231*> INFO is INTEGER
232*> On exit, if INFO is set to
233*> =0: Successful exit
234*> <0: If INFO = -i, input argument number i is illegal.
235*> >0: The matrix pairs (A, D) and (B, E) have common or very
236*> close eigenvalues.
237*> \endverbatim
238*
239* Authors:
240* ========
241*
242*> \author Univ. of Tennessee
243*> \author Univ. of California Berkeley
244*> \author Univ. of Colorado Denver
245*> \author NAG Ltd.
246*
247*> \ingroup tgsy2
248*
249*> \par Contributors:
250* ==================
251*>
252*> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
253*> Umea University, S-901 87 Umea, Sweden.
254*
255* =====================================================================
256 SUBROUTINE ctgsy2( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
257 \$ LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL,
258 \$ INFO )
259*
260* -- LAPACK auxiliary routine --
261* -- LAPACK is a software package provided by Univ. of Tennessee, --
262* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
263*
264* .. Scalar Arguments ..
265 CHARACTER TRANS
266 INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N
267 REAL RDSCAL, RDSUM, SCALE
268* ..
269* .. Array Arguments ..
270 COMPLEX A( LDA, * ), B( LDB, * ), C( LDC, * ),
271 \$ D( LDD, * ), E( LDE, * ), F( LDF, * )
272* ..
273*
274* =====================================================================
275*
276* .. Parameters ..
277 REAL ZERO, ONE
278 INTEGER LDZ
279 parameter( zero = 0.0e+0, one = 1.0e+0, ldz = 2 )
280* ..
281* .. Local Scalars ..
282 LOGICAL NOTRAN
283 INTEGER I, IERR, J, K
284 REAL SCALOC
285 COMPLEX ALPHA
286* ..
287* .. Local Arrays ..
288 INTEGER IPIV( LDZ ), JPIV( LDZ )
289 COMPLEX RHS( LDZ ), Z( LDZ, LDZ )
290* ..
291* .. External Functions ..
292 LOGICAL LSAME
293 EXTERNAL LSAME
294* ..
295* .. External Subroutines ..
296 EXTERNAL caxpy, cgesc2, cgetc2, cscal, clatdf, xerbla
297* ..
298* .. Intrinsic Functions ..
299 INTRINSIC cmplx, conjg, max
300* ..
301* .. Executable Statements ..
302*
303* Decode and test input parameters
304*
305 info = 0
306 ierr = 0
307 notran = lsame( trans, 'N' )
308 IF( .NOT.notran .AND. .NOT.lsame( trans, 'C' ) ) THEN
309 info = -1
310 ELSE IF( notran ) THEN
311 IF( ( ijob.LT.0 ) .OR. ( ijob.GT.2 ) ) THEN
312 info = -2
313 END IF
314 END IF
315 IF( info.EQ.0 ) THEN
316 IF( m.LE.0 ) THEN
317 info = -3
318 ELSE IF( n.LE.0 ) THEN
319 info = -4
320 ELSE IF( lda.LT.max( 1, m ) ) THEN
321 info = -6
322 ELSE IF( ldb.LT.max( 1, n ) ) THEN
323 info = -8
324 ELSE IF( ldc.LT.max( 1, m ) ) THEN
325 info = -10
326 ELSE IF( ldd.LT.max( 1, m ) ) THEN
327 info = -12
328 ELSE IF( lde.LT.max( 1, n ) ) THEN
329 info = -14
330 ELSE IF( ldf.LT.max( 1, m ) ) THEN
331 info = -16
332 END IF
333 END IF
334 IF( info.NE.0 ) THEN
335 CALL xerbla( 'CTGSY2', -info )
336 RETURN
337 END IF
338*
339 IF( notran ) THEN
340*
341* Solve (I, J) - system
342* A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
343* D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
344* for I = M, M - 1, ..., 1; J = 1, 2, ..., N
345*
346 scale = one
347 scaloc = one
348 DO 30 j = 1, n
349 DO 20 i = m, 1, -1
350*
351* Build 2 by 2 system
352*
353 z( 1, 1 ) = a( i, i )
354 z( 2, 1 ) = d( i, i )
355 z( 1, 2 ) = -b( j, j )
356 z( 2, 2 ) = -e( j, j )
357*
358* Set up right hand side(s)
359*
360 rhs( 1 ) = c( i, j )
361 rhs( 2 ) = f( i, j )
362*
363* Solve Z * x = RHS
364*
365 CALL cgetc2( ldz, z, ldz, ipiv, jpiv, ierr )
366 IF( ierr.GT.0 )
367 \$ info = ierr
368 IF( ijob.EQ.0 ) THEN
369 CALL cgesc2( ldz, z, ldz, rhs, ipiv, jpiv, scaloc )
370 IF( scaloc.NE.one ) THEN
371 DO 10 k = 1, n
372 CALL cscal( m, cmplx( scaloc, zero ), c( 1, k ),
373 \$ 1 )
374 CALL cscal( m, cmplx( scaloc, zero ), f( 1, k ),
375 \$ 1 )
376 10 CONTINUE
377 scale = scale*scaloc
378 END IF
379 ELSE
380 CALL clatdf( ijob, ldz, z, ldz, rhs, rdsum, rdscal,
381 \$ ipiv, jpiv )
382 END IF
383*
384* Unpack solution vector(s)
385*
386 c( i, j ) = rhs( 1 )
387 f( i, j ) = rhs( 2 )
388*
389* Substitute R(I, J) and L(I, J) into remaining equation.
390*
391 IF( i.GT.1 ) THEN
392 alpha = -rhs( 1 )
393 CALL caxpy( i-1, alpha, a( 1, i ), 1, c( 1, j ), 1 )
394 CALL caxpy( i-1, alpha, d( 1, i ), 1, f( 1, j ), 1 )
395 END IF
396 IF( j.LT.n ) THEN
397 CALL caxpy( n-j, rhs( 2 ), b( j, j+1 ), ldb,
398 \$ c( i, j+1 ), ldc )
399 CALL caxpy( n-j, rhs( 2 ), e( j, j+1 ), lde,
400 \$ f( i, j+1 ), ldf )
401 END IF
402*
403 20 CONTINUE
404 30 CONTINUE
405 ELSE
406*
407* Solve transposed (I, J) - system:
408* A(I, I)**H * R(I, J) + D(I, I)**H * L(J, J) = C(I, J)
409* R(I, I) * B(J, J) + L(I, J) * E(J, J) = -F(I, J)
410* for I = 1, 2, ..., M, J = N, N - 1, ..., 1
411*
412 scale = one
413 scaloc = one
414 DO 80 i = 1, m
415 DO 70 j = n, 1, -1
416*
417* Build 2 by 2 system Z**H
418*
419 z( 1, 1 ) = conjg( a( i, i ) )
420 z( 2, 1 ) = -conjg( b( j, j ) )
421 z( 1, 2 ) = conjg( d( i, i ) )
422 z( 2, 2 ) = -conjg( e( j, j ) )
423*
424*
425* Set up right hand side(s)
426*
427 rhs( 1 ) = c( i, j )
428 rhs( 2 ) = f( i, j )
429*
430* Solve Z**H * x = RHS
431*
432 CALL cgetc2( ldz, z, ldz, ipiv, jpiv, ierr )
433 IF( ierr.GT.0 )
434 \$ info = ierr
435 CALL cgesc2( ldz, z, ldz, rhs, ipiv, jpiv, scaloc )
436 IF( scaloc.NE.one ) THEN
437 DO 40 k = 1, n
438 CALL cscal( m, cmplx( scaloc, zero ), c( 1, k ),
439 \$ 1 )
440 CALL cscal( m, cmplx( scaloc, zero ), f( 1, k ),
441 \$ 1 )
442 40 CONTINUE
443 scale = scale*scaloc
444 END IF
445*
446* Unpack solution vector(s)
447*
448 c( i, j ) = rhs( 1 )
449 f( i, j ) = rhs( 2 )
450*
451* Substitute R(I, J) and L(I, J) into remaining equation.
452*
453 DO 50 k = 1, j - 1
454 f( i, k ) = f( i, k ) + rhs( 1 )*conjg( b( k, j ) ) +
455 \$ rhs( 2 )*conjg( e( k, j ) )
456 50 CONTINUE
457 DO 60 k = i + 1, m
458 c( k, j ) = c( k, j ) - conjg( a( i, k ) )*rhs( 1 ) -
459 \$ conjg( d( i, k ) )*rhs( 2 )
460 60 CONTINUE
461*
462 70 CONTINUE
463 80 CONTINUE
464 END IF
465 RETURN
466*
467* End of CTGSY2
468*
469 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine caxpy(n, ca, cx, incx, cy, incy)
CAXPY
Definition caxpy.f:88
subroutine cgesc2(n, a, lda, rhs, ipiv, jpiv, scale)
CGESC2 solves a system of linear equations using the LU factorization with complete pivoting computed...
Definition cgesc2.f:115
subroutine cgetc2(n, a, lda, ipiv, jpiv, info)
CGETC2 computes the LU factorization with complete pivoting of the general n-by-n matrix.
Definition cgetc2.f:111
subroutine clatdf(ijob, n, z, ldz, rhs, rdsum, rdscal, ipiv, jpiv)
CLATDF uses the LU factorization of the n-by-n matrix computed by sgetc2 and computes a contribution ...
Definition clatdf.f:169
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