LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
strsyl.f
Go to the documentation of this file.
1*> \brief \b STRSYL
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download STRSYL + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/strsyl.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/strsyl.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/strsyl.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE STRSYL( 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* REAL SCALE
26* ..
27* .. Array Arguments ..
28* REAL A( LDA, * ), B( LDB, * ), C( LDC, * )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> STRSYL solves the real 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**T, and A and B are both upper quasi-
43*> triangular. A is M-by-M and B is N-by-N; the right hand side C and
44*> the solution X are M-by-N; and scale is an output scale factor, set
45*> <= 1 to avoid overflow in X.
46*>
47*> A and B must be in Schur canonical form (as returned by SHSEQR), that
48*> is, block upper triangular with 1-by-1 and 2-by-2 diagonal blocks;
49*> each 2-by-2 diagonal block has its diagonal elements equal and its
50*> off-diagonal elements of opposite sign.
51*> \endverbatim
52*
53* Arguments:
54* ==========
55*
56*> \param[in] TRANA
57*> \verbatim
58*> TRANA is CHARACTER*1
59*> Specifies the option op(A):
60*> = 'N': op(A) = A (No transpose)
61*> = 'T': op(A) = A**T (Transpose)
62*> = 'C': op(A) = A**H (Conjugate transpose = Transpose)
63*> \endverbatim
64*>
65*> \param[in] TRANB
66*> \verbatim
67*> TRANB is CHARACTER*1
68*> Specifies the option op(B):
69*> = 'N': op(B) = B (No transpose)
70*> = 'T': op(B) = B**T (Transpose)
71*> = 'C': op(B) = B**H (Conjugate transpose = Transpose)
72*> \endverbatim
73*>
74*> \param[in] ISGN
75*> \verbatim
76*> ISGN is INTEGER
77*> Specifies the sign in the equation:
78*> = +1: solve op(A)*X + X*op(B) = scale*C
79*> = -1: solve op(A)*X - X*op(B) = scale*C
80*> \endverbatim
81*>
82*> \param[in] M
83*> \verbatim
84*> M is INTEGER
85*> The order of the matrix A, and the number of rows in the
86*> matrices X and C. M >= 0.
87*> \endverbatim
88*>
89*> \param[in] N
90*> \verbatim
91*> N is INTEGER
92*> The order of the matrix B, and the number of columns in the
93*> matrices X and C. N >= 0.
94*> \endverbatim
95*>
96*> \param[in] A
97*> \verbatim
98*> A is REAL array, dimension (LDA,M)
99*> The upper quasi-triangular matrix A, in Schur canonical form.
100*> \endverbatim
101*>
102*> \param[in] LDA
103*> \verbatim
104*> LDA is INTEGER
105*> The leading dimension of the array A. LDA >= max(1,M).
106*> \endverbatim
107*>
108*> \param[in] B
109*> \verbatim
110*> B is REAL array, dimension (LDB,N)
111*> The upper quasi-triangular matrix B, in Schur canonical form.
112*> \endverbatim
113*>
114*> \param[in] LDB
115*> \verbatim
116*> LDB is INTEGER
117*> The leading dimension of the array B. LDB >= max(1,N).
118*> \endverbatim
119*>
120*> \param[in,out] C
121*> \verbatim
122*> C is REAL array, dimension (LDC,N)
123*> On entry, the M-by-N right hand side matrix C.
124*> On exit, C is overwritten by the solution matrix X.
125*> \endverbatim
126*>
127*> \param[in] LDC
128*> \verbatim
129*> LDC is INTEGER
130*> The leading dimension of the array C. LDC >= max(1,M)
131*> \endverbatim
132*>
133*> \param[out] SCALE
134*> \verbatim
135*> SCALE is REAL
136*> The scale factor, scale, set <= 1 to avoid overflow in X.
137*> \endverbatim
138*>
139*> \param[out] INFO
140*> \verbatim
141*> INFO is INTEGER
142*> = 0: successful exit
143*> < 0: if INFO = -i, the i-th argument had an illegal value
144*> = 1: A and B have common or very close eigenvalues; perturbed
145*> values were used to solve the equation (but the matrices
146*> A and B are unchanged).
147*> \endverbatim
148*
149* Authors:
150* ========
151*
152*> \author Univ. of Tennessee
153*> \author Univ. of California Berkeley
154*> \author Univ. of Colorado Denver
155*> \author NAG Ltd.
156*
157*> \ingroup trsyl
158*
159* =====================================================================
160 SUBROUTINE strsyl( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C,
161 $ LDC, SCALE, INFO )
162*
163* -- LAPACK computational routine --
164* -- LAPACK is a software package provided by Univ. of Tennessee, --
165* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
166*
167* .. Scalar Arguments ..
168 CHARACTER TRANA, TRANB
169 INTEGER INFO, ISGN, LDA, LDB, LDC, M, N
170 REAL SCALE
171* ..
172* .. Array Arguments ..
173 REAL A( LDA, * ), B( LDB, * ), C( LDC, * )
174* ..
175*
176* =====================================================================
177*
178* .. Parameters ..
179 REAL ZERO, ONE
180 parameter( zero = 0.0e+0, one = 1.0e+0 )
181* ..
182* .. Local Scalars ..
183 LOGICAL NOTRNA, NOTRNB
184 INTEGER IERR, J, K, K1, K2, KNEXT, L, L1, L2, LNEXT
185 REAL A11, BIGNUM, DA11, DB, EPS, SCALOC, SGN, SMIN,
186 $ smlnum, suml, sumr, xnorm
187* ..
188* .. Local Arrays ..
189 REAL DUM( 1 ), VEC( 2, 2 ), X( 2, 2 )
190* ..
191* .. External Functions ..
192 LOGICAL LSAME
193 REAL SDOT, SLAMCH, SLANGE
194 EXTERNAL lsame, sdot, slamch, slange
195* ..
196* .. External Subroutines ..
197 EXTERNAL slaln2, slasy2, sscal, xerbla
198* ..
199* .. Intrinsic Functions ..
200 INTRINSIC abs, max, min, real
201* ..
202* .. Executable Statements ..
203*
204* Decode and Test input parameters
205*
206 notrna = lsame( trana, 'N' )
207 notrnb = lsame( tranb, 'N' )
208*
209 info = 0
210 IF( .NOT.notrna .AND. .NOT.lsame( trana, 'T' ) .AND. .NOT.
211 $ lsame( trana, 'C' ) ) THEN
212 info = -1
213 ELSE IF( .NOT.notrnb .AND. .NOT.lsame( tranb, 'T' ) .AND. .NOT.
214 $ lsame( tranb, 'C' ) ) THEN
215 info = -2
216 ELSE IF( isgn.NE.1 .AND. isgn.NE.-1 ) THEN
217 info = -3
218 ELSE IF( m.LT.0 ) THEN
219 info = -4
220 ELSE IF( n.LT.0 ) THEN
221 info = -5
222 ELSE IF( lda.LT.max( 1, m ) ) THEN
223 info = -7
224 ELSE IF( ldb.LT.max( 1, n ) ) THEN
225 info = -9
226 ELSE IF( ldc.LT.max( 1, m ) ) THEN
227 info = -11
228 END IF
229 IF( info.NE.0 ) THEN
230 CALL xerbla( 'STRSYL', -info )
231 RETURN
232 END IF
233*
234* Quick return if possible
235*
236 scale = one
237 IF( m.EQ.0 .OR. n.EQ.0 )
238 $ RETURN
239*
240* Set constants to control overflow
241*
242 eps = slamch( 'P' )
243 smlnum = slamch( 'S' )
244 bignum = one / smlnum
245 smlnum = smlnum*real( m*n ) / eps
246 bignum = one / smlnum
247*
248 smin = max( smlnum, eps*slange( 'M', m, m, a, lda, dum ),
249 $ eps*slange( 'M', n, n, b, ldb, dum ) )
250*
251 sgn = real( isgn )
252*
253 IF( notrna .AND. notrnb ) THEN
254*
255* Solve A*X + ISGN*X*B = scale*C.
256*
257* The (K,L)th block of X is determined starting from
258* bottom-left corner column by column by
259*
260* A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
261*
262* Where
263* M L-1
264* R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(J,L)].
265* I=K+1 J=1
266*
267* Start column loop (index = L)
268* L1 (L2) : column index of the first (first) row of X(K,L).
269*
270 lnext = 1
271 DO 70 l = 1, n
272 IF( l.LT.lnext )
273 $ GO TO 70
274 IF( l.EQ.n ) THEN
275 l1 = l
276 l2 = l
277 ELSE
278 IF( b( l+1, l ).NE.zero ) THEN
279 l1 = l
280 l2 = l + 1
281 lnext = l + 2
282 ELSE
283 l1 = l
284 l2 = l
285 lnext = l + 1
286 END IF
287 END IF
288*
289* Start row loop (index = K)
290* K1 (K2): row index of the first (last) row of X(K,L).
291*
292 knext = m
293 DO 60 k = m, 1, -1
294 IF( k.GT.knext )
295 $ GO TO 60
296 IF( k.EQ.1 ) THEN
297 k1 = k
298 k2 = k
299 ELSE
300 IF( a( k, k-1 ).NE.zero ) THEN
301 k1 = k - 1
302 k2 = k
303 knext = k - 2
304 ELSE
305 k1 = k
306 k2 = k
307 knext = k - 1
308 END IF
309 END IF
310*
311 IF( l1.EQ.l2 .AND. k1.EQ.k2 ) THEN
312 suml = sdot( m-k1, a( k1, min( k1+1, m ) ), lda,
313 $ c( min( k1+1, m ), l1 ), 1 )
314 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
315 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
316 scaloc = one
317*
318 a11 = a( k1, k1 ) + sgn*b( l1, l1 )
319 da11 = abs( a11 )
320 IF( da11.LE.smin ) THEN
321 a11 = smin
322 da11 = smin
323 info = 1
324 END IF
325 db = abs( vec( 1, 1 ) )
326 IF( da11.LT.one .AND. db.GT.one ) THEN
327 IF( db.GT.bignum*da11 )
328 $ scaloc = one / db
329 END IF
330 x( 1, 1 ) = ( vec( 1, 1 )*scaloc ) / a11
331*
332 IF( scaloc.NE.one ) THEN
333 DO 10 j = 1, n
334 CALL sscal( m, scaloc, c( 1, j ), 1 )
335 10 CONTINUE
336 scale = scale*scaloc
337 END IF
338 c( k1, l1 ) = x( 1, 1 )
339*
340 ELSE IF( l1.EQ.l2 .AND. k1.NE.k2 ) THEN
341*
342 suml = sdot( m-k2, a( k1, min( k2+1, m ) ), lda,
343 $ c( min( k2+1, m ), l1 ), 1 )
344 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
345 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
346*
347 suml = sdot( m-k2, a( k2, min( k2+1, m ) ), lda,
348 $ c( min( k2+1, m ), l1 ), 1 )
349 sumr = sdot( l1-1, c( k2, 1 ), ldc, b( 1, l1 ), 1 )
350 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
351*
352 CALL slaln2( .false., 2, 1, smin, one, a( k1, k1 ),
353 $ lda, one, one, vec, 2, -sgn*b( l1, l1 ),
354 $ zero, x, 2, scaloc, xnorm, ierr )
355 IF( ierr.NE.0 )
356 $ info = 1
357*
358 IF( scaloc.NE.one ) THEN
359 DO 20 j = 1, n
360 CALL sscal( m, scaloc, c( 1, j ), 1 )
361 20 CONTINUE
362 scale = scale*scaloc
363 END IF
364 c( k1, l1 ) = x( 1, 1 )
365 c( k2, l1 ) = x( 2, 1 )
366*
367 ELSE IF( l1.NE.l2 .AND. k1.EQ.k2 ) THEN
368*
369 suml = sdot( m-k1, a( k1, min( k1+1, m ) ), lda,
370 $ c( min( k1+1, m ), l1 ), 1 )
371 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
372 vec( 1, 1 ) = sgn*( c( k1, l1 )-( suml+sgn*sumr ) )
373*
374 suml = sdot( m-k1, a( k1, min( k1+1, m ) ), lda,
375 $ c( min( k1+1, m ), l2 ), 1 )
376 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l2 ), 1 )
377 vec( 2, 1 ) = sgn*( c( k1, l2 )-( suml+sgn*sumr ) )
378*
379 CALL slaln2( .true., 2, 1, smin, one, b( l1, l1 ),
380 $ ldb, one, one, vec, 2, -sgn*a( k1, k1 ),
381 $ zero, x, 2, scaloc, xnorm, ierr )
382 IF( ierr.NE.0 )
383 $ info = 1
384*
385 IF( scaloc.NE.one ) THEN
386 DO 40 j = 1, n
387 CALL sscal( m, scaloc, c( 1, j ), 1 )
388 40 CONTINUE
389 scale = scale*scaloc
390 END IF
391 c( k1, l1 ) = x( 1, 1 )
392 c( k1, l2 ) = x( 2, 1 )
393*
394 ELSE IF( l1.NE.l2 .AND. k1.NE.k2 ) THEN
395*
396 suml = sdot( m-k2, a( k1, min( k2+1, m ) ), lda,
397 $ c( min( k2+1, m ), l1 ), 1 )
398 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
399 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
400*
401 suml = sdot( m-k2, a( k1, min( k2+1, m ) ), lda,
402 $ c( min( k2+1, m ), l2 ), 1 )
403 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l2 ), 1 )
404 vec( 1, 2 ) = c( k1, l2 ) - ( suml+sgn*sumr )
405*
406 suml = sdot( m-k2, a( k2, min( k2+1, m ) ), lda,
407 $ c( min( k2+1, m ), l1 ), 1 )
408 sumr = sdot( l1-1, c( k2, 1 ), ldc, b( 1, l1 ), 1 )
409 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
410*
411 suml = sdot( m-k2, a( k2, min( k2+1, m ) ), lda,
412 $ c( min( k2+1, m ), l2 ), 1 )
413 sumr = sdot( l1-1, c( k2, 1 ), ldc, b( 1, l2 ), 1 )
414 vec( 2, 2 ) = c( k2, l2 ) - ( suml+sgn*sumr )
415*
416 CALL slasy2( .false., .false., isgn, 2, 2,
417 $ a( k1, k1 ), lda, b( l1, l1 ), ldb, vec,
418 $ 2, scaloc, x, 2, xnorm, ierr )
419 IF( ierr.NE.0 )
420 $ info = 1
421*
422 IF( scaloc.NE.one ) THEN
423 DO 50 j = 1, n
424 CALL sscal( m, scaloc, c( 1, j ), 1 )
425 50 CONTINUE
426 scale = scale*scaloc
427 END IF
428 c( k1, l1 ) = x( 1, 1 )
429 c( k1, l2 ) = x( 1, 2 )
430 c( k2, l1 ) = x( 2, 1 )
431 c( k2, l2 ) = x( 2, 2 )
432 END IF
433*
434 60 CONTINUE
435*
436 70 CONTINUE
437*
438 ELSE IF( .NOT.notrna .AND. notrnb ) THEN
439*
440* Solve A**T *X + ISGN*X*B = scale*C.
441*
442* The (K,L)th block of X is determined starting from
443* upper-left corner column by column by
444*
445* A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
446*
447* Where
448* K-1 L-1
449* R(K,L) = SUM [A(I,K)**T*X(I,L)] +ISGN*SUM [X(K,J)*B(J,L)]
450* I=1 J=1
451*
452* Start column loop (index = L)
453* L1 (L2): column index of the first (last) row of X(K,L)
454*
455 lnext = 1
456 DO 130 l = 1, n
457 IF( l.LT.lnext )
458 $ GO TO 130
459 IF( l.EQ.n ) THEN
460 l1 = l
461 l2 = l
462 ELSE
463 IF( b( l+1, l ).NE.zero ) THEN
464 l1 = l
465 l2 = l + 1
466 lnext = l + 2
467 ELSE
468 l1 = l
469 l2 = l
470 lnext = l + 1
471 END IF
472 END IF
473*
474* Start row loop (index = K)
475* K1 (K2): row index of the first (last) row of X(K,L)
476*
477 knext = 1
478 DO 120 k = 1, m
479 IF( k.LT.knext )
480 $ GO TO 120
481 IF( k.EQ.m ) THEN
482 k1 = k
483 k2 = k
484 ELSE
485 IF( a( k+1, k ).NE.zero ) THEN
486 k1 = k
487 k2 = k + 1
488 knext = k + 2
489 ELSE
490 k1 = k
491 k2 = k
492 knext = k + 1
493 END IF
494 END IF
495*
496 IF( l1.EQ.l2 .AND. k1.EQ.k2 ) THEN
497 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
498 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
499 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
500 scaloc = one
501*
502 a11 = a( k1, k1 ) + sgn*b( l1, l1 )
503 da11 = abs( a11 )
504 IF( da11.LE.smin ) THEN
505 a11 = smin
506 da11 = smin
507 info = 1
508 END IF
509 db = abs( vec( 1, 1 ) )
510 IF( da11.LT.one .AND. db.GT.one ) THEN
511 IF( db.GT.bignum*da11 )
512 $ scaloc = one / db
513 END IF
514 x( 1, 1 ) = ( vec( 1, 1 )*scaloc ) / a11
515*
516 IF( scaloc.NE.one ) THEN
517 DO 80 j = 1, n
518 CALL sscal( m, scaloc, c( 1, j ), 1 )
519 80 CONTINUE
520 scale = scale*scaloc
521 END IF
522 c( k1, l1 ) = x( 1, 1 )
523*
524 ELSE IF( l1.EQ.l2 .AND. k1.NE.k2 ) THEN
525*
526 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
527 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
528 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
529*
530 suml = sdot( k1-1, a( 1, k2 ), 1, c( 1, l1 ), 1 )
531 sumr = sdot( l1-1, c( k2, 1 ), ldc, b( 1, l1 ), 1 )
532 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
533*
534 CALL slaln2( .true., 2, 1, smin, one, a( k1, k1 ),
535 $ lda, one, one, vec, 2, -sgn*b( l1, l1 ),
536 $ zero, x, 2, scaloc, xnorm, ierr )
537 IF( ierr.NE.0 )
538 $ info = 1
539*
540 IF( scaloc.NE.one ) THEN
541 DO 90 j = 1, n
542 CALL sscal( m, scaloc, c( 1, j ), 1 )
543 90 CONTINUE
544 scale = scale*scaloc
545 END IF
546 c( k1, l1 ) = x( 1, 1 )
547 c( k2, l1 ) = x( 2, 1 )
548*
549 ELSE IF( l1.NE.l2 .AND. k1.EQ.k2 ) THEN
550*
551 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
552 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
553 vec( 1, 1 ) = sgn*( c( k1, l1 )-( suml+sgn*sumr ) )
554*
555 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l2 ), 1 )
556 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l2 ), 1 )
557 vec( 2, 1 ) = sgn*( c( k1, l2 )-( suml+sgn*sumr ) )
558*
559 CALL slaln2( .true., 2, 1, smin, one, b( l1, l1 ),
560 $ ldb, one, one, vec, 2, -sgn*a( k1, k1 ),
561 $ zero, x, 2, scaloc, xnorm, ierr )
562 IF( ierr.NE.0 )
563 $ info = 1
564*
565 IF( scaloc.NE.one ) THEN
566 DO 100 j = 1, n
567 CALL sscal( m, scaloc, c( 1, j ), 1 )
568 100 CONTINUE
569 scale = scale*scaloc
570 END IF
571 c( k1, l1 ) = x( 1, 1 )
572 c( k1, l2 ) = x( 2, 1 )
573*
574 ELSE IF( l1.NE.l2 .AND. k1.NE.k2 ) THEN
575*
576 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
577 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
578 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
579*
580 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l2 ), 1 )
581 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l2 ), 1 )
582 vec( 1, 2 ) = c( k1, l2 ) - ( suml+sgn*sumr )
583*
584 suml = sdot( k1-1, a( 1, k2 ), 1, c( 1, l1 ), 1 )
585 sumr = sdot( l1-1, c( k2, 1 ), ldc, b( 1, l1 ), 1 )
586 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
587*
588 suml = sdot( k1-1, a( 1, k2 ), 1, c( 1, l2 ), 1 )
589 sumr = sdot( l1-1, c( k2, 1 ), ldc, b( 1, l2 ), 1 )
590 vec( 2, 2 ) = c( k2, l2 ) - ( suml+sgn*sumr )
591*
592 CALL slasy2( .true., .false., isgn, 2, 2, a( k1,
593 $ k1 ),
594 $ lda, b( l1, l1 ), ldb, vec, 2, scaloc, x,
595 $ 2, xnorm, ierr )
596 IF( ierr.NE.0 )
597 $ info = 1
598*
599 IF( scaloc.NE.one ) THEN
600 DO 110 j = 1, n
601 CALL sscal( m, scaloc, c( 1, j ), 1 )
602 110 CONTINUE
603 scale = scale*scaloc
604 END IF
605 c( k1, l1 ) = x( 1, 1 )
606 c( k1, l2 ) = x( 1, 2 )
607 c( k2, l1 ) = x( 2, 1 )
608 c( k2, l2 ) = x( 2, 2 )
609 END IF
610*
611 120 CONTINUE
612 130 CONTINUE
613*
614 ELSE IF( .NOT.notrna .AND. .NOT.notrnb ) THEN
615*
616* Solve A**T*X + ISGN*X*B**T = scale*C.
617*
618* The (K,L)th block of X is determined starting from
619* top-right corner column by column by
620*
621* A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L)**T = C(K,L) - R(K,L)
622*
623* Where
624* K-1 N
625* R(K,L) = SUM [A(I,K)**T*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**T].
626* I=1 J=L+1
627*
628* Start column loop (index = L)
629* L1 (L2): column index of the first (last) row of X(K,L)
630*
631 lnext = n
632 DO 190 l = n, 1, -1
633 IF( l.GT.lnext )
634 $ GO TO 190
635 IF( l.EQ.1 ) THEN
636 l1 = l
637 l2 = l
638 ELSE
639 IF( b( l, l-1 ).NE.zero ) THEN
640 l1 = l - 1
641 l2 = l
642 lnext = l - 2
643 ELSE
644 l1 = l
645 l2 = l
646 lnext = l - 1
647 END IF
648 END IF
649*
650* Start row loop (index = K)
651* K1 (K2): row index of the first (last) row of X(K,L)
652*
653 knext = 1
654 DO 180 k = 1, m
655 IF( k.LT.knext )
656 $ GO TO 180
657 IF( k.EQ.m ) THEN
658 k1 = k
659 k2 = k
660 ELSE
661 IF( a( k+1, k ).NE.zero ) THEN
662 k1 = k
663 k2 = k + 1
664 knext = k + 2
665 ELSE
666 k1 = k
667 k2 = k
668 knext = k + 1
669 END IF
670 END IF
671*
672 IF( l1.EQ.l2 .AND. k1.EQ.k2 ) THEN
673 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
674 sumr = sdot( n-l1, c( k1, min( l1+1, n ) ), ldc,
675 $ b( l1, min( l1+1, n ) ), ldb )
676 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
677 scaloc = one
678*
679 a11 = a( k1, k1 ) + sgn*b( l1, l1 )
680 da11 = abs( a11 )
681 IF( da11.LE.smin ) THEN
682 a11 = smin
683 da11 = smin
684 info = 1
685 END IF
686 db = abs( vec( 1, 1 ) )
687 IF( da11.LT.one .AND. db.GT.one ) THEN
688 IF( db.GT.bignum*da11 )
689 $ scaloc = one / db
690 END IF
691 x( 1, 1 ) = ( vec( 1, 1 )*scaloc ) / a11
692*
693 IF( scaloc.NE.one ) THEN
694 DO 140 j = 1, n
695 CALL sscal( m, scaloc, c( 1, j ), 1 )
696 140 CONTINUE
697 scale = scale*scaloc
698 END IF
699 c( k1, l1 ) = x( 1, 1 )
700*
701 ELSE IF( l1.EQ.l2 .AND. k1.NE.k2 ) THEN
702*
703 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
704 sumr = sdot( n-l2, c( k1, min( l2+1, n ) ), ldc,
705 $ b( l1, min( l2+1, n ) ), ldb )
706 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
707*
708 suml = sdot( k1-1, a( 1, k2 ), 1, c( 1, l1 ), 1 )
709 sumr = sdot( n-l2, c( k2, min( l2+1, n ) ), ldc,
710 $ b( l1, min( l2+1, n ) ), ldb )
711 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
712*
713 CALL slaln2( .true., 2, 1, smin, one, a( k1, k1 ),
714 $ lda, one, one, vec, 2, -sgn*b( l1, l1 ),
715 $ zero, x, 2, scaloc, xnorm, ierr )
716 IF( ierr.NE.0 )
717 $ info = 1
718*
719 IF( scaloc.NE.one ) THEN
720 DO 150 j = 1, n
721 CALL sscal( m, scaloc, c( 1, j ), 1 )
722 150 CONTINUE
723 scale = scale*scaloc
724 END IF
725 c( k1, l1 ) = x( 1, 1 )
726 c( k2, l1 ) = x( 2, 1 )
727*
728 ELSE IF( l1.NE.l2 .AND. k1.EQ.k2 ) THEN
729*
730 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
731 sumr = sdot( n-l2, c( k1, min( l2+1, n ) ), ldc,
732 $ b( l1, min( l2+1, n ) ), ldb )
733 vec( 1, 1 ) = sgn*( c( k1, l1 )-( suml+sgn*sumr ) )
734*
735 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l2 ), 1 )
736 sumr = sdot( n-l2, c( k1, min( l2+1, n ) ), ldc,
737 $ b( l2, min( l2+1, n ) ), ldb )
738 vec( 2, 1 ) = sgn*( c( k1, l2 )-( suml+sgn*sumr ) )
739*
740 CALL slaln2( .false., 2, 1, smin, one, b( l1, l1 ),
741 $ ldb, one, one, vec, 2, -sgn*a( k1, k1 ),
742 $ zero, x, 2, scaloc, xnorm, ierr )
743 IF( ierr.NE.0 )
744 $ info = 1
745*
746 IF( scaloc.NE.one ) THEN
747 DO 160 j = 1, n
748 CALL sscal( m, scaloc, c( 1, j ), 1 )
749 160 CONTINUE
750 scale = scale*scaloc
751 END IF
752 c( k1, l1 ) = x( 1, 1 )
753 c( k1, l2 ) = x( 2, 1 )
754*
755 ELSE IF( l1.NE.l2 .AND. k1.NE.k2 ) THEN
756*
757 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
758 sumr = sdot( n-l2, c( k1, min( l2+1, n ) ), ldc,
759 $ b( l1, min( l2+1, n ) ), ldb )
760 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
761*
762 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l2 ), 1 )
763 sumr = sdot( n-l2, c( k1, min( l2+1, n ) ), ldc,
764 $ b( l2, min( l2+1, n ) ), ldb )
765 vec( 1, 2 ) = c( k1, l2 ) - ( suml+sgn*sumr )
766*
767 suml = sdot( k1-1, a( 1, k2 ), 1, c( 1, l1 ), 1 )
768 sumr = sdot( n-l2, c( k2, min( l2+1, n ) ), ldc,
769 $ b( l1, min( l2+1, n ) ), ldb )
770 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
771*
772 suml = sdot( k1-1, a( 1, k2 ), 1, c( 1, l2 ), 1 )
773 sumr = sdot( n-l2, c( k2, min( l2+1, n ) ), ldc,
774 $ b( l2, min(l2+1, n ) ), ldb )
775 vec( 2, 2 ) = c( k2, l2 ) - ( suml+sgn*sumr )
776*
777 CALL slasy2( .true., .true., isgn, 2, 2, a( k1,
778 $ k1 ),
779 $ lda, b( l1, l1 ), ldb, vec, 2, scaloc, x,
780 $ 2, xnorm, ierr )
781 IF( ierr.NE.0 )
782 $ info = 1
783*
784 IF( scaloc.NE.one ) THEN
785 DO 170 j = 1, n
786 CALL sscal( m, scaloc, c( 1, j ), 1 )
787 170 CONTINUE
788 scale = scale*scaloc
789 END IF
790 c( k1, l1 ) = x( 1, 1 )
791 c( k1, l2 ) = x( 1, 2 )
792 c( k2, l1 ) = x( 2, 1 )
793 c( k2, l2 ) = x( 2, 2 )
794 END IF
795*
796 180 CONTINUE
797 190 CONTINUE
798*
799 ELSE IF( notrna .AND. .NOT.notrnb ) THEN
800*
801* Solve A*X + ISGN*X*B**T = scale*C.
802*
803* The (K,L)th block of X is determined starting from
804* bottom-right corner column by column by
805*
806* A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L)**T = C(K,L) - R(K,L)
807*
808* Where
809* M N
810* R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**T].
811* I=K+1 J=L+1
812*
813* Start column loop (index = L)
814* L1 (L2): column index of the first (last) row of X(K,L)
815*
816 lnext = n
817 DO 250 l = n, 1, -1
818 IF( l.GT.lnext )
819 $ GO TO 250
820 IF( l.EQ.1 ) THEN
821 l1 = l
822 l2 = l
823 ELSE
824 IF( b( l, l-1 ).NE.zero ) THEN
825 l1 = l - 1
826 l2 = l
827 lnext = l - 2
828 ELSE
829 l1 = l
830 l2 = l
831 lnext = l - 1
832 END IF
833 END IF
834*
835* Start row loop (index = K)
836* K1 (K2): row index of the first (last) row of X(K,L)
837*
838 knext = m
839 DO 240 k = m, 1, -1
840 IF( k.GT.knext )
841 $ GO TO 240
842 IF( k.EQ.1 ) THEN
843 k1 = k
844 k2 = k
845 ELSE
846 IF( a( k, k-1 ).NE.zero ) THEN
847 k1 = k - 1
848 k2 = k
849 knext = k - 2
850 ELSE
851 k1 = k
852 k2 = k
853 knext = k - 1
854 END IF
855 END IF
856*
857 IF( l1.EQ.l2 .AND. k1.EQ.k2 ) THEN
858 suml = sdot( m-k1, a( k1, min(k1+1, m ) ), lda,
859 $ c( min( k1+1, m ), l1 ), 1 )
860 sumr = sdot( n-l1, c( k1, min( l1+1, n ) ), ldc,
861 $ b( l1, min( l1+1, n ) ), ldb )
862 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
863 scaloc = one
864*
865 a11 = a( k1, k1 ) + sgn*b( l1, l1 )
866 da11 = abs( a11 )
867 IF( da11.LE.smin ) THEN
868 a11 = smin
869 da11 = smin
870 info = 1
871 END IF
872 db = abs( vec( 1, 1 ) )
873 IF( da11.LT.one .AND. db.GT.one ) THEN
874 IF( db.GT.bignum*da11 )
875 $ scaloc = one / db
876 END IF
877 x( 1, 1 ) = ( vec( 1, 1 )*scaloc ) / a11
878*
879 IF( scaloc.NE.one ) THEN
880 DO 200 j = 1, n
881 CALL sscal( m, scaloc, c( 1, j ), 1 )
882 200 CONTINUE
883 scale = scale*scaloc
884 END IF
885 c( k1, l1 ) = x( 1, 1 )
886*
887 ELSE IF( l1.EQ.l2 .AND. k1.NE.k2 ) THEN
888*
889 suml = sdot( m-k2, a( k1, min( k2+1, m ) ), lda,
890 $ c( min( k2+1, m ), l1 ), 1 )
891 sumr = sdot( n-l2, c( k1, min( l2+1, n ) ), ldc,
892 $ b( l1, min( l2+1, n ) ), ldb )
893 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
894*
895 suml = sdot( m-k2, a( k2, min( k2+1, m ) ), lda,
896 $ c( min( k2+1, m ), l1 ), 1 )
897 sumr = sdot( n-l2, c( k2, min( l2+1, n ) ), ldc,
898 $ b( l1, min( l2+1, n ) ), ldb )
899 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
900*
901 CALL slaln2( .false., 2, 1, smin, one, a( k1, k1 ),
902 $ lda, one, one, vec, 2, -sgn*b( l1, l1 ),
903 $ zero, x, 2, scaloc, xnorm, ierr )
904 IF( ierr.NE.0 )
905 $ info = 1
906*
907 IF( scaloc.NE.one ) THEN
908 DO 210 j = 1, n
909 CALL sscal( m, scaloc, c( 1, j ), 1 )
910 210 CONTINUE
911 scale = scale*scaloc
912 END IF
913 c( k1, l1 ) = x( 1, 1 )
914 c( k2, l1 ) = x( 2, 1 )
915*
916 ELSE IF( l1.NE.l2 .AND. k1.EQ.k2 ) THEN
917*
918 suml = sdot( m-k1, a( k1, min( k1+1, m ) ), lda,
919 $ c( min( k1+1, m ), l1 ), 1 )
920 sumr = sdot( n-l2, c( k1, min( l2+1, n ) ), ldc,
921 $ b( l1, min( l2+1, n ) ), ldb )
922 vec( 1, 1 ) = sgn*( c( k1, l1 )-( suml+sgn*sumr ) )
923*
924 suml = sdot( m-k1, a( k1, min( k1+1, m ) ), lda,
925 $ c( min( k1+1, m ), l2 ), 1 )
926 sumr = sdot( n-l2, c( k1, min( l2+1, n ) ), ldc,
927 $ b( l2, min( l2+1, n ) ), ldb )
928 vec( 2, 1 ) = sgn*( c( k1, l2 )-( suml+sgn*sumr ) )
929*
930 CALL slaln2( .false., 2, 1, smin, one, b( l1, l1 ),
931 $ ldb, one, one, vec, 2, -sgn*a( k1, k1 ),
932 $ zero, x, 2, scaloc, xnorm, ierr )
933 IF( ierr.NE.0 )
934 $ info = 1
935*
936 IF( scaloc.NE.one ) THEN
937 DO 220 j = 1, n
938 CALL sscal( m, scaloc, c( 1, j ), 1 )
939 220 CONTINUE
940 scale = scale*scaloc
941 END IF
942 c( k1, l1 ) = x( 1, 1 )
943 c( k1, l2 ) = x( 2, 1 )
944*
945 ELSE IF( l1.NE.l2 .AND. k1.NE.k2 ) THEN
946*
947 suml = sdot( m-k2, a( k1, min( k2+1, m ) ), lda,
948 $ c( min( k2+1, m ), l1 ), 1 )
949 sumr = sdot( n-l2, c( k1, min( l2+1, n ) ), ldc,
950 $ b( l1, min( l2+1, n ) ), ldb )
951 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
952*
953 suml = sdot( m-k2, a( k1, min( k2+1, m ) ), lda,
954 $ c( min( k2+1, m ), l2 ), 1 )
955 sumr = sdot( n-l2, c( k1, min( l2+1, n ) ), ldc,
956 $ b( l2, min( l2+1, n ) ), ldb )
957 vec( 1, 2 ) = c( k1, l2 ) - ( suml+sgn*sumr )
958*
959 suml = sdot( m-k2, a( k2, min( k2+1, m ) ), lda,
960 $ c( min( k2+1, m ), l1 ), 1 )
961 sumr = sdot( n-l2, c( k2, min( l2+1, n ) ), ldc,
962 $ b( l1, min( l2+1, n ) ), ldb )
963 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
964*
965 suml = sdot( m-k2, a( k2, min( k2+1, m ) ), lda,
966 $ c( min( k2+1, m ), l2 ), 1 )
967 sumr = sdot( n-l2, c( k2, min( l2+1, n ) ), ldc,
968 $ b( l2, min( l2+1, n ) ), ldb )
969 vec( 2, 2 ) = c( k2, l2 ) - ( suml+sgn*sumr )
970*
971 CALL slasy2( .false., .true., isgn, 2, 2, a( k1,
972 $ k1 ),
973 $ lda, b( l1, l1 ), ldb, vec, 2, scaloc, x,
974 $ 2, xnorm, ierr )
975 IF( ierr.NE.0 )
976 $ info = 1
977*
978 IF( scaloc.NE.one ) THEN
979 DO 230 j = 1, n
980 CALL sscal( m, scaloc, c( 1, j ), 1 )
981 230 CONTINUE
982 scale = scale*scaloc
983 END IF
984 c( k1, l1 ) = x( 1, 1 )
985 c( k1, l2 ) = x( 1, 2 )
986 c( k2, l1 ) = x( 2, 1 )
987 c( k2, l2 ) = x( 2, 2 )
988 END IF
989*
990 240 CONTINUE
991 250 CONTINUE
992*
993 END IF
994*
995 RETURN
996*
997* End of STRSYL
998*
999 END
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
slaln2
subroutine slaln2(ltrans, na, nw, smin, ca, a, lda, d1, d2, b, ldb, wr, wi, x, ldx, scale, xnorm, info)
SLALN2 solves a 1-by-1 or 2-by-2 linear system of equations of the specified form.
Definition slaln2.f:216
slasy2
subroutine slasy2(ltranl, ltranr, isgn, n1, n2, tl, ldtl, tr, ldtr, b, ldb, scale, x, ldx, xnorm, info)
SLASY2 solves the Sylvester matrix equation where the matrices are of order 1 or 2.
Definition slasy2.f:172
sscal
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
strsyl
subroutine strsyl(trana, tranb, isgn, m, n, a, lda, b, ldb, c, ldc, scale, info)
STRSYL
Definition strsyl.f:162