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