LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ dtrsyl3()

 subroutine dtrsyl3 ( character trana, character tranb, integer isgn, integer m, integer n, double precision, dimension( lda, * ) a, integer lda, double precision, dimension( ldb, * ) b, integer ldb, double precision, dimension( ldc, * ) c, integer ldc, double precision scale, integer, dimension( * ) iwork, integer liwork, double precision, dimension( ldswork, * ) swork, integer ldswork, integer info )

DTRSYL3

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