LAPACK  3.6.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 communicaton 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 *> \date November 2015
263 *
264 *> \ingroup doubleSYauxiliary
265 *
266 *> \par Contributors:
267 * ==================
268 *>
269 *> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
270 *> Umea University, S-901 87 Umea, Sweden.
271 *
272 * =====================================================================
273  SUBROUTINE dtgsy2( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
274  $ ldd, e, lde, f, ldf, scale, rdsum, rdscal,
275  $ iwork, pq, info )
276 *
277 * -- LAPACK auxiliary routine (version 3.6.0) --
278 * -- LAPACK is a software package provided by Univ. of Tennessee, --
279 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
280 * November 2015
281 *
282 * .. Scalar Arguments ..
283  CHARACTER TRANS
284  INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N,
285  $ pq
286  DOUBLE PRECISION RDSCAL, RDSUM, SCALE
287 * ..
288 * .. Array Arguments ..
289  INTEGER IWORK( * )
290  DOUBLE PRECISION A( lda, * ), B( ldb, * ), C( ldc, * ),
291  $ d( ldd, * ), e( lde, * ), f( ldf, * )
292 * ..
293 *
294 * =====================================================================
295 * Replaced various illegal calls to DCOPY by calls to DLASET.
296 * Sven Hammarling, 27/5/02.
297 *
298 * .. Parameters ..
299  INTEGER LDZ
300  parameter ( ldz = 8 )
301  DOUBLE PRECISION ZERO, ONE
302  parameter ( zero = 0.0d+0, one = 1.0d+0 )
303 * ..
304 * .. Local Scalars ..
305  LOGICAL NOTRAN
306  INTEGER I, IE, IERR, II, IS, ISP1, J, JE, JJ, JS, JSP1,
307  $ k, mb, nb, p, q, zdim
308  DOUBLE PRECISION ALPHA, SCALOC
309 * ..
310 * .. Local Arrays ..
311  INTEGER IPIV( ldz ), JPIV( ldz )
312  DOUBLE PRECISION RHS( ldz ), Z( ldz, ldz )
313 * ..
314 * .. External Functions ..
315  LOGICAL LSAME
316  EXTERNAL lsame
317 * ..
318 * .. External Subroutines ..
319  EXTERNAL daxpy, dcopy, dgemm, dgemv, dger, dgesc2,
320  $ dgetc2, dlaset, dlatdf, dscal, xerbla
321 * ..
322 * .. Intrinsic Functions ..
323  INTRINSIC max
324 * ..
325 * .. Executable Statements ..
326 *
327 * Decode and test input parameters
328 *
329  info = 0
330  ierr = 0
331  notran = lsame( trans, 'N' )
332  IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) ) THEN
333  info = -1
334  ELSE IF( notran ) THEN
335  IF( ( ijob.LT.0 ) .OR. ( ijob.GT.2 ) ) THEN
336  info = -2
337  END IF
338  END IF
339  IF( info.EQ.0 ) THEN
340  IF( m.LE.0 ) THEN
341  info = -3
342  ELSE IF( n.LE.0 ) THEN
343  info = -4
344  ELSE IF( lda.LT.max( 1, m ) ) THEN
345  info = -6
346  ELSE IF( ldb.LT.max( 1, n ) ) THEN
347  info = -8
348  ELSE IF( ldc.LT.max( 1, m ) ) THEN
349  info = -10
350  ELSE IF( ldd.LT.max( 1, m ) ) THEN
351  info = -12
352  ELSE IF( lde.LT.max( 1, n ) ) THEN
353  info = -14
354  ELSE IF( ldf.LT.max( 1, m ) ) THEN
355  info = -16
356  END IF
357  END IF
358  IF( info.NE.0 ) THEN
359  CALL xerbla( 'DTGSY2', -info )
360  RETURN
361  END IF
362 *
363 * Determine block structure of A
364 *
365  pq = 0
366  p = 0
367  i = 1
368  10 CONTINUE
369  IF( i.GT.m )
370  $ GO TO 20
371  p = p + 1
372  iwork( p ) = i
373  IF( i.EQ.m )
374  $ GO TO 20
375  IF( a( i+1, i ).NE.zero ) THEN
376  i = i + 2
377  ELSE
378  i = i + 1
379  END IF
380  GO TO 10
381  20 CONTINUE
382  iwork( p+1 ) = m + 1
383 *
384 * Determine block structure of B
385 *
386  q = p + 1
387  j = 1
388  30 CONTINUE
389  IF( j.GT.n )
390  $ GO TO 40
391  q = q + 1
392  iwork( q ) = j
393  IF( j.EQ.n )
394  $ GO TO 40
395  IF( b( j+1, j ).NE.zero ) THEN
396  j = j + 2
397  ELSE
398  j = j + 1
399  END IF
400  GO TO 30
401  40 CONTINUE
402  iwork( q+1 ) = n + 1
403  pq = p*( q-p-1 )
404 *
405  IF( notran ) THEN
406 *
407 * Solve (I, J) - subsystem
408 * A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
409 * D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
410 * for I = P, P - 1, ..., 1; J = 1, 2, ..., Q
411 *
412  scale = one
413  scaloc = one
414  DO 120 j = p + 2, q
415  js = iwork( j )
416  jsp1 = js + 1
417  je = iwork( j+1 ) - 1
418  nb = je - js + 1
419  DO 110 i = p, 1, -1
420 *
421  is = iwork( i )
422  isp1 = is + 1
423  ie = iwork( i+1 ) - 1
424  mb = ie - is + 1
425  zdim = mb*nb*2
426 *
427  IF( ( mb.EQ.1 ) .AND. ( nb.EQ.1 ) ) THEN
428 *
429 * Build a 2-by-2 system Z * x = RHS
430 *
431  z( 1, 1 ) = a( is, is )
432  z( 2, 1 ) = d( is, is )
433  z( 1, 2 ) = -b( js, js )
434  z( 2, 2 ) = -e( js, js )
435 *
436 * Set up right hand side(s)
437 *
438  rhs( 1 ) = c( is, js )
439  rhs( 2 ) = f( is, js )
440 *
441 * Solve Z * x = RHS
442 *
443  CALL dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
444  IF( ierr.GT.0 )
445  $ info = ierr
446 *
447  IF( ijob.EQ.0 ) THEN
448  CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
449  $ scaloc )
450  IF( scaloc.NE.one ) THEN
451  DO 50 k = 1, n
452  CALL dscal( m, scaloc, c( 1, k ), 1 )
453  CALL dscal( m, scaloc, f( 1, k ), 1 )
454  50 CONTINUE
455  scale = scale*scaloc
456  END IF
457  ELSE
458  CALL dlatdf( ijob, zdim, z, ldz, rhs, rdsum,
459  $ rdscal, ipiv, jpiv )
460  END IF
461 *
462 * Unpack solution vector(s)
463 *
464  c( is, js ) = rhs( 1 )
465  f( is, js ) = rhs( 2 )
466 *
467 * Substitute R(I, J) and L(I, J) into remaining
468 * equation.
469 *
470  IF( i.GT.1 ) THEN
471  alpha = -rhs( 1 )
472  CALL daxpy( is-1, alpha, a( 1, is ), 1, c( 1, js ),
473  $ 1 )
474  CALL daxpy( is-1, alpha, d( 1, is ), 1, f( 1, js ),
475  $ 1 )
476  END IF
477  IF( j.LT.q ) THEN
478  CALL daxpy( n-je, rhs( 2 ), b( js, je+1 ), ldb,
479  $ c( is, je+1 ), ldc )
480  CALL daxpy( n-je, rhs( 2 ), e( js, je+1 ), lde,
481  $ f( is, je+1 ), ldf )
482  END IF
483 *
484  ELSE IF( ( mb.EQ.1 ) .AND. ( nb.EQ.2 ) ) THEN
485 *
486 * Build a 4-by-4 system Z * x = RHS
487 *
488  z( 1, 1 ) = a( is, is )
489  z( 2, 1 ) = zero
490  z( 3, 1 ) = d( is, is )
491  z( 4, 1 ) = zero
492 *
493  z( 1, 2 ) = zero
494  z( 2, 2 ) = a( is, is )
495  z( 3, 2 ) = zero
496  z( 4, 2 ) = d( is, is )
497 *
498  z( 1, 3 ) = -b( js, js )
499  z( 2, 3 ) = -b( js, jsp1 )
500  z( 3, 3 ) = -e( js, js )
501  z( 4, 3 ) = -e( js, jsp1 )
502 *
503  z( 1, 4 ) = -b( jsp1, js )
504  z( 2, 4 ) = -b( jsp1, jsp1 )
505  z( 3, 4 ) = zero
506  z( 4, 4 ) = -e( jsp1, jsp1 )
507 *
508 * Set up right hand side(s)
509 *
510  rhs( 1 ) = c( is, js )
511  rhs( 2 ) = c( is, jsp1 )
512  rhs( 3 ) = f( is, js )
513  rhs( 4 ) = f( is, jsp1 )
514 *
515 * Solve Z * x = RHS
516 *
517  CALL dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
518  IF( ierr.GT.0 )
519  $ info = ierr
520 *
521  IF( ijob.EQ.0 ) THEN
522  CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
523  $ scaloc )
524  IF( scaloc.NE.one ) THEN
525  DO 60 k = 1, n
526  CALL dscal( m, scaloc, c( 1, k ), 1 )
527  CALL dscal( m, scaloc, f( 1, k ), 1 )
528  60 CONTINUE
529  scale = scale*scaloc
530  END IF
531  ELSE
532  CALL dlatdf( ijob, zdim, z, ldz, rhs, rdsum,
533  $ rdscal, ipiv, jpiv )
534  END IF
535 *
536 * Unpack solution vector(s)
537 *
538  c( is, js ) = rhs( 1 )
539  c( is, jsp1 ) = rhs( 2 )
540  f( is, js ) = rhs( 3 )
541  f( is, jsp1 ) = rhs( 4 )
542 *
543 * Substitute R(I, J) and L(I, J) into remaining
544 * equation.
545 *
546  IF( i.GT.1 ) THEN
547  CALL dger( is-1, nb, -one, a( 1, is ), 1, rhs( 1 ),
548  $ 1, c( 1, js ), ldc )
549  CALL dger( is-1, nb, -one, d( 1, is ), 1, rhs( 1 ),
550  $ 1, f( 1, js ), ldf )
551  END IF
552  IF( j.LT.q ) THEN
553  CALL daxpy( n-je, rhs( 3 ), b( js, je+1 ), ldb,
554  $ c( is, je+1 ), ldc )
555  CALL daxpy( n-je, rhs( 3 ), e( js, je+1 ), lde,
556  $ f( is, je+1 ), ldf )
557  CALL daxpy( n-je, rhs( 4 ), b( jsp1, je+1 ), ldb,
558  $ c( is, je+1 ), ldc )
559  CALL daxpy( n-je, rhs( 4 ), e( jsp1, je+1 ), lde,
560  $ f( is, je+1 ), ldf )
561  END IF
562 *
563  ELSE IF( ( mb.EQ.2 ) .AND. ( nb.EQ.1 ) ) THEN
564 *
565 * Build a 4-by-4 system Z * x = RHS
566 *
567  z( 1, 1 ) = a( is, is )
568  z( 2, 1 ) = a( isp1, is )
569  z( 3, 1 ) = d( is, is )
570  z( 4, 1 ) = zero
571 *
572  z( 1, 2 ) = a( is, isp1 )
573  z( 2, 2 ) = a( isp1, isp1 )
574  z( 3, 2 ) = d( is, isp1 )
575  z( 4, 2 ) = d( isp1, isp1 )
576 *
577  z( 1, 3 ) = -b( js, js )
578  z( 2, 3 ) = zero
579  z( 3, 3 ) = -e( js, js )
580  z( 4, 3 ) = zero
581 *
582  z( 1, 4 ) = zero
583  z( 2, 4 ) = -b( js, js )
584  z( 3, 4 ) = zero
585  z( 4, 4 ) = -e( js, js )
586 *
587 * Set up right hand side(s)
588 *
589  rhs( 1 ) = c( is, js )
590  rhs( 2 ) = c( isp1, js )
591  rhs( 3 ) = f( is, js )
592  rhs( 4 ) = f( isp1, js )
593 *
594 * Solve Z * x = RHS
595 *
596  CALL dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
597  IF( ierr.GT.0 )
598  $ info = ierr
599  IF( ijob.EQ.0 ) THEN
600  CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
601  $ scaloc )
602  IF( scaloc.NE.one ) THEN
603  DO 70 k = 1, n
604  CALL dscal( m, scaloc, c( 1, k ), 1 )
605  CALL dscal( m, scaloc, f( 1, k ), 1 )
606  70 CONTINUE
607  scale = scale*scaloc
608  END IF
609  ELSE
610  CALL dlatdf( ijob, zdim, z, ldz, rhs, rdsum,
611  $ rdscal, ipiv, jpiv )
612  END IF
613 *
614 * Unpack solution vector(s)
615 *
616  c( is, js ) = rhs( 1 )
617  c( isp1, js ) = rhs( 2 )
618  f( is, js ) = rhs( 3 )
619  f( isp1, js ) = rhs( 4 )
620 *
621 * Substitute R(I, J) and L(I, J) into remaining
622 * equation.
623 *
624  IF( i.GT.1 ) THEN
625  CALL dgemv( 'N', is-1, mb, -one, a( 1, is ), lda,
626  $ rhs( 1 ), 1, one, c( 1, js ), 1 )
627  CALL dgemv( 'N', is-1, mb, -one, d( 1, is ), ldd,
628  $ rhs( 1 ), 1, one, f( 1, js ), 1 )
629  END IF
630  IF( j.LT.q ) THEN
631  CALL dger( mb, n-je, one, rhs( 3 ), 1,
632  $ b( js, je+1 ), ldb, c( is, je+1 ), ldc )
633  CALL dger( mb, n-je, one, rhs( 3 ), 1,
634  $ e( js, je+1 ), lde, f( is, je+1 ), ldf )
635  END IF
636 *
637  ELSE IF( ( mb.EQ.2 ) .AND. ( nb.EQ.2 ) ) THEN
638 *
639 * Build an 8-by-8 system Z * x = RHS
640 *
641  CALL dlaset( 'F', ldz, ldz, zero, zero, z, ldz )
642 *
643  z( 1, 1 ) = a( is, is )
644  z( 2, 1 ) = a( isp1, is )
645  z( 5, 1 ) = d( is, is )
646 *
647  z( 1, 2 ) = a( is, isp1 )
648  z( 2, 2 ) = a( isp1, isp1 )
649  z( 5, 2 ) = d( is, isp1 )
650  z( 6, 2 ) = d( isp1, isp1 )
651 *
652  z( 3, 3 ) = a( is, is )
653  z( 4, 3 ) = a( isp1, is )
654  z( 7, 3 ) = d( is, is )
655 *
656  z( 3, 4 ) = a( is, isp1 )
657  z( 4, 4 ) = a( isp1, isp1 )
658  z( 7, 4 ) = d( is, isp1 )
659  z( 8, 4 ) = d( isp1, isp1 )
660 *
661  z( 1, 5 ) = -b( js, js )
662  z( 3, 5 ) = -b( js, jsp1 )
663  z( 5, 5 ) = -e( js, js )
664  z( 7, 5 ) = -e( js, jsp1 )
665 *
666  z( 2, 6 ) = -b( js, js )
667  z( 4, 6 ) = -b( js, jsp1 )
668  z( 6, 6 ) = -e( js, js )
669  z( 8, 6 ) = -e( js, jsp1 )
670 *
671  z( 1, 7 ) = -b( jsp1, js )
672  z( 3, 7 ) = -b( jsp1, jsp1 )
673  z( 7, 7 ) = -e( jsp1, jsp1 )
674 *
675  z( 2, 8 ) = -b( jsp1, js )
676  z( 4, 8 ) = -b( jsp1, jsp1 )
677  z( 8, 8 ) = -e( jsp1, jsp1 )
678 *
679 * Set up right hand side(s)
680 *
681  k = 1
682  ii = mb*nb + 1
683  DO 80 jj = 0, nb - 1
684  CALL dcopy( mb, c( is, js+jj ), 1, rhs( k ), 1 )
685  CALL dcopy( mb, f( is, js+jj ), 1, rhs( ii ), 1 )
686  k = k + mb
687  ii = ii + mb
688  80 CONTINUE
689 *
690 * Solve Z * x = RHS
691 *
692  CALL dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
693  IF( ierr.GT.0 )
694  $ info = ierr
695  IF( ijob.EQ.0 ) THEN
696  CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
697  $ scaloc )
698  IF( scaloc.NE.one ) THEN
699  DO 90 k = 1, n
700  CALL dscal( m, scaloc, c( 1, k ), 1 )
701  CALL dscal( m, scaloc, f( 1, k ), 1 )
702  90 CONTINUE
703  scale = scale*scaloc
704  END IF
705  ELSE
706  CALL dlatdf( ijob, zdim, z, ldz, rhs, rdsum,
707  $ rdscal, ipiv, jpiv )
708  END IF
709 *
710 * Unpack solution vector(s)
711 *
712  k = 1
713  ii = mb*nb + 1
714  DO 100 jj = 0, nb - 1
715  CALL dcopy( mb, rhs( k ), 1, c( is, js+jj ), 1 )
716  CALL dcopy( mb, rhs( ii ), 1, f( is, js+jj ), 1 )
717  k = k + mb
718  ii = ii + mb
719  100 CONTINUE
720 *
721 * Substitute R(I, J) and L(I, J) into remaining
722 * equation.
723 *
724  IF( i.GT.1 ) THEN
725  CALL dgemm( 'N', 'N', is-1, nb, mb, -one,
726  $ a( 1, is ), lda, rhs( 1 ), mb, one,
727  $ c( 1, js ), ldc )
728  CALL dgemm( 'N', 'N', is-1, nb, mb, -one,
729  $ d( 1, is ), ldd, rhs( 1 ), mb, one,
730  $ f( 1, js ), ldf )
731  END IF
732  IF( j.LT.q ) THEN
733  k = mb*nb + 1
734  CALL dgemm( 'N', 'N', mb, n-je, nb, one, rhs( k ),
735  $ mb, b( js, je+1 ), ldb, one,
736  $ c( is, je+1 ), ldc )
737  CALL dgemm( 'N', 'N', mb, n-je, nb, one, rhs( k ),
738  $ mb, e( js, je+1 ), lde, one,
739  $ f( is, je+1 ), ldf )
740  END IF
741 *
742  END IF
743 *
744  110 CONTINUE
745  120 CONTINUE
746  ELSE
747 *
748 * Solve (I, J) - subsystem
749 * A(I, I)**T * R(I, J) + D(I, I)**T * L(J, J) = C(I, J)
750 * R(I, I) * B(J, J) + L(I, J) * E(J, J) = -F(I, J)
751 * for I = 1, 2, ..., P, J = Q, Q - 1, ..., 1
752 *
753  scale = one
754  scaloc = one
755  DO 200 i = 1, p
756 *
757  is = iwork( i )
758  isp1 = is + 1
759  ie = iwork( i+1 ) - 1
760  mb = ie - is + 1
761  DO 190 j = q, p + 2, -1
762 *
763  js = iwork( j )
764  jsp1 = js + 1
765  je = iwork( j+1 ) - 1
766  nb = je - js + 1
767  zdim = mb*nb*2
768  IF( ( mb.EQ.1 ) .AND. ( nb.EQ.1 ) ) THEN
769 *
770 * Build a 2-by-2 system Z**T * x = RHS
771 *
772  z( 1, 1 ) = a( is, is )
773  z( 2, 1 ) = -b( js, js )
774  z( 1, 2 ) = d( is, is )
775  z( 2, 2 ) = -e( js, js )
776 *
777 * Set up right hand side(s)
778 *
779  rhs( 1 ) = c( is, js )
780  rhs( 2 ) = f( is, js )
781 *
782 * Solve Z**T * x = RHS
783 *
784  CALL dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
785  IF( ierr.GT.0 )
786  $ info = ierr
787 *
788  CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
789  IF( scaloc.NE.one ) THEN
790  DO 130 k = 1, n
791  CALL dscal( m, scaloc, c( 1, k ), 1 )
792  CALL dscal( m, scaloc, f( 1, k ), 1 )
793  130 CONTINUE
794  scale = scale*scaloc
795  END IF
796 *
797 * Unpack solution vector(s)
798 *
799  c( is, js ) = rhs( 1 )
800  f( is, js ) = rhs( 2 )
801 *
802 * Substitute R(I, J) and L(I, J) into remaining
803 * equation.
804 *
805  IF( j.GT.p+2 ) THEN
806  alpha = rhs( 1 )
807  CALL daxpy( js-1, alpha, b( 1, js ), 1, f( is, 1 ),
808  $ ldf )
809  alpha = rhs( 2 )
810  CALL daxpy( js-1, alpha, e( 1, js ), 1, f( is, 1 ),
811  $ ldf )
812  END IF
813  IF( i.LT.p ) THEN
814  alpha = -rhs( 1 )
815  CALL daxpy( m-ie, alpha, a( is, ie+1 ), lda,
816  $ c( ie+1, js ), 1 )
817  alpha = -rhs( 2 )
818  CALL daxpy( m-ie, alpha, d( is, ie+1 ), ldd,
819  $ c( ie+1, js ), 1 )
820  END IF
821 *
822  ELSE IF( ( mb.EQ.1 ) .AND. ( nb.EQ.2 ) ) THEN
823 *
824 * Build a 4-by-4 system Z**T * x = RHS
825 *
826  z( 1, 1 ) = a( is, is )
827  z( 2, 1 ) = zero
828  z( 3, 1 ) = -b( js, js )
829  z( 4, 1 ) = -b( jsp1, js )
830 *
831  z( 1, 2 ) = zero
832  z( 2, 2 ) = a( is, is )
833  z( 3, 2 ) = -b( js, jsp1 )
834  z( 4, 2 ) = -b( jsp1, jsp1 )
835 *
836  z( 1, 3 ) = d( is, is )
837  z( 2, 3 ) = zero
838  z( 3, 3 ) = -e( js, js )
839  z( 4, 3 ) = zero
840 *
841  z( 1, 4 ) = zero
842  z( 2, 4 ) = d( is, is )
843  z( 3, 4 ) = -e( js, jsp1 )
844  z( 4, 4 ) = -e( jsp1, jsp1 )
845 *
846 * Set up right hand side(s)
847 *
848  rhs( 1 ) = c( is, js )
849  rhs( 2 ) = c( is, jsp1 )
850  rhs( 3 ) = f( is, js )
851  rhs( 4 ) = f( is, jsp1 )
852 *
853 * Solve Z**T * x = RHS
854 *
855  CALL dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
856  IF( ierr.GT.0 )
857  $ info = ierr
858  CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
859  IF( scaloc.NE.one ) THEN
860  DO 140 k = 1, n
861  CALL dscal( m, scaloc, c( 1, k ), 1 )
862  CALL dscal( m, scaloc, f( 1, k ), 1 )
863  140 CONTINUE
864  scale = scale*scaloc
865  END IF
866 *
867 * Unpack solution vector(s)
868 *
869  c( is, js ) = rhs( 1 )
870  c( is, jsp1 ) = rhs( 2 )
871  f( is, js ) = rhs( 3 )
872  f( is, jsp1 ) = rhs( 4 )
873 *
874 * Substitute R(I, J) and L(I, J) into remaining
875 * equation.
876 *
877  IF( j.GT.p+2 ) THEN
878  CALL daxpy( js-1, rhs( 1 ), b( 1, js ), 1,
879  $ f( is, 1 ), ldf )
880  CALL daxpy( js-1, rhs( 2 ), b( 1, jsp1 ), 1,
881  $ f( is, 1 ), ldf )
882  CALL daxpy( js-1, rhs( 3 ), e( 1, js ), 1,
883  $ f( is, 1 ), ldf )
884  CALL daxpy( js-1, rhs( 4 ), e( 1, jsp1 ), 1,
885  $ f( is, 1 ), ldf )
886  END IF
887  IF( i.LT.p ) THEN
888  CALL dger( m-ie, nb, -one, a( is, ie+1 ), lda,
889  $ rhs( 1 ), 1, c( ie+1, js ), ldc )
890  CALL dger( m-ie, nb, -one, d( is, ie+1 ), ldd,
891  $ rhs( 3 ), 1, c( ie+1, js ), ldc )
892  END IF
893 *
894  ELSE IF( ( mb.EQ.2 ) .AND. ( nb.EQ.1 ) ) THEN
895 *
896 * Build a 4-by-4 system Z**T * x = RHS
897 *
898  z( 1, 1 ) = a( is, is )
899  z( 2, 1 ) = a( is, isp1 )
900  z( 3, 1 ) = -b( js, js )
901  z( 4, 1 ) = zero
902 *
903  z( 1, 2 ) = a( isp1, is )
904  z( 2, 2 ) = a( isp1, isp1 )
905  z( 3, 2 ) = zero
906  z( 4, 2 ) = -b( js, js )
907 *
908  z( 1, 3 ) = d( is, is )
909  z( 2, 3 ) = d( is, isp1 )
910  z( 3, 3 ) = -e( js, js )
911  z( 4, 3 ) = zero
912 *
913  z( 1, 4 ) = zero
914  z( 2, 4 ) = d( isp1, isp1 )
915  z( 3, 4 ) = zero
916  z( 4, 4 ) = -e( js, js )
917 *
918 * Set up right hand side(s)
919 *
920  rhs( 1 ) = c( is, js )
921  rhs( 2 ) = c( isp1, js )
922  rhs( 3 ) = f( is, js )
923  rhs( 4 ) = f( isp1, js )
924 *
925 * Solve Z**T * x = RHS
926 *
927  CALL dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
928  IF( ierr.GT.0 )
929  $ info = ierr
930 *
931  CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
932  IF( scaloc.NE.one ) THEN
933  DO 150 k = 1, n
934  CALL dscal( m, scaloc, c( 1, k ), 1 )
935  CALL dscal( m, scaloc, f( 1, k ), 1 )
936  150 CONTINUE
937  scale = scale*scaloc
938  END IF
939 *
940 * Unpack solution vector(s)
941 *
942  c( is, js ) = rhs( 1 )
943  c( isp1, js ) = rhs( 2 )
944  f( is, js ) = rhs( 3 )
945  f( isp1, js ) = rhs( 4 )
946 *
947 * Substitute R(I, J) and L(I, J) into remaining
948 * equation.
949 *
950  IF( j.GT.p+2 ) THEN
951  CALL dger( mb, js-1, one, rhs( 1 ), 1, b( 1, js ),
952  $ 1, f( is, 1 ), ldf )
953  CALL dger( mb, js-1, one, rhs( 3 ), 1, e( 1, js ),
954  $ 1, f( is, 1 ), ldf )
955  END IF
956  IF( i.LT.p ) THEN
957  CALL dgemv( 'T', mb, m-ie, -one, a( is, ie+1 ),
958  $ lda, rhs( 1 ), 1, one, c( ie+1, js ),
959  $ 1 )
960  CALL dgemv( 'T', mb, m-ie, -one, d( is, ie+1 ),
961  $ ldd, rhs( 3 ), 1, one, c( ie+1, js ),
962  $ 1 )
963  END IF
964 *
965  ELSE IF( ( mb.EQ.2 ) .AND. ( nb.EQ.2 ) ) THEN
966 *
967 * Build an 8-by-8 system Z**T * x = RHS
968 *
969  CALL dlaset( 'F', ldz, ldz, zero, zero, z, ldz )
970 *
971  z( 1, 1 ) = a( is, is )
972  z( 2, 1 ) = a( is, isp1 )
973  z( 5, 1 ) = -b( js, js )
974  z( 7, 1 ) = -b( jsp1, js )
975 *
976  z( 1, 2 ) = a( isp1, is )
977  z( 2, 2 ) = a( isp1, isp1 )
978  z( 6, 2 ) = -b( js, js )
979  z( 8, 2 ) = -b( jsp1, js )
980 *
981  z( 3, 3 ) = a( is, is )
982  z( 4, 3 ) = a( is, isp1 )
983  z( 5, 3 ) = -b( js, jsp1 )
984  z( 7, 3 ) = -b( jsp1, jsp1 )
985 *
986  z( 3, 4 ) = a( isp1, is )
987  z( 4, 4 ) = a( isp1, isp1 )
988  z( 6, 4 ) = -b( js, jsp1 )
989  z( 8, 4 ) = -b( jsp1, jsp1 )
990 *
991  z( 1, 5 ) = d( is, is )
992  z( 2, 5 ) = d( is, isp1 )
993  z( 5, 5 ) = -e( js, js )
994 *
995  z( 2, 6 ) = d( isp1, isp1 )
996  z( 6, 6 ) = -e( js, js )
997 *
998  z( 3, 7 ) = d( is, is )
999  z( 4, 7 ) = d( is, isp1 )
1000  z( 5, 7 ) = -e( js, jsp1 )
1001  z( 7, 7 ) = -e( jsp1, jsp1 )
1002 *
1003  z( 4, 8 ) = d( isp1, isp1 )
1004  z( 6, 8 ) = -e( js, jsp1 )
1005  z( 8, 8 ) = -e( jsp1, jsp1 )
1006 *
1007 * Set up right hand side(s)
1008 *
1009  k = 1
1010  ii = mb*nb + 1
1011  DO 160 jj = 0, nb - 1
1012  CALL dcopy( mb, c( is, js+jj ), 1, rhs( k ), 1 )
1013  CALL dcopy( mb, f( is, js+jj ), 1, rhs( ii ), 1 )
1014  k = k + mb
1015  ii = ii + mb
1016  160 CONTINUE
1017 *
1018 *
1019 * Solve Z**T * x = RHS
1020 *
1021  CALL dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
1022  IF( ierr.GT.0 )
1023  $ info = ierr
1024 *
1025  CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
1026  IF( scaloc.NE.one ) THEN
1027  DO 170 k = 1, n
1028  CALL dscal( m, scaloc, c( 1, k ), 1 )
1029  CALL dscal( m, scaloc, f( 1, k ), 1 )
1030  170 CONTINUE
1031  scale = scale*scaloc
1032  END IF
1033 *
1034 * Unpack solution vector(s)
1035 *
1036  k = 1
1037  ii = mb*nb + 1
1038  DO 180 jj = 0, nb - 1
1039  CALL dcopy( mb, rhs( k ), 1, c( is, js+jj ), 1 )
1040  CALL dcopy( mb, rhs( ii ), 1, f( is, js+jj ), 1 )
1041  k = k + mb
1042  ii = ii + mb
1043  180 CONTINUE
1044 *
1045 * Substitute R(I, J) and L(I, J) into remaining
1046 * equation.
1047 *
1048  IF( j.GT.p+2 ) THEN
1049  CALL dgemm( 'N', 'T', mb, js-1, nb, one,
1050  $ c( is, js ), ldc, b( 1, js ), ldb, one,
1051  $ f( is, 1 ), ldf )
1052  CALL dgemm( 'N', 'T', mb, js-1, nb, one,
1053  $ f( is, js ), ldf, e( 1, js ), lde, one,
1054  $ f( is, 1 ), ldf )
1055  END IF
1056  IF( i.LT.p ) THEN
1057  CALL dgemm( 'T', 'N', m-ie, nb, mb, -one,
1058  $ a( is, ie+1 ), lda, c( is, js ), ldc,
1059  $ one, c( ie+1, js ), ldc )
1060  CALL dgemm( 'T', 'N', m-ie, nb, mb, -one,
1061  $ d( is, ie+1 ), ldd, f( is, js ), ldf,
1062  $ one, c( ie+1, js ), ldc )
1063  END IF
1064 *
1065  END IF
1066 *
1067  190 CONTINUE
1068  200 CONTINUE
1069 *
1070  END IF
1071  RETURN
1072 *
1073 * End of DTGSY2
1074 *
1075  END
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:112
dcopy
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
dgemv
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:158
lde
logical function lde(RI, RJ, LR)
Definition: dblat2.f:2945
daxpy
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:54
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:116
dgemm
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:189
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:173
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
dger
subroutine dger(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
DGER
Definition: dger.f:132
dscal
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55
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:276
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:113