LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ dtgsy2()

subroutine dtgsy2 ( character  TRANS,
integer  IJOB,
integer  M,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( ldb, * )  B,
integer  LDB,
double precision, dimension( ldc, * )  C,
integer  LDC,
double precision, dimension( ldd, * )  D,
integer  LDD,
double precision, dimension( lde, * )  E,
integer  LDE,
double precision, dimension( ldf, * )  F,
integer  LDF,
double precision  SCALE,
double precision  RDSUM,
double precision  RDSCAL,
integer, dimension( * )  IWORK,
integer  PQ,
integer  INFO 
)

DTGSY2 solves the generalized Sylvester equation (unblocked algorithm).

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

Purpose:
 DTGSY2 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 communicaton with DLACON.

 DTGSY2 also (IJOB >= 1) contributes to the computation in DTGSYL
 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
 DTGSYL. See DTGSYL 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. (DGECON 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
          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 DOUBLE PRECISION
          On entry, the sum of squares of computed contributions to
          the Dif-estimate under computation by DTGSYL, 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 DTGSY2 is called by DTGSYL.
[in,out]RDSCAL
          RDSCAL is DOUBLE PRECISION
          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 DTGSY2 is called by
                DTGSYL.
[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.
Date
December 2016
Contributors:
Bo Kagstrom and Peter Poromaa, Department of Computing Science, Umea University, S-901 87 Umea, Sweden.

Definition at line 276 of file dtgsy2.f.

276 *
277 * -- LAPACK auxiliary routine (version 3.7.0) --
278 * -- LAPACK is a software package provided by Univ. of Tennessee, --
279 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
280 * December 2016
281 *
282 * .. Scalar Arguments ..
283  CHARACTER trans
284  INTEGER ijob, info, lda, ldb, ldc, ldd, lde, ldf, m, n,
285  $ pq
286  DOUBLE PRECISION rdscal, rdsum, scale
287 * ..
288 * .. Array Arguments ..
289  INTEGER iwork( * )
290  DOUBLE PRECISION a( lda, * ), b( ldb, * ), c( ldc, * ),
291  $ d( ldd, * ), e( lde, * ), f( ldf, * )
292 * ..
293 *
294 * =====================================================================
295 * Replaced various illegal calls to DCOPY by calls to DLASET.
296 * Sven Hammarling, 27/5/02.
297 *
298 * .. Parameters ..
299  INTEGER ldz
300  parameter( ldz = 8 )
301  DOUBLE PRECISION zero, one
302  parameter( zero = 0.0d+0, one = 1.0d+0 )
303 * ..
304 * .. Local Scalars ..
305  LOGICAL notran
306  INTEGER i, ie, ierr, ii, is, isp1, j, je, jj, js, jsp1,
307  $ k, mb, nb, p, q, zdim
308  DOUBLE PRECISION alpha, scaloc
309 * ..
310 * .. Local Arrays ..
311  INTEGER ipiv( ldz ), jpiv( ldz )
312  DOUBLE PRECISION rhs( ldz ), z( ldz, ldz )
313 * ..
314 * .. External Functions ..
315  LOGICAL lsame
316  EXTERNAL lsame
317 * ..
318 * .. External Subroutines ..
319  EXTERNAL daxpy, dcopy, dgemm, dgemv, dger, dgesc2,
321 * ..
322 * .. Intrinsic Functions ..
323  INTRINSIC max
324 * ..
325 * .. Executable Statements ..
326 *
327 * Decode and test input parameters
328 *
329  info = 0
330  ierr = 0
331  notran = lsame( trans, 'N' )
332  IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) ) THEN
333  info = -1
334  ELSE IF( notran ) THEN
335  IF( ( ijob.LT.0 ) .OR. ( ijob.GT.2 ) ) THEN
336  info = -2
337  END IF
338  END IF
339  IF( info.EQ.0 ) THEN
340  IF( m.LE.0 ) THEN
341  info = -3
342  ELSE IF( n.LE.0 ) THEN
343  info = -4
344  ELSE IF( lda.LT.max( 1, m ) ) THEN
345  info = -6
346  ELSE IF( ldb.LT.max( 1, n ) ) THEN
347  info = -8
348  ELSE IF( ldc.LT.max( 1, m ) ) THEN
349  info = -10
350  ELSE IF( ldd.LT.max( 1, m ) ) THEN
351  info = -12
352  ELSE IF( lde.LT.max( 1, n ) ) THEN
353  info = -14
354  ELSE IF( ldf.LT.max( 1, m ) ) THEN
355  info = -16
356  END IF
357  END IF
358  IF( info.NE.0 ) THEN
359  CALL xerbla( 'DTGSY2', -info )
360  RETURN
361  END IF
362 *
363 * Determine block structure of A
364 *
365  pq = 0
366  p = 0
367  i = 1
368  10 CONTINUE
369  IF( i.GT.m )
370  $ GO TO 20
371  p = p + 1
372  iwork( p ) = i
373  IF( i.EQ.m )
374  $ GO TO 20
375  IF( a( i+1, i ).NE.zero ) THEN
376  i = i + 2
377  ELSE
378  i = i + 1
379  END IF
380  GO TO 10
381  20 CONTINUE
382  iwork( p+1 ) = m + 1
383 *
384 * Determine block structure of B
385 *
386  q = p + 1
387  j = 1
388  30 CONTINUE
389  IF( j.GT.n )
390  $ GO TO 40
391  q = q + 1
392  iwork( q ) = j
393  IF( j.EQ.n )
394  $ GO TO 40
395  IF( b( j+1, j ).NE.zero ) THEN
396  j = j + 2
397  ELSE
398  j = j + 1
399  END IF
400  GO TO 30
401  40 CONTINUE
402  iwork( q+1 ) = n + 1
403  pq = p*( q-p-1 )
404 *
405  IF( notran ) THEN
406 *
407 * Solve (I, J) - subsystem
408 * A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
409 * D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
410 * for I = P, P - 1, ..., 1; J = 1, 2, ..., Q
411 *
412  scale = one
413  scaloc = one
414  DO 120 j = p + 2, q
415  js = iwork( j )
416  jsp1 = js + 1
417  je = iwork( j+1 ) - 1
418  nb = je - js + 1
419  DO 110 i = p, 1, -1
420 *
421  is = iwork( i )
422  isp1 = is + 1
423  ie = iwork( i+1 ) - 1
424  mb = ie - is + 1
425  zdim = mb*nb*2
426 *
427  IF( ( mb.EQ.1 ) .AND. ( nb.EQ.1 ) ) THEN
428 *
429 * Build a 2-by-2 system Z * x = RHS
430 *
431  z( 1, 1 ) = a( is, is )
432  z( 2, 1 ) = d( is, is )
433  z( 1, 2 ) = -b( js, js )
434  z( 2, 2 ) = -e( js, js )
435 *
436 * Set up right hand side(s)
437 *
438  rhs( 1 ) = c( is, js )
439  rhs( 2 ) = f( is, js )
440 *
441 * Solve Z * x = RHS
442 *
443  CALL dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
444  IF( ierr.GT.0 )
445  $ info = ierr
446 *
447  IF( ijob.EQ.0 ) THEN
448  CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
449  $ scaloc )
450  IF( scaloc.NE.one ) THEN
451  DO 50 k = 1, n
452  CALL dscal( m, scaloc, c( 1, k ), 1 )
453  CALL dscal( m, scaloc, f( 1, k ), 1 )
454  50 CONTINUE
455  scale = scale*scaloc
456  END IF
457  ELSE
458  CALL dlatdf( ijob, zdim, z, ldz, rhs, rdsum,
459  $ rdscal, ipiv, jpiv )
460  END IF
461 *
462 * Unpack solution vector(s)
463 *
464  c( is, js ) = rhs( 1 )
465  f( is, js ) = rhs( 2 )
466 *
467 * Substitute R(I, J) and L(I, J) into remaining
468 * equation.
469 *
470  IF( i.GT.1 ) THEN
471  alpha = -rhs( 1 )
472  CALL daxpy( is-1, alpha, a( 1, is ), 1, c( 1, js ),
473  $ 1 )
474  CALL daxpy( is-1, alpha, d( 1, is ), 1, f( 1, js ),
475  $ 1 )
476  END IF
477  IF( j.LT.q ) THEN
478  CALL daxpy( n-je, rhs( 2 ), b( js, je+1 ), ldb,
479  $ c( is, je+1 ), ldc )
480  CALL daxpy( n-je, rhs( 2 ), e( js, je+1 ), lde,
481  $ f( is, je+1 ), ldf )
482  END IF
483 *
484  ELSE IF( ( mb.EQ.1 ) .AND. ( nb.EQ.2 ) ) THEN
485 *
486 * Build a 4-by-4 system Z * x = RHS
487 *
488  z( 1, 1 ) = a( is, is )
489  z( 2, 1 ) = zero
490  z( 3, 1 ) = d( is, is )
491  z( 4, 1 ) = zero
492 *
493  z( 1, 2 ) = zero
494  z( 2, 2 ) = a( is, is )
495  z( 3, 2 ) = zero
496  z( 4, 2 ) = d( is, is )
497 *
498  z( 1, 3 ) = -b( js, js )
499  z( 2, 3 ) = -b( js, jsp1 )
500  z( 3, 3 ) = -e( js, js )
501  z( 4, 3 ) = -e( js, jsp1 )
502 *
503  z( 1, 4 ) = -b( jsp1, js )
504  z( 2, 4 ) = -b( jsp1, jsp1 )
505  z( 3, 4 ) = zero
506  z( 4, 4 ) = -e( jsp1, jsp1 )
507 *
508 * Set up right hand side(s)
509 *
510  rhs( 1 ) = c( is, js )
511  rhs( 2 ) = c( is, jsp1 )
512  rhs( 3 ) = f( is, js )
513  rhs( 4 ) = f( is, jsp1 )
514 *
515 * Solve Z * x = RHS
516 *
517  CALL dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
518  IF( ierr.GT.0 )
519  $ info = ierr
520 *
521  IF( ijob.EQ.0 ) THEN
522  CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
523  $ scaloc )
524  IF( scaloc.NE.one ) THEN
525  DO 60 k = 1, n
526  CALL dscal( m, scaloc, c( 1, k ), 1 )
527  CALL dscal( m, scaloc, f( 1, k ), 1 )
528  60 CONTINUE
529  scale = scale*scaloc
530  END IF
531  ELSE
532  CALL dlatdf( ijob, zdim, z, ldz, rhs, rdsum,
533  $ rdscal, ipiv, jpiv )
534  END IF
535 *
536 * Unpack solution vector(s)
537 *
538  c( is, js ) = rhs( 1 )
539  c( is, jsp1 ) = rhs( 2 )
540  f( is, js ) = rhs( 3 )
541  f( is, jsp1 ) = rhs( 4 )
542 *
543 * Substitute R(I, J) and L(I, J) into remaining
544 * equation.
545 *
546  IF( i.GT.1 ) THEN
547  CALL dger( is-1, nb, -one, a( 1, is ), 1, rhs( 1 ),
548  $ 1, c( 1, js ), ldc )
549  CALL dger( is-1, nb, -one, d( 1, is ), 1, rhs( 1 ),
550  $ 1, f( 1, js ), ldf )
551  END IF
552  IF( j.LT.q ) THEN
553  CALL daxpy( n-je, rhs( 3 ), b( js, je+1 ), ldb,
554  $ c( is, je+1 ), ldc )
555  CALL daxpy( n-je, rhs( 3 ), e( js, je+1 ), lde,
556  $ f( is, je+1 ), ldf )
557  CALL daxpy( n-je, rhs( 4 ), b( jsp1, je+1 ), ldb,
558  $ c( is, je+1 ), ldc )
559  CALL daxpy( n-je, rhs( 4 ), e( jsp1, je+1 ), lde,
560  $ f( is, je+1 ), ldf )
561  END IF
562 *
563  ELSE IF( ( mb.EQ.2 ) .AND. ( nb.EQ.1 ) ) THEN
564 *
565 * Build a 4-by-4 system Z * x = RHS
566 *
567  z( 1, 1 ) = a( is, is )
568  z( 2, 1 ) = a( isp1, is )
569  z( 3, 1 ) = d( is, is )
570  z( 4, 1 ) = zero
571 *
572  z( 1, 2 ) = a( is, isp1 )
573  z( 2, 2 ) = a( isp1, isp1 )
574  z( 3, 2 ) = d( is, isp1 )
575  z( 4, 2 ) = d( isp1, isp1 )
576 *
577  z( 1, 3 ) = -b( js, js )
578  z( 2, 3 ) = zero
579  z( 3, 3 ) = -e( js, js )
580  z( 4, 3 ) = zero
581 *
582  z( 1, 4 ) = zero
583  z( 2, 4 ) = -b( js, js )
584  z( 3, 4 ) = zero
585  z( 4, 4 ) = -e( js, js )
586 *
587 * Set up right hand side(s)
588 *
589  rhs( 1 ) = c( is, js )
590  rhs( 2 ) = c( isp1, js )
591  rhs( 3 ) = f( is, js )
592  rhs( 4 ) = f( isp1, js )
593 *
594 * Solve Z * x = RHS
595 *
596  CALL dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
597  IF( ierr.GT.0 )
598  $ info = ierr
599  IF( ijob.EQ.0 ) THEN
600  CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
601  $ scaloc )
602  IF( scaloc.NE.one ) THEN
603  DO 70 k = 1, n
604  CALL dscal( m, scaloc, c( 1, k ), 1 )
605  CALL dscal( m, scaloc, f( 1, k ), 1 )
606  70 CONTINUE
607  scale = scale*scaloc
608  END IF
609  ELSE
610  CALL dlatdf( ijob, zdim, z, ldz, rhs, rdsum,
611  $ rdscal, ipiv, jpiv )
612  END IF
613 *
614 * Unpack solution vector(s)
615 *
616  c( is, js ) = rhs( 1 )
617  c( isp1, js ) = rhs( 2 )
618  f( is, js ) = rhs( 3 )
619  f( isp1, js ) = rhs( 4 )
620 *
621 * Substitute R(I, J) and L(I, J) into remaining
622 * equation.
623 *
624  IF( i.GT.1 ) THEN
625  CALL dgemv( 'N', is-1, mb, -one, a( 1, is ), lda,
626  $ rhs( 1 ), 1, one, c( 1, js ), 1 )
627  CALL dgemv( 'N', is-1, mb, -one, d( 1, is ), ldd,
628  $ rhs( 1 ), 1, one, f( 1, js ), 1 )
629  END IF
630  IF( j.LT.q ) THEN
631  CALL dger( mb, n-je, one, rhs( 3 ), 1,
632  $ b( js, je+1 ), ldb, c( is, je+1 ), ldc )
633  CALL dger( mb, n-je, one, rhs( 3 ), 1,
634  $ e( js, je+1 ), lde, f( is, je+1 ), ldf )
635  END IF
636 *
637  ELSE IF( ( mb.EQ.2 ) .AND. ( nb.EQ.2 ) ) THEN
638 *
639 * Build an 8-by-8 system Z * x = RHS
640 *
641  CALL dlaset( 'F', ldz, ldz, zero, zero, z, ldz )
642 *
643  z( 1, 1 ) = a( is, is )
644  z( 2, 1 ) = a( isp1, is )
645  z( 5, 1 ) = d( is, is )
646 *
647  z( 1, 2 ) = a( is, isp1 )
648  z( 2, 2 ) = a( isp1, isp1 )
649  z( 5, 2 ) = d( is, isp1 )
650  z( 6, 2 ) = d( isp1, isp1 )
651 *
652  z( 3, 3 ) = a( is, is )
653  z( 4, 3 ) = a( isp1, is )
654  z( 7, 3 ) = d( is, is )
655 *
656  z( 3, 4 ) = a( is, isp1 )
657  z( 4, 4 ) = a( isp1, isp1 )
658  z( 7, 4 ) = d( is, isp1 )
659  z( 8, 4 ) = d( isp1, isp1 )
660 *
661  z( 1, 5 ) = -b( js, js )
662  z( 3, 5 ) = -b( js, jsp1 )
663  z( 5, 5 ) = -e( js, js )
664  z( 7, 5 ) = -e( js, jsp1 )
665 *
666  z( 2, 6 ) = -b( js, js )
667  z( 4, 6 ) = -b( js, jsp1 )
668  z( 6, 6 ) = -e( js, js )
669  z( 8, 6 ) = -e( js, jsp1 )
670 *
671  z( 1, 7 ) = -b( jsp1, js )
672  z( 3, 7 ) = -b( jsp1, jsp1 )
673  z( 7, 7 ) = -e( jsp1, jsp1 )
674 *
675  z( 2, 8 ) = -b( jsp1, js )
676  z( 4, 8 ) = -b( jsp1, jsp1 )
677  z( 8, 8 ) = -e( jsp1, jsp1 )
678 *
679 * Set up right hand side(s)
680 *
681  k = 1
682  ii = mb*nb + 1
683  DO 80 jj = 0, nb - 1
684  CALL dcopy( mb, c( is, js+jj ), 1, rhs( k ), 1 )
685  CALL dcopy( mb, f( is, js+jj ), 1, rhs( ii ), 1 )
686  k = k + mb
687  ii = ii + mb
688  80 CONTINUE
689 *
690 * Solve Z * x = RHS
691 *
692  CALL dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
693  IF( ierr.GT.0 )
694  $ info = ierr
695  IF( ijob.EQ.0 ) THEN
696  CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
697  $ scaloc )
698  IF( scaloc.NE.one ) THEN
699  DO 90 k = 1, n
700  CALL dscal( m, scaloc, c( 1, k ), 1 )
701  CALL dscal( m, scaloc, f( 1, k ), 1 )
702  90 CONTINUE
703  scale = scale*scaloc
704  END IF
705  ELSE
706  CALL dlatdf( ijob, zdim, z, ldz, rhs, rdsum,
707  $ rdscal, ipiv, jpiv )
708  END IF
709 *
710 * Unpack solution vector(s)
711 *
712  k = 1
713  ii = mb*nb + 1
714  DO 100 jj = 0, nb - 1
715  CALL dcopy( mb, rhs( k ), 1, c( is, js+jj ), 1 )
716  CALL dcopy( mb, rhs( ii ), 1, f( is, js+jj ), 1 )
717  k = k + mb
718  ii = ii + mb
719  100 CONTINUE
720 *
721 * Substitute R(I, J) and L(I, J) into remaining
722 * equation.
723 *
724  IF( i.GT.1 ) THEN
725  CALL dgemm( 'N', 'N', is-1, nb, mb, -one,
726  $ a( 1, is ), lda, rhs( 1 ), mb, one,
727  $ c( 1, js ), ldc )
728  CALL dgemm( 'N', 'N', is-1, nb, mb, -one,
729  $ d( 1, is ), ldd, rhs( 1 ), mb, one,
730  $ f( 1, js ), ldf )
731  END IF
732  IF( j.LT.q ) THEN
733  k = mb*nb + 1
734  CALL dgemm( 'N', 'N', mb, n-je, nb, one, rhs( k ),
735  $ mb, b( js, je+1 ), ldb, one,
736  $ c( is, je+1 ), ldc )
737  CALL dgemm( 'N', 'N', mb, n-je, nb, one, rhs( k ),
738  $ mb, e( js, je+1 ), lde, one,
739  $ f( is, je+1 ), ldf )
740  END IF
741 *
742  END IF
743 *
744  110 CONTINUE
745  120 CONTINUE
746  ELSE
747 *
748 * Solve (I, J) - subsystem
749 * A(I, I)**T * R(I, J) + D(I, I)**T * L(J, J) = C(I, J)
750 * R(I, I) * B(J, J) + L(I, J) * E(J, J) = -F(I, J)
751 * for I = 1, 2, ..., P, J = Q, Q - 1, ..., 1
752 *
753  scale = one
754  scaloc = one
755  DO 200 i = 1, p
756 *
757  is = iwork( i )
758  isp1 = is + 1
759  ie = iwork( i+1 ) - 1
760  mb = ie - is + 1
761  DO 190 j = q, p + 2, -1
762 *
763  js = iwork( j )
764  jsp1 = js + 1
765  je = iwork( j+1 ) - 1
766  nb = je - js + 1
767  zdim = mb*nb*2
768  IF( ( mb.EQ.1 ) .AND. ( nb.EQ.1 ) ) THEN
769 *
770 * Build a 2-by-2 system Z**T * x = RHS
771 *
772  z( 1, 1 ) = a( is, is )
773  z( 2, 1 ) = -b( js, js )
774  z( 1, 2 ) = d( is, is )
775  z( 2, 2 ) = -e( js, js )
776 *
777 * Set up right hand side(s)
778 *
779  rhs( 1 ) = c( is, js )
780  rhs( 2 ) = f( is, js )
781 *
782 * Solve Z**T * x = RHS
783 *
784  CALL dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
785  IF( ierr.GT.0 )
786  $ info = ierr
787 *
788  CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
789  IF( scaloc.NE.one ) THEN
790  DO 130 k = 1, n
791  CALL dscal( m, scaloc, c( 1, k ), 1 )
792  CALL dscal( m, scaloc, f( 1, k ), 1 )
793  130 CONTINUE
794  scale = scale*scaloc
795  END IF
796 *
797 * Unpack solution vector(s)
798 *
799  c( is, js ) = rhs( 1 )
800  f( is, js ) = rhs( 2 )
801 *
802 * Substitute R(I, J) and L(I, J) into remaining
803 * equation.
804 *
805  IF( j.GT.p+2 ) THEN
806  alpha = rhs( 1 )
807  CALL daxpy( js-1, alpha, b( 1, js ), 1, f( is, 1 ),
808  $ ldf )
809  alpha = rhs( 2 )
810  CALL daxpy( js-1, alpha, e( 1, js ), 1, f( is, 1 ),
811  $ ldf )
812  END IF
813  IF( i.LT.p ) THEN
814  alpha = -rhs( 1 )
815  CALL daxpy( m-ie, alpha, a( is, ie+1 ), lda,
816  $ c( ie+1, js ), 1 )
817  alpha = -rhs( 2 )
818  CALL daxpy( m-ie, alpha, d( is, ie+1 ), ldd,
819  $ c( ie+1, js ), 1 )
820  END IF
821 *
822  ELSE IF( ( mb.EQ.1 ) .AND. ( nb.EQ.2 ) ) THEN
823 *
824 * Build a 4-by-4 system Z**T * x = RHS
825 *
826  z( 1, 1 ) = a( is, is )
827  z( 2, 1 ) = zero
828  z( 3, 1 ) = -b( js, js )
829  z( 4, 1 ) = -b( jsp1, js )
830 *
831  z( 1, 2 ) = zero
832  z( 2, 2 ) = a( is, is )
833  z( 3, 2 ) = -b( js, jsp1 )
834  z( 4, 2 ) = -b( jsp1, jsp1 )
835 *
836  z( 1, 3 ) = d( is, is )
837  z( 2, 3 ) = zero
838  z( 3, 3 ) = -e( js, js )
839  z( 4, 3 ) = zero
840 *
841  z( 1, 4 ) = zero
842  z( 2, 4 ) = d( is, is )
843  z( 3, 4 ) = -e( js, jsp1 )
844  z( 4, 4 ) = -e( jsp1, jsp1 )
845 *
846 * Set up right hand side(s)
847 *
848  rhs( 1 ) = c( is, js )
849  rhs( 2 ) = c( is, jsp1 )
850  rhs( 3 ) = f( is, js )
851  rhs( 4 ) = f( is, jsp1 )
852 *
853 * Solve Z**T * x = RHS
854 *
855  CALL dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
856  IF( ierr.GT.0 )
857  $ info = ierr
858  CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
859  IF( scaloc.NE.one ) THEN
860  DO 140 k = 1, n
861  CALL dscal( m, scaloc, c( 1, k ), 1 )
862  CALL dscal( m, scaloc, f( 1, k ), 1 )
863  140 CONTINUE
864  scale = scale*scaloc
865  END IF
866 *
867 * Unpack solution vector(s)
868 *
869  c( is, js ) = rhs( 1 )
870  c( is, jsp1 ) = rhs( 2 )
871  f( is, js ) = rhs( 3 )
872  f( is, jsp1 ) = rhs( 4 )
873 *
874 * Substitute R(I, J) and L(I, J) into remaining
875 * equation.
876 *
877  IF( j.GT.p+2 ) THEN
878  CALL daxpy( js-1, rhs( 1 ), b( 1, js ), 1,
879  $ f( is, 1 ), ldf )
880  CALL daxpy( js-1, rhs( 2 ), b( 1, jsp1 ), 1,
881  $ f( is, 1 ), ldf )
882  CALL daxpy( js-1, rhs( 3 ), e( 1, js ), 1,
883  $ f( is, 1 ), ldf )
884  CALL daxpy( js-1, rhs( 4 ), e( 1, jsp1 ), 1,
885  $ f( is, 1 ), ldf )
886  END IF
887  IF( i.LT.p ) THEN
888  CALL dger( m-ie, nb, -one, a( is, ie+1 ), lda,
889  $ rhs( 1 ), 1, c( ie+1, js ), ldc )
890  CALL dger( m-ie, nb, -one, d( is, ie+1 ), ldd,
891  $ rhs( 3 ), 1, c( ie+1, js ), ldc )
892  END IF
893 *
894  ELSE IF( ( mb.EQ.2 ) .AND. ( nb.EQ.1 ) ) THEN
895 *
896 * Build a 4-by-4 system Z**T * x = RHS
897 *
898  z( 1, 1 ) = a( is, is )
899  z( 2, 1 ) = a( is, isp1 )
900  z( 3, 1 ) = -b( js, js )
901  z( 4, 1 ) = zero
902 *
903  z( 1, 2 ) = a( isp1, is )
904  z( 2, 2 ) = a( isp1, isp1 )
905  z( 3, 2 ) = zero
906  z( 4, 2 ) = -b( js, js )
907 *
908  z( 1, 3 ) = d( is, is )
909  z( 2, 3 ) = d( is, isp1 )
910  z( 3, 3 ) = -e( js, js )
911  z( 4, 3 ) = zero
912 *
913  z( 1, 4 ) = zero
914  z( 2, 4 ) = d( isp1, isp1 )
915  z( 3, 4 ) = zero
916  z( 4, 4 ) = -e( js, js )
917 *
918 * Set up right hand side(s)
919 *
920  rhs( 1 ) = c( is, js )
921  rhs( 2 ) = c( isp1, js )
922  rhs( 3 ) = f( is, js )
923  rhs( 4 ) = f( isp1, js )
924 *
925 * Solve Z**T * x = RHS
926 *
927  CALL dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
928  IF( ierr.GT.0 )
929  $ info = ierr
930 *
931  CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
932  IF( scaloc.NE.one ) THEN
933  DO 150 k = 1, n
934  CALL dscal( m, scaloc, c( 1, k ), 1 )
935  CALL dscal( m, scaloc, f( 1, k ), 1 )
936  150 CONTINUE
937  scale = scale*scaloc
938  END IF
939 *
940 * Unpack solution vector(s)
941 *
942  c( is, js ) = rhs( 1 )
943  c( isp1, js ) = rhs( 2 )
944  f( is, js ) = rhs( 3 )
945  f( isp1, js ) = rhs( 4 )
946 *
947 * Substitute R(I, J) and L(I, J) into remaining
948 * equation.
949 *
950  IF( j.GT.p+2 ) THEN
951  CALL dger( mb, js-1, one, rhs( 1 ), 1, b( 1, js ),
952  $ 1, f( is, 1 ), ldf )
953  CALL dger( mb, js-1, one, rhs( 3 ), 1, e( 1, js ),
954  $ 1, f( is, 1 ), ldf )
955  END IF
956  IF( i.LT.p ) THEN
957  CALL dgemv( 'T', mb, m-ie, -one, a( is, ie+1 ),
958  $ lda, rhs( 1 ), 1, one, c( ie+1, js ),
959  $ 1 )
960  CALL dgemv( 'T', mb, m-ie, -one, d( is, ie+1 ),
961  $ ldd, rhs( 3 ), 1, one, c( ie+1, js ),
962  $ 1 )
963  END IF
964 *
965  ELSE IF( ( mb.EQ.2 ) .AND. ( nb.EQ.2 ) ) THEN
966 *
967 * Build an 8-by-8 system Z**T * x = RHS
968 *
969  CALL dlaset( 'F', ldz, ldz, zero, zero, z, ldz )
970 *
971  z( 1, 1 ) = a( is, is )
972  z( 2, 1 ) = a( is, isp1 )
973  z( 5, 1 ) = -b( js, js )
974  z( 7, 1 ) = -b( jsp1, js )
975 *
976  z( 1, 2 ) = a( isp1, is )
977  z( 2, 2 ) = a( isp1, isp1 )
978  z( 6, 2 ) = -b( js, js )
979  z( 8, 2 ) = -b( jsp1, js )
980 *
981  z( 3, 3 ) = a( is, is )
982  z( 4, 3 ) = a( is, isp1 )
983  z( 5, 3 ) = -b( js, jsp1 )
984  z( 7, 3 ) = -b( jsp1, jsp1 )
985 *
986  z( 3, 4 ) = a( isp1, is )
987  z( 4, 4 ) = a( isp1, isp1 )
988  z( 6, 4 ) = -b( js, jsp1 )
989  z( 8, 4 ) = -b( jsp1, jsp1 )
990 *
991  z( 1, 5 ) = d( is, is )
992  z( 2, 5 ) = d( is, isp1 )
993  z( 5, 5 ) = -e( js, js )
994 *
995  z( 2, 6 ) = d( isp1, isp1 )
996  z( 6, 6 ) = -e( js, js )
997 *
998  z( 3, 7 ) = d( is, is )
999  z( 4, 7 ) = d( is, isp1 )
1000  z( 5, 7 ) = -e( js, jsp1 )
1001  z( 7, 7 ) = -e( jsp1, jsp1 )
1002 *
1003  z( 4, 8 ) = d( isp1, isp1 )
1004  z( 6, 8 ) = -e( js, jsp1 )
1005  z( 8, 8 ) = -e( jsp1, jsp1 )
1006 *
1007 * Set up right hand side(s)
1008 *
1009  k = 1
1010  ii = mb*nb + 1
1011  DO 160 jj = 0, nb - 1
1012  CALL dcopy( mb, c( is, js+jj ), 1, rhs( k ), 1 )
1013  CALL dcopy( mb, f( is, js+jj ), 1, rhs( ii ), 1 )
1014  k = k + mb
1015  ii = ii + mb
1016  160 CONTINUE
1017 *
1018 *
1019 * Solve Z**T * x = RHS
1020 *
1021  CALL dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
1022  IF( ierr.GT.0 )
1023  $ info = ierr
1024 *
1025  CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
1026  IF( scaloc.NE.one ) THEN
1027  DO 170 k = 1, n
1028  CALL dscal( m, scaloc, c( 1, k ), 1 )
1029  CALL dscal( m, scaloc, f( 1, k ), 1 )
1030  170 CONTINUE
1031  scale = scale*scaloc
1032  END IF
1033 *
1034 * Unpack solution vector(s)
1035 *
1036  k = 1
1037  ii = mb*nb + 1
1038  DO 180 jj = 0, nb - 1
1039  CALL dcopy( mb, rhs( k ), 1, c( is, js+jj ), 1 )
1040  CALL dcopy( mb, rhs( ii ), 1, f( is, js+jj ), 1 )
1041  k = k + mb
1042  ii = ii + mb
1043  180 CONTINUE
1044 *
1045 * Substitute R(I, J) and L(I, J) into remaining
1046 * equation.
1047 *
1048  IF( j.GT.p+2 ) THEN
1049  CALL dgemm( 'N', 'T', mb, js-1, nb, one,
1050  $ c( is, js ), ldc, b( 1, js ), ldb, one,
1051  $ f( is, 1 ), ldf )
1052  CALL dgemm( 'N', 'T', mb, js-1, nb, one,
1053  $ f( is, js ), ldf, e( 1, js ), lde, one,
1054  $ f( is, 1 ), ldf )
1055  END IF
1056  IF( i.LT.p ) THEN
1057  CALL dgemm( 'T', 'N', m-ie, nb, mb, -one,
1058  $ a( is, ie+1 ), lda, c( is, js ), ldc,
1059  $ one, c( ie+1, js ), ldc )
1060  CALL dgemm( 'T', 'N', m-ie, nb, mb, -one,
1061  $ d( is, ie+1 ), ldd, f( is, js ), ldf,
1062  $ one, c( ie+1, js ), ldc )
1063  END IF
1064 *
1065  END IF
1066 *
1067  190 CONTINUE
1068  200 CONTINUE
1069 *
1070  END IF
1071  RETURN
1072 *
1073 * End of DTGSY2
1074 *
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:84
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:158
logical function lde(RI, RJ, LR)
Definition: dblat2.f:2945
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:91
subroutine dgesc2(N, A, LDA, RHS, IPIV, JPIV, SCALE)
DGESC2 solves a system of linear equations using the LU factorization with complete pivoting computed...
Definition: dgesc2.f:116
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: dlaset.f:112
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:189
subroutine dger(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
DGER
Definition: dger.f:132
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:81
subroutine dlatdf(IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV, JPIV)
DLATDF uses the LU factorization of the n-by-n matrix computed by sgetc2 and computes a contribution ...
Definition: dlatdf.f:173
subroutine dgetc2(N, A, LDA, IPIV, JPIV, INFO)
DGETC2 computes the LU factorization with complete pivoting of the general n-by-n matrix...
Definition: dgetc2.f:113
Here is the call graph for this function:
Here is the caller graph for this function: