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