LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ stgsy2()

subroutine stgsy2 ( character trans,
integer ijob,
integer m,
integer n,
real, dimension( lda, * ) a,
integer lda,
real, dimension( ldb, * ) b,
integer ldb,
real, dimension( ldc, * ) c,
integer ldc,
real, dimension( ldd, * ) d,
integer ldd,
real, dimension( lde, * ) e,
integer lde,
real, dimension( ldf, * ) f,
integer ldf,
real scale,
real rdsum,
real rdscal,
integer, dimension( * ) iwork,
integer pq,
integer info )

STGSY2 solves the generalized Sylvester equation (unblocked algorithm).

Download STGSY2 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> STGSY2 solves the generalized Sylvester equation:
!>
!>             A * R - L * B = scale * C                (1)
!>             D * R - L * E = scale * F,
!>
!> using Level 1 and 2 BLAS. where R and L are unknown M-by-N matrices,
!> (A, D), (B, E) and (C, F) are given matrix pairs of size M-by-M,
!> N-by-N and M-by-N, respectively, with real entries. (A, D) and (B, E)
!> must be in generalized Schur canonical form, i.e. A, B are upper
!> quasi triangular and D, E are upper triangular. The solution (R, L)
!> overwrites (C, F). 0 <= SCALE <= 1 is an output scaling factor
!> chosen to avoid overflow.
!>
!> In matrix notation solving equation (1) corresponds to solve
!> Z*x = scale*b, where Z is defined as
!>
!>        Z = [ kron(In, A)  -kron(B**T, Im) ]             (2)
!>            [ kron(In, D)  -kron(E**T, Im) ],
!>
!> Ik is the identity matrix of size k and X**T is the transpose of X.
!> kron(X, Y) is the Kronecker product between the matrices X and Y.
!> In the process of solving (1), we solve a number of such systems
!> where Dim(In), Dim(In) = 1 or 2.
!>
!> If TRANS = 'T', solve the transposed system Z**T*y = scale*b for y,
!> which is equivalent to solve for R and L in
!>
!>             A**T * R  + D**T * L   = scale * C           (3)
!>             R  * B**T + L  * E**T  = scale * -F
!>
!> This case is used to compute an estimate of Dif[(A, D), (B, E)] =
!> sigma_min(Z) using reverse communication with SLACON.
!>
!> STGSY2 also (IJOB >= 1) contributes to the computation in STGSYL
!> of an upper bound on the separation between to matrix pairs. Then
!> the input (A, D), (B, E) are sub-pencils of the matrix pair in
!> STGSYL. See STGSYL for details.
!>
Parameters
[in]TRANS
!>          TRANS is CHARACTER*1
!>          = 'N': solve the generalized Sylvester equation (1).
!>          = 'T': solve the 'transposed' system (3).
!>
[in]IJOB
!>          IJOB is INTEGER
!>          Specifies what kind of functionality to be performed.
!>          = 0: solve (1) only.
!>          = 1: A contribution from this subsystem to a Frobenius
!>               norm-based estimate of the separation between two matrix
!>               pairs is computed. (look ahead strategy is used).
!>          = 2: A contribution from this subsystem to a Frobenius
!>               norm-based estimate of the separation between two matrix
!>               pairs is computed. (SGECON on sub-systems is used.)
!>          Not referenced if TRANS = 'T'.
!>
[in]M
!>          M is INTEGER
!>          On entry, M specifies the order of A and D, and the row
!>          dimension of C, F, R and L.
!>
[in]N
!>          N is INTEGER
!>          On entry, N specifies the order of B and E, and the column
!>          dimension of C, F, R and L.
!>
[in]A
!>          A is REAL array, dimension (LDA, M)
!>          On entry, A contains an upper quasi triangular matrix.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the matrix A. LDA >= max(1, M).
!>
[in]B
!>          B is REAL array, dimension (LDB, N)
!>          On entry, B contains an upper quasi triangular matrix.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the matrix B. LDB >= max(1, N).
!>
[in,out]C
!>          C is REAL array, dimension (LDC, N)
!>          On entry, C contains the right-hand-side of the first matrix
!>          equation in (1).
!>          On exit, if IJOB = 0, C has been overwritten by the
!>          solution R.
!>
[in]LDC
!>          LDC is INTEGER
!>          The leading dimension of the matrix C. LDC >= max(1, M).
!>
[in]D
!>          D is REAL array, dimension (LDD, M)
!>          On entry, D contains an upper triangular matrix.
!>
[in]LDD
!>          LDD is INTEGER
!>          The leading dimension of the matrix D. LDD >= max(1, M).
!>
[in]E
!>          E is REAL array, dimension (LDE, N)
!>          On entry, E contains an upper triangular matrix.
!>
[in]LDE
!>          LDE is INTEGER
!>          The leading dimension of the matrix E. LDE >= max(1, N).
!>
[in,out]F
!>          F is REAL array, dimension (LDF, N)
!>          On entry, F contains the right-hand-side of the second matrix
!>          equation in (1).
!>          On exit, if IJOB = 0, F has been overwritten by the
!>          solution L.
!>
[in]LDF
!>          LDF is INTEGER
!>          The leading dimension of the matrix F. LDF >= max(1, M).
!>
[out]SCALE
!>          SCALE is REAL
!>          On exit, 0 <= SCALE <= 1. If 0 < SCALE < 1, the solutions
!>          R and L (C and F on entry) will hold the solutions to a
!>          slightly perturbed system but the input matrices A, B, D and
!>          E have not been changed. If SCALE = 0, R and L will hold the
!>          solutions to the homogeneous system with C = F = 0. Normally,
!>          SCALE = 1.
!>
[in,out]RDSUM
!>          RDSUM is REAL
!>          On entry, the sum of squares of computed contributions to
!>          the Dif-estimate under computation by STGSYL, where the
!>          scaling factor RDSCAL (see below) has been factored out.
!>          On exit, the corresponding sum of squares updated with the
!>          contributions from the current sub-system.
!>          If TRANS = 'T' RDSUM is not touched.
!>          NOTE: RDSUM only makes sense when STGSY2 is called by STGSYL.
!>
[in,out]RDSCAL
!>          RDSCAL is REAL
!>          On entry, scaling factor used to prevent overflow in RDSUM.
!>          On exit, RDSCAL is updated w.r.t. the current contributions
!>          in RDSUM.
!>          If TRANS = 'T', RDSCAL is not touched.
!>          NOTE: RDSCAL only makes sense when STGSY2 is called by
!>                STGSYL.
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (M+N+2)
!>
[out]PQ
!>          PQ is INTEGER
!>          On exit, the number of subsystems (of size 2-by-2, 4-by-4 and
!>          8-by-8) solved by this routine.
!>
[out]INFO
!>          INFO is INTEGER
!>          On exit, if INFO is set to
!>            =0: Successful exit
!>            <0: If INFO = -i, the i-th argument had an illegal value.
!>            >0: The matrix pairs (A, D) and (B, E) have common or very
!>                close eigenvalues.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Bo Kagstrom and Peter Poromaa, Department of Computing Science, Umea University, S-901 87 Umea, Sweden.

Definition at line 269 of file stgsy2.f.

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*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
lde
logical function lde(ri, rj, lr)
Definition dblat2.f:2970
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
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
sscal
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
Here is the call graph for this function:
Here is the caller graph for this function: