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