LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
ztrsyl.f
Go to the documentation of this file.
1*> \brief \b ZTRSYL
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download ZTRSYL + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztrsyl.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztrsyl.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztrsyl.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE ZTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C,
20* LDC, SCALE, INFO )
21*
22* .. Scalar Arguments ..
23* CHARACTER TRANA, TRANB
24* INTEGER INFO, ISGN, LDA, LDB, LDC, M, N
25* DOUBLE PRECISION SCALE
26* ..
27* .. Array Arguments ..
28* COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> ZTRSYL solves the complex Sylvester matrix equation:
38*>
39*> op(A)*X + X*op(B) = scale*C or
40*> op(A)*X - X*op(B) = scale*C,
41*>
42*> where op(A) = A or A**H, and A and B are both upper triangular. A is
43*> M-by-M and B is N-by-N; the right hand side C and the solution X are
44*> M-by-N; and scale is an output scale factor, set <= 1 to avoid
45*> overflow in X.
46*> \endverbatim
47*
48* Arguments:
49* ==========
50*
51*> \param[in] TRANA
52*> \verbatim
53*> TRANA is CHARACTER*1
54*> Specifies the option op(A):
55*> = 'N': op(A) = A (No transpose)
56*> = 'C': op(A) = A**H (Conjugate transpose)
57*> \endverbatim
58*>
59*> \param[in] TRANB
60*> \verbatim
61*> TRANB is CHARACTER*1
62*> Specifies the option op(B):
63*> = 'N': op(B) = B (No transpose)
64*> = 'C': op(B) = B**H (Conjugate transpose)
65*> \endverbatim
66*>
67*> \param[in] ISGN
68*> \verbatim
69*> ISGN is INTEGER
70*> Specifies the sign in the equation:
71*> = +1: solve op(A)*X + X*op(B) = scale*C
72*> = -1: solve op(A)*X - X*op(B) = scale*C
73*> \endverbatim
74*>
75*> \param[in] M
76*> \verbatim
77*> M is INTEGER
78*> The order of the matrix A, and the number of rows in the
79*> matrices X and C. M >= 0.
80*> \endverbatim
81*>
82*> \param[in] N
83*> \verbatim
84*> N is INTEGER
85*> The order of the matrix B, and the number of columns in the
86*> matrices X and C. N >= 0.
87*> \endverbatim
88*>
89*> \param[in] A
90*> \verbatim
91*> A is COMPLEX*16 array, dimension (LDA,M)
92*> The upper triangular matrix A.
93*> \endverbatim
94*>
95*> \param[in] LDA
96*> \verbatim
97*> LDA is INTEGER
98*> The leading dimension of the array A. LDA >= max(1,M).
99*> \endverbatim
100*>
101*> \param[in] B
102*> \verbatim
103*> B is COMPLEX*16 array, dimension (LDB,N)
104*> The upper triangular matrix B.
105*> \endverbatim
106*>
107*> \param[in] LDB
108*> \verbatim
109*> LDB is INTEGER
110*> The leading dimension of the array B. LDB >= max(1,N).
111*> \endverbatim
112*>
113*> \param[in,out] C
114*> \verbatim
115*> C is COMPLEX*16 array, dimension (LDC,N)
116*> On entry, the M-by-N right hand side matrix C.
117*> On exit, C is overwritten by the solution matrix X.
118*> \endverbatim
119*>
120*> \param[in] LDC
121*> \verbatim
122*> LDC is INTEGER
123*> The leading dimension of the array C. LDC >= max(1,M)
124*> \endverbatim
125*>
126*> \param[out] SCALE
127*> \verbatim
128*> SCALE is DOUBLE PRECISION
129*> The scale factor, scale, set <= 1 to avoid overflow in X.
130*> \endverbatim
131*>
132*> \param[out] INFO
133*> \verbatim
134*> INFO is INTEGER
135*> = 0: successful exit
136*> < 0: if INFO = -i, the i-th argument had an illegal value
137*> = 1: A and B have common or very close eigenvalues; perturbed
138*> values were used to solve the equation (but the matrices
139*> A and B are unchanged).
140*> \endverbatim
141*
142* Authors:
143* ========
144*
145*> \author Univ. of Tennessee
146*> \author Univ. of California Berkeley
147*> \author Univ. of Colorado Denver
148*> \author NAG Ltd.
149*
150*> \ingroup trsyl
151*
152* =====================================================================
153 SUBROUTINE ztrsyl( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C,
154 $ LDC, SCALE, INFO )
155*
156* -- LAPACK computational routine --
157* -- LAPACK is a software package provided by Univ. of Tennessee, --
158* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
159*
160* .. Scalar Arguments ..
161 CHARACTER TRANA, TRANB
162 INTEGER INFO, ISGN, LDA, LDB, LDC, M, N
163 DOUBLE PRECISION SCALE
164* ..
165* .. Array Arguments ..
166 COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * )
167* ..
168*
169* =====================================================================
170*
171* .. Parameters ..
172 DOUBLE PRECISION ONE
173 parameter( one = 1.0d+0 )
174* ..
175* .. Local Scalars ..
176 LOGICAL NOTRNA, NOTRNB
177 INTEGER J, K, L
178 DOUBLE PRECISION BIGNUM, DA11, DB, EPS, SCALOC, SGN, SMIN,
179 $ smlnum
180 COMPLEX*16 A11, SUML, SUMR, VEC, X11
181* ..
182* .. Local Arrays ..
183 DOUBLE PRECISION DUM( 1 )
184* ..
185* .. External Functions ..
186 LOGICAL LSAME
187 DOUBLE PRECISION DLAMCH, ZLANGE
188 COMPLEX*16 ZDOTC, ZDOTU, ZLADIV
189 EXTERNAL lsame, dlamch, zlange, zdotc, zdotu,
190 $ zladiv
191* ..
192* .. External Subroutines ..
193 EXTERNAL xerbla, zdscal
194* ..
195* .. Intrinsic Functions ..
196 INTRINSIC abs, dble, dcmplx, dconjg, dimag, max, min
197* ..
198* .. Executable Statements ..
199*
200* Decode and Test input parameters
201*
202 notrna = lsame( trana, 'N' )
203 notrnb = lsame( tranb, 'N' )
204*
205 info = 0
206 IF( .NOT.notrna .AND. .NOT.lsame( trana, 'C' ) ) THEN
207 info = -1
208 ELSE IF( .NOT.notrnb .AND. .NOT.lsame( tranb, 'C' ) ) THEN
209 info = -2
210 ELSE IF( isgn.NE.1 .AND. isgn.NE.-1 ) THEN
211 info = -3
212 ELSE IF( m.LT.0 ) THEN
213 info = -4
214 ELSE IF( n.LT.0 ) THEN
215 info = -5
216 ELSE IF( lda.LT.max( 1, m ) ) THEN
217 info = -7
218 ELSE IF( ldb.LT.max( 1, n ) ) THEN
219 info = -9
220 ELSE IF( ldc.LT.max( 1, m ) ) THEN
221 info = -11
222 END IF
223 IF( info.NE.0 ) THEN
224 CALL xerbla( 'ZTRSYL', -info )
225 RETURN
226 END IF
227*
228* Quick return if possible
229*
230 scale = one
231 IF( m.EQ.0 .OR. n.EQ.0 )
232 $ RETURN
233*
234* Set constants to control overflow
235*
236 eps = dlamch( 'P' )
237 smlnum = dlamch( 'S' )
238 bignum = one / smlnum
239 smlnum = smlnum*dble( m*n ) / eps
240 bignum = one / smlnum
241 smin = max( smlnum, eps*zlange( 'M', m, m, a, lda, dum ),
242 $ eps*zlange( 'M', n, n, b, ldb, dum ) )
243 sgn = isgn
244*
245 IF( notrna .AND. notrnb ) THEN
246*
247* Solve A*X + ISGN*X*B = scale*C.
248*
249* The (K,L)th block of X is determined starting from
250* bottom-left corner column by column by
251*
252* A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
253*
254* Where
255* M L-1
256* R(K,L) = SUM [A(K,I)*X(I,L)] +ISGN*SUM [X(K,J)*B(J,L)].
257* I=K+1 J=1
258*
259 DO 30 l = 1, n
260 DO 20 k = m, 1, -1
261*
262 suml = zdotu( m-k, a( k, min( k+1, m ) ), lda,
263 $ c( min( k+1, m ), l ), 1 )
264 sumr = zdotu( l-1, c( k, 1 ), ldc, b( 1, l ), 1 )
265 vec = c( k, l ) - ( suml+sgn*sumr )
266*
267 scaloc = one
268 a11 = a( k, k ) + sgn*b( l, l )
269 da11 = abs( dble( a11 ) ) + abs( dimag( a11 ) )
270 IF( da11.LE.smin ) THEN
271 a11 = smin
272 da11 = smin
273 info = 1
274 END IF
275 db = abs( dble( vec ) ) + abs( dimag( vec ) )
276 IF( da11.LT.one .AND. db.GT.one ) THEN
277 IF( db.GT.bignum*da11 )
278 $ scaloc = one / db
279 END IF
280 x11 = zladiv( vec*dcmplx( scaloc ), a11 )
281*
282 IF( scaloc.NE.one ) THEN
283 DO 10 j = 1, n
284 CALL zdscal( m, scaloc, c( 1, j ), 1 )
285 10 CONTINUE
286 scale = scale*scaloc
287 END IF
288 c( k, l ) = x11
289*
290 20 CONTINUE
291 30 CONTINUE
292*
293 ELSE IF( .NOT.notrna .AND. notrnb ) THEN
294*
295* Solve A**H *X + ISGN*X*B = scale*C.
296*
297* The (K,L)th block of X is determined starting from
298* upper-left corner column by column by
299*
300* A**H(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
301*
302* Where
303* K-1 L-1
304* R(K,L) = SUM [A**H(I,K)*X(I,L)] + ISGN*SUM [X(K,J)*B(J,L)]
305* I=1 J=1
306*
307 DO 60 l = 1, n
308 DO 50 k = 1, m
309*
310 suml = zdotc( k-1, a( 1, k ), 1, c( 1, l ), 1 )
311 sumr = zdotu( l-1, c( k, 1 ), ldc, b( 1, l ), 1 )
312 vec = c( k, l ) - ( suml+sgn*sumr )
313*
314 scaloc = one
315 a11 = dconjg( a( k, k ) ) + sgn*b( l, l )
316 da11 = abs( dble( a11 ) ) + abs( dimag( a11 ) )
317 IF( da11.LE.smin ) THEN
318 a11 = smin
319 da11 = smin
320 info = 1
321 END IF
322 db = abs( dble( vec ) ) + abs( dimag( vec ) )
323 IF( da11.LT.one .AND. db.GT.one ) THEN
324 IF( db.GT.bignum*da11 )
325 $ scaloc = one / db
326 END IF
327*
328 x11 = zladiv( vec*dcmplx( scaloc ), a11 )
329*
330 IF( scaloc.NE.one ) THEN
331 DO 40 j = 1, n
332 CALL zdscal( m, scaloc, c( 1, j ), 1 )
333 40 CONTINUE
334 scale = scale*scaloc
335 END IF
336 c( k, l ) = x11
337*
338 50 CONTINUE
339 60 CONTINUE
340*
341 ELSE IF( .NOT.notrna .AND. .NOT.notrnb ) THEN
342*
343* Solve A**H*X + ISGN*X*B**H = C.
344*
345* The (K,L)th block of X is determined starting from
346* upper-right corner column by column by
347*
348* A**H(K,K)*X(K,L) + ISGN*X(K,L)*B**H(L,L) = C(K,L) - R(K,L)
349*
350* Where
351* K-1
352* R(K,L) = SUM [A**H(I,K)*X(I,L)] +
353* I=1
354* N
355* ISGN*SUM [X(K,J)*B**H(L,J)].
356* J=L+1
357*
358 DO 90 l = n, 1, -1
359 DO 80 k = 1, m
360*
361 suml = zdotc( k-1, a( 1, k ), 1, c( 1, l ), 1 )
362 sumr = zdotc( n-l, c( k, min( l+1, n ) ), ldc,
363 $ b( l, min( l+1, n ) ), ldb )
364 vec = c( k, l ) - ( suml+sgn*dconjg( sumr ) )
365*
366 scaloc = one
367 a11 = dconjg( a( k, k )+sgn*b( l, l ) )
368 da11 = abs( dble( a11 ) ) + abs( dimag( a11 ) )
369 IF( da11.LE.smin ) THEN
370 a11 = smin
371 da11 = smin
372 info = 1
373 END IF
374 db = abs( dble( vec ) ) + abs( dimag( vec ) )
375 IF( da11.LT.one .AND. db.GT.one ) THEN
376 IF( db.GT.bignum*da11 )
377 $ scaloc = one / db
378 END IF
379*
380 x11 = zladiv( vec*dcmplx( scaloc ), a11 )
381*
382 IF( scaloc.NE.one ) THEN
383 DO 70 j = 1, n
384 CALL zdscal( m, scaloc, c( 1, j ), 1 )
385 70 CONTINUE
386 scale = scale*scaloc
387 END IF
388 c( k, l ) = x11
389*
390 80 CONTINUE
391 90 CONTINUE
392*
393 ELSE IF( notrna .AND. .NOT.notrnb ) THEN
394*
395* Solve A*X + ISGN*X*B**H = C.
396*
397* The (K,L)th block of X is determined starting from
398* bottom-left corner column by column by
399*
400* A(K,K)*X(K,L) + ISGN*X(K,L)*B**H(L,L) = C(K,L) - R(K,L)
401*
402* Where
403* M N
404* R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B**H(L,J)]
405* I=K+1 J=L+1
406*
407 DO 120 l = n, 1, -1
408 DO 110 k = m, 1, -1
409*
410 suml = zdotu( m-k, a( k, min( k+1, m ) ), lda,
411 $ c( min( k+1, m ), l ), 1 )
412 sumr = zdotc( n-l, c( k, min( l+1, n ) ), ldc,
413 $ b( l, min( l+1, n ) ), ldb )
414 vec = c( k, l ) - ( suml+sgn*dconjg( sumr ) )
415*
416 scaloc = one
417 a11 = a( k, k ) + sgn*dconjg( b( l, l ) )
418 da11 = abs( dble( a11 ) ) + abs( dimag( a11 ) )
419 IF( da11.LE.smin ) THEN
420 a11 = smin
421 da11 = smin
422 info = 1
423 END IF
424 db = abs( dble( vec ) ) + abs( dimag( vec ) )
425 IF( da11.LT.one .AND. db.GT.one ) THEN
426 IF( db.GT.bignum*da11 )
427 $ scaloc = one / db
428 END IF
429*
430 x11 = zladiv( vec*dcmplx( scaloc ), a11 )
431*
432 IF( scaloc.NE.one ) THEN
433 DO 100 j = 1, n
434 CALL zdscal( m, scaloc, c( 1, j ), 1 )
435 100 CONTINUE
436 scale = scale*scaloc
437 END IF
438 c( k, l ) = x11
439*
440 110 CONTINUE
441 120 CONTINUE
442*
443 END IF
444*
445 RETURN
446*
447* End of ZTRSYL
448*
449 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine zdscal(n, da, zx, incx)
ZDSCAL
Definition zdscal.f:78
subroutine ztrsyl(trana, tranb, isgn, m, n, a, lda, b, ldb, c, ldc, scale, info)
ZTRSYL
Definition ztrsyl.f:155