OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
dtrsyl.f
Go to the documentation of this file.
1*> \brief \b DTRSYL
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DTRSYL + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrsyl.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrsyl.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrsyl.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C,
22* LDC, SCALE, INFO )
23*
24* .. Scalar Arguments ..
25* CHARACTER TRANA, TRANB
26* INTEGER INFO, ISGN, LDA, LDB, LDC, M, N
27* DOUBLE PRECISION SCALE
28* ..
29* .. Array Arguments ..
30* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> DTRSYL solves the real Sylvester matrix equation:
40*>
41*> op(A)*X + X*op(B) = scale*C or
42*> op(A)*X - X*op(B) = scale*C,
43*>
44*> where op(A) = A or A**T, and A and B are both upper quasi-
45*> triangular. A is M-by-M and B is N-by-N; the right hand side C and
46*> the solution X are M-by-N; and scale is an output scale factor, set
47*> <= 1 to avoid overflow in X.
48*>
49*> A and B must be in Schur canonical form (as returned by DHSEQR), that
50*> is, block upper triangular with 1-by-1 and 2-by-2 diagonal blocks;
51*> each 2-by-2 diagonal block has its diagonal elements equal and its
52*> off-diagonal elements of opposite sign.
53*> \endverbatim
54*
55* Arguments:
56* ==========
57*
58*> \param[in] TRANA
59*> \verbatim
60*> TRANA is CHARACTER*1
61*> Specifies the option op(A):
62*> = 'N': op(A) = A (No transpose)
63*> = 'T': op(A) = A**T (Transpose)
64*> = 'C': op(A) = A**H (Conjugate transpose = Transpose)
65*> \endverbatim
66*>
67*> \param[in] TRANB
68*> \verbatim
69*> TRANB is CHARACTER*1
70*> Specifies the option op(B):
71*> = 'N': op(B) = B (No transpose)
72*> = 'T': op(B) = B**T (Transpose)
73*> = 'C': op(B) = B**H (Conjugate transpose = Transpose)
74*> \endverbatim
75*>
76*> \param[in] ISGN
77*> \verbatim
78*> ISGN is INTEGER
79*> Specifies the sign in the equation:
80*> = +1: solve op(A)*X + X*op(B) = scale*C
81*> = -1: solve op(A)*X - X*op(B) = scale*C
82*> \endverbatim
83*>
84*> \param[in] M
85*> \verbatim
86*> M is INTEGER
87*> The order of the matrix A, and the number of rows in the
88*> matrices X and C. M >= 0.
89*> \endverbatim
90*>
91*> \param[in] N
92*> \verbatim
93*> N is INTEGER
94*> The order of the matrix B, and the number of columns in the
95*> matrices X and C. N >= 0.
96*> \endverbatim
97*>
98*> \param[in] A
99*> \verbatim
100*> A is DOUBLE PRECISION array, dimension (LDA,M)
101*> The upper quasi-triangular matrix A, in Schur canonical form.
102*> \endverbatim
103*>
104*> \param[in] LDA
105*> \verbatim
106*> LDA is INTEGER
107*> The leading dimension of the array A. LDA >= max(1,M).
108*> \endverbatim
109*>
110*> \param[in] B
111*> \verbatim
112*> B is DOUBLE PRECISION array, dimension (LDB,N)
113*> The upper quasi-triangular matrix B, in Schur canonical form.
114*> \endverbatim
115*>
116*> \param[in] LDB
117*> \verbatim
118*> LDB is INTEGER
119*> The leading dimension of the array B. LDB >= max(1,N).
120*> \endverbatim
121*>
122*> \param[in,out] C
123*> \verbatim
124*> C is DOUBLE PRECISION array, dimension (LDC,N)
125*> On entry, the M-by-N right hand side matrix C.
126*> On exit, C is overwritten by the solution matrix X.
127*> \endverbatim
128*>
129*> \param[in] LDC
130*> \verbatim
131*> LDC is INTEGER
132*> The leading dimension of the array C. LDC >= max(1,M)
133*> \endverbatim
134*>
135*> \param[out] SCALE
136*> \verbatim
137*> SCALE is DOUBLE PRECISION
138*> The scale factor, scale, set <= 1 to avoid overflow in X.
139*> \endverbatim
140*>
141*> \param[out] INFO
142*> \verbatim
143*> INFO is INTEGER
144*> = 0: successful exit
145*> < 0: if INFO = -i, the i-th argument had an illegal value
146*> = 1: A and B have common or very close eigenvalues; perturbed
147*> values were used to solve the equation (but the matrices
148*> A and B are unchanged).
149*> \endverbatim
150*
151* Authors:
152* ========
153*
154*> \author Univ. of Tennessee
155*> \author Univ. of California Berkeley
156*> \author Univ. of Colorado Denver
157*> \author NAG Ltd.
158*
159*> \ingroup doubleSYcomputational
160*
161* =====================================================================
162 SUBROUTINE dtrsyl( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C,
163 $ LDC, SCALE, INFO )
164*
165* -- LAPACK computational routine --
166* -- LAPACK is a software package provided by Univ. of Tennessee, --
167* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
168*
169* .. Scalar Arguments ..
170 CHARACTER TRANA, TRANB
171 INTEGER INFO, ISGN, LDA, LDB, LDC, M, N
172 DOUBLE PRECISION SCALE
173* ..
174* .. Array Arguments ..
175 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * )
176* ..
177*
178* =====================================================================
179*
180* .. Parameters ..
181 DOUBLE PRECISION ZERO, ONE
182 parameter( zero = 0.0d+0, one = 1.0d+0 )
183* ..
184* .. Local Scalars ..
185 LOGICAL NOTRNA, NOTRNB
186 INTEGER IERR, J, K, K1, K2, KNEXT, L, L1, L2, LNEXT
187 DOUBLE PRECISION A11, BIGNUM, DA11, DB, EPS, SCALOC, SGN, SMIN,
188 $ smlnum, suml, sumr, xnorm
189* ..
190* .. Local Arrays ..
191 DOUBLE PRECISION DUM( 1 ), VEC( 2, 2 ), X( 2, 2 )
192* ..
193* .. External Functions ..
194 LOGICAL LSAME
195 DOUBLE PRECISION DDOT, DLAMCH, DLANGE
196 EXTERNAL lsame, ddot, dlamch, dlange
197* ..
198* .. External Subroutines ..
199 EXTERNAL dlabad, dlaln2, dlasy2, dscal, xerbla
200* ..
201* .. Intrinsic Functions ..
202 INTRINSIC abs, dble, max, min
203* ..
204* .. Executable Statements ..
205*
206* Decode and Test input parameters
207*
208 notrna = lsame( trana, 'N' )
209 notrnb = lsame( tranb, 'N' )
210*
211 info = 0
212 IF( .NOT.notrna .AND. .NOT.lsame( trana, 'T' ) .AND. .NOT.
213 $ lsame( trana, 'C' ) ) THEN
214 info = -1
215 ELSE IF( .NOT.notrnb .AND. .NOT.lsame( tranb, 'T' ) .AND. .NOT.
216 $ lsame( tranb, 'C' ) ) THEN
217 info = -2
218 ELSE IF( isgn.NE.1 .AND. isgn.NE.-1 ) THEN
219 info = -3
220 ELSE IF( m.LT.0 ) THEN
221 info = -4
222 ELSE IF( n.LT.0 ) THEN
223 info = -5
224 ELSE IF( lda.LT.max( 1, m ) ) THEN
225 info = -7
226 ELSE IF( ldb.LT.max( 1, n ) ) THEN
227 info = -9
228 ELSE IF( ldc.LT.max( 1, m ) ) THEN
229 info = -11
230 END IF
231 IF( info.NE.0 ) THEN
232 CALL xerbla( 'DTRSYL', -info )
233 RETURN
234 END IF
235*
236* Quick return if possible
237*
238 scale = one
239 IF( m.EQ.0 .OR. n.EQ.0 )
240 $ RETURN
241*
242* Set constants to control overflow
243*
244 eps = dlamch( 'P' )
245 smlnum = dlamch( 'S' )
246 bignum = one / smlnum
247 CALL dlabad( smlnum, bignum )
248 smlnum = smlnum*dble( m*n ) / eps
249 bignum = one / smlnum
250*
251 smin = max( smlnum, eps*dlange( 'M', m, m, a, lda, dum ),
252 $ eps*dlange( 'M', n, n, b, ldb, dum ) )
253*
254 sgn = isgn
255*
256 IF( notrna .AND. notrnb ) THEN
257*
258* Solve A*X + ISGN*X*B = scale*C.
259*
260* The (K,L)th block of X is determined starting from
261* bottom-left corner column by column by
262*
263* A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
264*
265* Where
266* M L-1
267* R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(J,L)].
268* I=K+1 J=1
269*
270* Start column loop (index = L)
271* L1 (L2) : column index of the first (first) row of X(K,L).
272*
273 lnext = 1
274 DO 60 l = 1, n
275 IF( l.LT.lnext )
276 $ GO TO 60
277 IF( l.EQ.n ) THEN
278 l1 = l
279 l2 = l
280 ELSE
281 IF( b( l+1, l ).NE.zero ) THEN
282 l1 = l
283 l2 = l + 1
284 lnext = l + 2
285 ELSE
286 l1 = l
287 l2 = l
288 lnext = l + 1
289 END IF
290 END IF
291*
292* Start row loop (index = K)
293* K1 (K2): row index of the first (last) row of X(K,L).
294*
295 knext = m
296 DO 50 k = m, 1, -1
297 IF( k.GT.knext )
298 $ GO TO 50
299 IF( k.EQ.1 ) THEN
300 k1 = k
301 k2 = k
302 ELSE
303 IF( a( k, k-1 ).NE.zero ) THEN
304 k1 = k - 1
305 k2 = k
306 knext = k - 2
307 ELSE
308 k1 = k
309 k2 = k
310 knext = k - 1
311 END IF
312 END IF
313*
314 IF( l1.EQ.l2 .AND. k1.EQ.k2 ) THEN
315 suml = ddot( m-k1, a( k1, min( k1+1, m ) ), lda,
316 $ c( min( k1+1, m ), l1 ), 1 )
317 sumr = ddot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
318 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
319 scaloc = one
320*
321 a11 = a( k1, k1 ) + sgn*b( l1, l1 )
322 da11 = abs( a11 )
323 IF( da11.LE.smin ) THEN
324 a11 = smin
325 da11 = smin
326 info = 1
327 END IF
328 db = abs( vec( 1, 1 ) )
329 IF( da11.LT.one .AND. db.GT.one ) THEN
330 IF( db.GT.bignum*da11 )
331 $ scaloc = one / db
332 END IF
333 x( 1, 1 ) = ( vec( 1, 1 )*scaloc ) / a11
334*
335 IF( scaloc.NE.one ) THEN
336 DO 10 j = 1, n
337 CALL dscal( m, scaloc, c( 1, j ), 1 )
338 10 CONTINUE
339 scale = scale*scaloc
340 END IF
341 c( k1, l1 ) = x( 1, 1 )
342*
343 ELSE IF( l1.EQ.l2 .AND. k1.NE.k2 ) THEN
344*
345 suml = ddot( m-k2, a( k1, min( k2+1, m ) ), lda,
346 $ c( min( k2+1, m ), l1 ), 1 )
347 sumr = ddot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
348 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
349*
350 suml = ddot( m-k2, a( k2, min( k2+1, m ) ), lda,
351 $ c( min( k2+1, m ), l1 ), 1 )
352 sumr = ddot( l1-1, c( k2, 1 ), ldc, b( 1, l1 ), 1 )
353 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
354*
355 CALL dlaln2( .false., 2, 1, smin, one, a( k1, k1 ),
356 $ lda, one, one, vec, 2, -sgn*b( l1, l1 ),
357 $ zero, x, 2, scaloc, xnorm, ierr )
358 IF( ierr.NE.0 )
359 $ info = 1
360*
361 IF( scaloc.NE.one ) THEN
362 DO 20 j = 1, n
363 CALL dscal( m, scaloc, c( 1, j ), 1 )
364 20 CONTINUE
365 scale = scale*scaloc
366 END IF
367 c( k1, l1 ) = x( 1, 1 )
368 c( k2, l1 ) = x( 2, 1 )
369*
370 ELSE IF( l1.NE.l2 .AND. k1.EQ.k2 ) THEN
371*
372 suml = ddot( m-k1, a( k1, min( k1+1, m ) ), lda,
373 $ c( min( k1+1, m ), l1 ), 1 )
374 sumr = ddot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
375 vec( 1, 1 ) = sgn*( c( k1, l1 )-( suml+sgn*sumr ) )
376*
377 suml = ddot( m-k1, a( k1, min( k1+1, m ) ), lda,
378 $ c( min( k1+1, m ), l2 ), 1 )
379 sumr = ddot( l1-1, c( k1, 1 ), ldc, b( 1, l2 ), 1 )
380 vec( 2, 1 ) = sgn*( c( k1, l2 )-( suml+sgn*sumr ) )
381*
382 CALL dlaln2( .true., 2, 1, smin, one, b( l1, l1 ),
383 $ ldb, one, one, vec, 2, -sgn*a( k1, k1 ),
384 $ zero, x, 2, scaloc, xnorm, ierr )
385 IF( ierr.NE.0 )
386 $ info = 1
387*
388 IF( scaloc.NE.one ) THEN
389 DO 30 j = 1, n
390 CALL dscal( m, scaloc, c( 1, j ), 1 )
391 30 CONTINUE
392 scale = scale*scaloc
393 END IF
394 c( k1, l1 ) = x( 1, 1 )
395 c( k1, l2 ) = x( 2, 1 )
396*
397 ELSE IF( l1.NE.l2 .AND. k1.NE.k2 ) THEN
398*
399 suml = ddot( m-k2, a( k1, min( k2+1, m ) ), lda,
400 $ c( min( k2+1, m ), l1 ), 1 )
401 sumr = ddot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
402 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
403*
404 suml = ddot( m-k2, a( k1, min( k2+1, m ) ), lda,
405 $ c( min( k2+1, m ), l2 ), 1 )
406 sumr = ddot( l1-1, c( k1, 1 ), ldc, b( 1, l2 ), 1 )
407 vec( 1, 2 ) = c( k1, l2 ) - ( suml+sgn*sumr )
408*
409 suml = ddot( m-k2, a( k2, min( k2+1, m ) ), lda,
410 $ c( min( k2+1, m ), l1 ), 1 )
411 sumr = ddot( l1-1, c( k2, 1 ), ldc, b( 1, l1 ), 1 )
412 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
413*
414 suml = ddot( m-k2, a( k2, min( k2+1, m ) ), lda,
415 $ c( min( k2+1, m ), l2 ), 1 )
416 sumr = ddot( l1-1, c( k2, 1 ), ldc, b( 1, l2 ), 1 )
417 vec( 2, 2 ) = c( k2, l2 ) - ( suml+sgn*sumr )
418*
419 CALL dlasy2( .false., .false., isgn, 2, 2,
420 $ a( k1, k1 ), lda, b( l1, l1 ), ldb, vec,
421 $ 2, scaloc, x, 2, xnorm, ierr )
422 IF( ierr.NE.0 )
423 $ info = 1
424*
425 IF( scaloc.NE.one ) THEN
426 DO 40 j = 1, n
427 CALL dscal( m, scaloc, c( 1, j ), 1 )
428 40 CONTINUE
429 scale = scale*scaloc
430 END IF
431 c( k1, l1 ) = x( 1, 1 )
432 c( k1, l2 ) = x( 1, 2 )
433 c( k2, l1 ) = x( 2, 1 )
434 c( k2, l2 ) = x( 2, 2 )
435 END IF
436*
437 50 CONTINUE
438*
439 60 CONTINUE
440*
441 ELSE IF( .NOT.notrna .AND. notrnb ) THEN
442*
443* Solve A**T *X + ISGN*X*B = scale*C.
444*
445* The (K,L)th block of X is determined starting from
446* upper-left corner column by column by
447*
448* A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
449*
450* Where
451* K-1 T L-1
452* R(K,L) = SUM [A(I,K)**T*X(I,L)] +ISGN*SUM [X(K,J)*B(J,L)]
453* I=1 J=1
454*
455* Start column loop (index = L)
456* L1 (L2): column index of the first (last) row of X(K,L)
457*
458 lnext = 1
459 DO 120 l = 1, n
460 IF( l.LT.lnext )
461 $ GO TO 120
462 IF( l.EQ.n ) THEN
463 l1 = l
464 l2 = l
465 ELSE
466 IF( b( l+1, l ).NE.zero ) THEN
467 l1 = l
468 l2 = l + 1
469 lnext = l + 2
470 ELSE
471 l1 = l
472 l2 = l
473 lnext = l + 1
474 END IF
475 END IF
476*
477* Start row loop (index = K)
478* K1 (K2): row index of the first (last) row of X(K,L)
479*
480 knext = 1
481 DO 110 k = 1, m
482 IF( k.LT.knext )
483 $ GO TO 110
484 IF( k.EQ.m ) THEN
485 k1 = k
486 k2 = k
487 ELSE
488 IF( a( k+1, k ).NE.zero ) THEN
489 k1 = k
490 k2 = k + 1
491 knext = k + 2
492 ELSE
493 k1 = k
494 k2 = k
495 knext = k + 1
496 END IF
497 END IF
498*
499 IF( l1.EQ.l2 .AND. k1.EQ.k2 ) THEN
500 suml = ddot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
501 sumr = ddot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
502 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
503 scaloc = one
504*
505 a11 = a( k1, k1 ) + sgn*b( l1, l1 )
506 da11 = abs( a11 )
507 IF( da11.LE.smin ) THEN
508 a11 = smin
509 da11 = smin
510 info = 1
511 END IF
512 db = abs( vec( 1, 1 ) )
513 IF( da11.LT.one .AND. db.GT.one ) THEN
514 IF( db.GT.bignum*da11 )
515 $ scaloc = one / db
516 END IF
517 x( 1, 1 ) = ( vec( 1, 1 )*scaloc ) / a11
518*
519 IF( scaloc.NE.one ) THEN
520 DO 70 j = 1, n
521 CALL dscal( m, scaloc, c( 1, j ), 1 )
522 70 CONTINUE
523 scale = scale*scaloc
524 END IF
525 c( k1, l1 ) = x( 1, 1 )
526*
527 ELSE IF( l1.EQ.l2 .AND. k1.NE.k2 ) THEN
528*
529 suml = ddot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
530 sumr = ddot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
531 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
532*
533 suml = ddot( k1-1, a( 1, k2 ), 1, c( 1, l1 ), 1 )
534 sumr = ddot( l1-1, c( k2, 1 ), ldc, b( 1, l1 ), 1 )
535 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
536*
537 CALL dlaln2( .true., 2, 1, smin, one, a( k1, k1 ),
538 $ lda, one, one, vec, 2, -sgn*b( l1, l1 ),
539 $ zero, x, 2, scaloc, xnorm, ierr )
540 IF( ierr.NE.0 )
541 $ info = 1
542*
543 IF( scaloc.NE.one ) THEN
544 DO 80 j = 1, n
545 CALL dscal( m, scaloc, c( 1, j ), 1 )
546 80 CONTINUE
547 scale = scale*scaloc
548 END IF
549 c( k1, l1 ) = x( 1, 1 )
550 c( k2, l1 ) = x( 2, 1 )
551*
552 ELSE IF( l1.NE.l2 .AND. k1.EQ.k2 ) THEN
553*
554 suml = ddot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
555 sumr = ddot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
556 vec( 1, 1 ) = sgn*( c( k1, l1 )-( suml+sgn*sumr ) )
557*
558 suml = ddot( k1-1, a( 1, k1 ), 1, c( 1, l2 ), 1 )
559 sumr = ddot( l1-1, c( k1, 1 ), ldc, b( 1, l2 ), 1 )
560 vec( 2, 1 ) = sgn*( c( k1, l2 )-( suml+sgn*sumr ) )
561*
562 CALL dlaln2( .true., 2, 1, smin, one, b( l1, l1 ),
563 $ ldb, one, one, vec, 2, -sgn*a( k1, k1 ),
564 $ zero, x, 2, scaloc, xnorm, ierr )
565 IF( ierr.NE.0 )
566 $ info = 1
567*
568 IF( scaloc.NE.one ) THEN
569 DO 90 j = 1, n
570 CALL dscal( m, scaloc, c( 1, j ), 1 )
571 90 CONTINUE
572 scale = scale*scaloc
573 END IF
574 c( k1, l1 ) = x( 1, 1 )
575 c( k1, l2 ) = x( 2, 1 )
576*
577 ELSE IF( l1.NE.l2 .AND. k1.NE.k2 ) THEN
578*
579 suml = ddot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
580 sumr = ddot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
581 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
582*
583 suml = ddot( k1-1, a( 1, k1 ), 1, c( 1, l2 ), 1 )
584 sumr = ddot( l1-1, c( k1, 1 ), ldc, b( 1, l2 ), 1 )
585 vec( 1, 2 ) = c( k1, l2 ) - ( suml+sgn*sumr )
586*
587 suml = ddot( k1-1, a( 1, k2 ), 1, c( 1, l1 ), 1 )
588 sumr = ddot( l1-1, c( k2, 1 ), ldc, b( 1, l1 ), 1 )
589 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
590*
591 suml = ddot( k1-1, a( 1, k2 ), 1, c( 1, l2 ), 1 )
592 sumr = ddot( l1-1, c( k2, 1 ), ldc, b( 1, l2 ), 1 )
593 vec( 2, 2 ) = c( k2, l2 ) - ( suml+sgn*sumr )
594*
595 CALL dlasy2( .true., .false., isgn, 2, 2, a( k1, k1 ),
596 $ lda, b( l1, l1 ), ldb, vec, 2, scaloc, x,
597 $ 2, xnorm, ierr )
598 IF( ierr.NE.0 )
599 $ info = 1
600*
601 IF( scaloc.NE.one ) THEN
602 DO 100 j = 1, n
603 CALL dscal( m, scaloc, c( 1, j ), 1 )
604 100 CONTINUE
605 scale = scale*scaloc
606 END IF
607 c( k1, l1 ) = x( 1, 1 )
608 c( k1, l2 ) = x( 1, 2 )
609 c( k2, l1 ) = x( 2, 1 )
610 c( k2, l2 ) = x( 2, 2 )
611 END IF
612*
613 110 CONTINUE
614 120 CONTINUE
615*
616 ELSE IF( .NOT.notrna .AND. .NOT.notrnb ) THEN
617*
618* Solve A**T*X + ISGN*X*B**T = scale*C.
619*
620* The (K,L)th block of X is determined starting from
621* top-right corner column by column by
622*
623* A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L)**T = C(K,L) - R(K,L)
624*
625* Where
626* K-1 N
627* R(K,L) = SUM [A(I,K)**T*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**T].
628* I=1 J=L+1
629*
630* Start column loop (index = L)
631* L1 (L2): column index of the first (last) row of X(K,L)
632*
633 lnext = n
634 DO 180 l = n, 1, -1
635 IF( l.GT.lnext )
636 $ GO TO 180
637 IF( l.EQ.1 ) THEN
638 l1 = l
639 l2 = l
640 ELSE
641 IF( b( l, l-1 ).NE.zero ) THEN
642 l1 = l - 1
643 l2 = l
644 lnext = l - 2
645 ELSE
646 l1 = l
647 l2 = l
648 lnext = l - 1
649 END IF
650 END IF
651*
652* Start row loop (index = K)
653* K1 (K2): row index of the first (last) row of X(K,L)
654*
655 knext = 1
656 DO 170 k = 1, m
657 IF( k.LT.knext )
658 $ GO TO 170
659 IF( k.EQ.m ) THEN
660 k1 = k
661 k2 = k
662 ELSE
663 IF( a( k+1, k ).NE.zero ) THEN
664 k1 = k
665 k2 = k + 1
666 knext = k + 2
667 ELSE
668 k1 = k
669 k2 = k
670 knext = k + 1
671 END IF
672 END IF
673*
674 IF( l1.EQ.l2 .AND. k1.EQ.k2 ) THEN
675 suml = ddot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
676 sumr = ddot( n-l1, c( k1, min( l1+1, n ) ), ldc,
677 $ b( l1, min( l1+1, n ) ), ldb )
678 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
679 scaloc = one
680*
681 a11 = a( k1, k1 ) + sgn*b( l1, l1 )
682 da11 = abs( a11 )
683 IF( da11.LE.smin ) THEN
684 a11 = smin
685 da11 = smin
686 info = 1
687 END IF
688 db = abs( vec( 1, 1 ) )
689 IF( da11.LT.one .AND. db.GT.one ) THEN
690 IF( db.GT.bignum*da11 )
691 $ scaloc = one / db
692 END IF
693 x( 1, 1 ) = ( vec( 1, 1 )*scaloc ) / a11
694*
695 IF( scaloc.NE.one ) THEN
696 DO 130 j = 1, n
697 CALL dscal( m, scaloc, c( 1, j ), 1 )
698 130 CONTINUE
699 scale = scale*scaloc
700 END IF
701 c( k1, l1 ) = x( 1, 1 )
702*
703 ELSE IF( l1.EQ.l2 .AND. k1.NE.k2 ) THEN
704*
705 suml = ddot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
706 sumr = ddot( n-l2, c( k1, min( l2+1, n ) ), ldc,
707 $ b( l1, min( l2+1, n ) ), ldb )
708 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
709*
710 suml = ddot( k1-1, a( 1, k2 ), 1, c( 1, l1 ), 1 )
711 sumr = ddot( n-l2, c( k2, min( l2+1, n ) ), ldc,
712 $ b( l1, min( l2+1, n ) ), ldb )
713 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
714*
715 CALL dlaln2( .true., 2, 1, smin, one, a( k1, k1 ),
716 $ lda, one, one, vec, 2, -sgn*b( l1, l1 ),
717 $ zero, x, 2, scaloc, xnorm, ierr )
718 IF( ierr.NE.0 )
719 $ info = 1
720*
721 IF( scaloc.NE.one ) THEN
722 DO 140 j = 1, n
723 CALL dscal( m, scaloc, c( 1, j ), 1 )
724 140 CONTINUE
725 scale = scale*scaloc
726 END IF
727 c( k1, l1 ) = x( 1, 1 )
728 c( k2, l1 ) = x( 2, 1 )
729*
730 ELSE IF( l1.NE.l2 .AND. k1.EQ.k2 ) THEN
731*
732 suml = ddot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
733 sumr = ddot( n-l2, c( k1, min( l2+1, n ) ), ldc,
734 $ b( l1, min( l2+1, n ) ), ldb )
735 vec( 1, 1 ) = sgn*( c( k1, l1 )-( suml+sgn*sumr ) )
736*
737 suml = ddot( k1-1, a( 1, k1 ), 1, c( 1, l2 ), 1 )
738 sumr = ddot( n-l2, c( k1, min( l2+1, n ) ), ldc,
739 $ b( l2, min( l2+1, n ) ), ldb )
740 vec( 2, 1 ) = sgn*( c( k1, l2 )-( suml+sgn*sumr ) )
741*
742 CALL dlaln2( .false., 2, 1, smin, one, b( l1, l1 ),
743 $ ldb, one, one, vec, 2, -sgn*a( k1, k1 ),
744 $ zero, x, 2, scaloc, xnorm, ierr )
745 IF( ierr.NE.0 )
746 $ info = 1
747*
748 IF( scaloc.NE.one ) THEN
749 DO 150 j = 1, n
750 CALL dscal( m, scaloc, c( 1, j ), 1 )
751 150 CONTINUE
752 scale = scale*scaloc
753 END IF
754 c( k1, l1 ) = x( 1, 1 )
755 c( k1, l2 ) = x( 2, 1 )
756*
757 ELSE IF( l1.NE.l2 .AND. k1.NE.k2 ) THEN
758*
759 suml = ddot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
760 sumr = ddot( n-l2, c( k1, min( l2+1, n ) ), ldc,
761 $ b( l1, min( l2+1, n ) ), ldb )
762 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
763*
764 suml = ddot( k1-1, a( 1, k1 ), 1, c( 1, l2 ), 1 )
765 sumr = ddot( n-l2, c( k1, min( l2+1, n ) ), ldc,
766 $ b( l2, min( l2+1, n ) ), ldb )
767 vec( 1, 2 ) = c( k1, l2 ) - ( suml+sgn*sumr )
768*
769 suml = ddot( k1-1, a( 1, k2 ), 1, c( 1, l1 ), 1 )
770 sumr = ddot( n-l2, c( k2, min( l2+1, n ) ), ldc,
771 $ b( l1, min( l2+1, n ) ), ldb )
772 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
773*
774 suml = ddot( k1-1, a( 1, k2 ), 1, c( 1, l2 ), 1 )
775 sumr = ddot( n-l2, c( k2, min( l2+1, n ) ), ldc,
776 $ b( l2, min( l2+1, n ) ), ldb )
777 vec( 2, 2 ) = c( k2, l2 ) - ( suml+sgn*sumr )
778*
779 CALL dlasy2( .true., .true., isgn, 2, 2, a( k1, k1 ),
780 $ lda, b( l1, l1 ), ldb, vec, 2, scaloc, x,
781 $ 2, xnorm, ierr )
782 IF( ierr.NE.0 )
783 $ info = 1
784*
785 IF( scaloc.NE.one ) THEN
786 DO 160 j = 1, n
787 CALL dscal( m, scaloc, c( 1, j ), 1 )
788 160 CONTINUE
789 scale = scale*scaloc
790 END IF
791 c( k1, l1 ) = x( 1, 1 )
792 c( k1, l2 ) = x( 1, 2 )
793 c( k2, l1 ) = x( 2, 1 )
794 c( k2, l2 ) = x( 2, 2 )
795 END IF
796*
797 170 CONTINUE
798 180 CONTINUE
799*
800 ELSE IF( notrna .AND. .NOT.notrnb ) THEN
801*
802* Solve A*X + ISGN*X*B**T = scale*C.
803*
804* The (K,L)th block of X is determined starting from
805* bottom-right corner column by column by
806*
807* A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L)**T = C(K,L) - R(K,L)
808*
809* Where
810* M N
811* R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**T].
812* I=K+1 J=L+1
813*
814* Start column loop (index = L)
815* L1 (L2): column index of the first (last) row of X(K,L)
816*
817 lnext = n
818 DO 240 l = n, 1, -1
819 IF( l.GT.lnext )
820 $ GO TO 240
821 IF( l.EQ.1 ) THEN
822 l1 = l
823 l2 = l
824 ELSE
825 IF( b( l, l-1 ).NE.zero ) THEN
826 l1 = l - 1
827 l2 = l
828 lnext = l - 2
829 ELSE
830 l1 = l
831 l2 = l
832 lnext = l - 1
833 END IF
834 END IF
835*
836* Start row loop (index = K)
837* K1 (K2): row index of the first (last) row of X(K,L)
838*
839 knext = m
840 DO 230 k = m, 1, -1
841 IF( k.GT.knext )
842 $ GO TO 230
843 IF( k.EQ.1 ) THEN
844 k1 = k
845 k2 = k
846 ELSE
847 IF( a( k, k-1 ).NE.zero ) THEN
848 k1 = k - 1
849 k2 = k
850 knext = k - 2
851 ELSE
852 k1 = k
853 k2 = k
854 knext = k - 1
855 END IF
856 END IF
857*
858 IF( l1.EQ.l2 .AND. k1.EQ.k2 ) THEN
859 suml = ddot( m-k1, a( k1, min( k1+1, m ) ), lda,
860 $ c( min( k1+1, m ), l1 ), 1 )
861 sumr = ddot( n-l1, c( k1, min( l1+1, n ) ), ldc,
862 $ b( l1, min( l1+1, n ) ), ldb )
863 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
864 scaloc = one
865*
866 a11 = a( k1, k1 ) + sgn*b( l1, l1 )
867 da11 = abs( a11 )
868 IF( da11.LE.smin ) THEN
869 a11 = smin
870 da11 = smin
871 info = 1
872 END IF
873 db = abs( vec( 1, 1 ) )
874 IF( da11.LT.one .AND. db.GT.one ) THEN
875 IF( db.GT.bignum*da11 )
876 $ scaloc = one / db
877 END IF
878 x( 1, 1 ) = ( vec( 1, 1 )*scaloc ) / a11
879*
880 IF( scaloc.NE.one ) THEN
881 DO 190 j = 1, n
882 CALL dscal( m, scaloc, c( 1, j ), 1 )
883 190 CONTINUE
884 scale = scale*scaloc
885 END IF
886 c( k1, l1 ) = x( 1, 1 )
887*
888 ELSE IF( l1.EQ.l2 .AND. k1.NE.k2 ) THEN
889*
890 suml = ddot( m-k2, a( k1, min( k2+1, m ) ), lda,
891 $ c( min( k2+1, m ), l1 ), 1 )
892 sumr = ddot( n-l2, c( k1, min( l2+1, n ) ), ldc,
893 $ b( l1, min( l2+1, n ) ), ldb )
894 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
895*
896 suml = ddot( m-k2, a( k2, min( k2+1, m ) ), lda,
897 $ c( min( k2+1, m ), l1 ), 1 )
898 sumr = ddot( n-l2, c( k2, min( l2+1, n ) ), ldc,
899 $ b( l1, min( l2+1, n ) ), ldb )
900 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
901*
902 CALL dlaln2( .false., 2, 1, smin, one, a( k1, k1 ),
903 $ lda, one, one, vec, 2, -sgn*b( l1, l1 ),
904 $ zero, x, 2, scaloc, xnorm, ierr )
905 IF( ierr.NE.0 )
906 $ info = 1
907*
908 IF( scaloc.NE.one ) THEN
909 DO 200 j = 1, n
910 CALL dscal( m, scaloc, c( 1, j ), 1 )
911 200 CONTINUE
912 scale = scale*scaloc
913 END IF
914 c( k1, l1 ) = x( 1, 1 )
915 c( k2, l1 ) = x( 2, 1 )
916*
917 ELSE IF( l1.NE.l2 .AND. k1.EQ.k2 ) THEN
918*
919 suml = ddot( m-k1, a( k1, min( k1+1, m ) ), lda,
920 $ c( min( k1+1, m ), l1 ), 1 )
921 sumr = ddot( n-l2, c( k1, min( l2+1, n ) ), ldc,
922 $ b( l1, min( l2+1, n ) ), ldb )
923 vec( 1, 1 ) = sgn*( c( k1, l1 )-( suml+sgn*sumr ) )
924*
925 suml = ddot( m-k1, a( k1, min( k1+1, m ) ), lda,
926 $ c( min( k1+1, m ), l2 ), 1 )
927 sumr = ddot( n-l2, c( k1, min( l2+1, n ) ), ldc,
928 $ b( l2, min( l2+1, n ) ), ldb )
929 vec( 2, 1 ) = sgn*( c( k1, l2 )-( suml+sgn*sumr ) )
930*
931 CALL dlaln2( .false., 2, 1, smin, one, b( l1, l1 ),
932 $ ldb, one, one, vec, 2, -sgn*a( k1, k1 ),
933 $ zero, x, 2, scaloc, xnorm, ierr )
934 IF( ierr.NE.0 )
935 $ info = 1
936*
937 IF( scaloc.NE.one ) THEN
938 DO 210 j = 1, n
939 CALL dscal( m, scaloc, c( 1, j ), 1 )
940 210 CONTINUE
941 scale = scale*scaloc
942 END IF
943 c( k1, l1 ) = x( 1, 1 )
944 c( k1, l2 ) = x( 2, 1 )
945*
946 ELSE IF( l1.NE.l2 .AND. k1.NE.k2 ) THEN
947*
948 suml = ddot( m-k2, a( k1, min( k2+1, m ) ), lda,
949 $ c( min( k2+1, m ), l1 ), 1 )
950 sumr = ddot( n-l2, c( k1, min( l2+1, n ) ), ldc,
951 $ b( l1, min( l2+1, n ) ), ldb )
952 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
953*
954 suml = ddot( m-k2, a( k1, min( k2+1, m ) ), lda,
955 $ c( min( k2+1, m ), l2 ), 1 )
956 sumr = ddot( n-l2, c( k1, min( l2+1, n ) ), ldc,
957 $ b( l2, min( l2+1, n ) ), ldb )
958 vec( 1, 2 ) = c( k1, l2 ) - ( suml+sgn*sumr )
959*
960 suml = ddot( m-k2, a( k2, min( k2+1, m ) ), lda,
961 $ c( min( k2+1, m ), l1 ), 1 )
962 sumr = ddot( n-l2, c( k2, min( l2+1, n ) ), ldc,
963 $ b( l1, min( l2+1, n ) ), ldb )
964 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
965*
966 suml = ddot( m-k2, a( k2, min( k2+1, m ) ), lda,
967 $ c( min( k2+1, m ), l2 ), 1 )
968 sumr = ddot( n-l2, c( k2, min( l2+1, n ) ), ldc,
969 $ b( l2, min( l2+1, n ) ), ldb )
970 vec( 2, 2 ) = c( k2, l2 ) - ( suml+sgn*sumr )
971*
972 CALL dlasy2( .false., .true., isgn, 2, 2, a( k1, k1 ),
973 $ lda, b( l1, l1 ), ldb, vec, 2, scaloc, x,
974 $ 2, xnorm, ierr )
975 IF( ierr.NE.0 )
976 $ info = 1
977*
978 IF( scaloc.NE.one ) THEN
979 DO 220 j = 1, n
980 CALL dscal( m, scaloc, c( 1, j ), 1 )
981 220 CONTINUE
982 scale = scale*scaloc
983 END IF
984 c( k1, l1 ) = x( 1, 1 )
985 c( k1, l2 ) = x( 1, 2 )
986 c( k2, l1 ) = x( 2, 1 )
987 c( k2, l2 ) = x( 2, 2 )
988 END IF
989*
990 230 CONTINUE
991 240 CONTINUE
992*
993 END IF
994*
995 RETURN
996*
997* End of DTRSYL
998*
999 END
dlabad
subroutine dlabad(small, large)
DLABAD
Definition dlabad.f:74
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
dlaln2
subroutine dlaln2(ltrans, na, nw, smin, ca, a, lda, d1, d2, b, ldb, wr, wi, x, ldx, scale, xnorm, info)
DLALN2 solves a 1-by-1 or 2-by-2 linear system of equations of the specified form.
Definition dlaln2.f:218
dlasy2
subroutine dlasy2(ltranl, ltranr, isgn, n1, n2, tl, ldtl, tr, ldtr, b, ldb, scale, x, ldx, xnorm, info)
DLASY2 solves the Sylvester matrix equation where the matrices are of order 1 or 2.
Definition dlasy2.f:174
dtrsyl
subroutine dtrsyl(trana, tranb, isgn, m, n, a, lda, b, ldb, c, ldc, scale, info)
DTRSYL
Definition dtrsyl.f:164
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21