LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
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 doubleSYauxiliary
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,
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
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
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:89
subroutine dger(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
DGER
Definition: dger.f:130
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:156
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187
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
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
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
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