LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ strsyl()

subroutine strsyl ( character  TRANA,
character  TRANB,
integer  ISGN,
integer  M,
integer  N,
real, dimension( lda, * )  A,
integer  LDA,
real, dimension( ldb, * )  B,
integer  LDB,
real, dimension( ldc, * )  C,
integer  LDC,
real  SCALE,
integer  INFO 
)

STRSYL

Download STRSYL + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 STRSYL solves the real Sylvester matrix equation:

    op(A)*X + X*op(B) = scale*C or
    op(A)*X - X*op(B) = scale*C,

 where op(A) = A or A**T, and  A and B are both upper quasi-
 triangular. A is M-by-M and B is N-by-N; the right hand side C and
 the solution X are M-by-N; and scale is an output scale factor, set
 <= 1 to avoid overflow in X.

 A and B must be in Schur canonical form (as returned by SHSEQR), that
 is, block upper triangular with 1-by-1 and 2-by-2 diagonal blocks;
 each 2-by-2 diagonal block has its diagonal elements equal and its
 off-diagonal elements of opposite sign.
Parameters
[in]TRANA
          TRANA is CHARACTER*1
          Specifies the option op(A):
          = 'N': op(A) = A    (No transpose)
          = 'T': op(A) = A**T (Transpose)
          = 'C': op(A) = A**H (Conjugate transpose = Transpose)
[in]TRANB
          TRANB is CHARACTER*1
          Specifies the option op(B):
          = 'N': op(B) = B    (No transpose)
          = 'T': op(B) = B**T (Transpose)
          = 'C': op(B) = B**H (Conjugate transpose = Transpose)
[in]ISGN
          ISGN is INTEGER
          Specifies the sign in the equation:
          = +1: solve op(A)*X + X*op(B) = scale*C
          = -1: solve op(A)*X - X*op(B) = scale*C
[in]M
          M is INTEGER
          The order of the matrix A, and the number of rows in the
          matrices X and C. M >= 0.
[in]N
          N is INTEGER
          The order of the matrix B, and the number of columns in the
          matrices X and C. N >= 0.
[in]A
          A is REAL array, dimension (LDA,M)
          The upper quasi-triangular matrix A, in Schur canonical form.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A. LDA >= max(1,M).
[in]B
          B is REAL array, dimension (LDB,N)
          The upper quasi-triangular matrix B, in Schur canonical form.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B. LDB >= max(1,N).
[in,out]C
          C is REAL array, dimension (LDC,N)
          On entry, the M-by-N right hand side matrix C.
          On exit, C is overwritten by the solution matrix X.
[in]LDC
          LDC is INTEGER
          The leading dimension of the array C. LDC >= max(1,M)
[out]SCALE
          SCALE is REAL
          The scale factor, scale, set <= 1 to avoid overflow in X.
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value
          = 1: A and B have common or very close eigenvalues; perturbed
               values were used to solve the equation (but the matrices
               A and B are unchanged).
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 162 of file strsyl.f.

164*
165* -- LAPACK computational routine --
166* -- LAPACK is a software package provided by Univ. of Tennessee, --
167* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
168*
169* .. Scalar Arguments ..
170 CHARACTER TRANA, TRANB
171 INTEGER INFO, ISGN, LDA, LDB, LDC, M, N
172 REAL SCALE
173* ..
174* .. Array Arguments ..
175 REAL A( LDA, * ), B( LDB, * ), C( LDC, * )
176* ..
177*
178* =====================================================================
179*
180* .. Parameters ..
181 REAL ZERO, ONE
182 parameter( zero = 0.0e+0, one = 1.0e+0 )
183* ..
184* .. Local Scalars ..
185 LOGICAL NOTRNA, NOTRNB
186 INTEGER IERR, J, K, K1, K2, KNEXT, L, L1, L2, LNEXT
187 REAL A11, BIGNUM, DA11, DB, EPS, SCALOC, SGN, SMIN,
188 $ SMLNUM, SUML, SUMR, XNORM
189* ..
190* .. Local Arrays ..
191 REAL DUM( 1 ), VEC( 2, 2 ), X( 2, 2 )
192* ..
193* .. External Functions ..
194 LOGICAL LSAME
195 REAL SDOT, SLAMCH, SLANGE
196 EXTERNAL lsame, sdot, slamch, slange
197* ..
198* .. External Subroutines ..
199 EXTERNAL slabad, slaln2, slasy2, sscal, xerbla
200* ..
201* .. Intrinsic Functions ..
202 INTRINSIC abs, max, min, real
203* ..
204* .. Executable Statements ..
205*
206* Decode and Test input parameters
207*
208 notrna = lsame( trana, 'N' )
209 notrnb = lsame( tranb, 'N' )
210*
211 info = 0
212 IF( .NOT.notrna .AND. .NOT.lsame( trana, 'T' ) .AND. .NOT.
213 $ lsame( trana, 'C' ) ) THEN
214 info = -1
215 ELSE IF( .NOT.notrnb .AND. .NOT.lsame( tranb, 'T' ) .AND. .NOT.
216 $ lsame( tranb, 'C' ) ) THEN
217 info = -2
218 ELSE IF( isgn.NE.1 .AND. isgn.NE.-1 ) THEN
219 info = -3
220 ELSE IF( m.LT.0 ) THEN
221 info = -4
222 ELSE IF( n.LT.0 ) THEN
223 info = -5
224 ELSE IF( lda.LT.max( 1, m ) ) THEN
225 info = -7
226 ELSE IF( ldb.LT.max( 1, n ) ) THEN
227 info = -9
228 ELSE IF( ldc.LT.max( 1, m ) ) THEN
229 info = -11
230 END IF
231 IF( info.NE.0 ) THEN
232 CALL xerbla( 'STRSYL', -info )
233 RETURN
234 END IF
235*
236* Quick return if possible
237*
238 scale = one
239 IF( m.EQ.0 .OR. n.EQ.0 )
240 $ RETURN
241*
242* Set constants to control overflow
243*
244 eps = slamch( 'P' )
245 smlnum = slamch( 'S' )
246 bignum = one / smlnum
247 CALL slabad( smlnum, bignum )
248 smlnum = smlnum*real( m*n ) / eps
249 bignum = one / smlnum
250*
251 smin = max( smlnum, eps*slange( 'M', m, m, a, lda, dum ),
252 $ eps*slange( 'M', n, n, b, ldb, dum ) )
253*
254 sgn = isgn
255*
256 IF( notrna .AND. notrnb ) THEN
257*
258* Solve A*X + ISGN*X*B = scale*C.
259*
260* The (K,L)th block of X is determined starting from
261* bottom-left corner column by column by
262*
263* A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
264*
265* Where
266* M L-1
267* R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(J,L)].
268* I=K+1 J=1
269*
270* Start column loop (index = L)
271* L1 (L2) : column index of the first (first) row of X(K,L).
272*
273 lnext = 1
274 DO 70 l = 1, n
275 IF( l.LT.lnext )
276 $ GO TO 70
277 IF( l.EQ.n ) THEN
278 l1 = l
279 l2 = l
280 ELSE
281 IF( b( l+1, l ).NE.zero ) THEN
282 l1 = l
283 l2 = l + 1
284 lnext = l + 2
285 ELSE
286 l1 = l
287 l2 = l
288 lnext = l + 1
289 END IF
290 END IF
291*
292* Start row loop (index = K)
293* K1 (K2): row index of the first (last) row of X(K,L).
294*
295 knext = m
296 DO 60 k = m, 1, -1
297 IF( k.GT.knext )
298 $ GO TO 60
299 IF( k.EQ.1 ) THEN
300 k1 = k
301 k2 = k
302 ELSE
303 IF( a( k, k-1 ).NE.zero ) THEN
304 k1 = k - 1
305 k2 = k
306 knext = k - 2
307 ELSE
308 k1 = k
309 k2 = k
310 knext = k - 1
311 END IF
312 END IF
313*
314 IF( l1.EQ.l2 .AND. k1.EQ.k2 ) THEN
315 suml = sdot( m-k1, a( k1, min( k1+1, m ) ), lda,
316 $ c( min( k1+1, m ), l1 ), 1 )
317 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
318 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
319 scaloc = one
320*
321 a11 = a( k1, k1 ) + sgn*b( l1, l1 )
322 da11 = abs( a11 )
323 IF( da11.LE.smin ) THEN
324 a11 = smin
325 da11 = smin
326 info = 1
327 END IF
328 db = abs( vec( 1, 1 ) )
329 IF( da11.LT.one .AND. db.GT.one ) THEN
330 IF( db.GT.bignum*da11 )
331 $ scaloc = one / db
332 END IF
333 x( 1, 1 ) = ( vec( 1, 1 )*scaloc ) / a11
334*
335 IF( scaloc.NE.one ) THEN
336 DO 10 j = 1, n
337 CALL sscal( m, scaloc, c( 1, j ), 1 )
338 10 CONTINUE
339 scale = scale*scaloc
340 END IF
341 c( k1, l1 ) = x( 1, 1 )
342*
343 ELSE IF( l1.EQ.l2 .AND. k1.NE.k2 ) THEN
344*
345 suml = sdot( m-k2, a( k1, min( k2+1, m ) ), lda,
346 $ c( min( k2+1, m ), l1 ), 1 )
347 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
348 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
349*
350 suml = sdot( m-k2, a( k2, min( k2+1, m ) ), lda,
351 $ c( min( k2+1, m ), l1 ), 1 )
352 sumr = sdot( l1-1, c( k2, 1 ), ldc, b( 1, l1 ), 1 )
353 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
354*
355 CALL slaln2( .false., 2, 1, smin, one, a( k1, k1 ),
356 $ lda, one, one, vec, 2, -sgn*b( l1, l1 ),
357 $ zero, x, 2, scaloc, xnorm, ierr )
358 IF( ierr.NE.0 )
359 $ info = 1
360*
361 IF( scaloc.NE.one ) THEN
362 DO 20 j = 1, n
363 CALL sscal( m, scaloc, c( 1, j ), 1 )
364 20 CONTINUE
365 scale = scale*scaloc
366 END IF
367 c( k1, l1 ) = x( 1, 1 )
368 c( k2, l1 ) = x( 2, 1 )
369*
370 ELSE IF( l1.NE.l2 .AND. k1.EQ.k2 ) THEN
371*
372 suml = sdot( m-k1, a( k1, min( k1+1, m ) ), lda,
373 $ c( min( k1+1, m ), l1 ), 1 )
374 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
375 vec( 1, 1 ) = sgn*( c( k1, l1 )-( suml+sgn*sumr ) )
376*
377 suml = sdot( m-k1, a( k1, min( k1+1, m ) ), lda,
378 $ c( min( k1+1, m ), l2 ), 1 )
379 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l2 ), 1 )
380 vec( 2, 1 ) = sgn*( c( k1, l2 )-( suml+sgn*sumr ) )
381*
382 CALL slaln2( .true., 2, 1, smin, one, b( l1, l1 ),
383 $ ldb, one, one, vec, 2, -sgn*a( k1, k1 ),
384 $ zero, x, 2, scaloc, xnorm, ierr )
385 IF( ierr.NE.0 )
386 $ info = 1
387*
388 IF( scaloc.NE.one ) THEN
389 DO 40 j = 1, n
390 CALL sscal( m, scaloc, c( 1, j ), 1 )
391 40 CONTINUE
392 scale = scale*scaloc
393 END IF
394 c( k1, l1 ) = x( 1, 1 )
395 c( k1, l2 ) = x( 2, 1 )
396*
397 ELSE IF( l1.NE.l2 .AND. k1.NE.k2 ) THEN
398*
399 suml = sdot( m-k2, a( k1, min( k2+1, m ) ), lda,
400 $ c( min( k2+1, m ), l1 ), 1 )
401 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
402 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
403*
404 suml = sdot( m-k2, a( k1, min( k2+1, m ) ), lda,
405 $ c( min( k2+1, m ), l2 ), 1 )
406 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l2 ), 1 )
407 vec( 1, 2 ) = c( k1, l2 ) - ( suml+sgn*sumr )
408*
409 suml = sdot( m-k2, a( k2, min( k2+1, m ) ), lda,
410 $ c( min( k2+1, m ), l1 ), 1 )
411 sumr = sdot( l1-1, c( k2, 1 ), ldc, b( 1, l1 ), 1 )
412 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
413*
414 suml = sdot( m-k2, a( k2, min( k2+1, m ) ), lda,
415 $ c( min( k2+1, m ), l2 ), 1 )
416 sumr = sdot( l1-1, c( k2, 1 ), ldc, b( 1, l2 ), 1 )
417 vec( 2, 2 ) = c( k2, l2 ) - ( suml+sgn*sumr )
418*
419 CALL slasy2( .false., .false., isgn, 2, 2,
420 $ a( k1, k1 ), lda, b( l1, l1 ), ldb, vec,
421 $ 2, scaloc, x, 2, xnorm, ierr )
422 IF( ierr.NE.0 )
423 $ info = 1
424*
425 IF( scaloc.NE.one ) THEN
426 DO 50 j = 1, n
427 CALL sscal( m, scaloc, c( 1, j ), 1 )
428 50 CONTINUE
429 scale = scale*scaloc
430 END IF
431 c( k1, l1 ) = x( 1, 1 )
432 c( k1, l2 ) = x( 1, 2 )
433 c( k2, l1 ) = x( 2, 1 )
434 c( k2, l2 ) = x( 2, 2 )
435 END IF
436*
437 60 CONTINUE
438*
439 70 CONTINUE
440*
441 ELSE IF( .NOT.notrna .AND. notrnb ) THEN
442*
443* Solve A**T *X + ISGN*X*B = scale*C.
444*
445* The (K,L)th block of X is determined starting from
446* upper-left corner column by column by
447*
448* A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
449*
450* Where
451* K-1 L-1
452* R(K,L) = SUM [A(I,K)**T*X(I,L)] +ISGN*SUM [X(K,J)*B(J,L)]
453* I=1 J=1
454*
455* Start column loop (index = L)
456* L1 (L2): column index of the first (last) row of X(K,L)
457*
458 lnext = 1
459 DO 130 l = 1, n
460 IF( l.LT.lnext )
461 $ GO TO 130
462 IF( l.EQ.n ) THEN
463 l1 = l
464 l2 = l
465 ELSE
466 IF( b( l+1, l ).NE.zero ) THEN
467 l1 = l
468 l2 = l + 1
469 lnext = l + 2
470 ELSE
471 l1 = l
472 l2 = l
473 lnext = l + 1
474 END IF
475 END IF
476*
477* Start row loop (index = K)
478* K1 (K2): row index of the first (last) row of X(K,L)
479*
480 knext = 1
481 DO 120 k = 1, m
482 IF( k.LT.knext )
483 $ GO TO 120
484 IF( k.EQ.m ) THEN
485 k1 = k
486 k2 = k
487 ELSE
488 IF( a( k+1, k ).NE.zero ) THEN
489 k1 = k
490 k2 = k + 1
491 knext = k + 2
492 ELSE
493 k1 = k
494 k2 = k
495 knext = k + 1
496 END IF
497 END IF
498*
499 IF( l1.EQ.l2 .AND. k1.EQ.k2 ) THEN
500 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
501 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
502 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
503 scaloc = one
504*
505 a11 = a( k1, k1 ) + sgn*b( l1, l1 )
506 da11 = abs( a11 )
507 IF( da11.LE.smin ) THEN
508 a11 = smin
509 da11 = smin
510 info = 1
511 END IF
512 db = abs( vec( 1, 1 ) )
513 IF( da11.LT.one .AND. db.GT.one ) THEN
514 IF( db.GT.bignum*da11 )
515 $ scaloc = one / db
516 END IF
517 x( 1, 1 ) = ( vec( 1, 1 )*scaloc ) / a11
518*
519 IF( scaloc.NE.one ) THEN
520 DO 80 j = 1, n
521 CALL sscal( m, scaloc, c( 1, j ), 1 )
522 80 CONTINUE
523 scale = scale*scaloc
524 END IF
525 c( k1, l1 ) = x( 1, 1 )
526*
527 ELSE IF( l1.EQ.l2 .AND. k1.NE.k2 ) THEN
528*
529 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
530 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
531 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
532*
533 suml = sdot( k1-1, a( 1, k2 ), 1, c( 1, l1 ), 1 )
534 sumr = sdot( l1-1, c( k2, 1 ), ldc, b( 1, l1 ), 1 )
535 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
536*
537 CALL slaln2( .true., 2, 1, smin, one, a( k1, k1 ),
538 $ lda, one, one, vec, 2, -sgn*b( l1, l1 ),
539 $ zero, x, 2, scaloc, xnorm, ierr )
540 IF( ierr.NE.0 )
541 $ info = 1
542*
543 IF( scaloc.NE.one ) THEN
544 DO 90 j = 1, n
545 CALL sscal( m, scaloc, c( 1, j ), 1 )
546 90 CONTINUE
547 scale = scale*scaloc
548 END IF
549 c( k1, l1 ) = x( 1, 1 )
550 c( k2, l1 ) = x( 2, 1 )
551*
552 ELSE IF( l1.NE.l2 .AND. k1.EQ.k2 ) THEN
553*
554 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
555 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
556 vec( 1, 1 ) = sgn*( c( k1, l1 )-( suml+sgn*sumr ) )
557*
558 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l2 ), 1 )
559 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l2 ), 1 )
560 vec( 2, 1 ) = sgn*( c( k1, l2 )-( suml+sgn*sumr ) )
561*
562 CALL slaln2( .true., 2, 1, smin, one, b( l1, l1 ),
563 $ ldb, one, one, vec, 2, -sgn*a( k1, k1 ),
564 $ zero, x, 2, scaloc, xnorm, ierr )
565 IF( ierr.NE.0 )
566 $ info = 1
567*
568 IF( scaloc.NE.one ) THEN
569 DO 100 j = 1, n
570 CALL sscal( m, scaloc, c( 1, j ), 1 )
571 100 CONTINUE
572 scale = scale*scaloc
573 END IF
574 c( k1, l1 ) = x( 1, 1 )
575 c( k1, l2 ) = x( 2, 1 )
576*
577 ELSE IF( l1.NE.l2 .AND. k1.NE.k2 ) THEN
578*
579 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
580 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
581 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
582*
583 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l2 ), 1 )
584 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l2 ), 1 )
585 vec( 1, 2 ) = c( k1, l2 ) - ( suml+sgn*sumr )
586*
587 suml = sdot( k1-1, a( 1, k2 ), 1, c( 1, l1 ), 1 )
588 sumr = sdot( l1-1, c( k2, 1 ), ldc, b( 1, l1 ), 1 )
589 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
590*
591 suml = sdot( k1-1, a( 1, k2 ), 1, c( 1, l2 ), 1 )
592 sumr = sdot( l1-1, c( k2, 1 ), ldc, b( 1, l2 ), 1 )
593 vec( 2, 2 ) = c( k2, l2 ) - ( suml+sgn*sumr )
594*
595 CALL slasy2( .true., .false., isgn, 2, 2, a( k1, k1 ),
596 $ lda, b( l1, l1 ), ldb, vec, 2, scaloc, x,
597 $ 2, xnorm, ierr )
598 IF( ierr.NE.0 )
599 $ info = 1
600*
601 IF( scaloc.NE.one ) THEN
602 DO 110 j = 1, n
603 CALL sscal( m, scaloc, c( 1, j ), 1 )
604 110 CONTINUE
605 scale = scale*scaloc
606 END IF
607 c( k1, l1 ) = x( 1, 1 )
608 c( k1, l2 ) = x( 1, 2 )
609 c( k2, l1 ) = x( 2, 1 )
610 c( k2, l2 ) = x( 2, 2 )
611 END IF
612*
613 120 CONTINUE
614 130 CONTINUE
615*
616 ELSE IF( .NOT.notrna .AND. .NOT.notrnb ) THEN
617*
618* Solve A**T*X + ISGN*X*B**T = scale*C.
619*
620* The (K,L)th block of X is determined starting from
621* top-right corner column by column by
622*
623* A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L)**T = C(K,L) - R(K,L)
624*
625* Where
626* K-1 N
627* R(K,L) = SUM [A(I,K)**T*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**T].
628* I=1 J=L+1
629*
630* Start column loop (index = L)
631* L1 (L2): column index of the first (last) row of X(K,L)
632*
633 lnext = n
634 DO 190 l = n, 1, -1
635 IF( l.GT.lnext )
636 $ GO TO 190
637 IF( l.EQ.1 ) THEN
638 l1 = l
639 l2 = l
640 ELSE
641 IF( b( l, l-1 ).NE.zero ) THEN
642 l1 = l - 1
643 l2 = l
644 lnext = l - 2
645 ELSE
646 l1 = l
647 l2 = l
648 lnext = l - 1
649 END IF
650 END IF
651*
652* Start row loop (index = K)
653* K1 (K2): row index of the first (last) row of X(K,L)
654*
655 knext = 1
656 DO 180 k = 1, m
657 IF( k.LT.knext )
658 $ GO TO 180
659 IF( k.EQ.m ) THEN
660 k1 = k
661 k2 = k
662 ELSE
663 IF( a( k+1, k ).NE.zero ) THEN
664 k1 = k
665 k2 = k + 1
666 knext = k + 2
667 ELSE
668 k1 = k
669 k2 = k
670 knext = k + 1
671 END IF
672 END IF
673*
674 IF( l1.EQ.l2 .AND. k1.EQ.k2 ) THEN
675 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
676 sumr = sdot( n-l1, c( k1, min( l1+1, n ) ), ldc,
677 $ b( l1, min( l1+1, n ) ), ldb )
678 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
679 scaloc = one
680*
681 a11 = a( k1, k1 ) + sgn*b( l1, l1 )
682 da11 = abs( a11 )
683 IF( da11.LE.smin ) THEN
684 a11 = smin
685 da11 = smin
686 info = 1
687 END IF
688 db = abs( vec( 1, 1 ) )
689 IF( da11.LT.one .AND. db.GT.one ) THEN
690 IF( db.GT.bignum*da11 )
691 $ scaloc = one / db
692 END IF
693 x( 1, 1 ) = ( vec( 1, 1 )*scaloc ) / a11
694*
695 IF( scaloc.NE.one ) THEN
696 DO 140 j = 1, n
697 CALL sscal( m, scaloc, c( 1, j ), 1 )
698 140 CONTINUE
699 scale = scale*scaloc
700 END IF
701 c( k1, l1 ) = x( 1, 1 )
702*
703 ELSE IF( l1.EQ.l2 .AND. k1.NE.k2 ) THEN
704*
705 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
706 sumr = sdot( n-l2, c( k1, min( l2+1, n ) ), ldc,
707 $ b( l1, min( l2+1, n ) ), ldb )
708 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
709*
710 suml = sdot( k1-1, a( 1, k2 ), 1, c( 1, l1 ), 1 )
711 sumr = sdot( n-l2, c( k2, min( l2+1, n ) ), ldc,
712 $ b( l1, min( l2+1, n ) ), ldb )
713 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
714*
715 CALL slaln2( .true., 2, 1, smin, one, a( k1, k1 ),
716 $ lda, one, one, vec, 2, -sgn*b( l1, l1 ),
717 $ zero, x, 2, scaloc, xnorm, ierr )
718 IF( ierr.NE.0 )
719 $ info = 1
720*
721 IF( scaloc.NE.one ) THEN
722 DO 150 j = 1, n
723 CALL sscal( m, scaloc, c( 1, j ), 1 )
724 150 CONTINUE
725 scale = scale*scaloc
726 END IF
727 c( k1, l1 ) = x( 1, 1 )
728 c( k2, l1 ) = x( 2, 1 )
729*
730 ELSE IF( l1.NE.l2 .AND. k1.EQ.k2 ) THEN
731*
732 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
733 sumr = sdot( n-l2, c( k1, min( l2+1, n ) ), ldc,
734 $ b( l1, min( l2+1, n ) ), ldb )
735 vec( 1, 1 ) = sgn*( c( k1, l1 )-( suml+sgn*sumr ) )
736*
737 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l2 ), 1 )
738 sumr = sdot( n-l2, c( k1, min( l2+1, n ) ), ldc,
739 $ b( l2, min( l2+1, n ) ), ldb )
740 vec( 2, 1 ) = sgn*( c( k1, l2 )-( suml+sgn*sumr ) )
741*
742 CALL slaln2( .false., 2, 1, smin, one, b( l1, l1 ),
743 $ ldb, one, one, vec, 2, -sgn*a( k1, k1 ),
744 $ zero, x, 2, scaloc, xnorm, ierr )
745 IF( ierr.NE.0 )
746 $ info = 1
747*
748 IF( scaloc.NE.one ) THEN
749 DO 160 j = 1, n
750 CALL sscal( m, scaloc, c( 1, j ), 1 )
751 160 CONTINUE
752 scale = scale*scaloc
753 END IF
754 c( k1, l1 ) = x( 1, 1 )
755 c( k1, l2 ) = x( 2, 1 )
756*
757 ELSE IF( l1.NE.l2 .AND. k1.NE.k2 ) THEN
758*
759 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
760 sumr = sdot( n-l2, c( k1, min( l2+1, n ) ), ldc,
761 $ b( l1, min( l2+1, n ) ), ldb )
762 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
763*
764 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l2 ), 1 )
765 sumr = sdot( n-l2, c( k1, min( l2+1, n ) ), ldc,
766 $ b( l2, min( l2+1, n ) ), ldb )
767 vec( 1, 2 ) = c( k1, l2 ) - ( suml+sgn*sumr )
768*
769 suml = sdot( k1-1, a( 1, k2 ), 1, c( 1, l1 ), 1 )
770 sumr = sdot( n-l2, c( k2, min( l2+1, n ) ), ldc,
771 $ b( l1, min( l2+1, n ) ), ldb )
772 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
773*
774 suml = sdot( k1-1, a( 1, k2 ), 1, c( 1, l2 ), 1 )
775 sumr = sdot( n-l2, c( k2, min( l2+1, n ) ), ldc,
776 $ b( l2, min(l2+1, n ) ), ldb )
777 vec( 2, 2 ) = c( k2, l2 ) - ( suml+sgn*sumr )
778*
779 CALL slasy2( .true., .true., isgn, 2, 2, a( k1, k1 ),
780 $ lda, b( l1, l1 ), ldb, vec, 2, scaloc, x,
781 $ 2, xnorm, ierr )
782 IF( ierr.NE.0 )
783 $ info = 1
784*
785 IF( scaloc.NE.one ) THEN
786 DO 170 j = 1, n
787 CALL sscal( m, scaloc, c( 1, j ), 1 )
788 170 CONTINUE
789 scale = scale*scaloc
790 END IF
791 c( k1, l1 ) = x( 1, 1 )
792 c( k1, l2 ) = x( 1, 2 )
793 c( k2, l1 ) = x( 2, 1 )
794 c( k2, l2 ) = x( 2, 2 )
795 END IF
796*
797 180 CONTINUE
798 190 CONTINUE
799*
800 ELSE IF( notrna .AND. .NOT.notrnb ) THEN
801*
802* Solve A*X + ISGN*X*B**T = scale*C.
803*
804* The (K,L)th block of X is determined starting from
805* bottom-right corner column by column by
806*
807* A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L)**T = C(K,L) - R(K,L)
808*
809* Where
810* M N
811* R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**T].
812* I=K+1 J=L+1
813*
814* Start column loop (index = L)
815* L1 (L2): column index of the first (last) row of X(K,L)
816*
817 lnext = n
818 DO 250 l = n, 1, -1
819 IF( l.GT.lnext )
820 $ GO TO 250
821 IF( l.EQ.1 ) THEN
822 l1 = l
823 l2 = l
824 ELSE
825 IF( b( l, l-1 ).NE.zero ) THEN
826 l1 = l - 1
827 l2 = l
828 lnext = l - 2
829 ELSE
830 l1 = l
831 l2 = l
832 lnext = l - 1
833 END IF
834 END IF
835*
836* Start row loop (index = K)
837* K1 (K2): row index of the first (last) row of X(K,L)
838*
839 knext = m
840 DO 240 k = m, 1, -1
841 IF( k.GT.knext )
842 $ GO TO 240
843 IF( k.EQ.1 ) THEN
844 k1 = k
845 k2 = k
846 ELSE
847 IF( a( k, k-1 ).NE.zero ) THEN
848 k1 = k - 1
849 k2 = k
850 knext = k - 2
851 ELSE
852 k1 = k
853 k2 = k
854 knext = k - 1
855 END IF
856 END IF
857*
858 IF( l1.EQ.l2 .AND. k1.EQ.k2 ) THEN
859 suml = sdot( m-k1, a( k1, min(k1+1, m ) ), lda,
860 $ c( min( k1+1, m ), l1 ), 1 )
861 sumr = sdot( n-l1, c( k1, min( l1+1, n ) ), ldc,
862 $ b( l1, min( l1+1, n ) ), ldb )
863 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
864 scaloc = one
865*
866 a11 = a( k1, k1 ) + sgn*b( l1, l1 )
867 da11 = abs( a11 )
868 IF( da11.LE.smin ) THEN
869 a11 = smin
870 da11 = smin
871 info = 1
872 END IF
873 db = abs( vec( 1, 1 ) )
874 IF( da11.LT.one .AND. db.GT.one ) THEN
875 IF( db.GT.bignum*da11 )
876 $ scaloc = one / db
877 END IF
878 x( 1, 1 ) = ( vec( 1, 1 )*scaloc ) / a11
879*
880 IF( scaloc.NE.one ) THEN
881 DO 200 j = 1, n
882 CALL sscal( m, scaloc, c( 1, j ), 1 )
883 200 CONTINUE
884 scale = scale*scaloc
885 END IF
886 c( k1, l1 ) = x( 1, 1 )
887*
888 ELSE IF( l1.EQ.l2 .AND. k1.NE.k2 ) THEN
889*
890 suml = sdot( m-k2, a( k1, min( k2+1, m ) ), lda,
891 $ c( min( k2+1, m ), l1 ), 1 )
892 sumr = sdot( n-l2, c( k1, min( l2+1, n ) ), ldc,
893 $ b( l1, min( l2+1, n ) ), ldb )
894 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
895*
896 suml = sdot( m-k2, a( k2, min( k2+1, m ) ), lda,
897 $ c( min( k2+1, m ), l1 ), 1 )
898 sumr = sdot( n-l2, c( k2, min( l2+1, n ) ), ldc,
899 $ b( l1, min( l2+1, n ) ), ldb )
900 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
901*
902 CALL slaln2( .false., 2, 1, smin, one, a( k1, k1 ),
903 $ lda, one, one, vec, 2, -sgn*b( l1, l1 ),
904 $ zero, x, 2, scaloc, xnorm, ierr )
905 IF( ierr.NE.0 )
906 $ info = 1
907*
908 IF( scaloc.NE.one ) THEN
909 DO 210 j = 1, n
910 CALL sscal( m, scaloc, c( 1, j ), 1 )
911 210 CONTINUE
912 scale = scale*scaloc
913 END IF
914 c( k1, l1 ) = x( 1, 1 )
915 c( k2, l1 ) = x( 2, 1 )
916*
917 ELSE IF( l1.NE.l2 .AND. k1.EQ.k2 ) THEN
918*
919 suml = sdot( m-k1, a( k1, min( k1+1, m ) ), lda,
920 $ c( min( k1+1, m ), l1 ), 1 )
921 sumr = sdot( n-l2, c( k1, min( l2+1, n ) ), ldc,
922 $ b( l1, min( l2+1, n ) ), ldb )
923 vec( 1, 1 ) = sgn*( c( k1, l1 )-( suml+sgn*sumr ) )
924*
925 suml = sdot( m-k1, a( k1, min( k1+1, m ) ), lda,
926 $ c( min( k1+1, m ), l2 ), 1 )
927 sumr = sdot( n-l2, c( k1, min( l2+1, n ) ), ldc,
928 $ b( l2, min( l2+1, n ) ), ldb )
929 vec( 2, 1 ) = sgn*( c( k1, l2 )-( suml+sgn*sumr ) )
930*
931 CALL slaln2( .false., 2, 1, smin, one, b( l1, l1 ),
932 $ ldb, one, one, vec, 2, -sgn*a( k1, k1 ),
933 $ zero, x, 2, scaloc, xnorm, ierr )
934 IF( ierr.NE.0 )
935 $ info = 1
936*
937 IF( scaloc.NE.one ) THEN
938 DO 220 j = 1, n
939 CALL sscal( m, scaloc, c( 1, j ), 1 )
940 220 CONTINUE
941 scale = scale*scaloc
942 END IF
943 c( k1, l1 ) = x( 1, 1 )
944 c( k1, l2 ) = x( 2, 1 )
945*
946 ELSE IF( l1.NE.l2 .AND. k1.NE.k2 ) THEN
947*
948 suml = sdot( m-k2, a( k1, min( k2+1, m ) ), lda,
949 $ c( min( k2+1, m ), l1 ), 1 )
950 sumr = sdot( n-l2, c( k1, min( l2+1, n ) ), ldc,
951 $ b( l1, min( l2+1, n ) ), ldb )
952 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
953*
954 suml = sdot( m-k2, a( k1, min( k2+1, m ) ), lda,
955 $ c( min( k2+1, m ), l2 ), 1 )
956 sumr = sdot( n-l2, c( k1, min( l2+1, n ) ), ldc,
957 $ b( l2, min( l2+1, n ) ), ldb )
958 vec( 1, 2 ) = c( k1, l2 ) - ( suml+sgn*sumr )
959*
960 suml = sdot( m-k2, a( k2, min( k2+1, m ) ), lda,
961 $ c( min( k2+1, m ), l1 ), 1 )
962 sumr = sdot( n-l2, c( k2, min( l2+1, n ) ), ldc,
963 $ b( l1, min( l2+1, n ) ), ldb )
964 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
965*
966 suml = sdot( m-k2, a( k2, min( k2+1, m ) ), lda,
967 $ c( min( k2+1, m ), l2 ), 1 )
968 sumr = sdot( n-l2, c( k2, min( l2+1, n ) ), ldc,
969 $ b( l2, min( l2+1, n ) ), ldb )
970 vec( 2, 2 ) = c( k2, l2 ) - ( suml+sgn*sumr )
971*
972 CALL slasy2( .false., .true., isgn, 2, 2, a( k1, 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*
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
real function slange(NORM, M, N, A, LDA, WORK)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: slange.f:114
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:218
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:174
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
real function sdot(N, SX, INCX, SY, INCY)
SDOT
Definition: sdot.f:82
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: