LAPACK 3.12.0
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 271 of file stgsy2.f.

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*
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: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
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: