LAPACK 3.12.0
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*> \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)
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: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
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
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