LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dtgsy2.f
Go to the documentation of this file.
1*> \brief \b DTGSY2 solves the generalized Sylvester equation (unblocked algorithm).
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download DTGSY2 + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtgsy2.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtgsy2.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtgsy2.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE DTGSY2( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
20* LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL,
21* IWORK, PQ, INFO )
22*
23* .. Scalar Arguments ..
24* CHARACTER TRANS
25* INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N,
26* $ PQ
27* DOUBLE PRECISION RDSCAL, RDSUM, SCALE
28* ..
29* .. Array Arguments ..
30* INTEGER IWORK( * )
31* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ),
32* $ D( LDD, * ), E( LDE, * ), F( LDF, * )
33* ..
34*
35*
36*> \par Purpose:
37* =============
38*>
39*> \verbatim
40*>
41*> DTGSY2 solves the generalized Sylvester equation:
42*>
43*> A * R - L * B = scale * C (1)
44*> D * R - L * E = scale * F,
45*>
46*> using Level 1 and 2 BLAS. where R and L are unknown M-by-N matrices,
47*> (A, D), (B, E) and (C, F) are given matrix pairs of size M-by-M,
48*> N-by-N and M-by-N, respectively, with real entries. (A, D) and (B, E)
49*> must be in generalized Schur canonical form, i.e. A, B are upper
50*> quasi triangular and D, E are upper triangular. The solution (R, L)
51*> overwrites (C, F). 0 <= SCALE <= 1 is an output scaling factor
52*> chosen to avoid overflow.
53*>
54*> In matrix notation solving equation (1) corresponds to solve
55*> Z*x = scale*b, where Z is defined as
56*>
57*> Z = [ kron(In, A) -kron(B**T, Im) ] (2)
58*> [ kron(In, D) -kron(E**T, Im) ],
59*>
60*> Ik is the identity matrix of size k and X**T is the transpose of X.
61*> kron(X, Y) is the Kronecker product between the matrices X and Y.
62*> In the process of solving (1), we solve a number of such systems
63*> where Dim(In), Dim(In) = 1 or 2.
64*>
65*> If TRANS = 'T', solve the transposed system Z**T*y = scale*b for y,
66*> which is equivalent to solve for R and L in
67*>
68*> A**T * R + D**T * L = scale * C (3)
69*> R * B**T + L * E**T = scale * -F
70*>
71*> This case is used to compute an estimate of Dif[(A, D), (B, E)] =
72*> sigma_min(Z) using reverse communication with DLACON.
73*>
74*> DTGSY2 also (IJOB >= 1) contributes to the computation in DTGSYL
75*> of an upper bound on the separation between to matrix pairs. Then
76*> the input (A, D), (B, E) are sub-pencils of the matrix pair in
77*> DTGSYL. See DTGSYL for details.
78*> \endverbatim
79*
80* Arguments:
81* ==========
82*
83*> \param[in] TRANS
84*> \verbatim
85*> TRANS is CHARACTER*1
86*> = 'N': solve the generalized Sylvester equation (1).
87*> = 'T': solve the 'transposed' system (3).
88*> \endverbatim
89*>
90*> \param[in] IJOB
91*> \verbatim
92*> IJOB is INTEGER
93*> Specifies what kind of functionality to be performed.
94*> = 0: solve (1) only.
95*> = 1: A contribution from this subsystem to a Frobenius
96*> norm-based estimate of the separation between two matrix
97*> pairs is computed. (look ahead strategy is used).
98*> = 2: A contribution from this subsystem to a Frobenius
99*> norm-based estimate of the separation between two matrix
100*> pairs is computed. (DGECON on sub-systems is used.)
101*> Not referenced if TRANS = 'T'.
102*> \endverbatim
103*>
104*> \param[in] M
105*> \verbatim
106*> M is INTEGER
107*> On entry, M specifies the order of A and D, and the row
108*> dimension of C, F, R and L.
109*> \endverbatim
110*>
111*> \param[in] N
112*> \verbatim
113*> N is INTEGER
114*> On entry, N specifies the order of B and E, and the column
115*> dimension of C, F, R and L.
116*> \endverbatim
117*>
118*> \param[in] A
119*> \verbatim
120*> A is DOUBLE PRECISION array, dimension (LDA, M)
121*> On entry, A contains an upper quasi triangular matrix.
122*> \endverbatim
123*>
124*> \param[in] LDA
125*> \verbatim
126*> LDA is INTEGER
127*> The leading dimension of the matrix A. LDA >= max(1, M).
128*> \endverbatim
129*>
130*> \param[in] B
131*> \verbatim
132*> B is DOUBLE PRECISION array, dimension (LDB, N)
133*> On entry, B contains an upper quasi triangular matrix.
134*> \endverbatim
135*>
136*> \param[in] LDB
137*> \verbatim
138*> LDB is INTEGER
139*> The leading dimension of the matrix B. LDB >= max(1, N).
140*> \endverbatim
141*>
142*> \param[in,out] C
143*> \verbatim
144*> C is DOUBLE PRECISION array, dimension (LDC, N)
145*> On entry, C contains the right-hand-side of the first matrix
146*> equation in (1).
147*> On exit, if IJOB = 0, C has been overwritten by the
148*> solution R.
149*> \endverbatim
150*>
151*> \param[in] LDC
152*> \verbatim
153*> LDC is INTEGER
154*> The leading dimension of the matrix C. LDC >= max(1, M).
155*> \endverbatim
156*>
157*> \param[in] D
158*> \verbatim
159*> D is DOUBLE PRECISION array, dimension (LDD, M)
160*> On entry, D contains an upper triangular matrix.
161*> \endverbatim
162*>
163*> \param[in] LDD
164*> \verbatim
165*> LDD is INTEGER
166*> The leading dimension of the matrix D. LDD >= max(1, M).
167*> \endverbatim
168*>
169*> \param[in] E
170*> \verbatim
171*> E is DOUBLE PRECISION array, dimension (LDE, N)
172*> On entry, E contains an upper triangular matrix.
173*> \endverbatim
174*>
175*> \param[in] LDE
176*> \verbatim
177*> LDE is INTEGER
178*> The leading dimension of the matrix E. LDE >= max(1, N).
179*> \endverbatim
180*>
181*> \param[in,out] F
182*> \verbatim
183*> F is DOUBLE PRECISION array, dimension (LDF, N)
184*> On entry, F contains the right-hand-side of the second matrix
185*> equation in (1).
186*> On exit, if IJOB = 0, F has been overwritten by the
187*> solution L.
188*> \endverbatim
189*>
190*> \param[in] LDF
191*> \verbatim
192*> LDF is INTEGER
193*> The leading dimension of the matrix F. LDF >= max(1, M).
194*> \endverbatim
195*>
196*> \param[out] SCALE
197*> \verbatim
198*> SCALE is DOUBLE PRECISION
199*> On exit, 0 <= SCALE <= 1. If 0 < SCALE < 1, the solutions
200*> R and L (C and F on entry) will hold the solutions to a
201*> slightly perturbed system but the input matrices A, B, D and
202*> E have not been changed. If SCALE = 0, R and L will hold the
203*> solutions to the homogeneous system with C = F = 0. Normally,
204*> SCALE = 1.
205*> \endverbatim
206*>
207*> \param[in,out] RDSUM
208*> \verbatim
209*> RDSUM is DOUBLE PRECISION
210*> On entry, the sum of squares of computed contributions to
211*> the Dif-estimate under computation by DTGSYL, where the
212*> scaling factor RDSCAL (see below) has been factored out.
213*> On exit, the corresponding sum of squares updated with the
214*> contributions from the current sub-system.
215*> If TRANS = 'T' RDSUM is not touched.
216*> NOTE: RDSUM only makes sense when DTGSY2 is called by DTGSYL.
217*> \endverbatim
218*>
219*> \param[in,out] RDSCAL
220*> \verbatim
221*> RDSCAL is DOUBLE PRECISION
222*> On entry, scaling factor used to prevent overflow in RDSUM.
223*> On exit, RDSCAL is updated w.r.t. the current contributions
224*> in RDSUM.
225*> If TRANS = 'T', RDSCAL is not touched.
226*> NOTE: RDSCAL only makes sense when DTGSY2 is called by
227*> DTGSYL.
228*> \endverbatim
229*>
230*> \param[out] IWORK
231*> \verbatim
232*> IWORK is INTEGER array, dimension (M+N+2)
233*> \endverbatim
234*>
235*> \param[out] PQ
236*> \verbatim
237*> PQ is INTEGER
238*> On exit, the number of subsystems (of size 2-by-2, 4-by-4 and
239*> 8-by-8) solved by this routine.
240*> \endverbatim
241*>
242*> \param[out] INFO
243*> \verbatim
244*> INFO is INTEGER
245*> On exit, if INFO is set to
246*> =0: Successful exit
247*> <0: If INFO = -i, the i-th argument had an illegal value.
248*> >0: The matrix pairs (A, D) and (B, E) have common or very
249*> close eigenvalues.
250*> \endverbatim
251*
252* Authors:
253* ========
254*
255*> \author Univ. of Tennessee
256*> \author Univ. of California Berkeley
257*> \author Univ. of Colorado Denver
258*> \author NAG Ltd.
259*
260*> \ingroup tgsy2
261*
262*> \par Contributors:
263* ==================
264*>
265*> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
266*> Umea University, S-901 87 Umea, Sweden.
267*
268* =====================================================================
269 SUBROUTINE dtgsy2( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC,
270 $ D,
271 $ LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL,
272 $ IWORK, PQ, INFO )
273*
274* -- LAPACK auxiliary routine --
275* -- LAPACK is a software package provided by Univ. of Tennessee, --
276* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
277*
278* .. Scalar Arguments ..
279 CHARACTER TRANS
280 INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N,
281 $ PQ
282 DOUBLE PRECISION RDSCAL, RDSUM, SCALE
283* ..
284* .. Array Arguments ..
285 INTEGER IWORK( * )
286 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ),
287 $ D( LDD, * ), E( LDE, * ), F( LDF, * )
288* ..
289*
290* =====================================================================
291* Replaced various illegal calls to DCOPY by calls to DLASET.
292* Sven Hammarling, 27/5/02.
293*
294* .. Parameters ..
295 INTEGER LDZ
296 PARAMETER ( LDZ = 8 )
297 double precision zero, one
298 parameter( zero = 0.0d+0, one = 1.0d+0 )
299* ..
300* .. Local Scalars ..
301 LOGICAL NOTRAN
302 INTEGER I, IE, IERR, II, IS, ISP1, J, JE, JJ, JS, JSP1,
303 $ K, MB, NB, P, Q, ZDIM
304 DOUBLE PRECISION ALPHA, SCALOC
305* ..
306* .. Local Arrays ..
307 INTEGER IPIV( LDZ ), JPIV( LDZ )
308 DOUBLE PRECISION RHS( LDZ ), Z( LDZ, LDZ )
309* ..
310* .. External Functions ..
311 LOGICAL LSAME
312 EXTERNAL LSAME
313* ..
314* .. External Subroutines ..
315 EXTERNAL daxpy, dcopy, dgemm, dgemv, dger,
316 $ dgesc2,
317 $ dgetc2, dlaset, dlatdf, dscal, xerbla
318* ..
319* .. Intrinsic Functions ..
320 INTRINSIC max
321* ..
322* .. Executable Statements ..
323*
324* Decode and test input parameters
325*
326 info = 0
327 ierr = 0
328 notran = lsame( trans, 'N' )
329 IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) ) THEN
330 info = -1
331 ELSE IF( notran ) THEN
332 IF( ( ijob.LT.0 ) .OR. ( ijob.GT.2 ) ) THEN
333 info = -2
334 END IF
335 END IF
336 IF( info.EQ.0 ) THEN
337 IF( m.LE.0 ) THEN
338 info = -3
339 ELSE IF( n.LE.0 ) THEN
340 info = -4
341 ELSE IF( lda.LT.max( 1, m ) ) THEN
342 info = -6
343 ELSE IF( ldb.LT.max( 1, n ) ) THEN
344 info = -8
345 ELSE IF( ldc.LT.max( 1, m ) ) THEN
346 info = -10
347 ELSE IF( ldd.LT.max( 1, m ) ) THEN
348 info = -12
349 ELSE IF( lde.LT.max( 1, n ) ) THEN
350 info = -14
351 ELSE IF( ldf.LT.max( 1, m ) ) THEN
352 info = -16
353 END IF
354 END IF
355 IF( info.NE.0 ) THEN
356 CALL xerbla( 'DTGSY2', -info )
357 RETURN
358 END IF
359*
360* Determine block structure of A
361*
362 pq = 0
363 p = 0
364 i = 1
365 10 CONTINUE
366 IF( i.GT.m )
367 $ GO TO 20
368 p = p + 1
369 iwork( p ) = i
370 IF( i.EQ.m )
371 $ GO TO 20
372 IF( a( i+1, i ).NE.zero ) THEN
373 i = i + 2
374 ELSE
375 i = i + 1
376 END IF
377 GO TO 10
378 20 CONTINUE
379 iwork( p+1 ) = m + 1
380*
381* Determine block structure of B
382*
383 q = p + 1
384 j = 1
385 30 CONTINUE
386 IF( j.GT.n )
387 $ GO TO 40
388 q = q + 1
389 iwork( q ) = j
390 IF( j.EQ.n )
391 $ GO TO 40
392 IF( b( j+1, j ).NE.zero ) THEN
393 j = j + 2
394 ELSE
395 j = j + 1
396 END IF
397 GO TO 30
398 40 CONTINUE
399 iwork( q+1 ) = n + 1
400 pq = p*( q-p-1 )
401*
402 IF( notran ) THEN
403*
404* Solve (I, J) - subsystem
405* A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
406* D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
407* for I = P, P - 1, ..., 1; J = 1, 2, ..., Q
408*
409 scale = one
410 scaloc = one
411 DO 120 j = p + 2, q
412 js = iwork( j )
413 jsp1 = js + 1
414 je = iwork( j+1 ) - 1
415 nb = je - js + 1
416 DO 110 i = p, 1, -1
417*
418 is = iwork( i )
419 isp1 = is + 1
420 ie = iwork( i+1 ) - 1
421 mb = ie - is + 1
422 zdim = mb*nb*2
423*
424 IF( ( mb.EQ.1 ) .AND. ( nb.EQ.1 ) ) THEN
425*
426* Build a 2-by-2 system Z * x = RHS
427*
428 z( 1, 1 ) = a( is, is )
429 z( 2, 1 ) = d( is, is )
430 z( 1, 2 ) = -b( js, js )
431 z( 2, 2 ) = -e( js, js )
432*
433* Set up right hand side(s)
434*
435 rhs( 1 ) = c( is, js )
436 rhs( 2 ) = f( is, js )
437*
438* Solve Z * x = RHS
439*
440 CALL dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
441 IF( ierr.GT.0 )
442 $ info = ierr
443*
444 IF( ijob.EQ.0 ) THEN
445 CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
446 $ scaloc )
447 IF( scaloc.NE.one ) THEN
448 DO 50 k = 1, n
449 CALL dscal( m, scaloc, c( 1, k ), 1 )
450 CALL dscal( m, scaloc, f( 1, k ), 1 )
451 50 CONTINUE
452 scale = scale*scaloc
453 END IF
454 ELSE
455 CALL dlatdf( ijob, zdim, z, ldz, rhs, rdsum,
456 $ rdscal, ipiv, jpiv )
457 END IF
458*
459* Unpack solution vector(s)
460*
461 c( is, js ) = rhs( 1 )
462 f( is, js ) = rhs( 2 )
463*
464* Substitute R(I, J) and L(I, J) into remaining
465* equation.
466*
467 IF( i.GT.1 ) THEN
468 alpha = -rhs( 1 )
469 CALL daxpy( is-1, alpha, a( 1, is ), 1, c( 1,
470 $ js ),
471 $ 1 )
472 CALL daxpy( is-1, alpha, d( 1, is ), 1, f( 1,
473 $ js ),
474 $ 1 )
475 END IF
476 IF( j.LT.q ) THEN
477 CALL daxpy( n-je, rhs( 2 ), b( js, je+1 ), ldb,
478 $ c( is, je+1 ), ldc )
479 CALL daxpy( n-je, rhs( 2 ), e( js, je+1 ), lde,
480 $ f( is, je+1 ), ldf )
481 END IF
482*
483 ELSE IF( ( mb.EQ.1 ) .AND. ( nb.EQ.2 ) ) THEN
484*
485* Build a 4-by-4 system Z * x = RHS
486*
487 z( 1, 1 ) = a( is, is )
488 z( 2, 1 ) = zero
489 z( 3, 1 ) = d( is, is )
490 z( 4, 1 ) = zero
491*
492 z( 1, 2 ) = zero
493 z( 2, 2 ) = a( is, is )
494 z( 3, 2 ) = zero
495 z( 4, 2 ) = d( is, is )
496*
497 z( 1, 3 ) = -b( js, js )
498 z( 2, 3 ) = -b( js, jsp1 )
499 z( 3, 3 ) = -e( js, js )
500 z( 4, 3 ) = -e( js, jsp1 )
501*
502 z( 1, 4 ) = -b( jsp1, js )
503 z( 2, 4 ) = -b( jsp1, jsp1 )
504 z( 3, 4 ) = zero
505 z( 4, 4 ) = -e( jsp1, jsp1 )
506*
507* Set up right hand side(s)
508*
509 rhs( 1 ) = c( is, js )
510 rhs( 2 ) = c( is, jsp1 )
511 rhs( 3 ) = f( is, js )
512 rhs( 4 ) = f( is, jsp1 )
513*
514* Solve Z * x = RHS
515*
516 CALL dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
517 IF( ierr.GT.0 )
518 $ info = ierr
519*
520 IF( ijob.EQ.0 ) THEN
521 CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
522 $ scaloc )
523 IF( scaloc.NE.one ) THEN
524 DO 60 k = 1, n
525 CALL dscal( m, scaloc, c( 1, k ), 1 )
526 CALL dscal( m, scaloc, f( 1, k ), 1 )
527 60 CONTINUE
528 scale = scale*scaloc
529 END IF
530 ELSE
531 CALL dlatdf( ijob, zdim, z, ldz, rhs, rdsum,
532 $ rdscal, ipiv, jpiv )
533 END IF
534*
535* Unpack solution vector(s)
536*
537 c( is, js ) = rhs( 1 )
538 c( is, jsp1 ) = rhs( 2 )
539 f( is, js ) = rhs( 3 )
540 f( is, jsp1 ) = rhs( 4 )
541*
542* Substitute R(I, J) and L(I, J) into remaining
543* equation.
544*
545 IF( i.GT.1 ) THEN
546 CALL dger( is-1, nb, -one, a( 1, is ), 1,
547 $ rhs( 1 ),
548 $ 1, c( 1, js ), ldc )
549 CALL dger( is-1, nb, -one, d( 1, is ), 1,
550 $ rhs( 1 ),
551 $ 1, f( 1, js ), ldf )
552 END IF
553 IF( j.LT.q ) THEN
554 CALL daxpy( n-je, rhs( 3 ), b( js, je+1 ), ldb,
555 $ c( is, je+1 ), ldc )
556 CALL daxpy( n-je, rhs( 3 ), e( js, je+1 ), lde,
557 $ f( is, je+1 ), ldf )
558 CALL daxpy( n-je, rhs( 4 ), b( jsp1, je+1 ),
559 $ ldb,
560 $ c( is, je+1 ), ldc )
561 CALL daxpy( n-je, rhs( 4 ), e( jsp1, je+1 ),
562 $ lde,
563 $ f( is, je+1 ), ldf )
564 END IF
565*
566 ELSE IF( ( mb.EQ.2 ) .AND. ( nb.EQ.1 ) ) THEN
567*
568* Build a 4-by-4 system Z * x = RHS
569*
570 z( 1, 1 ) = a( is, is )
571 z( 2, 1 ) = a( isp1, is )
572 z( 3, 1 ) = d( is, is )
573 z( 4, 1 ) = zero
574*
575 z( 1, 2 ) = a( is, isp1 )
576 z( 2, 2 ) = a( isp1, isp1 )
577 z( 3, 2 ) = d( is, isp1 )
578 z( 4, 2 ) = d( isp1, isp1 )
579*
580 z( 1, 3 ) = -b( js, js )
581 z( 2, 3 ) = zero
582 z( 3, 3 ) = -e( js, js )
583 z( 4, 3 ) = zero
584*
585 z( 1, 4 ) = zero
586 z( 2, 4 ) = -b( js, js )
587 z( 3, 4 ) = zero
588 z( 4, 4 ) = -e( js, js )
589*
590* Set up right hand side(s)
591*
592 rhs( 1 ) = c( is, js )
593 rhs( 2 ) = c( isp1, js )
594 rhs( 3 ) = f( is, js )
595 rhs( 4 ) = f( isp1, js )
596*
597* Solve Z * x = RHS
598*
599 CALL dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
600 IF( ierr.GT.0 )
601 $ info = ierr
602 IF( ijob.EQ.0 ) THEN
603 CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
604 $ scaloc )
605 IF( scaloc.NE.one ) THEN
606 DO 70 k = 1, n
607 CALL dscal( m, scaloc, c( 1, k ), 1 )
608 CALL dscal( m, scaloc, f( 1, k ), 1 )
609 70 CONTINUE
610 scale = scale*scaloc
611 END IF
612 ELSE
613 CALL dlatdf( ijob, zdim, z, ldz, rhs, rdsum,
614 $ rdscal, ipiv, jpiv )
615 END IF
616*
617* Unpack solution vector(s)
618*
619 c( is, js ) = rhs( 1 )
620 c( isp1, js ) = rhs( 2 )
621 f( is, js ) = rhs( 3 )
622 f( isp1, js ) = rhs( 4 )
623*
624* Substitute R(I, J) and L(I, J) into remaining
625* equation.
626*
627 IF( i.GT.1 ) THEN
628 CALL dgemv( 'N', is-1, mb, -one, a( 1, is ),
629 $ lda,
630 $ rhs( 1 ), 1, one, c( 1, js ), 1 )
631 CALL dgemv( 'N', is-1, mb, -one, d( 1, is ),
632 $ ldd,
633 $ rhs( 1 ), 1, one, f( 1, js ), 1 )
634 END IF
635 IF( j.LT.q ) THEN
636 CALL dger( mb, n-je, one, rhs( 3 ), 1,
637 $ b( js, je+1 ), ldb, c( is, je+1 ), ldc )
638 CALL dger( mb, n-je, one, rhs( 3 ), 1,
639 $ e( js, je+1 ), lde, f( is, je+1 ), ldf )
640 END IF
641*
642 ELSE IF( ( mb.EQ.2 ) .AND. ( nb.EQ.2 ) ) THEN
643*
644* Build an 8-by-8 system Z * x = RHS
645*
646 CALL dlaset( 'F', ldz, ldz, zero, zero, z, ldz )
647*
648 z( 1, 1 ) = a( is, is )
649 z( 2, 1 ) = a( isp1, is )
650 z( 5, 1 ) = d( is, is )
651*
652 z( 1, 2 ) = a( is, isp1 )
653 z( 2, 2 ) = a( isp1, isp1 )
654 z( 5, 2 ) = d( is, isp1 )
655 z( 6, 2 ) = d( isp1, isp1 )
656*
657 z( 3, 3 ) = a( is, is )
658 z( 4, 3 ) = a( isp1, is )
659 z( 7, 3 ) = d( is, is )
660*
661 z( 3, 4 ) = a( is, isp1 )
662 z( 4, 4 ) = a( isp1, isp1 )
663 z( 7, 4 ) = d( is, isp1 )
664 z( 8, 4 ) = d( isp1, isp1 )
665*
666 z( 1, 5 ) = -b( js, js )
667 z( 3, 5 ) = -b( js, jsp1 )
668 z( 5, 5 ) = -e( js, js )
669 z( 7, 5 ) = -e( js, jsp1 )
670*
671 z( 2, 6 ) = -b( js, js )
672 z( 4, 6 ) = -b( js, jsp1 )
673 z( 6, 6 ) = -e( js, js )
674 z( 8, 6 ) = -e( js, jsp1 )
675*
676 z( 1, 7 ) = -b( jsp1, js )
677 z( 3, 7 ) = -b( jsp1, jsp1 )
678 z( 7, 7 ) = -e( jsp1, jsp1 )
679*
680 z( 2, 8 ) = -b( jsp1, js )
681 z( 4, 8 ) = -b( jsp1, jsp1 )
682 z( 8, 8 ) = -e( jsp1, jsp1 )
683*
684* Set up right hand side(s)
685*
686 k = 1
687 ii = mb*nb + 1
688 DO 80 jj = 0, nb - 1
689 CALL dcopy( mb, c( is, js+jj ), 1, rhs( k ), 1 )
690 CALL dcopy( mb, f( is, js+jj ), 1, rhs( ii ),
691 $ 1 )
692 k = k + mb
693 ii = ii + mb
694 80 CONTINUE
695*
696* Solve Z * x = RHS
697*
698 CALL dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
699 IF( ierr.GT.0 )
700 $ info = ierr
701 IF( ijob.EQ.0 ) THEN
702 CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
703 $ scaloc )
704 IF( scaloc.NE.one ) THEN
705 DO 90 k = 1, n
706 CALL dscal( m, scaloc, c( 1, k ), 1 )
707 CALL dscal( m, scaloc, f( 1, k ), 1 )
708 90 CONTINUE
709 scale = scale*scaloc
710 END IF
711 ELSE
712 CALL dlatdf( ijob, zdim, z, ldz, rhs, rdsum,
713 $ rdscal, ipiv, jpiv )
714 END IF
715*
716* Unpack solution vector(s)
717*
718 k = 1
719 ii = mb*nb + 1
720 DO 100 jj = 0, nb - 1
721 CALL dcopy( mb, rhs( k ), 1, c( is, js+jj ), 1 )
722 CALL dcopy( mb, rhs( ii ), 1, f( is, js+jj ),
723 $ 1 )
724 k = k + mb
725 ii = ii + mb
726 100 CONTINUE
727*
728* Substitute R(I, J) and L(I, J) into remaining
729* equation.
730*
731 IF( i.GT.1 ) THEN
732 CALL dgemm( 'N', 'N', is-1, nb, mb, -one,
733 $ a( 1, is ), lda, rhs( 1 ), mb, one,
734 $ c( 1, js ), ldc )
735 CALL dgemm( 'N', 'N', is-1, nb, mb, -one,
736 $ d( 1, is ), ldd, rhs( 1 ), mb, one,
737 $ f( 1, js ), ldf )
738 END IF
739 IF( j.LT.q ) THEN
740 k = mb*nb + 1
741 CALL dgemm( 'N', 'N', mb, n-je, nb, one,
742 $ rhs( k ),
743 $ mb, b( js, je+1 ), ldb, one,
744 $ c( is, je+1 ), ldc )
745 CALL dgemm( 'N', 'N', mb, n-je, nb, one,
746 $ rhs( k ),
747 $ mb, e( js, je+1 ), lde, one,
748 $ f( is, je+1 ), ldf )
749 END IF
750*
751 END IF
752*
753 110 CONTINUE
754 120 CONTINUE
755 ELSE
756*
757* Solve (I, J) - subsystem
758* A(I, I)**T * R(I, J) + D(I, I)**T * L(J, J) = C(I, J)
759* R(I, I) * B(J, J) + L(I, J) * E(J, J) = -F(I, J)
760* for I = 1, 2, ..., P, J = Q, Q - 1, ..., 1
761*
762 scale = one
763 scaloc = one
764 DO 200 i = 1, p
765*
766 is = iwork( i )
767 isp1 = is + 1
768 ie = iwork( i+1 ) - 1
769 mb = ie - is + 1
770 DO 190 j = q, p + 2, -1
771*
772 js = iwork( j )
773 jsp1 = js + 1
774 je = iwork( j+1 ) - 1
775 nb = je - js + 1
776 zdim = mb*nb*2
777 IF( ( mb.EQ.1 ) .AND. ( nb.EQ.1 ) ) THEN
778*
779* Build a 2-by-2 system Z**T * x = RHS
780*
781 z( 1, 1 ) = a( is, is )
782 z( 2, 1 ) = -b( js, js )
783 z( 1, 2 ) = d( is, is )
784 z( 2, 2 ) = -e( js, js )
785*
786* Set up right hand side(s)
787*
788 rhs( 1 ) = c( is, js )
789 rhs( 2 ) = f( is, js )
790*
791* Solve Z**T * x = RHS
792*
793 CALL dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
794 IF( ierr.GT.0 )
795 $ info = ierr
796*
797 CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
798 $ scaloc )
799 IF( scaloc.NE.one ) THEN
800 DO 130 k = 1, n
801 CALL dscal( m, scaloc, c( 1, k ), 1 )
802 CALL dscal( m, scaloc, f( 1, k ), 1 )
803 130 CONTINUE
804 scale = scale*scaloc
805 END IF
806*
807* Unpack solution vector(s)
808*
809 c( is, js ) = rhs( 1 )
810 f( is, js ) = rhs( 2 )
811*
812* Substitute R(I, J) and L(I, J) into remaining
813* equation.
814*
815 IF( j.GT.p+2 ) THEN
816 alpha = rhs( 1 )
817 CALL daxpy( js-1, alpha, b( 1, js ), 1, f( is,
818 $ 1 ),
819 $ ldf )
820 alpha = rhs( 2 )
821 CALL daxpy( js-1, alpha, e( 1, js ), 1, f( is,
822 $ 1 ),
823 $ ldf )
824 END IF
825 IF( i.LT.p ) THEN
826 alpha = -rhs( 1 )
827 CALL daxpy( m-ie, alpha, a( is, ie+1 ), lda,
828 $ c( ie+1, js ), 1 )
829 alpha = -rhs( 2 )
830 CALL daxpy( m-ie, alpha, d( is, ie+1 ), ldd,
831 $ c( ie+1, js ), 1 )
832 END IF
833*
834 ELSE IF( ( mb.EQ.1 ) .AND. ( nb.EQ.2 ) ) THEN
835*
836* Build a 4-by-4 system Z**T * x = RHS
837*
838 z( 1, 1 ) = a( is, is )
839 z( 2, 1 ) = zero
840 z( 3, 1 ) = -b( js, js )
841 z( 4, 1 ) = -b( jsp1, js )
842*
843 z( 1, 2 ) = zero
844 z( 2, 2 ) = a( is, is )
845 z( 3, 2 ) = -b( js, jsp1 )
846 z( 4, 2 ) = -b( jsp1, jsp1 )
847*
848 z( 1, 3 ) = d( is, is )
849 z( 2, 3 ) = zero
850 z( 3, 3 ) = -e( js, js )
851 z( 4, 3 ) = zero
852*
853 z( 1, 4 ) = zero
854 z( 2, 4 ) = d( is, is )
855 z( 3, 4 ) = -e( js, jsp1 )
856 z( 4, 4 ) = -e( jsp1, jsp1 )
857*
858* Set up right hand side(s)
859*
860 rhs( 1 ) = c( is, js )
861 rhs( 2 ) = c( is, jsp1 )
862 rhs( 3 ) = f( is, js )
863 rhs( 4 ) = f( is, jsp1 )
864*
865* Solve Z**T * x = RHS
866*
867 CALL dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
868 IF( ierr.GT.0 )
869 $ info = ierr
870 CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
871 $ scaloc )
872 IF( scaloc.NE.one ) THEN
873 DO 140 k = 1, n
874 CALL dscal( m, scaloc, c( 1, k ), 1 )
875 CALL dscal( m, scaloc, f( 1, k ), 1 )
876 140 CONTINUE
877 scale = scale*scaloc
878 END IF
879*
880* Unpack solution vector(s)
881*
882 c( is, js ) = rhs( 1 )
883 c( is, jsp1 ) = rhs( 2 )
884 f( is, js ) = rhs( 3 )
885 f( is, jsp1 ) = rhs( 4 )
886*
887* Substitute R(I, J) and L(I, J) into remaining
888* equation.
889*
890 IF( j.GT.p+2 ) THEN
891 CALL daxpy( js-1, rhs( 1 ), b( 1, js ), 1,
892 $ f( is, 1 ), ldf )
893 CALL daxpy( js-1, rhs( 2 ), b( 1, jsp1 ), 1,
894 $ f( is, 1 ), ldf )
895 CALL daxpy( js-1, rhs( 3 ), e( 1, js ), 1,
896 $ f( is, 1 ), ldf )
897 CALL daxpy( js-1, rhs( 4 ), e( 1, jsp1 ), 1,
898 $ f( is, 1 ), ldf )
899 END IF
900 IF( i.LT.p ) THEN
901 CALL dger( m-ie, nb, -one, a( is, ie+1 ), lda,
902 $ rhs( 1 ), 1, c( ie+1, js ), ldc )
903 CALL dger( m-ie, nb, -one, d( is, ie+1 ), ldd,
904 $ rhs( 3 ), 1, c( ie+1, js ), ldc )
905 END IF
906*
907 ELSE IF( ( mb.EQ.2 ) .AND. ( nb.EQ.1 ) ) THEN
908*
909* Build a 4-by-4 system Z**T * x = RHS
910*
911 z( 1, 1 ) = a( is, is )
912 z( 2, 1 ) = a( is, isp1 )
913 z( 3, 1 ) = -b( js, js )
914 z( 4, 1 ) = zero
915*
916 z( 1, 2 ) = a( isp1, is )
917 z( 2, 2 ) = a( isp1, isp1 )
918 z( 3, 2 ) = zero
919 z( 4, 2 ) = -b( js, js )
920*
921 z( 1, 3 ) = d( is, is )
922 z( 2, 3 ) = d( is, isp1 )
923 z( 3, 3 ) = -e( js, js )
924 z( 4, 3 ) = zero
925*
926 z( 1, 4 ) = zero
927 z( 2, 4 ) = d( isp1, isp1 )
928 z( 3, 4 ) = zero
929 z( 4, 4 ) = -e( js, js )
930*
931* Set up right hand side(s)
932*
933 rhs( 1 ) = c( is, js )
934 rhs( 2 ) = c( isp1, js )
935 rhs( 3 ) = f( is, js )
936 rhs( 4 ) = f( isp1, js )
937*
938* Solve Z**T * x = RHS
939*
940 CALL dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
941 IF( ierr.GT.0 )
942 $ info = ierr
943*
944 CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
945 $ scaloc )
946 IF( scaloc.NE.one ) THEN
947 DO 150 k = 1, n
948 CALL dscal( m, scaloc, c( 1, k ), 1 )
949 CALL dscal( m, scaloc, f( 1, k ), 1 )
950 150 CONTINUE
951 scale = scale*scaloc
952 END IF
953*
954* Unpack solution vector(s)
955*
956 c( is, js ) = rhs( 1 )
957 c( isp1, js ) = rhs( 2 )
958 f( is, js ) = rhs( 3 )
959 f( isp1, js ) = rhs( 4 )
960*
961* Substitute R(I, J) and L(I, J) into remaining
962* equation.
963*
964 IF( j.GT.p+2 ) THEN
965 CALL dger( mb, js-1, one, rhs( 1 ), 1, b( 1,
966 $ js ),
967 $ 1, f( is, 1 ), ldf )
968 CALL dger( mb, js-1, one, rhs( 3 ), 1, e( 1,
969 $ js ),
970 $ 1, f( is, 1 ), ldf )
971 END IF
972 IF( i.LT.p ) THEN
973 CALL dgemv( 'T', mb, m-ie, -one, a( is, ie+1 ),
974 $ lda, rhs( 1 ), 1, one, c( ie+1, js ),
975 $ 1 )
976 CALL dgemv( 'T', mb, m-ie, -one, d( is, ie+1 ),
977 $ ldd, rhs( 3 ), 1, one, c( ie+1, js ),
978 $ 1 )
979 END IF
980*
981 ELSE IF( ( mb.EQ.2 ) .AND. ( nb.EQ.2 ) ) THEN
982*
983* Build an 8-by-8 system Z**T * x = RHS
984*
985 CALL dlaset( 'F', ldz, ldz, zero, zero, z, ldz )
986*
987 z( 1, 1 ) = a( is, is )
988 z( 2, 1 ) = a( is, isp1 )
989 z( 5, 1 ) = -b( js, js )
990 z( 7, 1 ) = -b( jsp1, js )
991*
992 z( 1, 2 ) = a( isp1, is )
993 z( 2, 2 ) = a( isp1, isp1 )
994 z( 6, 2 ) = -b( js, js )
995 z( 8, 2 ) = -b( jsp1, js )
996*
997 z( 3, 3 ) = a( is, is )
998 z( 4, 3 ) = a( is, isp1 )
999 z( 5, 3 ) = -b( js, jsp1 )
1000 z( 7, 3 ) = -b( jsp1, jsp1 )
1001*
1002 z( 3, 4 ) = a( isp1, is )
1003 z( 4, 4 ) = a( isp1, isp1 )
1004 z( 6, 4 ) = -b( js, jsp1 )
1005 z( 8, 4 ) = -b( jsp1, jsp1 )
1006*
1007 z( 1, 5 ) = d( is, is )
1008 z( 2, 5 ) = d( is, isp1 )
1009 z( 5, 5 ) = -e( js, js )
1010*
1011 z( 2, 6 ) = d( isp1, isp1 )
1012 z( 6, 6 ) = -e( js, js )
1013*
1014 z( 3, 7 ) = d( is, is )
1015 z( 4, 7 ) = d( is, isp1 )
1016 z( 5, 7 ) = -e( js, jsp1 )
1017 z( 7, 7 ) = -e( jsp1, jsp1 )
1018*
1019 z( 4, 8 ) = d( isp1, isp1 )
1020 z( 6, 8 ) = -e( js, jsp1 )
1021 z( 8, 8 ) = -e( jsp1, jsp1 )
1022*
1023* Set up right hand side(s)
1024*
1025 k = 1
1026 ii = mb*nb + 1
1027 DO 160 jj = 0, nb - 1
1028 CALL dcopy( mb, c( is, js+jj ), 1, rhs( k ), 1 )
1029 CALL dcopy( mb, f( is, js+jj ), 1, rhs( ii ),
1030 $ 1 )
1031 k = k + mb
1032 ii = ii + mb
1033 160 CONTINUE
1034*
1035*
1036* Solve Z**T * x = RHS
1037*
1038 CALL dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
1039 IF( ierr.GT.0 )
1040 $ info = ierr
1041*
1042 CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
1043 $ scaloc )
1044 IF( scaloc.NE.one ) THEN
1045 DO 170 k = 1, n
1046 CALL dscal( m, scaloc, c( 1, k ), 1 )
1047 CALL dscal( m, scaloc, f( 1, k ), 1 )
1048 170 CONTINUE
1049 scale = scale*scaloc
1050 END IF
1051*
1052* Unpack solution vector(s)
1053*
1054 k = 1
1055 ii = mb*nb + 1
1056 DO 180 jj = 0, nb - 1
1057 CALL dcopy( mb, rhs( k ), 1, c( is, js+jj ), 1 )
1058 CALL dcopy( mb, rhs( ii ), 1, f( is, js+jj ),
1059 $ 1 )
1060 k = k + mb
1061 ii = ii + mb
1062 180 CONTINUE
1063*
1064* Substitute R(I, J) and L(I, J) into remaining
1065* equation.
1066*
1067 IF( j.GT.p+2 ) THEN
1068 CALL dgemm( 'N', 'T', mb, js-1, nb, one,
1069 $ c( is, js ), ldc, b( 1, js ), ldb, one,
1070 $ f( is, 1 ), ldf )
1071 CALL dgemm( 'N', 'T', mb, js-1, nb, one,
1072 $ f( is, js ), ldf, e( 1, js ), lde, one,
1073 $ f( is, 1 ), ldf )
1074 END IF
1075 IF( i.LT.p ) THEN
1076 CALL dgemm( 'T', 'N', m-ie, nb, mb, -one,
1077 $ a( is, ie+1 ), lda, c( is, js ), ldc,
1078 $ one, c( ie+1, js ), ldc )
1079 CALL dgemm( 'T', 'N', m-ie, nb, mb, -one,
1080 $ d( is, ie+1 ), ldd, f( is, js ), ldf,
1081 $ one, c( ie+1, js ), ldc )
1082 END IF
1083*
1084 END IF
1085*
1086 190 CONTINUE
1087 200 CONTINUE
1088*
1089 END IF
1090 RETURN
1091*
1092* End of DTGSY2
1093*
1094 END
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
daxpy
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
dgemv
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:158
dger
subroutine dger(m, n, alpha, x, incx, y, incy, a, lda)
DGER
Definition dger.f:130
dgesc2
subroutine dgesc2(n, a, lda, rhs, ipiv, jpiv, scale)
DGESC2 solves a system of linear equations using the LU factorization with complete pivoting computed...
Definition dgesc2.f:112
dgetc2
subroutine dgetc2(n, a, lda, ipiv, jpiv, info)
DGETC2 computes the LU factorization with complete pivoting of the general n-by-n matrix.
Definition dgetc2.f:109
dlaset
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:108
dlatdf
subroutine dlatdf(ijob, n, z, ldz, rhs, rdsum, rdscal, ipiv, jpiv)
DLATDF uses the LU factorization of the n-by-n matrix computed by sgetc2 and computes a contribution ...
Definition dlatdf.f:169
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
dtgsy2
subroutine dtgsy2(trans, ijob, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, rdsum, rdscal, iwork, pq, info)
DTGSY2 solves the generalized Sylvester equation (unblocked algorithm).
Definition dtgsy2.f:273