LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dtrsyl3.f
Go to the documentation of this file.
1*> \brief \b DTRSYL3
2*
3* Definition:
4* ===========
5*
6* SUBROUTINE DTRSYL3( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB,
7* C,
8* LDC, SCALE, IWORK, LIWORK, SWORK, LDSWORK,
9* INFO )
10*
11* .. Scalar Arguments ..
12* CHARACTER TRANA, TRANB
13* INTEGER INFO, ISGN, LDA, LDB, LDC, M, N,
14* LIWORK, LDSWORK
15* DOUBLE PRECISION SCALE
16* ..
17* .. Array Arguments ..
18* INTEGER IWORK( * )
19* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ),
20* SWORK( LDSWORK, * )
21* ..
22*
23*
24*> \par Purpose:
25* =============
26*>
27*> \verbatim
28*>
29*> DTRSYL3 solves the real Sylvester matrix equation:
30*>
31*> op(A)*X + X*op(B) = scale*C or
32*> op(A)*X - X*op(B) = scale*C,
33*>
34*> where op(A) = A or A**T, and A and B are both upper quasi-
35*> triangular. A is M-by-M and B is N-by-N; the right hand side C and
36*> the solution X are M-by-N; and scale is an output scale factor, set
37*> <= 1 to avoid overflow in X.
38*>
39*> A and B must be in Schur canonical form (as returned by DHSEQR), that
40*> is, block upper triangular with 1-by-1 and 2-by-2 diagonal blocks;
41*> each 2-by-2 diagonal block has its diagonal elements equal and its
42*> off-diagonal elements of opposite sign.
43*>
44*> This is the block version of the algorithm.
45*> \endverbatim
46*
47* Arguments:
48* ==========
49*
50*> \param[in] TRANA
51*> \verbatim
52*> TRANA is CHARACTER*1
53*> Specifies the option op(A):
54*> = 'N': op(A) = A (No transpose)
55*> = 'T': op(A) = A**T (Transpose)
56*> = 'C': op(A) = A**H (Conjugate transpose = Transpose)
57*> \endverbatim
58*>
59*> \param[in] TRANB
60*> \verbatim
61*> TRANB is CHARACTER*1
62*> Specifies the option op(B):
63*> = 'N': op(B) = B (No transpose)
64*> = 'T': op(B) = B**T (Transpose)
65*> = 'C': op(B) = B**H (Conjugate transpose = Transpose)
66*> \endverbatim
67*>
68*> \param[in] ISGN
69*> \verbatim
70*> ISGN is INTEGER
71*> Specifies the sign in the equation:
72*> = +1: solve op(A)*X + X*op(B) = scale*C
73*> = -1: solve op(A)*X - X*op(B) = scale*C
74*> \endverbatim
75*>
76*> \param[in] M
77*> \verbatim
78*> M is INTEGER
79*> The order of the matrix A, and the number of rows in the
80*> matrices X and C. M >= 0.
81*> \endverbatim
82*>
83*> \param[in] N
84*> \verbatim
85*> N is INTEGER
86*> The order of the matrix B, and the number of columns in the
87*> matrices X and C. N >= 0.
88*> \endverbatim
89*>
90*> \param[in] A
91*> \verbatim
92*> A is DOUBLE PRECISION array, dimension (LDA,M)
93*> The upper quasi-triangular matrix A, in Schur canonical form.
94*> \endverbatim
95*>
96*> \param[in] LDA
97*> \verbatim
98*> LDA is INTEGER
99*> The leading dimension of the array A. LDA >= max(1,M).
100*> \endverbatim
101*>
102*> \param[in] B
103*> \verbatim
104*> B is DOUBLE PRECISION array, dimension (LDB,N)
105*> The upper quasi-triangular matrix B, in Schur canonical form.
106*> \endverbatim
107*>
108*> \param[in] LDB
109*> \verbatim
110*> LDB is INTEGER
111*> The leading dimension of the array B. LDB >= max(1,N).
112*> \endverbatim
113*>
114*> \param[in,out] C
115*> \verbatim
116*> C is DOUBLE PRECISION array, dimension (LDC,N)
117*> On entry, the M-by-N right hand side matrix C.
118*> On exit, C is overwritten by the solution matrix X.
119*> \endverbatim
120*>
121*> \param[in] LDC
122*> \verbatim
123*> LDC is INTEGER
124*> The leading dimension of the array C. LDC >= max(1,M)
125*> \endverbatim
126*>
127*> \param[out] SCALE
128*> \verbatim
129*> SCALE is DOUBLE PRECISION
130*> The scale factor, scale, set <= 1 to avoid overflow in X.
131*> \endverbatim
132*>
133*> \param[out] IWORK
134*> \verbatim
135*> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
136*> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
137*> \endverbatim
138*>
139*> \param[in] LIWORK
140*> \verbatim
141*> IWORK is INTEGER
142*> The dimension of the array IWORK. LIWORK >= ((M + NB - 1) / NB + 1)
143*> + ((N + NB - 1) / NB + 1), where NB is the optimal block size.
144*>
145*> If LIWORK = -1, then a workspace query is assumed; the routine
146*> only calculates the optimal dimension of the IWORK array,
147*> returns this value as the first entry of the IWORK array, and
148*> no error message related to LIWORK is issued by XERBLA.
149*> \endverbatim
150*>
151*> \param[out] SWORK
152*> \verbatim
153*> SWORK is DOUBLE PRECISION array, dimension (MAX(2, ROWS),
154*> MAX(1,COLS)).
155*> On exit, if INFO = 0, SWORK(1) returns the optimal value ROWS
156*> and SWORK(2) returns the optimal COLS.
157*> \endverbatim
158*>
159*> \param[in] LDSWORK
160*> \verbatim
161*> LDSWORK is INTEGER
162*> LDSWORK >= MAX(2,ROWS), where ROWS = ((M + NB - 1) / NB + 1)
163*> and NB is the optimal block size.
164*>
165*> If LDSWORK = -1, then a workspace query is assumed; the routine
166*> only calculates the optimal dimensions of the SWORK matrix,
167*> returns these values as the first and second entry of the SWORK
168*> matrix, and no error message related LWORK is issued by XERBLA.
169*> \endverbatim
170*>
171*> \param[out] INFO
172*> \verbatim
173*> INFO is INTEGER
174*> = 0: successful exit
175*> < 0: if INFO = -i, the i-th argument had an illegal value
176*> = 1: A and B have common or very close eigenvalues; perturbed
177*> values were used to solve the equation (but the matrices
178*> A and B are unchanged).
179*> \endverbatim
180*
181*> \ingroup trsyl3
182*
183* =====================================================================
184* References:
185* E. S. Quintana-Orti and R. A. Van De Geijn (2003). Formal derivation of
186* algorithms: The triangular Sylvester equation, ACM Transactions
187* on Mathematical Software (TOMS), volume 29, pages 218--243.
188*
189* A. Schwarz and C. C. Kjelgaard Mikkelsen (2020). Robust Task-Parallel
190* Solution of the Triangular Sylvester Equation. Lecture Notes in
191* Computer Science, vol 12043, pages 82--92, Springer.
192*
193* Contributor:
194* Angelika Schwarz, Umea University, Sweden.
195*
196* =====================================================================
197 SUBROUTINE dtrsyl3( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB,
198 $ C, LDC, SCALE, IWORK, LIWORK, SWORK,
199 $ LDSWORK, INFO )
200 IMPLICIT NONE
201*
202* .. Scalar Arguments ..
203 CHARACTER TRANA, TRANB
204 INTEGER INFO, ISGN, LDA, LDB, LDC, M, N,
205 $ liwork, ldswork
206 DOUBLE PRECISION SCALE
207* ..
208* .. Array Arguments ..
209 INTEGER IWORK( * )
210 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ),
211 $ swork( ldswork, * )
212* ..
213* .. Parameters ..
214 DOUBLE PRECISION ZERO, ONE
215 PARAMETER ( ZERO = 0.0d+0, one = 1.0d+0 )
216* ..
217* .. Local Scalars ..
218 LOGICAL NOTRNA, NOTRNB, LQUERY, SKIP
219 INTEGER AWRK, BWRK, I, I1, I2, IINFO, J, J1, J2, JJ,
220 $ k, k1, k2, l, l1, l2, ll, nba, nb, nbb, pc
221 DOUBLE PRECISION ANRM, BIGNUM, BNRM, CNRM, SCAL, SCALOC,
222 $ SCAMIN, SGN, XNRM, BUF, SMLNUM
223* ..
224* .. Local Arrays ..
225 DOUBLE PRECISION WNRM( MAX( M, N ) )
226* ..
227* .. External Functions ..
228 LOGICAL LSAME
229 INTEGER ILAENV
230 DOUBLE PRECISION DLANGE, DLAMCH, DLARMM
231 EXTERNAL dlange, dlamch, dlarmm, ilaenv,
232 $ lsame
233* ..
234* .. External Subroutines ..
235 EXTERNAL dgemm, dlascl, dscal, dtrsyl,
236 $ xerbla
237* ..
238* .. Intrinsic Functions ..
239 INTRINSIC abs, dble, exponent, max, min
240* ..
241* .. Executable Statements ..
242*
243* Decode and Test input parameters
244*
245 notrna = lsame( trana, 'N' )
246 notrnb = lsame( tranb, 'N' )
247*
248* Use the same block size for all matrices.
249*
250 nb = max(8, ilaenv( 1, 'DTRSYL', '', m, n, -1, -1) )
251*
252* Compute number of blocks in A and B
253*
254 nba = max( 1, (m + nb - 1) / nb )
255 nbb = max( 1, (n + nb - 1) / nb )
256*
257* Compute workspace
258*
259 info = 0
260 lquery = ( liwork.EQ.-1 .OR. ldswork.EQ.-1 )
261 iwork( 1 ) = nba + nbb + 2
262 IF( lquery ) THEN
263 ldswork = 2
264 swork( 1, 1 ) = max( nba, nbb )
265 swork( 2, 1 ) = 2 * nbb + nba
266 END IF
267*
268* Test the input arguments
269*
270 IF( .NOT.notrna .AND. .NOT.lsame( trana, 'T' ) .AND. .NOT.
271 $ lsame( trana, 'C' ) ) THEN
272 info = -1
273 ELSE IF( .NOT.notrnb .AND. .NOT.lsame( tranb, 'T' ) .AND. .NOT.
274 $ lsame( tranb, 'C' ) ) THEN
275 info = -2
276 ELSE IF( isgn.NE.1 .AND. isgn.NE.-1 ) THEN
277 info = -3
278 ELSE IF( m.LT.0 ) THEN
279 info = -4
280 ELSE IF( n.LT.0 ) THEN
281 info = -5
282 ELSE IF( lda.LT.max( 1, m ) ) THEN
283 info = -7
284 ELSE IF( ldb.LT.max( 1, n ) ) THEN
285 info = -9
286 ELSE IF( ldc.LT.max( 1, m ) ) THEN
287 info = -11
288 END IF
289 IF( info.NE.0 ) THEN
290 CALL xerbla( 'DTRSYL3', -info )
291 RETURN
292 ELSE IF( lquery ) THEN
293 RETURN
294 END IF
295*
296* Quick return if possible
297*
298 scale = one
299 IF( m.EQ.0 .OR. n.EQ.0 )
300 $ RETURN
301*
302* Use unblocked code for small problems or if insufficient
303* workspaces are provided
304*
305 IF( min( nba, nbb ).EQ.1 .OR. ldswork.LT.max( nba, nbb ) .OR.
306 $ liwork.LT.iwork(1) ) THEN
307 CALL dtrsyl( trana, tranb, isgn, m, n, a, lda, b, ldb,
308 $ c, ldc, scale, info )
309 RETURN
310 END IF
311*
312* Set constants to control overflow
313*
314 smlnum = dlamch( 'S' )
315 bignum = one / smlnum
316*
317* Partition A such that 2-by-2 blocks on the diagonal are not split
318*
319 skip = .false.
320 DO i = 1, nba
321 iwork( i ) = ( i - 1 ) * nb + 1
322 END DO
323 iwork( nba + 1 ) = m + 1
324 DO k = 1, nba
325 l1 = iwork( k )
326 l2 = iwork( k + 1 ) - 1
327 DO l = l1, l2
328 IF( skip ) THEN
329 skip = .false.
330 cycle
331 END IF
332 IF( l.GE.m ) THEN
333* A( M, M ) is a 1-by-1 block
334 cycle
335 END IF
336 IF( a( l, l+1 ).NE.zero .AND. a( l+1, l ).NE.zero ) THEN
337* Check if 2-by-2 block is split
338 IF( l + 1 .EQ. iwork( k + 1 ) ) THEN
339 iwork( k + 1 ) = iwork( k + 1 ) + 1
340 cycle
341 END IF
342 skip = .true.
343 END IF
344 END DO
345 END DO
346 iwork( nba + 1 ) = m + 1
347 IF( iwork( nba ).GE.iwork( nba + 1 ) ) THEN
348 iwork( nba ) = iwork( nba + 1 )
349 nba = nba - 1
350 END IF
351*
352* Partition B such that 2-by-2 blocks on the diagonal are not split
353*
354 pc = nba + 1
355 skip = .false.
356 DO i = 1, nbb
357 iwork( pc + i ) = ( i - 1 ) * nb + 1
358 END DO
359 iwork( pc + nbb + 1 ) = n + 1
360 DO k = 1, nbb
361 l1 = iwork( pc + k )
362 l2 = iwork( pc + k + 1 ) - 1
363 DO l = l1, l2
364 IF( skip ) THEN
365 skip = .false.
366 cycle
367 END IF
368 IF( l.GE.n ) THEN
369* B( N, N ) is a 1-by-1 block
370 cycle
371 END IF
372 IF( b( l, l+1 ).NE.zero .AND. b( l+1, l ).NE.zero ) THEN
373* Check if 2-by-2 block is split
374 IF( l + 1 .EQ. iwork( pc + k + 1 ) ) THEN
375 iwork( pc + k + 1 ) = iwork( pc + k + 1 ) + 1
376 cycle
377 END IF
378 skip = .true.
379 END IF
380 END DO
381 END DO
382 iwork( pc + nbb + 1 ) = n + 1
383 IF( iwork( pc + nbb ).GE.iwork( pc + nbb + 1 ) ) THEN
384 iwork( pc + nbb ) = iwork( pc + nbb + 1 )
385 nbb = nbb - 1
386 END IF
387*
388* Set local scaling factors - must never attain zero.
389*
390 DO l = 1, nbb
391 DO k = 1, nba
392 swork( k, l ) = one
393 END DO
394 END DO
395*
396* Fallback scaling factor to prevent flushing of SWORK( K, L ) to zero.
397* This scaling is to ensure compatibility with TRSYL and may get flushed.
398*
399 buf = one
400*
401* Compute upper bounds of blocks of A and B
402*
403 awrk = nbb
404 DO k = 1, nba
405 k1 = iwork( k )
406 k2 = iwork( k + 1 )
407 DO l = k, nba
408 l1 = iwork( l )
409 l2 = iwork( l + 1 )
410 IF( notrna ) THEN
411 swork( k, awrk + l ) = dlange( 'I', k2-k1, l2-l1,
412 $ a( k1, l1 ), lda, wnrm )
413 ELSE
414 swork( l, awrk + k ) = dlange( '1', k2-k1, l2-l1,
415 $ a( k1, l1 ), lda, wnrm )
416 END IF
417 END DO
418 END DO
419 bwrk = nbb + nba
420 DO k = 1, nbb
421 k1 = iwork( pc + k )
422 k2 = iwork( pc + k + 1 )
423 DO l = k, nbb
424 l1 = iwork( pc + l )
425 l2 = iwork( pc + l + 1 )
426 IF( notrnb ) THEN
427 swork( k, bwrk + l ) = dlange( 'I', k2-k1, l2-l1,
428 $ b( k1, l1 ), ldb, wnrm )
429 ELSE
430 swork( l, bwrk + k ) = dlange( '1', k2-k1, l2-l1,
431 $ b( k1, l1 ), ldb, wnrm )
432 END IF
433 END DO
434 END DO
435*
436 sgn = dble( isgn )
437*
438 IF( notrna .AND. notrnb ) THEN
439*
440* Solve A*X + ISGN*X*B = scale*C.
441*
442* The (K,L)th block of X is determined starting from
443* bottom-left corner column by column by
444*
445* A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
446*
447* Where
448* M L-1
449* R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(J,L)].
450* I=K+1 J=1
451*
452* Start loop over block rows (index = K) and block columns (index = L)
453*
454 DO k = nba, 1, -1
455*
456* K1: row index of the first row in X( K, L )
457* K2: row index of the first row in X( K+1, L )
458* so the K2 - K1 is the column count of the block X( K, L )
459*
460 k1 = iwork( k )
461 k2 = iwork( k + 1 )
462 DO l = 1, nbb
463*
464* L1: column index of the first column in X( K, L )
465* L2: column index of the first column in X( K, L + 1)
466* so that L2 - L1 is the row count of the block X( K, L )
467*
468 l1 = iwork( pc + l )
469 l2 = iwork( pc + l + 1 )
470*
471 CALL dtrsyl( trana, tranb, isgn, k2-k1, l2-l1,
472 $ a( k1, k1 ), lda,
473 $ b( l1, l1 ), ldb,
474 $ c( k1, l1 ), ldc, scaloc, iinfo )
475 info = max( info, iinfo )
476*
477 IF ( scaloc * swork( k, l ) .EQ. zero ) THEN
478 IF( scaloc .EQ. zero ) THEN
479* The magnitude of the largest entry of X(K1:K2-1, L1:L2-1)
480* is larger than the product of BIGNUM**2 and cannot be
481* represented in the form (1/SCALE)*X(K1:K2-1, L1:L2-1).
482* Mark the computation as pointless.
483 buf = zero
484 ELSE
485* Use second scaling factor to prevent flushing to zero.
486 buf = buf*2.d0**exponent( scaloc )
487 END IF
488 DO jj = 1, nbb
489 DO ll = 1, nba
490* Bound by BIGNUM to not introduce Inf. The value
491* is irrelevant; corresponding entries of the
492* solution will be flushed in consistency scaling.
493 swork( ll, jj ) = min( bignum,
494 $ swork( ll, jj ) / 2.d0**exponent( scaloc ) )
495 END DO
496 END DO
497 END IF
498 swork( k, l ) = scaloc * swork( k, l )
499 xnrm = dlange( 'I', k2-k1, l2-l1, c( k1, l1 ), ldc,
500 $ wnrm )
501*
502 DO i = k - 1, 1, -1
503*
504* C( I, L ) := C( I, L ) - A( I, K ) * C( K, L )
505*
506 i1 = iwork( i )
507 i2 = iwork( i + 1 )
508*
509* Compute scaling factor to survive the linear update
510* simulating consistent scaling.
511*
512 cnrm = dlange( 'I', i2-i1, l2-l1, c( i1, l1 ),
513 $ ldc, wnrm )
514 scamin = min( swork( i, l ), swork( k, l ) )
515 cnrm = cnrm * ( scamin / swork( i, l ) )
516 xnrm = xnrm * ( scamin / swork( k, l ) )
517 anrm = swork( i, awrk + k )
518 scaloc = dlarmm( anrm, xnrm, cnrm )
519 IF( scaloc * scamin .EQ. zero ) THEN
520* Use second scaling factor to prevent flushing to zero.
521 buf = buf*2.d0**exponent( scaloc )
522 DO jj = 1, nbb
523 DO ll = 1, nba
524 swork( ll, jj ) = min( bignum,
525 $ swork( ll, jj ) / 2.d0**exponent( scaloc ) )
526 END DO
527 END DO
528 scamin = scamin / 2.d0**exponent( scaloc )
529 scaloc = scaloc / 2.d0**exponent( scaloc )
530 END IF
531 cnrm = cnrm * scaloc
532 xnrm = xnrm * scaloc
533*
534* Simultaneously apply the robust update factor and the
535* consistency scaling factor to C( I, L ) and C( K, L ).
536*
537 scal = ( scamin / swork( k, l ) ) * scaloc
538 IF (scal .NE. one) THEN
539 DO jj = l1, l2-1
540 CALL dscal( k2-k1, scal, c( k1, jj ), 1)
541 END DO
542 ENDIF
543*
544 scal = ( scamin / swork( i, l ) ) * scaloc
545 IF (scal .NE. one) THEN
546 DO ll = l1, l2-1
547 CALL dscal( i2-i1, scal, c( i1, ll ), 1)
548 END DO
549 ENDIF
550*
551* Record current scaling factor
552*
553 swork( k, l ) = scamin * scaloc
554 swork( i, l ) = scamin * scaloc
555*
556 CALL dgemm( 'N', 'N', i2-i1, l2-l1, k2-k1, -one,
557 $ a( i1, k1 ), lda, c( k1, l1 ), ldc,
558 $ one, c( i1, l1 ), ldc )
559*
560 END DO
561*
562 DO j = l + 1, nbb
563*
564* C( K, J ) := C( K, J ) - SGN * C( K, L ) * B( L, J )
565*
566 j1 = iwork( pc + j )
567 j2 = iwork( pc + j + 1 )
568*
569* Compute scaling factor to survive the linear update
570* simulating consistent scaling.
571*
572 cnrm = dlange( 'I', k2-k1, j2-j1, c( k1, j1 ),
573 $ ldc, wnrm )
574 scamin = min( swork( k, j ), swork( k, l ) )
575 cnrm = cnrm * ( scamin / swork( k, j ) )
576 xnrm = xnrm * ( scamin / swork( k, l ) )
577 bnrm = swork(l, bwrk + j)
578 scaloc = dlarmm( bnrm, xnrm, cnrm )
579 IF( scaloc * scamin .EQ. zero ) THEN
580* Use second scaling factor to prevent flushing to zero.
581 buf = buf*2.d0**exponent( scaloc )
582 DO jj = 1, nbb
583 DO ll = 1, nba
584 swork( ll, jj ) = min( bignum,
585 $ swork( ll, jj ) / 2.d0**exponent( scaloc ) )
586 END DO
587 END DO
588 scamin = scamin / 2.d0**exponent( scaloc )
589 scaloc = scaloc / 2.d0**exponent( scaloc )
590 END IF
591 cnrm = cnrm * scaloc
592 xnrm = xnrm * scaloc
593*
594* Simultaneously apply the robust update factor and the
595* consistency scaling factor to C( K, J ) and C( K, L).
596*
597 scal = ( scamin / swork( k, l ) ) * scaloc
598 IF( scal .NE. one ) THEN
599 DO ll = l1, l2-1
600 CALL dscal( k2-k1, scal, c( k1, ll ), 1 )
601 END DO
602 ENDIF
603*
604 scal = ( scamin / swork( k, j ) ) * scaloc
605 IF( scal .NE. one ) THEN
606 DO jj = j1, j2-1
607 CALL dscal( k2-k1, scal, c( k1, jj ), 1 )
608 END DO
609 ENDIF
610*
611* Record current scaling factor
612*
613 swork( k, l ) = scamin * scaloc
614 swork( k, j ) = scamin * scaloc
615*
616 CALL dgemm( 'N', 'N', k2-k1, j2-j1, l2-l1, -sgn,
617 $ c( k1, l1 ), ldc, b( l1, j1 ), ldb,
618 $ one, c( k1, j1 ), ldc )
619 END DO
620 END DO
621 END DO
622 ELSE IF( .NOT.notrna .AND. notrnb ) THEN
623*
624* Solve A**T*X + ISGN*X*B = scale*C.
625*
626* The (K,L)th block of X is determined starting from
627* upper-left corner column by column by
628*
629* A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
630*
631* Where
632* K-1 L-1
633* R(K,L) = SUM [A(I,K)**T*X(I,L)] +ISGN*SUM [X(K,J)*B(J,L)]
634* I=1 J=1
635*
636* Start loop over block rows (index = K) and block columns (index = L)
637*
638 DO k = 1, nba
639*
640* K1: row index of the first row in X( K, L )
641* K2: row index of the first row in X( K+1, L )
642* so the K2 - K1 is the column count of the block X( K, L )
643*
644 k1 = iwork( k )
645 k2 = iwork( k + 1 )
646 DO l = 1, nbb
647*
648* L1: column index of the first column in X( K, L )
649* L2: column index of the first column in X( K, L + 1)
650* so that L2 - L1 is the row count of the block X( K, L )
651*
652 l1 = iwork( pc + l )
653 l2 = iwork( pc + l + 1 )
654*
655 CALL dtrsyl( trana, tranb, isgn, k2-k1, l2-l1,
656 $ a( k1, k1 ), lda,
657 $ b( l1, l1 ), ldb,
658 $ c( k1, l1 ), ldc, scaloc, iinfo )
659 info = max( info, iinfo )
660*
661 IF( scaloc * swork( k, l ) .EQ. zero ) THEN
662 IF( scaloc .EQ. zero ) THEN
663* The magnitude of the largest entry of X(K1:K2-1, L1:L2-1)
664* is larger than the product of BIGNUM**2 and cannot be
665* represented in the form (1/SCALE)*X(K1:K2-1, L1:L2-1).
666* Mark the computation as pointless.
667 buf = zero
668 ELSE
669* Use second scaling factor to prevent flushing to zero.
670 buf = buf*2.d0**exponent( scaloc )
671 END IF
672 DO jj = 1, nbb
673 DO ll = 1, nba
674* Bound by BIGNUM to not introduce Inf. The value
675* is irrelevant; corresponding entries of the
676* solution will be flushed in consistency scaling.
677 swork( ll, jj ) = min( bignum,
678 $ swork( ll, jj ) / 2.d0**exponent( scaloc ) )
679 END DO
680 END DO
681 END IF
682 swork( k, l ) = scaloc * swork( k, l )
683 xnrm = dlange( 'I', k2-k1, l2-l1, c( k1, l1 ), ldc,
684 $ wnrm )
685*
686 DO i = k + 1, nba
687*
688* C( I, L ) := C( I, L ) - A( K, I )**T * C( K, L )
689*
690 i1 = iwork( i )
691 i2 = iwork( i + 1 )
692*
693* Compute scaling factor to survive the linear update
694* simulating consistent scaling.
695*
696 cnrm = dlange( 'I', i2-i1, l2-l1, c( i1, l1 ),
697 $ ldc, wnrm )
698 scamin = min( swork( i, l ), swork( k, l ) )
699 cnrm = cnrm * ( scamin / swork( i, l ) )
700 xnrm = xnrm * ( scamin / swork( k, l ) )
701 anrm = swork( i, awrk + k )
702 scaloc = dlarmm( anrm, xnrm, cnrm )
703 IF( scaloc * scamin .EQ. zero ) THEN
704* Use second scaling factor to prevent flushing to zero.
705 buf = buf*2.d0**exponent( scaloc )
706 DO jj = 1, nbb
707 DO ll = 1, nba
708 swork( ll, jj ) = min( bignum,
709 $ swork( ll, jj ) / 2.d0**exponent( scaloc ) )
710 END DO
711 END DO
712 scamin = scamin / 2.d0**exponent( scaloc )
713 scaloc = scaloc / 2.d0**exponent( scaloc )
714 END IF
715 cnrm = cnrm * scaloc
716 xnrm = xnrm * scaloc
717*
718* Simultaneously apply the robust update factor and the
719* consistency scaling factor to to C( I, L ) and C( K, L ).
720*
721 scal = ( scamin / swork( k, l ) ) * scaloc
722 IF (scal .NE. one) THEN
723 DO ll = l1, l2-1
724 CALL dscal( k2-k1, scal, c( k1, ll ), 1 )
725 END DO
726 ENDIF
727*
728 scal = ( scamin / swork( i, l ) ) * scaloc
729 IF (scal .NE. one) THEN
730 DO ll = l1, l2-1
731 CALL dscal( i2-i1, scal, c( i1, ll ), 1 )
732 END DO
733 ENDIF
734*
735* Record current scaling factor
736*
737 swork( k, l ) = scamin * scaloc
738 swork( i, l ) = scamin * scaloc
739*
740 CALL dgemm( 'T', 'N', i2-i1, l2-l1, k2-k1, -one,
741 $ a( k1, i1 ), lda, c( k1, l1 ), ldc,
742 $ one, c( i1, l1 ), ldc )
743 END DO
744*
745 DO j = l + 1, nbb
746*
747* C( K, J ) := C( K, J ) - SGN * C( K, L ) * B( L, J )
748*
749 j1 = iwork( pc + j )
750 j2 = iwork( pc + j + 1 )
751*
752* Compute scaling factor to survive the linear update
753* simulating consistent scaling.
754*
755 cnrm = dlange( 'I', k2-k1, j2-j1, c( k1, j1 ),
756 $ ldc, wnrm )
757 scamin = min( swork( k, j ), swork( k, l ) )
758 cnrm = cnrm * ( scamin / swork( k, j ) )
759 xnrm = xnrm * ( scamin / swork( k, l ) )
760 bnrm = swork( l, bwrk + j )
761 scaloc = dlarmm( bnrm, xnrm, cnrm )
762 IF( scaloc * scamin .EQ. zero ) THEN
763* Use second scaling factor to prevent flushing to zero.
764 buf = buf*2.d0**exponent( scaloc )
765 DO jj = 1, nbb
766 DO ll = 1, nba
767 swork( ll, jj ) = min( bignum,
768 $ swork( ll, jj ) / 2.d0**exponent( scaloc ) )
769 END DO
770 END DO
771 scamin = scamin / 2.d0**exponent( scaloc )
772 scaloc = scaloc / 2.d0**exponent( scaloc )
773 END IF
774 cnrm = cnrm * scaloc
775 xnrm = xnrm * scaloc
776*
777* Simultaneously apply the robust update factor and the
778* consistency scaling factor to to C( K, J ) and C( K, L ).
779*
780 scal = ( scamin / swork( k, l ) ) * scaloc
781 IF( scal .NE. one ) THEN
782 DO ll = l1, l2-1
783 CALL dscal( k2-k1, scal, c( k1, ll ), 1 )
784 END DO
785 ENDIF
786*
787 scal = ( scamin / swork( k, j ) ) * scaloc
788 IF( scal .NE. one ) THEN
789 DO jj = j1, j2-1
790 CALL dscal( k2-k1, scal, c( k1, jj ), 1 )
791 END DO
792 ENDIF
793*
794* Record current scaling factor
795*
796 swork( k, l ) = scamin * scaloc
797 swork( k, j ) = scamin * scaloc
798*
799 CALL dgemm( 'N', 'N', k2-k1, j2-j1, l2-l1, -sgn,
800 $ c( k1, l1 ), ldc, b( l1, j1 ), ldb,
801 $ one, c( k1, j1 ), ldc )
802 END DO
803 END DO
804 END DO
805 ELSE IF( .NOT.notrna .AND. .NOT.notrnb ) THEN
806*
807* Solve A**T*X + ISGN*X*B**T = scale*C.
808*
809* The (K,L)th block of X is determined starting from
810* top-right corner column by column by
811*
812* A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L)**T = C(K,L) - R(K,L)
813*
814* Where
815* K-1 N
816* R(K,L) = SUM [A(I,K)**T*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**T].
817* I=1 J=L+1
818*
819* Start loop over block rows (index = K) and block columns (index = L)
820*
821 DO k = 1, nba
822*
823* K1: row index of the first row in X( K, L )
824* K2: row index of the first row in X( K+1, L )
825* so the K2 - K1 is the column count of the block X( K, L )
826*
827 k1 = iwork( k )
828 k2 = iwork( k + 1 )
829 DO l = nbb, 1, -1
830*
831* L1: column index of the first column in X( K, L )
832* L2: column index of the first column in X( K, L + 1)
833* so that L2 - L1 is the row count of the block X( K, L )
834*
835 l1 = iwork( pc + l )
836 l2 = iwork( pc + l + 1 )
837*
838 CALL dtrsyl( trana, tranb, isgn, k2-k1, l2-l1,
839 $ a( k1, k1 ), lda,
840 $ b( l1, l1 ), ldb,
841 $ c( k1, l1 ), ldc, scaloc, iinfo )
842 info = max( info, iinfo )
843*
844 swork( k, l ) = scaloc * swork( k, l )
845 IF( scaloc * swork( k, l ) .EQ. zero ) THEN
846 IF( scaloc .EQ. zero ) THEN
847* The magnitude of the largest entry of X(K1:K2-1, L1:L2-1)
848* is larger than the product of BIGNUM**2 and cannot be
849* represented in the form (1/SCALE)*X(K1:K2-1, L1:L2-1).
850* Mark the computation as pointless.
851 buf = zero
852 ELSE
853* Use second scaling factor to prevent flushing to zero.
854 buf = buf*2.d0**exponent( scaloc )
855 END IF
856 DO jj = 1, nbb
857 DO ll = 1, nba
858* Bound by BIGNUM to not introduce Inf. The value
859* is irrelevant; corresponding entries of the
860* solution will be flushed in consistency scaling.
861 swork( ll, jj ) = min( bignum,
862 $ swork( ll, jj ) / 2.d0**exponent( scaloc ) )
863 END DO
864 END DO
865 END IF
866 xnrm = dlange( 'I', k2-k1, l2-l1, c( k1, l1 ), ldc,
867 $ wnrm )
868*
869 DO i = k + 1, nba
870*
871* C( I, L ) := C( I, L ) - A( K, I )**T * C( K, L )
872*
873 i1 = iwork( i )
874 i2 = iwork( i + 1 )
875*
876* Compute scaling factor to survive the linear update
877* simulating consistent scaling.
878*
879 cnrm = dlange( 'I', i2-i1, l2-l1, c( i1, l1 ),
880 $ ldc, wnrm )
881 scamin = min( swork( i, l ), swork( k, l ) )
882 cnrm = cnrm * ( scamin / swork( i, l ) )
883 xnrm = xnrm * ( scamin / swork( k, l ) )
884 anrm = swork( i, awrk + k )
885 scaloc = dlarmm( anrm, xnrm, cnrm )
886 IF( scaloc * scamin .EQ. zero ) THEN
887* Use second scaling factor to prevent flushing to zero.
888 buf = buf*2.d0**exponent( scaloc )
889 DO jj = 1, nbb
890 DO ll = 1, nba
891 swork( ll, jj ) = min( bignum,
892 $ swork( ll, jj ) / 2.d0**exponent( scaloc ) )
893 END DO
894 END DO
895 scamin = scamin / 2.d0**exponent( scaloc )
896 scaloc = scaloc / 2.d0**exponent( scaloc )
897 END IF
898 cnrm = cnrm * scaloc
899 xnrm = xnrm * scaloc
900*
901* Simultaneously apply the robust update factor and the
902* consistency scaling factor to C( I, L ) and C( K, L ).
903*
904 scal = ( scamin / swork( k, l ) ) * scaloc
905 IF (scal .NE. one) THEN
906 DO ll = l1, l2-1
907 CALL dscal( k2-k1, scal, c( k1, ll ), 1 )
908 END DO
909 ENDIF
910*
911 scal = ( scamin / swork( i, l ) ) * scaloc
912 IF (scal .NE. one) THEN
913 DO ll = l1, l2-1
914 CALL dscal( i2-i1, scal, c( i1, ll ), 1 )
915 END DO
916 ENDIF
917*
918* Record current scaling factor
919*
920 swork( k, l ) = scamin * scaloc
921 swork( i, l ) = scamin * scaloc
922*
923 CALL dgemm( 'T', 'N', i2-i1, l2-l1, k2-k1, -one,
924 $ a( k1, i1 ), lda, c( k1, l1 ), ldc,
925 $ one, c( i1, l1 ), ldc )
926 END DO
927*
928 DO j = 1, l - 1
929*
930* C( K, J ) := C( K, J ) - SGN * C( K, L ) * B( J, L )**T
931*
932 j1 = iwork( pc + j )
933 j2 = iwork( pc + j + 1 )
934*
935* Compute scaling factor to survive the linear update
936* simulating consistent scaling.
937*
938 cnrm = dlange( 'I', k2-k1, j2-j1, c( k1, j1 ),
939 $ ldc, wnrm )
940 scamin = min( swork( k, j ), swork( k, l ) )
941 cnrm = cnrm * ( scamin / swork( k, j ) )
942 xnrm = xnrm * ( scamin / swork( k, l ) )
943 bnrm = swork( l, bwrk + j )
944 scaloc = dlarmm( bnrm, xnrm, cnrm )
945 IF( scaloc * scamin .EQ. zero ) THEN
946* Use second scaling factor to prevent flushing to zero.
947 buf = buf*2.d0**exponent( scaloc )
948 DO jj = 1, nbb
949 DO ll = 1, nba
950 swork( ll, jj ) = min( bignum,
951 $ swork( ll, jj ) / 2.d0**exponent( scaloc ) )
952 END DO
953 END DO
954 scamin = scamin / 2.d0**exponent( scaloc )
955 scaloc = scaloc / 2.d0**exponent( scaloc )
956 END IF
957 cnrm = cnrm * scaloc
958 xnrm = xnrm * scaloc
959*
960* Simultaneously apply the robust update factor and the
961* consistency scaling factor to C( K, J ) and C( K, L ).
962*
963 scal = ( scamin / swork( k, l ) ) * scaloc
964 IF( scal .NE. one ) THEN
965 DO ll = l1, l2-1
966 CALL dscal( k2-k1, scal, c( k1, ll ), 1)
967 END DO
968 ENDIF
969*
970 scal = ( scamin / swork( k, j ) ) * scaloc
971 IF( scal .NE. one ) THEN
972 DO jj = j1, j2-1
973 CALL dscal( k2-k1, scal, c( k1, jj ), 1 )
974 END DO
975 ENDIF
976*
977* Record current scaling factor
978*
979 swork( k, l ) = scamin * scaloc
980 swork( k, j ) = scamin * scaloc
981*
982 CALL dgemm( 'N', 'T', k2-k1, j2-j1, l2-l1, -sgn,
983 $ c( k1, l1 ), ldc, b( j1, l1 ), ldb,
984 $ one, c( k1, j1 ), ldc )
985 END DO
986 END DO
987 END DO
988 ELSE IF( notrna .AND. .NOT.notrnb ) THEN
989*
990* Solve A*X + ISGN*X*B**T = scale*C.
991*
992* The (K,L)th block of X is determined starting from
993* bottom-right corner column by column by
994*
995* A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L)**T = C(K,L) - R(K,L)
996*
997* Where
998* M N
999* R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**T].
1000* I=K+1 J=L+1
1001*
1002* Start loop over block rows (index = K) and block columns (index = L)
1003*
1004 DO k = nba, 1, -1
1005*
1006* K1: row index of the first row in X( K, L )
1007* K2: row index of the first row in X( K+1, L )
1008* so the K2 - K1 is the column count of the block X( K, L )
1009*
1010 k1 = iwork( k )
1011 k2 = iwork( k + 1 )
1012 DO l = nbb, 1, -1
1013*
1014* L1: column index of the first column in X( K, L )
1015* L2: column index of the first column in X( K, L + 1)
1016* so that L2 - L1 is the row count of the block X( K, L )
1017*
1018 l1 = iwork( pc + l )
1019 l2 = iwork( pc + l + 1 )
1020*
1021 CALL dtrsyl( trana, tranb, isgn, k2-k1, l2-l1,
1022 $ a( k1, k1 ), lda,
1023 $ b( l1, l1 ), ldb,
1024 $ c( k1, l1 ), ldc, scaloc, iinfo )
1025 info = max( info, iinfo )
1026*
1027 IF( scaloc * swork( k, l ) .EQ. zero ) THEN
1028 IF( scaloc .EQ. zero ) THEN
1029* The magnitude of the largest entry of X(K1:K2-1, L1:L2-1)
1030* is larger than the product of BIGNUM**2 and cannot be
1031* represented in the form (1/SCALE)*X(K1:K2-1, L1:L2-1).
1032* Mark the computation as pointless.
1033 buf = zero
1034 ELSE
1035* Use second scaling factor to prevent flushing to zero.
1036 buf = buf*2.d0**exponent( scaloc )
1037 END IF
1038 DO jj = 1, nbb
1039 DO ll = 1, nba
1040* Bound by BIGNUM to not introduce Inf. The value
1041* is irrelevant; corresponding entries of the
1042* solution will be flushed in consistency scaling.
1043 swork( ll, jj ) = min( bignum,
1044 $ swork( ll, jj ) / 2.d0**exponent( scaloc ) )
1045 END DO
1046 END DO
1047 END IF
1048 swork( k, l ) = scaloc * swork( k, l )
1049 xnrm = dlange( 'I', k2-k1, l2-l1, c( k1, l1 ), ldc,
1050 $ wnrm )
1051*
1052 DO i = 1, k - 1
1053*
1054* C( I, L ) := C( I, L ) - A( I, K ) * C( K, L )
1055*
1056 i1 = iwork( i )
1057 i2 = iwork( i + 1 )
1058*
1059* Compute scaling factor to survive the linear update
1060* simulating consistent scaling.
1061*
1062 cnrm = dlange( 'I', i2-i1, l2-l1, c( i1, l1 ),
1063 $ ldc, wnrm )
1064 scamin = min( swork( i, l ), swork( k, l ) )
1065 cnrm = cnrm * ( scamin / swork( i, l ) )
1066 xnrm = xnrm * ( scamin / swork( k, l ) )
1067 anrm = swork( i, awrk + k )
1068 scaloc = dlarmm( anrm, xnrm, cnrm )
1069 IF( scaloc * scamin .EQ. zero ) THEN
1070* Use second scaling factor to prevent flushing to zero.
1071 buf = buf*2.d0**exponent( scaloc )
1072 DO jj = 1, nbb
1073 DO ll = 1, nba
1074 swork( ll, jj ) = min( bignum,
1075 $ swork( ll, jj ) / 2.d0**exponent( scaloc ) )
1076 END DO
1077 END DO
1078 scamin = scamin / 2.d0**exponent( scaloc )
1079 scaloc = scaloc / 2.d0**exponent( scaloc )
1080 END IF
1081 cnrm = cnrm * scaloc
1082 xnrm = xnrm * scaloc
1083*
1084* Simultaneously apply the robust update factor and the
1085* consistency scaling factor to C( I, L ) and C( K, L ).
1086*
1087 scal = ( scamin / swork( k, l ) ) * scaloc
1088 IF (scal .NE. one) THEN
1089 DO ll = l1, l2-1
1090 CALL dscal( k2-k1, scal, c( k1, ll ), 1 )
1091 END DO
1092 ENDIF
1093*
1094 scal = ( scamin / swork( i, l ) ) * scaloc
1095 IF (scal .NE. one) THEN
1096 DO ll = l1, l2-1
1097 CALL dscal( i2-i1, scal, c( i1, ll ), 1 )
1098 END DO
1099 ENDIF
1100*
1101* Record current scaling factor
1102*
1103 swork( k, l ) = scamin * scaloc
1104 swork( i, l ) = scamin * scaloc
1105*
1106 CALL dgemm( 'N', 'N', i2-i1, l2-l1, k2-k1, -one,
1107 $ a( i1, k1 ), lda, c( k1, l1 ), ldc,
1108 $ one, c( i1, l1 ), ldc )
1109*
1110 END DO
1111*
1112 DO j = 1, l - 1
1113*
1114* C( K, J ) := C( K, J ) - SGN * C( K, L ) * B( J, L )**T
1115*
1116 j1 = iwork( pc + j )
1117 j2 = iwork( pc + j + 1 )
1118*
1119* Compute scaling factor to survive the linear update
1120* simulating consistent scaling.
1121*
1122 cnrm = dlange( 'I', k2-k1, j2-j1, c( k1, j1 ),
1123 $ ldc, wnrm )
1124 scamin = min( swork( k, j ), swork( k, l ) )
1125 cnrm = cnrm * ( scamin / swork( k, j ) )
1126 xnrm = xnrm * ( scamin / swork( k, l ) )
1127 bnrm = swork( l, bwrk + j )
1128 scaloc = dlarmm( bnrm, xnrm, cnrm )
1129 IF( scaloc * scamin .EQ. zero ) THEN
1130* Use second scaling factor to prevent flushing to zero.
1131 buf = buf*2.d0**exponent( scaloc )
1132 DO jj = 1, nbb
1133 DO ll = 1, nba
1134 swork( ll, jj ) = min( bignum,
1135 $ swork( ll, jj ) / 2.d0**exponent( scaloc ) )
1136 END DO
1137 END DO
1138 scamin = scamin / 2.d0**exponent( scaloc )
1139 scaloc = scaloc / 2.d0**exponent( scaloc )
1140 END IF
1141 cnrm = cnrm * scaloc
1142 xnrm = xnrm * scaloc
1143*
1144* Simultaneously apply the robust update factor and the
1145* consistency scaling factor to C( K, J ) and C( K, L ).
1146*
1147 scal = ( scamin / swork( k, l ) ) * scaloc
1148 IF( scal .NE. one ) THEN
1149 DO jj = l1, l2-1
1150 CALL dscal( k2-k1, scal, c( k1, jj ), 1 )
1151 END DO
1152 ENDIF
1153*
1154 scal = ( scamin / swork( k, j ) ) * scaloc
1155 IF( scal .NE. one ) THEN
1156 DO jj = j1, j2-1
1157 CALL dscal( k2-k1, scal, c( k1, jj ), 1 )
1158 END DO
1159 ENDIF
1160*
1161* Record current scaling factor
1162*
1163 swork( k, l ) = scamin * scaloc
1164 swork( k, j ) = scamin * scaloc
1165*
1166 CALL dgemm( 'N', 'T', k2-k1, j2-j1, l2-l1, -sgn,
1167 $ c( k1, l1 ), ldc, b( j1, l1 ), ldb,
1168 $ one, c( k1, j1 ), ldc )
1169 END DO
1170 END DO
1171 END DO
1172*
1173 END IF
1174*
1175* Reduce local scaling factors
1176*
1177 scale = swork( 1, 1 )
1178 DO k = 1, nba
1179 DO l = 1, nbb
1180 scale = min( scale, swork( k, l ) )
1181 END DO
1182 END DO
1183*
1184 IF( scale .EQ. zero ) THEN
1185*
1186* The magnitude of the largest entry of the solution is larger
1187* than the product of BIGNUM**2 and cannot be represented in the
1188* form (1/SCALE)*X if SCALE is DOUBLE PRECISION. Set SCALE to
1189* zero and give up.
1190*
1191 iwork(1) = nba + nbb + 2
1192 swork(1,1) = max( nba, nbb )
1193 swork(2,1) = 2 * nbb + nba
1194 RETURN
1195 END IF
1196*
1197* Realize consistent scaling
1198*
1199 DO k = 1, nba
1200 k1 = iwork( k )
1201 k2 = iwork( k + 1 )
1202 DO l = 1, nbb
1203 l1 = iwork( pc + l )
1204 l2 = iwork( pc + l + 1 )
1205 scal = scale / swork( k, l )
1206 IF( scal .NE. one ) THEN
1207 DO ll = l1, l2-1
1208 CALL dscal( k2-k1, scal, c( k1, ll ), 1 )
1209 END DO
1210 ENDIF
1211 END DO
1212 END DO
1213*
1214 IF( buf .NE. one .AND. buf.GT.zero ) THEN
1215*
1216* Decrease SCALE as much as possible.
1217*
1218 scaloc = min( scale / smlnum, one / buf )
1219 buf = buf * scaloc
1220 scale = scale / scaloc
1221 END IF
1222
1223 IF( buf.NE.one .AND. buf.GT.zero ) THEN
1224*
1225* In case of overly aggressive scaling during the computation,
1226* flushing of the global scale factor may be prevented by
1227* undoing some of the scaling. This step is to ensure that
1228* this routine flushes only scale factors that TRSYL also
1229* flushes and be usable as a drop-in replacement.
1230*
1231* How much can the normwise largest entry be upscaled?
1232*
1233 scal = c( 1, 1 )
1234 DO k = 1, m
1235 DO l = 1, n
1236 scal = max( scal, abs( c( k, l ) ) )
1237 END DO
1238 END DO
1239*
1240* Increase BUF as close to 1 as possible and apply scaling.
1241*
1242 scaloc = min( bignum / scal, one / buf )
1243 buf = buf * scaloc
1244 CALL dlascl( 'G', -1, -1, one, scaloc, m, n, c, ldc,
1245 $ iwork(1) )
1246 END IF
1247*
1248* Combine with buffer scaling factor. SCALE will be flushed if
1249* BUF is less than one here.
1250*
1251 scale = scale * buf
1252*
1253* Restore workspace dimensions
1254*
1255 iwork(1) = nba + nbb + 2
1256 swork(1,1) = max( nba, nbb )
1257 swork(2,1) = 2 * nbb + nba
1258*
1259 RETURN
1260*
1261* End of DTRSYL3
1262*
1263 END
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
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:142
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
dtrsyl3
subroutine dtrsyl3(trana, tranb, isgn, m, n, a, lda, b, ldb, c, ldc, scale, iwork, liwork, swork, ldswork, info)
DTRSYL3
Definition dtrsyl3.f:200
dtrsyl
subroutine dtrsyl(trana, tranb, isgn, m, n, a, lda, b, ldb, c, ldc, scale, info)
DTRSYL
Definition dtrsyl.f:162