LAPACK  3.9.0
LAPACK: Linear Algebra PACKage
stgsy2.f
Go to the documentation of this file.
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 *> \date December 2016
263 *
264 *> \ingroup realSYauxiliary
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 stgsy2( 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.7.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 * December 2016
281 *
282 * .. Scalar Arguments ..
283  CHARACTER TRANS
284  INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N,
285  $ pq
286  REAL RDSCAL, RDSUM, SCALE
287 * ..
288 * .. Array Arguments ..
289  INTEGER IWORK( * )
290  REAL A( LDA, * ), B( LDB, * ), C( LDC, * ),
291  $ d( ldd, * ), e( lde, * ), f( ldf, * )
292 * ..
293 *
294 * =====================================================================
295 * Replaced various illegal calls to SCOPY by calls to SLASET.
296 * Sven Hammarling, 27/5/02.
297 *
298 * .. Parameters ..
299  INTEGER LDZ
300  PARAMETER ( LDZ = 8 )
301  REAL ZERO, ONE
302  parameter( zero = 0.0e+0, one = 1.0e+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  REAL ALPHA, SCALOC
309 * ..
310 * .. Local Arrays ..
311  INTEGER IPIV( LDZ ), JPIV( LDZ )
312  REAL RHS( LDZ ), Z( LDZ, LDZ )
313 * ..
314 * .. External Functions ..
315  LOGICAL LSAME
316  EXTERNAL LSAME
317 * ..
318 * .. External Subroutines ..
319  EXTERNAL saxpy, scopy, sgemm, sgemv, sger, sgesc2,
320  $ sgetc2, sscal, slaset, slatdf, 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( 'STGSY2', -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 sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
444  IF( ierr.GT.0 )
445  $ info = ierr
446 *
447  IF( ijob.EQ.0 ) THEN
448  CALL sgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
449  $ scaloc )
450  IF( scaloc.NE.one ) THEN
451  DO 50 k = 1, n
452  CALL sscal( m, scaloc, c( 1, k ), 1 )
453  CALL sscal( m, scaloc, f( 1, k ), 1 )
454  50 CONTINUE
455  scale = scale*scaloc
456  END IF
457  ELSE
458  CALL slatdf( 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 saxpy( is-1, alpha, a( 1, is ), 1, c( 1, js ),
473  $ 1 )
474  CALL saxpy( is-1, alpha, d( 1, is ), 1, f( 1, js ),
475  $ 1 )
476  END IF
477  IF( j.LT.q ) THEN
478  CALL saxpy( n-je, rhs( 2 ), b( js, je+1 ), ldb,
479  $ c( is, je+1 ), ldc )
480  CALL saxpy( 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 sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
518  IF( ierr.GT.0 )
519  $ info = ierr
520 *
521  IF( ijob.EQ.0 ) THEN
522  CALL sgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
523  $ scaloc )
524  IF( scaloc.NE.one ) THEN
525  DO 60 k = 1, n
526  CALL sscal( m, scaloc, c( 1, k ), 1 )
527  CALL sscal( m, scaloc, f( 1, k ), 1 )
528  60 CONTINUE
529  scale = scale*scaloc
530  END IF
531  ELSE
532  CALL slatdf( 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 sger( is-1, nb, -one, a( 1, is ), 1, rhs( 1 ),
548  $ 1, c( 1, js ), ldc )
549  CALL sger( 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 saxpy( n-je, rhs( 3 ), b( js, je+1 ), ldb,
554  $ c( is, je+1 ), ldc )
555  CALL saxpy( n-je, rhs( 3 ), e( js, je+1 ), lde,
556  $ f( is, je+1 ), ldf )
557  CALL saxpy( n-je, rhs( 4 ), b( jsp1, je+1 ), ldb,
558  $ c( is, je+1 ), ldc )
559  CALL saxpy( 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 sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
597  IF( ierr.GT.0 )
598  $ info = ierr
599  IF( ijob.EQ.0 ) THEN
600  CALL sgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
601  $ scaloc )
602  IF( scaloc.NE.one ) THEN
603  DO 70 k = 1, n
604  CALL sscal( m, scaloc, c( 1, k ), 1 )
605  CALL sscal( m, scaloc, f( 1, k ), 1 )
606  70 CONTINUE
607  scale = scale*scaloc
608  END IF
609  ELSE
610  CALL slatdf( 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 sgemv( 'N', is-1, mb, -one, a( 1, is ), lda,
626  $ rhs( 1 ), 1, one, c( 1, js ), 1 )
627  CALL sgemv( '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 sger( mb, n-je, one, rhs( 3 ), 1,
632  $ b( js, je+1 ), ldb, c( is, je+1 ), ldc )
633  CALL sger( 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 slaset( '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 scopy( mb, c( is, js+jj ), 1, rhs( k ), 1 )
685  CALL scopy( 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 sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
693  IF( ierr.GT.0 )
694  $ info = ierr
695  IF( ijob.EQ.0 ) THEN
696  CALL sgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
697  $ scaloc )
698  IF( scaloc.NE.one ) THEN
699  DO 90 k = 1, n
700  CALL sscal( m, scaloc, c( 1, k ), 1 )
701  CALL sscal( m, scaloc, f( 1, k ), 1 )
702  90 CONTINUE
703  scale = scale*scaloc
704  END IF
705  ELSE
706  CALL slatdf( 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 scopy( mb, rhs( k ), 1, c( is, js+jj ), 1 )
716  CALL scopy( 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 sgemm( 'N', 'N', is-1, nb, mb, -one,
726  $ a( 1, is ), lda, rhs( 1 ), mb, one,
727  $ c( 1, js ), ldc )
728  CALL sgemm( '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 sgemm( '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 sgemm( '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 sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
785  IF( ierr.GT.0 )
786  $ info = ierr
787 *
788  CALL sgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
789  IF( scaloc.NE.one ) THEN
790  DO 130 k = 1, n
791  CALL sscal( m, scaloc, c( 1, k ), 1 )
792  CALL sscal( 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 saxpy( js-1, alpha, b( 1, js ), 1, f( is, 1 ),
808  $ ldf )
809  alpha = rhs( 2 )
810  CALL saxpy( 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 saxpy( m-ie, alpha, a( is, ie+1 ), lda,
816  $ c( ie+1, js ), 1 )
817  alpha = -rhs( 2 )
818  CALL saxpy( 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 sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
856  IF( ierr.GT.0 )
857  $ info = ierr
858  CALL sgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
859  IF( scaloc.NE.one ) THEN
860  DO 140 k = 1, n
861  CALL sscal( m, scaloc, c( 1, k ), 1 )
862  CALL sscal( 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 saxpy( js-1, rhs( 1 ), b( 1, js ), 1,
879  $ f( is, 1 ), ldf )
880  CALL saxpy( js-1, rhs( 2 ), b( 1, jsp1 ), 1,
881  $ f( is, 1 ), ldf )
882  CALL saxpy( js-1, rhs( 3 ), e( 1, js ), 1,
883  $ f( is, 1 ), ldf )
884  CALL saxpy( js-1, rhs( 4 ), e( 1, jsp1 ), 1,
885  $ f( is, 1 ), ldf )
886  END IF
887  IF( i.LT.p ) THEN
888  CALL sger( m-ie, nb, -one, a( is, ie+1 ), lda,
889  $ rhs( 1 ), 1, c( ie+1, js ), ldc )
890  CALL sger( 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 sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
928  IF( ierr.GT.0 )
929  $ info = ierr
930 *
931  CALL sgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
932  IF( scaloc.NE.one ) THEN
933  DO 150 k = 1, n
934  CALL sscal( m, scaloc, c( 1, k ), 1 )
935  CALL sscal( 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 sger( mb, js-1, one, rhs( 1 ), 1, b( 1, js ),
952  $ 1, f( is, 1 ), ldf )
953  CALL sger( 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 sgemv( 'T', mb, m-ie, -one, a( is, ie+1 ),
958  $ lda, rhs( 1 ), 1, one, c( ie+1, js ),
959  $ 1 )
960  CALL sgemv( '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 slaset( '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 scopy( mb, c( is, js+jj ), 1, rhs( k ), 1 )
1013  CALL scopy( 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 sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
1022  IF( ierr.GT.0 )
1023  $ info = ierr
1024 *
1025  CALL sgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
1026  IF( scaloc.NE.one ) THEN
1027  DO 170 k = 1, n
1028  CALL sscal( m, scaloc, c( 1, k ), 1 )
1029  CALL sscal( 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 scopy( mb, rhs( k ), 1, c( is, js+jj ), 1 )
1040  CALL scopy( 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 sgemm( 'N', 'T', mb, js-1, nb, one,
1050  $ c( is, js ), ldc, b( 1, js ), ldb, one,
1051  $ f( is, 1 ), ldf )
1052  CALL sgemm( '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 sgemm( '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 sgemm( '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 STGSY2
1074 *
1075  END
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:113
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:173
sgemm
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:189
sgemv
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:158
scopy
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:84
sscal
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:81
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
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:112
sger
subroutine sger(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
SGER
Definition: sger.f:132
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:276
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:116
saxpy
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
Definition: saxpy.f:91