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