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