LAPACK 3.12.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 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 smlnum = smlnum*real( m*n ) / eps
248 bignum = one / smlnum
249*
250 smin = max( smlnum, eps*slange( 'M', m, m, a, lda, dum ),
251 $ eps*slange( 'M', n, n, b, ldb, dum ) )
252*
253 sgn = isgn
254*
255 IF( notrna .AND. notrnb ) THEN
256*
257* Solve A*X + ISGN*X*B = scale*C.
258*
259* The (K,L)th block of X is determined starting from
260* bottom-left corner column by column by
261*
262* A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
263*
264* Where
265* M L-1
266* R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(J,L)].
267* I=K+1 J=1
268*
269* Start column loop (index = L)
270* L1 (L2) : column index of the first (first) row of X(K,L).
271*
272 lnext = 1
273 DO 70 l = 1, n
274 IF( l.LT.lnext )
275 $ GO TO 70
276 IF( l.EQ.n ) THEN
277 l1 = l
278 l2 = l
279 ELSE
280 IF( b( l+1, l ).NE.zero ) THEN
281 l1 = l
282 l2 = l + 1
283 lnext = l + 2
284 ELSE
285 l1 = l
286 l2 = l
287 lnext = l + 1
288 END IF
289 END IF
290*
291* Start row loop (index = K)
292* K1 (K2): row index of the first (last) row of X(K,L).
293*
294 knext = m
295 DO 60 k = m, 1, -1
296 IF( k.GT.knext )
297 $ GO TO 60
298 IF( k.EQ.1 ) THEN
299 k1 = k
300 k2 = k
301 ELSE
302 IF( a( k, k-1 ).NE.zero ) THEN
303 k1 = k - 1
304 k2 = k
305 knext = k - 2
306 ELSE
307 k1 = k
308 k2 = k
309 knext = k - 1
310 END IF
311 END IF
312*
313 IF( l1.EQ.l2 .AND. k1.EQ.k2 ) THEN
314 suml = sdot( m-k1, a( k1, min( k1+1, m ) ), lda,
315 $ c( min( k1+1, m ), l1 ), 1 )
316 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
317 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
318 scaloc = one
319*
320 a11 = a( k1, k1 ) + sgn*b( l1, l1 )
321 da11 = abs( a11 )
322 IF( da11.LE.smin ) THEN
323 a11 = smin
324 da11 = smin
325 info = 1
326 END IF
327 db = abs( vec( 1, 1 ) )
328 IF( da11.LT.one .AND. db.GT.one ) THEN
329 IF( db.GT.bignum*da11 )
330 $ scaloc = one / db
331 END IF
332 x( 1, 1 ) = ( vec( 1, 1 )*scaloc ) / a11
333*
334 IF( scaloc.NE.one ) THEN
335 DO 10 j = 1, n
336 CALL sscal( m, scaloc, c( 1, j ), 1 )
337 10 CONTINUE
338 scale = scale*scaloc
339 END IF
340 c( k1, l1 ) = x( 1, 1 )
341*
342 ELSE IF( l1.EQ.l2 .AND. k1.NE.k2 ) THEN
343*
344 suml = sdot( m-k2, a( k1, min( k2+1, m ) ), lda,
345 $ c( min( k2+1, m ), l1 ), 1 )
346 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
347 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
348*
349 suml = sdot( m-k2, a( k2, min( k2+1, m ) ), lda,
350 $ c( min( k2+1, m ), l1 ), 1 )
351 sumr = sdot( l1-1, c( k2, 1 ), ldc, b( 1, l1 ), 1 )
352 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
353*
354 CALL slaln2( .false., 2, 1, smin, one, a( k1, k1 ),
355 $ lda, one, one, vec, 2, -sgn*b( l1, l1 ),
356 $ zero, x, 2, scaloc, xnorm, ierr )
357 IF( ierr.NE.0 )
358 $ info = 1
359*
360 IF( scaloc.NE.one ) THEN
361 DO 20 j = 1, n
362 CALL sscal( m, scaloc, c( 1, j ), 1 )
363 20 CONTINUE
364 scale = scale*scaloc
365 END IF
366 c( k1, l1 ) = x( 1, 1 )
367 c( k2, l1 ) = x( 2, 1 )
368*
369 ELSE IF( l1.NE.l2 .AND. k1.EQ.k2 ) THEN
370*
371 suml = sdot( m-k1, a( k1, min( k1+1, m ) ), lda,
372 $ c( min( k1+1, m ), l1 ), 1 )
373 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
374 vec( 1, 1 ) = sgn*( c( k1, l1 )-( suml+sgn*sumr ) )
375*
376 suml = sdot( m-k1, a( k1, min( k1+1, m ) ), lda,
377 $ c( min( k1+1, m ), l2 ), 1 )
378 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l2 ), 1 )
379 vec( 2, 1 ) = sgn*( c( k1, l2 )-( suml+sgn*sumr ) )
380*
381 CALL slaln2( .true., 2, 1, smin, one, b( l1, l1 ),
382 $ ldb, one, one, vec, 2, -sgn*a( k1, k1 ),
383 $ zero, x, 2, scaloc, xnorm, ierr )
384 IF( ierr.NE.0 )
385 $ info = 1
386*
387 IF( scaloc.NE.one ) THEN
388 DO 40 j = 1, n
389 CALL sscal( m, scaloc, c( 1, j ), 1 )
390 40 CONTINUE
391 scale = scale*scaloc
392 END IF
393 c( k1, l1 ) = x( 1, 1 )
394 c( k1, l2 ) = x( 2, 1 )
395*
396 ELSE IF( l1.NE.l2 .AND. k1.NE.k2 ) THEN
397*
398 suml = sdot( m-k2, a( k1, min( k2+1, m ) ), lda,
399 $ c( min( k2+1, m ), l1 ), 1 )
400 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
401 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
402*
403 suml = sdot( m-k2, a( k1, min( k2+1, m ) ), lda,
404 $ c( min( k2+1, m ), l2 ), 1 )
405 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l2 ), 1 )
406 vec( 1, 2 ) = c( k1, l2 ) - ( suml+sgn*sumr )
407*
408 suml = sdot( m-k2, a( k2, min( k2+1, m ) ), lda,
409 $ c( min( k2+1, m ), l1 ), 1 )
410 sumr = sdot( l1-1, c( k2, 1 ), ldc, b( 1, l1 ), 1 )
411 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
412*
413 suml = sdot( m-k2, a( k2, min( k2+1, m ) ), lda,
414 $ c( min( k2+1, m ), l2 ), 1 )
415 sumr = sdot( l1-1, c( k2, 1 ), ldc, b( 1, l2 ), 1 )
416 vec( 2, 2 ) = c( k2, l2 ) - ( suml+sgn*sumr )
417*
418 CALL slasy2( .false., .false., isgn, 2, 2,
419 $ a( k1, k1 ), lda, b( l1, l1 ), ldb, vec,
420 $ 2, scaloc, x, 2, xnorm, ierr )
421 IF( ierr.NE.0 )
422 $ info = 1
423*
424 IF( scaloc.NE.one ) THEN
425 DO 50 j = 1, n
426 CALL sscal( m, scaloc, c( 1, j ), 1 )
427 50 CONTINUE
428 scale = scale*scaloc
429 END IF
430 c( k1, l1 ) = x( 1, 1 )
431 c( k1, l2 ) = x( 1, 2 )
432 c( k2, l1 ) = x( 2, 1 )
433 c( k2, l2 ) = x( 2, 2 )
434 END IF
435*
436 60 CONTINUE
437*
438 70 CONTINUE
439*
440 ELSE IF( .NOT.notrna .AND. notrnb ) THEN
441*
442* Solve A**T *X + ISGN*X*B = scale*C.
443*
444* The (K,L)th block of X is determined starting from
445* upper-left corner column by column by
446*
447* A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
448*
449* Where
450* K-1 L-1
451* R(K,L) = SUM [A(I,K)**T*X(I,L)] +ISGN*SUM [X(K,J)*B(J,L)]
452* I=1 J=1
453*
454* Start column loop (index = L)
455* L1 (L2): column index of the first (last) row of X(K,L)
456*
457 lnext = 1
458 DO 130 l = 1, n
459 IF( l.LT.lnext )
460 $ GO TO 130
461 IF( l.EQ.n ) THEN
462 l1 = l
463 l2 = l
464 ELSE
465 IF( b( l+1, l ).NE.zero ) THEN
466 l1 = l
467 l2 = l + 1
468 lnext = l + 2
469 ELSE
470 l1 = l
471 l2 = l
472 lnext = l + 1
473 END IF
474 END IF
475*
476* Start row loop (index = K)
477* K1 (K2): row index of the first (last) row of X(K,L)
478*
479 knext = 1
480 DO 120 k = 1, m
481 IF( k.LT.knext )
482 $ GO TO 120
483 IF( k.EQ.m ) THEN
484 k1 = k
485 k2 = k
486 ELSE
487 IF( a( k+1, k ).NE.zero ) THEN
488 k1 = k
489 k2 = k + 1
490 knext = k + 2
491 ELSE
492 k1 = k
493 k2 = k
494 knext = k + 1
495 END IF
496 END IF
497*
498 IF( l1.EQ.l2 .AND. k1.EQ.k2 ) THEN
499 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
500 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
501 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
502 scaloc = one
503*
504 a11 = a( k1, k1 ) + sgn*b( l1, l1 )
505 da11 = abs( a11 )
506 IF( da11.LE.smin ) THEN
507 a11 = smin
508 da11 = smin
509 info = 1
510 END IF
511 db = abs( vec( 1, 1 ) )
512 IF( da11.LT.one .AND. db.GT.one ) THEN
513 IF( db.GT.bignum*da11 )
514 $ scaloc = one / db
515 END IF
516 x( 1, 1 ) = ( vec( 1, 1 )*scaloc ) / a11
517*
518 IF( scaloc.NE.one ) THEN
519 DO 80 j = 1, n
520 CALL sscal( m, scaloc, c( 1, j ), 1 )
521 80 CONTINUE
522 scale = scale*scaloc
523 END IF
524 c( k1, l1 ) = x( 1, 1 )
525*
526 ELSE IF( l1.EQ.l2 .AND. k1.NE.k2 ) THEN
527*
528 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
529 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
530 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
531*
532 suml = sdot( k1-1, a( 1, k2 ), 1, c( 1, l1 ), 1 )
533 sumr = sdot( l1-1, c( k2, 1 ), ldc, b( 1, l1 ), 1 )
534 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
535*
536 CALL slaln2( .true., 2, 1, smin, one, a( k1, k1 ),
537 $ lda, one, one, vec, 2, -sgn*b( l1, l1 ),
538 $ zero, x, 2, scaloc, xnorm, ierr )
539 IF( ierr.NE.0 )
540 $ info = 1
541*
542 IF( scaloc.NE.one ) THEN
543 DO 90 j = 1, n
544 CALL sscal( m, scaloc, c( 1, j ), 1 )
545 90 CONTINUE
546 scale = scale*scaloc
547 END IF
548 c( k1, l1 ) = x( 1, 1 )
549 c( k2, l1 ) = x( 2, 1 )
550*
551 ELSE IF( l1.NE.l2 .AND. k1.EQ.k2 ) THEN
552*
553 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
554 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
555 vec( 1, 1 ) = sgn*( c( k1, l1 )-( suml+sgn*sumr ) )
556*
557 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l2 ), 1 )
558 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l2 ), 1 )
559 vec( 2, 1 ) = sgn*( c( k1, l2 )-( suml+sgn*sumr ) )
560*
561 CALL slaln2( .true., 2, 1, smin, one, b( l1, l1 ),
562 $ ldb, one, one, vec, 2, -sgn*a( k1, k1 ),
563 $ zero, x, 2, scaloc, xnorm, ierr )
564 IF( ierr.NE.0 )
565 $ info = 1
566*
567 IF( scaloc.NE.one ) THEN
568 DO 100 j = 1, n
569 CALL sscal( m, scaloc, c( 1, j ), 1 )
570 100 CONTINUE
571 scale = scale*scaloc
572 END IF
573 c( k1, l1 ) = x( 1, 1 )
574 c( k1, l2 ) = x( 2, 1 )
575*
576 ELSE IF( l1.NE.l2 .AND. k1.NE.k2 ) THEN
577*
578 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
579 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
580 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
581*
582 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l2 ), 1 )
583 sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l2 ), 1 )
584 vec( 1, 2 ) = c( k1, l2 ) - ( suml+sgn*sumr )
585*
586 suml = sdot( k1-1, a( 1, k2 ), 1, c( 1, l1 ), 1 )
587 sumr = sdot( l1-1, c( k2, 1 ), ldc, b( 1, l1 ), 1 )
588 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
589*
590 suml = sdot( k1-1, a( 1, k2 ), 1, c( 1, l2 ), 1 )
591 sumr = sdot( l1-1, c( k2, 1 ), ldc, b( 1, l2 ), 1 )
592 vec( 2, 2 ) = c( k2, l2 ) - ( suml+sgn*sumr )
593*
594 CALL slasy2( .true., .false., isgn, 2, 2, a( k1, k1 ),
595 $ lda, b( l1, l1 ), ldb, vec, 2, scaloc, x,
596 $ 2, xnorm, ierr )
597 IF( ierr.NE.0 )
598 $ info = 1
599*
600 IF( scaloc.NE.one ) THEN
601 DO 110 j = 1, n
602 CALL sscal( m, scaloc, c( 1, j ), 1 )
603 110 CONTINUE
604 scale = scale*scaloc
605 END IF
606 c( k1, l1 ) = x( 1, 1 )
607 c( k1, l2 ) = x( 1, 2 )
608 c( k2, l1 ) = x( 2, 1 )
609 c( k2, l2 ) = x( 2, 2 )
610 END IF
611*
612 120 CONTINUE
613 130 CONTINUE
614*
615 ELSE IF( .NOT.notrna .AND. .NOT.notrnb ) THEN
616*
617* Solve A**T*X + ISGN*X*B**T = scale*C.
618*
619* The (K,L)th block of X is determined starting from
620* top-right corner column by column by
621*
622* A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L)**T = C(K,L) - R(K,L)
623*
624* Where
625* K-1 N
626* R(K,L) = SUM [A(I,K)**T*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**T].
627* I=1 J=L+1
628*
629* Start column loop (index = L)
630* L1 (L2): column index of the first (last) row of X(K,L)
631*
632 lnext = n
633 DO 190 l = n, 1, -1
634 IF( l.GT.lnext )
635 $ GO TO 190
636 IF( l.EQ.1 ) THEN
637 l1 = l
638 l2 = l
639 ELSE
640 IF( b( l, l-1 ).NE.zero ) THEN
641 l1 = l - 1
642 l2 = l
643 lnext = l - 2
644 ELSE
645 l1 = l
646 l2 = l
647 lnext = l - 1
648 END IF
649 END IF
650*
651* Start row loop (index = K)
652* K1 (K2): row index of the first (last) row of X(K,L)
653*
654 knext = 1
655 DO 180 k = 1, m
656 IF( k.LT.knext )
657 $ GO TO 180
658 IF( k.EQ.m ) THEN
659 k1 = k
660 k2 = k
661 ELSE
662 IF( a( k+1, k ).NE.zero ) THEN
663 k1 = k
664 k2 = k + 1
665 knext = k + 2
666 ELSE
667 k1 = k
668 k2 = k
669 knext = k + 1
670 END IF
671 END IF
672*
673 IF( l1.EQ.l2 .AND. k1.EQ.k2 ) THEN
674 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
675 sumr = sdot( n-l1, c( k1, min( l1+1, n ) ), ldc,
676 $ b( l1, min( l1+1, n ) ), ldb )
677 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
678 scaloc = one
679*
680 a11 = a( k1, k1 ) + sgn*b( l1, l1 )
681 da11 = abs( a11 )
682 IF( da11.LE.smin ) THEN
683 a11 = smin
684 da11 = smin
685 info = 1
686 END IF
687 db = abs( vec( 1, 1 ) )
688 IF( da11.LT.one .AND. db.GT.one ) THEN
689 IF( db.GT.bignum*da11 )
690 $ scaloc = one / db
691 END IF
692 x( 1, 1 ) = ( vec( 1, 1 )*scaloc ) / a11
693*
694 IF( scaloc.NE.one ) THEN
695 DO 140 j = 1, n
696 CALL sscal( m, scaloc, c( 1, j ), 1 )
697 140 CONTINUE
698 scale = scale*scaloc
699 END IF
700 c( k1, l1 ) = x( 1, 1 )
701*
702 ELSE IF( l1.EQ.l2 .AND. k1.NE.k2 ) THEN
703*
704 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
705 sumr = sdot( n-l2, c( k1, min( l2+1, n ) ), ldc,
706 $ b( l1, min( l2+1, n ) ), ldb )
707 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
708*
709 suml = sdot( k1-1, a( 1, k2 ), 1, c( 1, l1 ), 1 )
710 sumr = sdot( n-l2, c( k2, min( l2+1, n ) ), ldc,
711 $ b( l1, min( l2+1, n ) ), ldb )
712 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
713*
714 CALL slaln2( .true., 2, 1, smin, one, a( k1, k1 ),
715 $ lda, one, one, vec, 2, -sgn*b( l1, l1 ),
716 $ zero, x, 2, scaloc, xnorm, ierr )
717 IF( ierr.NE.0 )
718 $ info = 1
719*
720 IF( scaloc.NE.one ) THEN
721 DO 150 j = 1, n
722 CALL sscal( m, scaloc, c( 1, j ), 1 )
723 150 CONTINUE
724 scale = scale*scaloc
725 END IF
726 c( k1, l1 ) = x( 1, 1 )
727 c( k2, l1 ) = x( 2, 1 )
728*
729 ELSE IF( l1.NE.l2 .AND. k1.EQ.k2 ) THEN
730*
731 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
732 sumr = sdot( n-l2, c( k1, min( l2+1, n ) ), ldc,
733 $ b( l1, min( l2+1, n ) ), ldb )
734 vec( 1, 1 ) = sgn*( c( k1, l1 )-( suml+sgn*sumr ) )
735*
736 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l2 ), 1 )
737 sumr = sdot( n-l2, c( k1, min( l2+1, n ) ), ldc,
738 $ b( l2, min( l2+1, n ) ), ldb )
739 vec( 2, 1 ) = sgn*( c( k1, l2 )-( suml+sgn*sumr ) )
740*
741 CALL slaln2( .false., 2, 1, smin, one, b( l1, l1 ),
742 $ ldb, one, one, vec, 2, -sgn*a( k1, k1 ),
743 $ zero, x, 2, scaloc, xnorm, ierr )
744 IF( ierr.NE.0 )
745 $ info = 1
746*
747 IF( scaloc.NE.one ) THEN
748 DO 160 j = 1, n
749 CALL sscal( m, scaloc, c( 1, j ), 1 )
750 160 CONTINUE
751 scale = scale*scaloc
752 END IF
753 c( k1, l1 ) = x( 1, 1 )
754 c( k1, l2 ) = x( 2, 1 )
755*
756 ELSE IF( l1.NE.l2 .AND. k1.NE.k2 ) THEN
757*
758 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
759 sumr = sdot( n-l2, c( k1, min( l2+1, n ) ), ldc,
760 $ b( l1, min( l2+1, n ) ), ldb )
761 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
762*
763 suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l2 ), 1 )
764 sumr = sdot( n-l2, c( k1, min( l2+1, n ) ), ldc,
765 $ b( l2, min( l2+1, n ) ), ldb )
766 vec( 1, 2 ) = c( k1, l2 ) - ( suml+sgn*sumr )
767*
768 suml = sdot( k1-1, a( 1, k2 ), 1, c( 1, l1 ), 1 )
769 sumr = sdot( n-l2, c( k2, min( l2+1, n ) ), ldc,
770 $ b( l1, min( l2+1, n ) ), ldb )
771 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
772*
773 suml = sdot( k1-1, a( 1, k2 ), 1, c( 1, l2 ), 1 )
774 sumr = sdot( n-l2, c( k2, min( l2+1, n ) ), ldc,
775 $ b( l2, min(l2+1, n ) ), ldb )
776 vec( 2, 2 ) = c( k2, l2 ) - ( suml+sgn*sumr )
777*
778 CALL slasy2( .true., .true., isgn, 2, 2, a( k1, 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, k1 ),
972 $ lda, b( l1, l1 ), ldb, vec, 2, scaloc, x,
973 $ 2, xnorm, ierr )
974 IF( ierr.NE.0 )
975 $ info = 1
976*
977 IF( scaloc.NE.one ) THEN
978 DO 230 j = 1, n
979 CALL sscal( m, scaloc, c( 1, j ), 1 )
980 230 CONTINUE
981 scale = scale*scaloc
982 END IF
983 c( k1, l1 ) = x( 1, 1 )
984 c( k1, l2 ) = x( 1, 2 )
985 c( k2, l1 ) = x( 2, 1 )
986 c( k2, l2 ) = x( 2, 2 )
987 END IF
988*
989 240 CONTINUE
990 250 CONTINUE
991*
992 END IF
993*
994 RETURN
995*
996* End of STRSYL
997*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
real function sdot(n, sx, incx, sy, incy)
SDOT
Definition sdot.f:82
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
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
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 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
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
Here is the call graph for this function:
Here is the caller graph for this function: