LAPACK  3.11.0
LAPACK: Linear Algebra PACKage
stgsy2.f
1 *> \brief \b STGSY2 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 STGSY2 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/stgsy2.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/stgsy2.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/stgsy2.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE STGSY2( 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 * REAL RDSCAL, RDSUM, SCALE
30 * ..
31 * .. Array Arguments ..
32 * INTEGER IWORK( * )
33 * REAL 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 *> STGSY2 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 SLACON.
75 *>
76 *> STGSY2 also (IJOB >= 1) contributes to the computation in STGSYL
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 *> STGSYL. See STGSYL 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. (SGECON 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL
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 REAL
212 *> On entry, the sum of squares of computed contributions to
213 *> the Dif-estimate under computation by STGSYL, 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 STGSY2 is called by STGSYL.
219 *> \endverbatim
220 *>
221 *> \param[in,out] RDSCAL
222 *> \verbatim
223 *> RDSCAL is REAL
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 STGSY2 is called by
229 *> STGSYL.
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 stgsy2( 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  REAL RDSCAL, RDSUM, SCALE
284 * ..
285 * .. Array Arguments ..
286  INTEGER IWORK( * )
287  REAL A( lda, * ), B( ldb, * ), C( ldc, * ),
288  $ d( ldd, * ), e( lde, * ), f( ldf, * )
289 * ..
290 *
291 * =====================================================================
292 * Replaced various illegal calls to SCOPY by calls to SLASET.
293 * Sven Hammarling, 27/5/02.
294 *
295 * .. Parameters ..
296  INTEGER LDZ
297  parameter( ldz = 8 )
298  REAL ZERO, ONE
299  parameter( zero = 0.0e+0, one = 1.0e+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  REAL ALPHA, SCALOC
306 * ..
307 * .. Local Arrays ..
308  INTEGER IPIV( ldz ), JPIV( ldz )
309  REAL RHS( ldz ), Z( ldz, ldz )
310 * ..
311 * .. External Functions ..
312  LOGICAL LSAME
313  EXTERNAL lsame
314 * ..
315 * .. External Subroutines ..
316  EXTERNAL saxpy, scopy, sgemm, sgemv, sger, sgesc2,
317  $ sgetc2, sscal, slaset, slatdf, 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( 'STGSY2', -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 sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
441  IF( ierr.GT.0 )
442  $ info = ierr
443 *
444  IF( ijob.EQ.0 ) THEN
445  CALL sgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
446  $ scaloc )
447  IF( scaloc.NE.one ) THEN
448  DO 50 k = 1, n
449  CALL sscal( m, scaloc, c( 1, k ), 1 )
450  CALL sscal( m, scaloc, f( 1, k ), 1 )
451  50 CONTINUE
452  scale = scale*scaloc
453  END IF
454  ELSE
455  CALL slatdf( 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 saxpy( is-1, alpha, a( 1, is ), 1, c( 1, js ),
470  $ 1 )
471  CALL saxpy( is-1, alpha, d( 1, is ), 1, f( 1, js ),
472  $ 1 )
473  END IF
474  IF( j.LT.q ) THEN
475  CALL saxpy( n-je, rhs( 2 ), b( js, je+1 ), ldb,
476  $ c( is, je+1 ), ldc )
477  CALL saxpy( 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 sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
515  IF( ierr.GT.0 )
516  $ info = ierr
517 *
518  IF( ijob.EQ.0 ) THEN
519  CALL sgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
520  $ scaloc )
521  IF( scaloc.NE.one ) THEN
522  DO 60 k = 1, n
523  CALL sscal( m, scaloc, c( 1, k ), 1 )
524  CALL sscal( m, scaloc, f( 1, k ), 1 )
525  60 CONTINUE
526  scale = scale*scaloc
527  END IF
528  ELSE
529  CALL slatdf( 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 sger( is-1, nb, -one, a( 1, is ), 1, rhs( 1 ),
545  $ 1, c( 1, js ), ldc )
546  CALL sger( 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 saxpy( n-je, rhs( 3 ), b( js, je+1 ), ldb,
551  $ c( is, je+1 ), ldc )
552  CALL saxpy( n-je, rhs( 3 ), e( js, je+1 ), lde,
553  $ f( is, je+1 ), ldf )
554  CALL saxpy( n-je, rhs( 4 ), b( jsp1, je+1 ), ldb,
555  $ c( is, je+1 ), ldc )
556  CALL saxpy( 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 sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
594  IF( ierr.GT.0 )
595  $ info = ierr
596  IF( ijob.EQ.0 ) THEN
597  CALL sgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
598  $ scaloc )
599  IF( scaloc.NE.one ) THEN
600  DO 70 k = 1, n
601  CALL sscal( m, scaloc, c( 1, k ), 1 )
602  CALL sscal( m, scaloc, f( 1, k ), 1 )
603  70 CONTINUE
604  scale = scale*scaloc
605  END IF
606  ELSE
607  CALL slatdf( 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 sgemv( 'N', is-1, mb, -one, a( 1, is ), lda,
623  $ rhs( 1 ), 1, one, c( 1, js ), 1 )
624  CALL sgemv( '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 sger( mb, n-je, one, rhs( 3 ), 1,
629  $ b( js, je+1 ), ldb, c( is, je+1 ), ldc )
630  CALL sger( 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 slaset( '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 scopy( mb, c( is, js+jj ), 1, rhs( k ), 1 )
682  CALL scopy( 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 sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
690  IF( ierr.GT.0 )
691  $ info = ierr
692  IF( ijob.EQ.0 ) THEN
693  CALL sgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
694  $ scaloc )
695  IF( scaloc.NE.one ) THEN
696  DO 90 k = 1, n
697  CALL sscal( m, scaloc, c( 1, k ), 1 )
698  CALL sscal( m, scaloc, f( 1, k ), 1 )
699  90 CONTINUE
700  scale = scale*scaloc
701  END IF
702  ELSE
703  CALL slatdf( 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 scopy( mb, rhs( k ), 1, c( is, js+jj ), 1 )
713  CALL scopy( 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 sgemm( 'N', 'N', is-1, nb, mb, -one,
723  $ a( 1, is ), lda, rhs( 1 ), mb, one,
724  $ c( 1, js ), ldc )
725  CALL sgemm( '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 sgemm( '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 sgemm( '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 sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
782  IF( ierr.GT.0 )
783  $ info = ierr
784 *
785  CALL sgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
786  IF( scaloc.NE.one ) THEN
787  DO 130 k = 1, n
788  CALL sscal( m, scaloc, c( 1, k ), 1 )
789  CALL sscal( 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 saxpy( js-1, alpha, b( 1, js ), 1, f( is, 1 ),
805  $ ldf )
806  alpha = rhs( 2 )
807  CALL saxpy( 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 saxpy( m-ie, alpha, a( is, ie+1 ), lda,
813  $ c( ie+1, js ), 1 )
814  alpha = -rhs( 2 )
815  CALL saxpy( 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 sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
853  IF( ierr.GT.0 )
854  $ info = ierr
855  CALL sgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
856  IF( scaloc.NE.one ) THEN
857  DO 140 k = 1, n
858  CALL sscal( m, scaloc, c( 1, k ), 1 )
859  CALL sscal( 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 saxpy( js-1, rhs( 1 ), b( 1, js ), 1,
876  $ f( is, 1 ), ldf )
877  CALL saxpy( js-1, rhs( 2 ), b( 1, jsp1 ), 1,
878  $ f( is, 1 ), ldf )
879  CALL saxpy( js-1, rhs( 3 ), e( 1, js ), 1,
880  $ f( is, 1 ), ldf )
881  CALL saxpy( js-1, rhs( 4 ), e( 1, jsp1 ), 1,
882  $ f( is, 1 ), ldf )
883  END IF
884  IF( i.LT.p ) THEN
885  CALL sger( m-ie, nb, -one, a( is, ie+1 ), lda,
886  $ rhs( 1 ), 1, c( ie+1, js ), ldc )
887  CALL sger( 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 sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
925  IF( ierr.GT.0 )
926  $ info = ierr
927 *
928  CALL sgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
929  IF( scaloc.NE.one ) THEN
930  DO 150 k = 1, n
931  CALL sscal( m, scaloc, c( 1, k ), 1 )
932  CALL sscal( 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 sger( mb, js-1, one, rhs( 1 ), 1, b( 1, js ),
949  $ 1, f( is, 1 ), ldf )
950  CALL sger( 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 sgemv( 'T', mb, m-ie, -one, a( is, ie+1 ),
955  $ lda, rhs( 1 ), 1, one, c( ie+1, js ),
956  $ 1 )
957  CALL sgemv( '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 slaset( '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 scopy( mb, c( is, js+jj ), 1, rhs( k ), 1 )
1010  CALL scopy( 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 sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
1019  IF( ierr.GT.0 )
1020  $ info = ierr
1021 *
1022  CALL sgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
1023  IF( scaloc.NE.one ) THEN
1024  DO 170 k = 1, n
1025  CALL sscal( m, scaloc, c( 1, k ), 1 )
1026  CALL sscal( 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 scopy( mb, rhs( k ), 1, c( is, js+jj ), 1 )
1037  CALL scopy( 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 sgemm( 'N', 'T', mb, js-1, nb, one,
1047  $ c( is, js ), ldc, b( 1, js ), ldb, one,
1048  $ f( is, 1 ), ldf )
1049  CALL sgemm( '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 sgemm( '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 sgemm( '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 STGSY2
1071 *
1072  END
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
stgsy2
subroutine stgsy2(TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL, IWORK, PQ, INFO)
STGSY2 solves the generalized Sylvester equation (unblocked algorithm).
Definition: stgsy2.f:274
sgemm
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187
sgesc2
subroutine sgesc2(N, A, LDA, RHS, IPIV, JPIV, SCALE)
SGESC2 solves a system of linear equations using the LU factorization with complete pivoting computed...
Definition: sgesc2.f:114
sgetc2
subroutine sgetc2(N, A, LDA, IPIV, JPIV, INFO)
SGETC2 computes the LU factorization with complete pivoting of the general n-by-n matrix...
Definition: sgetc2.f:111
slaset
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: slaset.f:110
scopy
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
saxpy
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
Definition: saxpy.f:89
slatdf
subroutine slatdf(IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV, JPIV)
SLATDF uses the LU factorization of the n-by-n matrix computed by sgetc2 and computes a contribution ...
Definition: slatdf.f:171
sscal
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
sgemv
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:156
sger
subroutine sger(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
SGER
Definition: sger.f:130