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