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