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

◆ strsyl3()

subroutine strsyl3 ( 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, dimension( * )  IWORK,
integer  LIWORK,
real, dimension( ldswork, * )  SWORK,
integer  LDSWORK,
integer  INFO 
)

STRSYL3

Purpose
  STRSYL3 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.

  This is the block version of the algorithm.
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]IWORK
          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
[in]LIWORK
          IWORK is INTEGER
          The dimension of the array IWORK. LIWORK >=  ((M + NB - 1) / NB + 1)
          + ((N + NB - 1) / NB + 1), where NB is the optimal block size.

          If LIWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal dimension of the IWORK array,
          returns this value as the first entry of the IWORK array, and
          no error message related to LIWORK is issued by XERBLA.
[out]SWORK
          SWORK is REAL array, dimension (MAX(2, ROWS),
          MAX(1,COLS)).
          On exit, if INFO = 0, SWORK(1) returns the optimal value ROWS
          and SWORK(2) returns the optimal COLS.
[in]LDSWORK
          LDSWORK is INTEGER
          LDSWORK >= MAX(2,ROWS), where ROWS = ((M + NB - 1) / NB + 1)
          and NB is the optimal block size.

          If LDSWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal dimensions of the SWORK matrix,
          returns these values as the first and second entry of the SWORK
          matrix, and no error message related LWORK is issued by XERBLA.
[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).

Definition at line 178 of file strsyl3.f.

181 IMPLICIT NONE
182*
183* .. Scalar Arguments ..
184 CHARACTER TRANA, TRANB
185 INTEGER INFO, ISGN, LDA, LDB, LDC, M, N,
186 $ LIWORK, LDSWORK
187 REAL SCALE
188* ..
189* .. Array Arguments ..
190 INTEGER IWORK( * )
191 REAL A( LDA, * ), B( LDB, * ), C( LDC, * ),
192 $ SWORK( LDSWORK, * )
193* ..
194* .. Parameters ..
195 REAL ZERO, ONE
196 parameter( zero = 0.0e+0, one = 1.0e+0 )
197* ..
198* .. Local Scalars ..
199 LOGICAL NOTRNA, NOTRNB, LQUERY, SKIP
200 INTEGER AWRK, BWRK, I, I1, I2, IINFO, J, J1, J2, JJ,
201 $ K, K1, K2, L, L1, L2, LL, NBA, NB, NBB, PC
202 REAL ANRM, BIGNUM, BNRM, CNRM, SCAL, SCALOC,
203 $ SCAMIN, SGN, XNRM, BUF, SMLNUM
204* ..
205* .. Local Arrays ..
206 REAL WNRM( MAX( M, N ) )
207* ..
208* .. External Functions ..
209 LOGICAL LSAME
210 INTEGER ILAENV
211 REAL SLANGE, SLAMCH, SLARMM
212 EXTERNAL slange, slamch, slarmm, ilaenv, lsame
213* ..
214* .. External Subroutines ..
215 EXTERNAL sgemm, slascl, sscal, strsyl, xerbla
216* ..
217* .. Intrinsic Functions ..
218 INTRINSIC abs, exponent, max, min, real
219* ..
220* .. Executable Statements ..
221*
222* Decode and Test input parameters
223*
224 notrna = lsame( trana, 'N' )
225 notrnb = lsame( tranb, 'N' )
226*
227* Use the same block size for all matrices.
228*
229 nb = max(8, ilaenv( 1, 'STRSYL', '', m, n, -1, -1) )
230*
231* Compute number of blocks in A and B
232*
233 nba = max( 1, (m + nb - 1) / nb )
234 nbb = max( 1, (n + nb - 1) / nb )
235*
236* Compute workspace
237*
238 info = 0
239 lquery = ( liwork.EQ.-1 .OR. ldswork.EQ.-1 )
240 iwork( 1 ) = nba + nbb + 2
241 IF( lquery ) THEN
242 ldswork = 2
243 swork( 1, 1 ) = max( nba, nbb )
244 swork( 2, 1 ) = 2 * nbb + nba
245 END IF
246*
247* Test the input arguments
248*
249 IF( .NOT.notrna .AND. .NOT.lsame( trana, 'T' ) .AND. .NOT.
250 $ lsame( trana, 'C' ) ) THEN
251 info = -1
252 ELSE IF( .NOT.notrnb .AND. .NOT.lsame( tranb, 'T' ) .AND. .NOT.
253 $ lsame( tranb, 'C' ) ) THEN
254 info = -2
255 ELSE IF( isgn.NE.1 .AND. isgn.NE.-1 ) THEN
256 info = -3
257 ELSE IF( m.LT.0 ) THEN
258 info = -4
259 ELSE IF( n.LT.0 ) THEN
260 info = -5
261 ELSE IF( lda.LT.max( 1, m ) ) THEN
262 info = -7
263 ELSE IF( ldb.LT.max( 1, n ) ) THEN
264 info = -9
265 ELSE IF( ldc.LT.max( 1, m ) ) THEN
266 info = -11
267 ELSE IF( .NOT.lquery .AND. liwork.LT.iwork(1) ) THEN
268 info = -14
269 ELSE IF( .NOT.lquery .AND. ldswork.LT.max( nba, nbb ) ) THEN
270 info = -16
271 END IF
272 IF( info.NE.0 ) THEN
273 CALL xerbla( 'STRSYL3', -info )
274 RETURN
275 ELSE IF( lquery ) THEN
276 RETURN
277 END IF
278*
279* Quick return if possible
280*
281 scale = one
282 IF( m.EQ.0 .OR. n.EQ.0 )
283 $ RETURN
284*
285* Use unblocked code for small problems or if insufficient
286* workspaces are provided
287*
288 IF( min( nba, nbb ).EQ.1 .OR. ldswork.LT.max( nba, nbb ) .OR.
289 $ liwork.LT.iwork(1) ) THEN
290 CALL strsyl( trana, tranb, isgn, m, n, a, lda, b, ldb,
291 $ c, ldc, scale, info )
292 RETURN
293 END IF
294*
295* Set constants to control overflow
296*
297 smlnum = slamch( 'S' )
298 bignum = one / smlnum
299*
300* Partition A such that 2-by-2 blocks on the diagonal are not split
301*
302 skip = .false.
303 DO i = 1, nba
304 iwork( i ) = ( i - 1 ) * nb + 1
305 END DO
306 iwork( nba + 1 ) = m + 1
307 DO k = 1, nba
308 l1 = iwork( k )
309 l2 = iwork( k + 1 ) - 1
310 DO l = l1, l2
311 IF( skip ) THEN
312 skip = .false.
313 cycle
314 END IF
315 IF( l.GE.m ) THEN
316* A( M, M ) is a 1-by-1 block
317 cycle
318 END IF
319 IF( a( l, l+1 ).NE.zero .AND. a( l+1, l ).NE.zero ) THEN
320* Check if 2-by-2 block is split
321 IF( l + 1 .EQ. iwork( k + 1 ) ) THEN
322 iwork( k + 1 ) = iwork( k + 1 ) + 1
323 cycle
324 END IF
325 skip = .true.
326 END IF
327 END DO
328 END DO
329 iwork( nba + 1 ) = m + 1
330 IF( iwork( nba ).GE.iwork( nba + 1 ) ) THEN
331 iwork( nba ) = iwork( nba + 1 )
332 nba = nba - 1
333 END IF
334*
335* Partition B such that 2-by-2 blocks on the diagonal are not split
336*
337 pc = nba + 1
338 skip = .false.
339 DO i = 1, nbb
340 iwork( pc + i ) = ( i - 1 ) * nb + 1
341 END DO
342 iwork( pc + nbb + 1 ) = n + 1
343 DO k = 1, nbb
344 l1 = iwork( pc + k )
345 l2 = iwork( pc + k + 1 ) - 1
346 DO l = l1, l2
347 IF( skip ) THEN
348 skip = .false.
349 cycle
350 END IF
351 IF( l.GE.n ) THEN
352* B( N, N ) is a 1-by-1 block
353 cycle
354 END IF
355 IF( b( l, l+1 ).NE.zero .AND. b( l+1, l ).NE.zero ) THEN
356* Check if 2-by-2 block is split
357 IF( l + 1 .EQ. iwork( pc + k + 1 ) ) THEN
358 iwork( pc + k + 1 ) = iwork( pc + k + 1 ) + 1
359 cycle
360 END IF
361 skip = .true.
362 END IF
363 END DO
364 END DO
365 iwork( pc + nbb + 1 ) = n + 1
366 IF( iwork( pc + nbb ).GE.iwork( pc + nbb + 1 ) ) THEN
367 iwork( pc + nbb ) = iwork( pc + nbb + 1 )
368 nbb = nbb - 1
369 END IF
370*
371* Set local scaling factors - must never attain zero.
372*
373 DO l = 1, nbb
374 DO k = 1, nba
375 swork( k, l ) = one
376 END DO
377 END DO
378*
379* Fallback scaling factor to prevent flushing of SWORK( K, L ) to zero.
380* This scaling is to ensure compatibility with TRSYL and may get flushed.
381*
382 buf = one
383*
384* Compute upper bounds of blocks of A and B
385*
386 awrk = nbb
387 DO k = 1, nba
388 k1 = iwork( k )
389 k2 = iwork( k + 1 )
390 DO l = k, nba
391 l1 = iwork( l )
392 l2 = iwork( l + 1 )
393 IF( notrna ) THEN
394 swork( k, awrk + l ) = slange( 'I', k2-k1, l2-l1,
395 $ a( k1, l1 ), lda, wnrm )
396 ELSE
397 swork( l, awrk + k ) = slange( '1', k2-k1, l2-l1,
398 $ a( k1, l1 ), lda, wnrm )
399 END IF
400 END DO
401 END DO
402 bwrk = nbb + nba
403 DO k = 1, nbb
404 k1 = iwork( pc + k )
405 k2 = iwork( pc + k + 1 )
406 DO l = k, nbb
407 l1 = iwork( pc + l )
408 l2 = iwork( pc + l + 1 )
409 IF( notrnb ) THEN
410 swork( k, bwrk + l ) = slange( 'I', k2-k1, l2-l1,
411 $ b( k1, l1 ), ldb, wnrm )
412 ELSE
413 swork( l, bwrk + k ) = slange( '1', k2-k1, l2-l1,
414 $ b( k1, l1 ), ldb, wnrm )
415 END IF
416 END DO
417 END DO
418*
419 sgn = real( isgn )
420*
421 IF( notrna .AND. notrnb ) THEN
422*
423* Solve A*X + ISGN*X*B = scale*C.
424*
425* The (K,L)th block of X is determined starting from
426* bottom-left corner column by column by
427*
428* A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
429*
430* Where
431* M L-1
432* R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(J,L)].
433* I=K+1 J=1
434*
435* Start loop over block rows (index = K) and block columns (index = L)
436*
437 DO k = nba, 1, -1
438*
439* K1: row index of the first row in X( K, L )
440* K2: row index of the first row in X( K+1, L )
441* so the K2 - K1 is the column count of the block X( K, L )
442*
443 k1 = iwork( k )
444 k2 = iwork( k + 1 )
445 DO l = 1, nbb
446*
447* L1: column index of the first column in X( K, L )
448* L2: column index of the first column in X( K, L + 1)
449* so that L2 - L1 is the row count of the block X( K, L )
450*
451 l1 = iwork( pc + l )
452 l2 = iwork( pc + l + 1 )
453*
454 CALL strsyl( trana, tranb, isgn, k2-k1, l2-l1,
455 $ a( k1, k1 ), lda,
456 $ b( l1, l1 ), ldb,
457 $ c( k1, l1 ), ldc, scaloc, iinfo )
458 info = max( info, iinfo )
459*
460 IF ( scaloc * swork( k, l ) .EQ. zero ) THEN
461 IF( scaloc .EQ. zero ) THEN
462* The magnitude of the largest entry of X(K1:K2-1, L1:L2-1)
463* is larger than the product of BIGNUM**2 and cannot be
464* represented in the form (1/SCALE)*X(K1:K2-1, L1:L2-1).
465* Mark the computation as pointless.
466 buf = zero
467 ELSE
468* Use second scaling factor to prevent flushing to zero.
469 buf = buf*2.e0**exponent( scaloc )
470 END IF
471 DO jj = 1, nbb
472 DO ll = 1, nba
473* Bound by BIGNUM to not introduce Inf. The value
474* is irrelevant; corresponding entries of the
475* solution will be flushed in consistency scaling.
476 swork( ll, jj ) = min( bignum,
477 $ swork( ll, jj ) / 2.e0**exponent( scaloc ) )
478 END DO
479 END DO
480 END IF
481 swork( k, l ) = scaloc * swork( k, l )
482 xnrm = slange( 'I', k2-k1, l2-l1, c( k1, l1 ), ldc,
483 $ wnrm )
484*
485 DO i = k - 1, 1, -1
486*
487* C( I, L ) := C( I, L ) - A( I, K ) * C( K, L )
488*
489 i1 = iwork( i )
490 i2 = iwork( i + 1 )
491*
492* Compute scaling factor to survive the linear update
493* simulating consistent scaling.
494*
495 cnrm = slange( 'I', i2-i1, l2-l1, c( i1, l1 ),
496 $ ldc, wnrm )
497 scamin = min( swork( i, l ), swork( k, l ) )
498 cnrm = cnrm * ( scamin / swork( i, l ) )
499 xnrm = xnrm * ( scamin / swork( k, l ) )
500 anrm = swork( i, awrk + k )
501 scaloc = slarmm( anrm, xnrm, cnrm )
502 IF( scaloc * scamin .EQ. zero ) THEN
503* Use second scaling factor to prevent flushing to zero.
504 buf = buf*2.e0**exponent( scaloc )
505 DO jj = 1, nbb
506 DO ll = 1, nba
507 swork( ll, jj ) = min( bignum,
508 $ swork( ll, jj ) / 2.e0**exponent( scaloc ) )
509 END DO
510 END DO
511 scamin = scamin / 2.e0**exponent( scaloc )
512 scaloc = scaloc / 2.e0**exponent( scaloc )
513 END IF
514 cnrm = cnrm * scaloc
515 xnrm = xnrm * scaloc
516*
517* Simultaneously apply the robust update factor and the
518* consistency scaling factor to C( I, L ) and C( K, L ).
519*
520 scal = ( scamin / swork( k, l ) ) * scaloc
521 IF (scal .NE. one) THEN
522 DO jj = l1, l2-1
523 CALL sscal( k2-k1, scal, c( k1, jj ), 1)
524 END DO
525 ENDIF
526*
527 scal = ( scamin / swork( i, l ) ) * scaloc
528 IF (scal .NE. one) THEN
529 DO ll = l1, l2-1
530 CALL sscal( i2-i1, scal, c( i1, ll ), 1)
531 END DO
532 ENDIF
533*
534* Record current scaling factor
535*
536 swork( k, l ) = scamin * scaloc
537 swork( i, l ) = scamin * scaloc
538*
539 CALL sgemm( 'N', 'N', i2-i1, l2-l1, k2-k1, -one,
540 $ a( i1, k1 ), lda, c( k1, l1 ), ldc,
541 $ one, c( i1, l1 ), ldc )
542*
543 END DO
544*
545 DO j = l + 1, nbb
546*
547* C( K, J ) := C( K, J ) - SGN * C( K, L ) * B( L, J )
548*
549 j1 = iwork( pc + j )
550 j2 = iwork( pc + j + 1 )
551*
552* Compute scaling factor to survive the linear update
553* simulating consistent scaling.
554*
555 cnrm = slange( 'I', k2-k1, j2-j1, c( k1, j1 ),
556 $ ldc, wnrm )
557 scamin = min( swork( k, j ), swork( k, l ) )
558 cnrm = cnrm * ( scamin / swork( k, j ) )
559 xnrm = xnrm * ( scamin / swork( k, l ) )
560 bnrm = swork(l, bwrk + j)
561 scaloc = slarmm( bnrm, xnrm, cnrm )
562 IF( scaloc * scamin .EQ. zero ) THEN
563* Use second scaling factor to prevent flushing to zero.
564 buf = buf*2.e0**exponent( scaloc )
565 DO jj = 1, nbb
566 DO ll = 1, nba
567 swork( ll, jj ) = min( bignum,
568 $ swork( ll, jj ) / 2.e0**exponent( scaloc ) )
569 END DO
570 END DO
571 scamin = scamin / 2.e0**exponent( scaloc )
572 scaloc = scaloc / 2.e0**exponent( scaloc )
573 END IF
574 cnrm = cnrm * scaloc
575 xnrm = xnrm * scaloc
576*
577* Simultaneously apply the robust update factor and the
578* consistency scaling factor to C( K, J ) and C( K, L).
579*
580 scal = ( scamin / swork( k, l ) ) * scaloc
581 IF( scal .NE. one ) THEN
582 DO ll = l1, l2-1
583 CALL sscal( k2-k1, scal, c( k1, ll ), 1 )
584 END DO
585 ENDIF
586*
587 scal = ( scamin / swork( k, j ) ) * scaloc
588 IF( scal .NE. one ) THEN
589 DO jj = j1, j2-1
590 CALL sscal( k2-k1, scal, c( k1, jj ), 1 )
591 END DO
592 ENDIF
593*
594* Record current scaling factor
595*
596 swork( k, l ) = scamin * scaloc
597 swork( k, j ) = scamin * scaloc
598*
599 CALL sgemm( 'N', 'N', k2-k1, j2-j1, l2-l1, -sgn,
600 $ c( k1, l1 ), ldc, b( l1, j1 ), ldb,
601 $ one, c( k1, j1 ), ldc )
602 END DO
603 END DO
604 END DO
605 ELSE IF( .NOT.notrna .AND. notrnb ) THEN
606*
607* Solve A**T*X + ISGN*X*B = scale*C.
608*
609* The (K,L)th block of X is determined starting from
610* upper-left corner column by column by
611*
612* A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
613*
614* Where
615* K-1 L-1
616* R(K,L) = SUM [A(I,K)**T*X(I,L)] +ISGN*SUM [X(K,J)*B(J,L)]
617* I=1 J=1
618*
619* Start loop over block rows (index = K) and block columns (index = L)
620*
621 DO k = 1, nba
622*
623* K1: row index of the first row in X( K, L )
624* K2: row index of the first row in X( K+1, L )
625* so the K2 - K1 is the column count of the block X( K, L )
626*
627 k1 = iwork( k )
628 k2 = iwork( k + 1 )
629 DO l = 1, nbb
630*
631* L1: column index of the first column in X( K, L )
632* L2: column index of the first column in X( K, L + 1)
633* so that L2 - L1 is the row count of the block X( K, L )
634*
635 l1 = iwork( pc + l )
636 l2 = iwork( pc + l + 1 )
637*
638 CALL strsyl( trana, tranb, isgn, k2-k1, l2-l1,
639 $ a( k1, k1 ), lda,
640 $ b( l1, l1 ), ldb,
641 $ c( k1, l1 ), ldc, scaloc, iinfo )
642 info = max( info, iinfo )
643*
644 IF( scaloc * swork( k, l ) .EQ. zero ) THEN
645 IF( scaloc .EQ. zero ) THEN
646* The magnitude of the largest entry of X(K1:K2-1, L1:L2-1)
647* is larger than the product of BIGNUM**2 and cannot be
648* represented in the form (1/SCALE)*X(K1:K2-1, L1:L2-1).
649* Mark the computation as pointless.
650 buf = zero
651 ELSE
652* Use second scaling factor to prevent flushing to zero.
653 buf = buf*2.e0**exponent( scaloc )
654 END IF
655 DO jj = 1, nbb
656 DO ll = 1, nba
657* Bound by BIGNUM to not introduce Inf. The value
658* is irrelevant; corresponding entries of the
659* solution will be flushed in consistency scaling.
660 swork( ll, jj ) = min( bignum,
661 $ swork( ll, jj ) / 2.e0**exponent( scaloc ) )
662 END DO
663 END DO
664 END IF
665 swork( k, l ) = scaloc * swork( k, l )
666 xnrm = slange( 'I', k2-k1, l2-l1, c( k1, l1 ), ldc,
667 $ wnrm )
668*
669 DO i = k + 1, nba
670*
671* C( I, L ) := C( I, L ) - A( K, I )**T * C( K, L )
672*
673 i1 = iwork( i )
674 i2 = iwork( i + 1 )
675*
676* Compute scaling factor to survive the linear update
677* simulating consistent scaling.
678*
679 cnrm = slange( 'I', i2-i1, l2-l1, c( i1, l1 ),
680 $ ldc, wnrm )
681 scamin = min( swork( i, l ), swork( k, l ) )
682 cnrm = cnrm * ( scamin / swork( i, l ) )
683 xnrm = xnrm * ( scamin / swork( k, l ) )
684 anrm = swork( i, awrk + k )
685 scaloc = slarmm( anrm, xnrm, cnrm )
686 IF( scaloc * scamin .EQ. zero ) THEN
687* Use second scaling factor to prevent flushing to zero.
688 buf = buf*2.e0**exponent( scaloc )
689 DO jj = 1, nbb
690 DO ll = 1, nba
691 swork( ll, jj ) = min( bignum,
692 $ swork( ll, jj ) / 2.e0**exponent( scaloc ) )
693 END DO
694 END DO
695 scamin = scamin / 2.e0**exponent( scaloc )
696 scaloc = scaloc / 2.e0**exponent( scaloc )
697 END IF
698 cnrm = cnrm * scaloc
699 xnrm = xnrm * scaloc
700*
701* Simultaneously apply the robust update factor and the
702* consistency scaling factor to to C( I, L ) and C( K, L ).
703*
704 scal = ( scamin / swork( k, l ) ) * scaloc
705 IF (scal .NE. one) THEN
706 DO ll = l1, l2-1
707 CALL sscal( k2-k1, scal, c( k1, ll ), 1 )
708 END DO
709 ENDIF
710*
711 scal = ( scamin / swork( i, l ) ) * scaloc
712 IF (scal .NE. one) THEN
713 DO ll = l1, l2-1
714 CALL sscal( i2-i1, scal, c( i1, ll ), 1 )
715 END DO
716 ENDIF
717*
718* Record current scaling factor
719*
720 swork( k, l ) = scamin * scaloc
721 swork( i, l ) = scamin * scaloc
722*
723 CALL sgemm( 'T', 'N', i2-i1, l2-l1, k2-k1, -one,
724 $ a( k1, i1 ), lda, c( k1, l1 ), ldc,
725 $ one, c( i1, l1 ), ldc )
726 END DO
727*
728 DO j = l + 1, nbb
729*
730* C( K, J ) := C( K, J ) - SGN * C( K, L ) * B( L, J )
731*
732 j1 = iwork( pc + j )
733 j2 = iwork( pc + j + 1 )
734*
735* Compute scaling factor to survive the linear update
736* simulating consistent scaling.
737*
738 cnrm = slange( 'I', k2-k1, j2-j1, c( k1, j1 ),
739 $ ldc, wnrm )
740 scamin = min( swork( k, j ), swork( k, l ) )
741 cnrm = cnrm * ( scamin / swork( k, j ) )
742 xnrm = xnrm * ( scamin / swork( k, l ) )
743 bnrm = swork( l, bwrk + j )
744 scaloc = slarmm( bnrm, xnrm, cnrm )
745 IF( scaloc * scamin .EQ. zero ) THEN
746* Use second scaling factor to prevent flushing to zero.
747 buf = buf*2.e0**exponent( scaloc )
748 DO jj = 1, nbb
749 DO ll = 1, nba
750 swork( ll, jj ) = min( bignum,
751 $ swork( ll, jj ) / 2.e0**exponent( scaloc ) )
752 END DO
753 END DO
754 scamin = scamin / 2.e0**exponent( scaloc )
755 scaloc = scaloc / 2.e0**exponent( scaloc )
756 END IF
757 cnrm = cnrm * scaloc
758 xnrm = xnrm * scaloc
759*
760* Simultaneously apply the robust update factor and the
761* consistency scaling factor to to C( K, J ) and C( K, L ).
762*
763 scal = ( scamin / swork( k, l ) ) * scaloc
764 IF( scal .NE. one ) THEN
765 DO ll = l1, l2-1
766 CALL sscal( k2-k1, scal, c( k1, ll ), 1 )
767 END DO
768 ENDIF
769*
770 scal = ( scamin / swork( k, j ) ) * scaloc
771 IF( scal .NE. one ) THEN
772 DO jj = j1, j2-1
773 CALL sscal( k2-k1, scal, c( k1, jj ), 1 )
774 END DO
775 ENDIF
776*
777* Record current scaling factor
778*
779 swork( k, l ) = scamin * scaloc
780 swork( k, j ) = scamin * scaloc
781*
782 CALL sgemm( 'N', 'N', k2-k1, j2-j1, l2-l1, -sgn,
783 $ c( k1, l1 ), ldc, b( l1, j1 ), ldb,
784 $ one, c( k1, j1 ), ldc )
785 END DO
786 END DO
787 END DO
788 ELSE IF( .NOT.notrna .AND. .NOT.notrnb ) THEN
789*
790* Solve A**T*X + ISGN*X*B**T = scale*C.
791*
792* The (K,L)th block of X is determined starting from
793* top-right corner column by column by
794*
795* A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L)**T = C(K,L) - R(K,L)
796*
797* Where
798* K-1 N
799* R(K,L) = SUM [A(I,K)**T*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**T].
800* I=1 J=L+1
801*
802* Start loop over block rows (index = K) and block columns (index = L)
803*
804 DO k = 1, nba
805*
806* K1: row index of the first row in X( K, L )
807* K2: row index of the first row in X( K+1, L )
808* so the K2 - K1 is the column count of the block X( K, L )
809*
810 k1 = iwork( k )
811 k2 = iwork( k + 1 )
812 DO l = nbb, 1, -1
813*
814* L1: column index of the first column in X( K, L )
815* L2: column index of the first column in X( K, L + 1)
816* so that L2 - L1 is the row count of the block X( K, L )
817*
818 l1 = iwork( pc + l )
819 l2 = iwork( pc + l + 1 )
820*
821 CALL strsyl( trana, tranb, isgn, k2-k1, l2-l1,
822 $ a( k1, k1 ), lda,
823 $ b( l1, l1 ), ldb,
824 $ c( k1, l1 ), ldc, scaloc, iinfo )
825 info = max( info, iinfo )
826*
827 IF( scaloc * swork( k, l ) .EQ. zero ) THEN
828 IF( scaloc .EQ. zero ) THEN
829* The magnitude of the largest entry of X(K1:K2-1, L1:L2-1)
830* is larger than the product of BIGNUM**2 and cannot be
831* represented in the form (1/SCALE)*X(K1:K2-1, L1:L2-1).
832* Mark the computation as pointless.
833 buf = zero
834 ELSE
835* Use second scaling factor to prevent flushing to zero.
836 buf = buf*2.e0**exponent( scaloc )
837 END IF
838 DO jj = 1, nbb
839 DO ll = 1, nba
840* Bound by BIGNUM to not introduce Inf. The value
841* is irrelevant; corresponding entries of the
842* solution will be flushed in consistency scaling.
843 swork( ll, jj ) = min( bignum,
844 $ swork( ll, jj ) / 2.e0**exponent( scaloc ) )
845 END DO
846 END DO
847 END IF
848 swork( k, l ) = scaloc * swork( k, l )
849 xnrm = slange( 'I', k2-k1, l2-l1, c( k1, l1 ), ldc,
850 $ wnrm )
851*
852 DO i = k + 1, nba
853*
854* C( I, L ) := C( I, L ) - A( K, I )**T * C( K, L )
855*
856 i1 = iwork( i )
857 i2 = iwork( i + 1 )
858*
859* Compute scaling factor to survive the linear update
860* simulating consistent scaling.
861*
862 cnrm = slange( 'I', i2-i1, l2-l1, c( i1, l1 ),
863 $ ldc, wnrm )
864 scamin = min( swork( i, l ), swork( k, l ) )
865 cnrm = cnrm * ( scamin / swork( i, l ) )
866 xnrm = xnrm * ( scamin / swork( k, l ) )
867 anrm = swork( i, awrk + k )
868 scaloc = slarmm( anrm, xnrm, cnrm )
869 IF( scaloc * scamin .EQ. zero ) THEN
870* Use second scaling factor to prevent flushing to zero.
871 buf = buf*2.e0**exponent( scaloc )
872 DO jj = 1, nbb
873 DO ll = 1, nba
874 swork( ll, jj ) = min( bignum,
875 $ swork( ll, jj ) / 2.e0**exponent( scaloc ) )
876 END DO
877 END DO
878 scamin = scamin / 2.e0**exponent( scaloc )
879 scaloc = scaloc / 2.e0**exponent( scaloc )
880 END IF
881 cnrm = cnrm * scaloc
882 xnrm = xnrm * scaloc
883*
884* Simultaneously apply the robust update factor and the
885* consistency scaling factor to C( I, L ) and C( K, L ).
886*
887 scal = ( scamin / swork( k, l ) ) * scaloc
888 IF (scal .NE. one) THEN
889 DO ll = l1, l2-1
890 CALL sscal( k2-k1, scal, c( k1, ll ), 1 )
891 END DO
892 ENDIF
893*
894 scal = ( scamin / swork( i, l ) ) * scaloc
895 IF (scal .NE. one) THEN
896 DO ll = l1, l2-1
897 CALL sscal( i2-i1, scal, c( i1, ll ), 1 )
898 END DO
899 ENDIF
900*
901* Record current scaling factor
902*
903 swork( k, l ) = scamin * scaloc
904 swork( i, l ) = scamin * scaloc
905*
906 CALL sgemm( 'T', 'N', i2-i1, l2-l1, k2-k1, -one,
907 $ a( k1, i1 ), lda, c( k1, l1 ), ldc,
908 $ one, c( i1, l1 ), ldc )
909 END DO
910*
911 DO j = 1, l - 1
912*
913* C( K, J ) := C( K, J ) - SGN * C( K, L ) * B( J, L )**T
914*
915 j1 = iwork( pc + j )
916 j2 = iwork( pc + j + 1 )
917*
918* Compute scaling factor to survive the linear update
919* simulating consistent scaling.
920*
921 cnrm = slange( 'I', k2-k1, j2-j1, c( k1, j1 ),
922 $ ldc, wnrm )
923 scamin = min( swork( k, j ), swork( k, l ) )
924 cnrm = cnrm * ( scamin / swork( k, j ) )
925 xnrm = xnrm * ( scamin / swork( k, l ) )
926 bnrm = swork( l, bwrk + j )
927 scaloc = slarmm( bnrm, xnrm, cnrm )
928 IF( scaloc * scamin .EQ. zero ) THEN
929* Use second scaling factor to prevent flushing to zero.
930 buf = buf*2.e0**exponent( scaloc )
931 DO jj = 1, nbb
932 DO ll = 1, nba
933 swork( ll, jj ) = min( bignum,
934 $ swork( ll, jj ) / 2.e0**exponent( scaloc ) )
935 END DO
936 END DO
937 scamin = scamin / 2.e0**exponent( scaloc )
938 scaloc = scaloc / 2.e0**exponent( scaloc )
939 END IF
940 cnrm = cnrm * scaloc
941 xnrm = xnrm * scaloc
942*
943* Simultaneously apply the robust update factor and the
944* consistency scaling factor to C( K, J ) and C( K, L ).
945*
946 scal = ( scamin / swork( k, l ) ) * scaloc
947 IF( scal .NE. one ) THEN
948 DO ll = l1, l2-1
949 CALL sscal( k2-k1, scal, c( k1, ll ), 1)
950 END DO
951 ENDIF
952*
953 scal = ( scamin / swork( k, j ) ) * scaloc
954 IF( scal .NE. one ) THEN
955 DO jj = j1, j2-1
956 CALL sscal( k2-k1, scal, c( k1, jj ), 1 )
957 END DO
958 ENDIF
959*
960* Record current scaling factor
961*
962 swork( k, l ) = scamin * scaloc
963 swork( k, j ) = scamin * scaloc
964*
965 CALL sgemm( 'N', 'T', k2-k1, j2-j1, l2-l1, -sgn,
966 $ c( k1, l1 ), ldc, b( j1, l1 ), ldb,
967 $ one, c( k1, j1 ), ldc )
968 END DO
969 END DO
970 END DO
971 ELSE IF( notrna .AND. .NOT.notrnb ) THEN
972*
973* Solve A*X + ISGN*X*B**T = scale*C.
974*
975* The (K,L)th block of X is determined starting from
976* bottom-right corner column by column by
977*
978* A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L)**T = C(K,L) - R(K,L)
979*
980* Where
981* M N
982* R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**T].
983* I=K+1 J=L+1
984*
985* Start loop over block rows (index = K) and block columns (index = L)
986*
987 DO k = nba, 1, -1
988*
989* K1: row index of the first row in X( K, L )
990* K2: row index of the first row in X( K+1, L )
991* so the K2 - K1 is the column count of the block X( K, L )
992*
993 k1 = iwork( k )
994 k2 = iwork( k + 1 )
995 DO l = nbb, 1, -1
996*
997* L1: column index of the first column in X( K, L )
998* L2: column index of the first column in X( K, L + 1)
999* so that L2 - L1 is the row count of the block X( K, L )
1000*
1001 l1 = iwork( pc + l )
1002 l2 = iwork( pc + l + 1 )
1003*
1004 CALL strsyl( trana, tranb, isgn, k2-k1, l2-l1,
1005 $ a( k1, k1 ), lda,
1006 $ b( l1, l1 ), ldb,
1007 $ c( k1, l1 ), ldc, scaloc, iinfo )
1008 info = max( info, iinfo )
1009*
1010 IF( scaloc * swork( k, l ) .EQ. zero ) THEN
1011 IF( scaloc .EQ. zero ) THEN
1012* The magnitude of the largest entry of X(K1:K2-1, L1:L2-1)
1013* is larger than the product of BIGNUM**2 and cannot be
1014* represented in the form (1/SCALE)*X(K1:K2-1, L1:L2-1).
1015* Mark the computation as pointless.
1016 buf = zero
1017 ELSE
1018* Use second scaling factor to prevent flushing to zero.
1019 buf = buf*2.e0**exponent( scaloc )
1020 END IF
1021 DO jj = 1, nbb
1022 DO ll = 1, nba
1023* Bound by BIGNUM to not introduce Inf. The value
1024* is irrelevant; corresponding entries of the
1025* solution will be flushed in consistency scaling.
1026 swork( ll, jj ) = min( bignum,
1027 $ swork( ll, jj ) / 2.e0**exponent( scaloc ) )
1028 END DO
1029 END DO
1030 END IF
1031 swork( k, l ) = scaloc * swork( k, l )
1032 xnrm = slange( 'I', k2-k1, l2-l1, c( k1, l1 ), ldc,
1033 $ wnrm )
1034*
1035 DO i = 1, k - 1
1036*
1037* C( I, L ) := C( I, L ) - A( I, K ) * C( K, L )
1038*
1039 i1 = iwork( i )
1040 i2 = iwork( i + 1 )
1041*
1042* Compute scaling factor to survive the linear update
1043* simulating consistent scaling.
1044*
1045 cnrm = slange( 'I', i2-i1, l2-l1, c( i1, l1 ),
1046 $ ldc, wnrm )
1047 scamin = min( swork( i, l ), swork( k, l ) )
1048 cnrm = cnrm * ( scamin / swork( i, l ) )
1049 xnrm = xnrm * ( scamin / swork( k, l ) )
1050 anrm = swork( i, awrk + k )
1051 scaloc = slarmm( anrm, xnrm, cnrm )
1052 IF( scaloc * scamin .EQ. zero ) THEN
1053* Use second scaling factor to prevent flushing to zero.
1054 buf = buf*2.e0**exponent( scaloc )
1055 DO jj = 1, nbb
1056 DO ll = 1, nba
1057 swork( ll, jj ) = min( bignum,
1058 $ swork( ll, jj ) / 2.e0**exponent( scaloc ) )
1059 END DO
1060 END DO
1061 scamin = scamin / 2.e0**exponent( scaloc )
1062 scaloc = scaloc / 2.e0**exponent( scaloc )
1063 END IF
1064 cnrm = cnrm * scaloc
1065 xnrm = xnrm * scaloc
1066*
1067* Simultaneously apply the robust update factor and the
1068* consistency scaling factor to C( I, L ) and C( K, L ).
1069*
1070 scal = ( scamin / swork( k, l ) ) * scaloc
1071 IF (scal .NE. one) THEN
1072 DO ll = l1, l2-1
1073 CALL sscal( k2-k1, scal, c( k1, ll ), 1 )
1074 END DO
1075 ENDIF
1076*
1077 scal = ( scamin / swork( i, l ) ) * scaloc
1078 IF (scal .NE. one) THEN
1079 DO ll = l1, l2-1
1080 CALL sscal( i2-i1, scal, c( i1, ll ), 1 )
1081 END DO
1082 ENDIF
1083*
1084* Record current scaling factor
1085*
1086 swork( k, l ) = scamin * scaloc
1087 swork( i, l ) = scamin * scaloc
1088*
1089 CALL sgemm( 'N', 'N', i2-i1, l2-l1, k2-k1, -one,
1090 $ a( i1, k1 ), lda, c( k1, l1 ), ldc,
1091 $ one, c( i1, l1 ), ldc )
1092*
1093 END DO
1094*
1095 DO j = 1, l - 1
1096*
1097* C( K, J ) := C( K, J ) - SGN * C( K, L ) * B( J, L )**T
1098*
1099 j1 = iwork( pc + j )
1100 j2 = iwork( pc + j + 1 )
1101*
1102* Compute scaling factor to survive the linear update
1103* simulating consistent scaling.
1104*
1105 cnrm = slange( 'I', k2-k1, j2-j1, c( k1, j1 ),
1106 $ ldc, wnrm )
1107 scamin = min( swork( k, j ), swork( k, l ) )
1108 cnrm = cnrm * ( scamin / swork( k, j ) )
1109 xnrm = xnrm * ( scamin / swork( k, l ) )
1110 bnrm = swork( l, bwrk + j )
1111 scaloc = slarmm( bnrm, xnrm, cnrm )
1112 IF( scaloc * scamin .EQ. zero ) THEN
1113* Use second scaling factor to prevent flushing to zero.
1114 buf = buf*2.e0**exponent( scaloc )
1115 DO jj = 1, nbb
1116 DO ll = 1, nba
1117 swork( ll, jj ) = min( bignum,
1118 $ swork( ll, jj ) / 2.e0**exponent( scaloc ) )
1119 END DO
1120 END DO
1121 scamin = scamin / 2.e0**exponent( scaloc )
1122 scaloc = scaloc / 2.e0**exponent( scaloc )
1123 END IF
1124 cnrm = cnrm * scaloc
1125 xnrm = xnrm * scaloc
1126*
1127* Simultaneously apply the robust update factor and the
1128* consistency scaling factor to C( K, J ) and C( K, L ).
1129*
1130 scal = ( scamin / swork( k, l ) ) * scaloc
1131 IF( scal .NE. one ) THEN
1132 DO jj = l1, l2-1
1133 CALL sscal( k2-k1, scal, c( k1, jj ), 1 )
1134 END DO
1135 ENDIF
1136*
1137 scal = ( scamin / swork( k, j ) ) * scaloc
1138 IF( scal .NE. one ) THEN
1139 DO jj = j1, j2-1
1140 CALL sscal( k2-k1, scal, c( k1, jj ), 1 )
1141 END DO
1142 ENDIF
1143*
1144* Record current scaling factor
1145*
1146 swork( k, l ) = scamin * scaloc
1147 swork( k, j ) = scamin * scaloc
1148*
1149 CALL sgemm( 'N', 'T', k2-k1, j2-j1, l2-l1, -sgn,
1150 $ c( k1, l1 ), ldc, b( j1, l1 ), ldb,
1151 $ one, c( k1, j1 ), ldc )
1152 END DO
1153 END DO
1154 END DO
1155*
1156 END IF
1157*
1158* Reduce local scaling factors
1159*
1160 scale = swork( 1, 1 )
1161 DO k = 1, nba
1162 DO l = 1, nbb
1163 scale = min( scale, swork( k, l ) )
1164 END DO
1165 END DO
1166*
1167 IF( scale .EQ. zero ) THEN
1168*
1169* The magnitude of the largest entry of the solution is larger
1170* than the product of BIGNUM**2 and cannot be represented in the
1171* form (1/SCALE)*X if SCALE is REAL. Set SCALE to zero and give up.
1172*
1173 iwork(1) = nba + nbb + 2
1174 swork(1,1) = max( nba, nbb )
1175 swork(2,1) = 2 * nbb + nba
1176 RETURN
1177 END IF
1178*
1179* Realize consistent scaling
1180*
1181 DO k = 1, nba
1182 k1 = iwork( k )
1183 k2 = iwork( k + 1 )
1184 DO l = 1, nbb
1185 l1 = iwork( pc + l )
1186 l2 = iwork( pc + l + 1 )
1187 scal = scale / swork( k, l )
1188 IF( scal .NE. one ) THEN
1189 DO ll = l1, l2-1
1190 CALL sscal( k2-k1, scal, c( k1, ll ), 1 )
1191 END DO
1192 ENDIF
1193 END DO
1194 END DO
1195*
1196 IF( buf .NE. one .AND. buf.GT.zero ) THEN
1197*
1198* Decrease SCALE as much as possible.
1199*
1200 scaloc = min( scale / smlnum, one / buf )
1201 buf = buf * scaloc
1202 scale = scale / scaloc
1203 END IF
1204
1205 IF( buf.NE.one .AND. buf.GT.zero ) THEN
1206*
1207* In case of overly aggressive scaling during the computation,
1208* flushing of the global scale factor may be prevented by
1209* undoing some of the scaling. This step is to ensure that
1210* this routine flushes only scale factors that TRSYL also
1211* flushes and be usable as a drop-in replacement.
1212*
1213* How much can the normwise largest entry be upscaled?
1214*
1215 scal = c( 1, 1 )
1216 DO k = 1, m
1217 DO l = 1, n
1218 scal = max( scal, abs( c( k, l ) ) )
1219 END DO
1220 END DO
1221*
1222* Increase BUF as close to 1 as possible and apply scaling.
1223*
1224 scaloc = min( bignum / scal, one / buf )
1225 buf = buf * scaloc
1226 CALL slascl( 'G', -1, -1, one, scaloc, m, n, c, ldc, iwork )
1227 END IF
1228*
1229* Combine with buffer scaling factor. SCALE will be flushed if
1230* BUF is less than one here.
1231*
1232 scale = scale * buf
1233*
1234* Restore workspace dimensions
1235*
1236 iwork(1) = nba + nbb + 2
1237 swork(1,1) = max( nba, nbb )
1238 swork(2,1) = 2 * nbb + nba
1239*
1240 RETURN
1241*
1242* End of STRSYL3
1243*
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: slascl.f:143
real function slarmm(ANORM, BNORM, CNORM)
SLARMM
Definition: slarmm.f:61
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
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 strsyl(TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C, LDC, SCALE, INFO)
STRSYL
Definition: strsyl.f:164
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187
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: