LAPACK  3.11.0
LAPACK: Linear Algebra PACKage
dtrsyl3.f
1 *> \brief \b DTRSYL3
2 *
3 * Definition:
4 * ===========
5 *
6 *
7 *> \par Purpose
8 * =============
9 *>
10 *> \verbatim
11 *>
12 *> DTRSYL3 solves the real Sylvester matrix equation:
13 *>
14 *> op(A)*X + X*op(B) = scale*C or
15 *> op(A)*X - X*op(B) = scale*C,
16 *>
17 *> where op(A) = A or A**T, and A and B are both upper quasi-
18 *> triangular. A is M-by-M and B is N-by-N; the right hand side C and
19 *> the solution X are M-by-N; and scale is an output scale factor, set
20 *> <= 1 to avoid overflow in X.
21 *>
22 *> A and B must be in Schur canonical form (as returned by DHSEQR), that
23 *> is, block upper triangular with 1-by-1 and 2-by-2 diagonal blocks;
24 *> each 2-by-2 diagonal block has its diagonal elements equal and its
25 *> off-diagonal elements of opposite sign.
26 *>
27 *> This is the block version of the algorithm.
28 *> \endverbatim
29 *
30 * Arguments
31 * =========
32 *
33 *> \param[in] TRANA
34 *> \verbatim
35 *> TRANA is CHARACTER*1
36 *> Specifies the option op(A):
37 *> = 'N': op(A) = A (No transpose)
38 *> = 'T': op(A) = A**T (Transpose)
39 *> = 'C': op(A) = A**H (Conjugate transpose = Transpose)
40 *> \endverbatim
41 *>
42 *> \param[in] TRANB
43 *> \verbatim
44 *> TRANB is CHARACTER*1
45 *> Specifies the option op(B):
46 *> = 'N': op(B) = B (No transpose)
47 *> = 'T': op(B) = B**T (Transpose)
48 *> = 'C': op(B) = B**H (Conjugate transpose = Transpose)
49 *> \endverbatim
50 *>
51 *> \param[in] ISGN
52 *> \verbatim
53 *> ISGN is INTEGER
54 *> Specifies the sign in the equation:
55 *> = +1: solve op(A)*X + X*op(B) = scale*C
56 *> = -1: solve op(A)*X - X*op(B) = scale*C
57 *> \endverbatim
58 *>
59 *> \param[in] M
60 *> \verbatim
61 *> M is INTEGER
62 *> The order of the matrix A, and the number of rows in the
63 *> matrices X and C. M >= 0.
64 *> \endverbatim
65 *>
66 *> \param[in] N
67 *> \verbatim
68 *> N is INTEGER
69 *> The order of the matrix B, and the number of columns in the
70 *> matrices X and C. N >= 0.
71 *> \endverbatim
72 *>
73 *> \param[in] A
74 *> \verbatim
75 *> A is DOUBLE PRECISION array, dimension (LDA,M)
76 *> The upper quasi-triangular matrix A, in Schur canonical form.
77 *> \endverbatim
78 *>
79 *> \param[in] LDA
80 *> \verbatim
81 *> LDA is INTEGER
82 *> The leading dimension of the array A. LDA >= max(1,M).
83 *> \endverbatim
84 *>
85 *> \param[in] B
86 *> \verbatim
87 *> B is DOUBLE PRECISION array, dimension (LDB,N)
88 *> The upper quasi-triangular matrix B, in Schur canonical form.
89 *> \endverbatim
90 *>
91 *> \param[in] LDB
92 *> \verbatim
93 *> LDB is INTEGER
94 *> The leading dimension of the array B. LDB >= max(1,N).
95 *> \endverbatim
96 *>
97 *> \param[in,out] C
98 *> \verbatim
99 *> C is DOUBLE PRECISION array, dimension (LDC,N)
100 *> On entry, the M-by-N right hand side matrix C.
101 *> On exit, C is overwritten by the solution matrix X.
102 *> \endverbatim
103 *>
104 *> \param[in] LDC
105 *> \verbatim
106 *> LDC is INTEGER
107 *> The leading dimension of the array C. LDC >= max(1,M)
108 *> \endverbatim
109 *>
110 *> \param[out] SCALE
111 *> \verbatim
112 *> SCALE is DOUBLE PRECISION
113 *> The scale factor, scale, set <= 1 to avoid overflow in X.
114 *> \endverbatim
115 *>
116 *> \param[out] IWORK
117 *> \verbatim
118 *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
119 *> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
120 *> \endverbatim
121 *>
122 *> \param[in] LIWORK
123 *> \verbatim
124 *> IWORK is INTEGER
125 *> The dimension of the array IWORK. LIWORK >= ((M + NB - 1) / NB + 1)
126 *> + ((N + NB - 1) / NB + 1), where NB is the optimal block size.
127 *>
128 *> If LIWORK = -1, then a workspace query is assumed; the routine
129 *> only calculates the optimal dimension of the IWORK array,
130 *> returns this value as the first entry of the IWORK array, and
131 *> no error message related to LIWORK is issued by XERBLA.
132 *> \endverbatim
133 *>
134 *> \param[out] SWORK
135 *> \verbatim
136 *> SWORK is DOUBLE PRECISION array, dimension (MAX(2, ROWS),
137 *> MAX(1,COLS)).
138 *> On exit, if INFO = 0, SWORK(1) returns the optimal value ROWS
139 *> and SWORK(2) returns the optimal COLS.
140 *> \endverbatim
141 *>
142 *> \param[in] LDSWORK
143 *> \verbatim
144 *> LDSWORK is INTEGER
145 *> LDSWORK >= MAX(2,ROWS), where ROWS = ((M + NB - 1) / NB + 1)
146 *> and NB is the optimal block size.
147 *>
148 *> If LDSWORK = -1, then a workspace query is assumed; the routine
149 *> only calculates the optimal dimensions of the SWORK matrix,
150 *> returns these values as the first and second entry of the SWORK
151 *> matrix, and no error message related LWORK is issued by XERBLA.
152 *> \endverbatim
153 *>
154 *> \param[out] INFO
155 *> \verbatim
156 *> INFO is INTEGER
157 *> = 0: successful exit
158 *> < 0: if INFO = -i, the i-th argument had an illegal value
159 *> = 1: A and B have common or very close eigenvalues; perturbed
160 *> values were used to solve the equation (but the matrices
161 *> A and B are unchanged).
162 *> \endverbatim
163 *
164 * =====================================================================
165 * References:
166 * E. S. Quintana-Orti and R. A. Van De Geijn (2003). Formal derivation of
167 * algorithms: The triangular Sylvester equation, ACM Transactions
168 * on Mathematical Software (TOMS), volume 29, pages 218--243.
169 *
170 * A. Schwarz and C. C. Kjelgaard Mikkelsen (2020). Robust Task-Parallel
171 * Solution of the Triangular Sylvester Equation. Lecture Notes in
172 * Computer Science, vol 12043, pages 82--92, Springer.
173 *
174 * Contributor:
175 * Angelika Schwarz, Umea University, Sweden.
176 *
177 * =====================================================================
178  SUBROUTINE dtrsyl3( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C,
179  $ LDC, SCALE, IWORK, LIWORK, SWORK, LDSWORK,
180  $ INFO )
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 )
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 *
1241  END
ilaenv
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
dlarmm
double precision function dlarmm(ANORM, BNORM, CNORM)
DLARMM
Definition: dlarmm.f:61
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
dgemm
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
dscal
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
dtrsyl
subroutine dtrsyl(TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C, LDC, SCALE, INFO)
DTRSYL
Definition: dtrsyl.f:164
dlange
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
dlascl
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