LAPACK  3.10.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 communication 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.
Contributors:
Bo Kagstrom and Peter Poromaa, Department of Computing Science, Umea University, S-901 87 Umea, Sweden.

Definition at line 271 of file dtgsy2.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  DOUBLE PRECISION RDSCAL, RDSUM, SCALE
284 * ..
285 * .. Array Arguments ..
286  INTEGER IWORK( * )
287  DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ),
288  $ D( LDD, * ), E( LDE, * ), F( LDF, * )
289 * ..
290 *
291 * =====================================================================
292 * Replaced various illegal calls to DCOPY by calls to DLASET.
293 * Sven Hammarling, 27/5/02.
294 *
295 * .. Parameters ..
296  INTEGER LDZ
297  parameter( ldz = 8 )
298  DOUBLE PRECISION ZERO, ONE
299  parameter( zero = 0.0d+0, one = 1.0d+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  DOUBLE PRECISION ALPHA, SCALOC
306 * ..
307 * .. Local Arrays ..
308  INTEGER IPIV( LDZ ), JPIV( LDZ )
309  DOUBLE PRECISION RHS( LDZ ), Z( LDZ, LDZ )
310 * ..
311 * .. External Functions ..
312  LOGICAL LSAME
313  EXTERNAL lsame
314 * ..
315 * .. External Subroutines ..
316  EXTERNAL daxpy, dcopy, dgemm, dgemv, dger, dgesc2,
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( 'DTGSY2', -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 dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
441  IF( ierr.GT.0 )
442  $ info = ierr
443 *
444  IF( ijob.EQ.0 ) THEN
445  CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
446  $ scaloc )
447  IF( scaloc.NE.one ) THEN
448  DO 50 k = 1, n
449  CALL dscal( m, scaloc, c( 1, k ), 1 )
450  CALL dscal( m, scaloc, f( 1, k ), 1 )
451  50 CONTINUE
452  scale = scale*scaloc
453  END IF
454  ELSE
455  CALL dlatdf( 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 daxpy( is-1, alpha, a( 1, is ), 1, c( 1, js ),
470  $ 1 )
471  CALL daxpy( is-1, alpha, d( 1, is ), 1, f( 1, js ),
472  $ 1 )
473  END IF
474  IF( j.LT.q ) THEN
475  CALL daxpy( n-je, rhs( 2 ), b( js, je+1 ), ldb,
476  $ c( is, je+1 ), ldc )
477  CALL daxpy( 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 dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
515  IF( ierr.GT.0 )
516  $ info = ierr
517 *
518  IF( ijob.EQ.0 ) THEN
519  CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
520  $ scaloc )
521  IF( scaloc.NE.one ) THEN
522  DO 60 k = 1, n
523  CALL dscal( m, scaloc, c( 1, k ), 1 )
524  CALL dscal( m, scaloc, f( 1, k ), 1 )
525  60 CONTINUE
526  scale = scale*scaloc
527  END IF
528  ELSE
529  CALL dlatdf( 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 dger( is-1, nb, -one, a( 1, is ), 1, rhs( 1 ),
545  $ 1, c( 1, js ), ldc )
546  CALL dger( 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 daxpy( n-je, rhs( 3 ), b( js, je+1 ), ldb,
551  $ c( is, je+1 ), ldc )
552  CALL daxpy( n-je, rhs( 3 ), e( js, je+1 ), lde,
553  $ f( is, je+1 ), ldf )
554  CALL daxpy( n-je, rhs( 4 ), b( jsp1, je+1 ), ldb,
555  $ c( is, je+1 ), ldc )
556  CALL daxpy( 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 dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
594  IF( ierr.GT.0 )
595  $ info = ierr
596  IF( ijob.EQ.0 ) THEN
597  CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
598  $ scaloc )
599  IF( scaloc.NE.one ) THEN
600  DO 70 k = 1, n
601  CALL dscal( m, scaloc, c( 1, k ), 1 )
602  CALL dscal( m, scaloc, f( 1, k ), 1 )
603  70 CONTINUE
604  scale = scale*scaloc
605  END IF
606  ELSE
607  CALL dlatdf( 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 dgemv( 'N', is-1, mb, -one, a( 1, is ), lda,
623  $ rhs( 1 ), 1, one, c( 1, js ), 1 )
624  CALL dgemv( '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 dger( mb, n-je, one, rhs( 3 ), 1,
629  $ b( js, je+1 ), ldb, c( is, je+1 ), ldc )
630  CALL dger( 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 dlaset( '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 dcopy( mb, c( is, js+jj ), 1, rhs( k ), 1 )
682  CALL dcopy( 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 dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
690  IF( ierr.GT.0 )
691  $ info = ierr
692  IF( ijob.EQ.0 ) THEN
693  CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
694  $ scaloc )
695  IF( scaloc.NE.one ) THEN
696  DO 90 k = 1, n
697  CALL dscal( m, scaloc, c( 1, k ), 1 )
698  CALL dscal( m, scaloc, f( 1, k ), 1 )
699  90 CONTINUE
700  scale = scale*scaloc
701  END IF
702  ELSE
703  CALL dlatdf( 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 dcopy( mb, rhs( k ), 1, c( is, js+jj ), 1 )
713  CALL dcopy( 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 dgemm( 'N', 'N', is-1, nb, mb, -one,
723  $ a( 1, is ), lda, rhs( 1 ), mb, one,
724  $ c( 1, js ), ldc )
725  CALL dgemm( '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 dgemm( '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 dgemm( '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 dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
782  IF( ierr.GT.0 )
783  $ info = ierr
784 *
785  CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
786  IF( scaloc.NE.one ) THEN
787  DO 130 k = 1, n
788  CALL dscal( m, scaloc, c( 1, k ), 1 )
789  CALL dscal( 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 daxpy( js-1, alpha, b( 1, js ), 1, f( is, 1 ),
805  $ ldf )
806  alpha = rhs( 2 )
807  CALL daxpy( 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 daxpy( m-ie, alpha, a( is, ie+1 ), lda,
813  $ c( ie+1, js ), 1 )
814  alpha = -rhs( 2 )
815  CALL daxpy( 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 dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
853  IF( ierr.GT.0 )
854  $ info = ierr
855  CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
856  IF( scaloc.NE.one ) THEN
857  DO 140 k = 1, n
858  CALL dscal( m, scaloc, c( 1, k ), 1 )
859  CALL dscal( 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 daxpy( js-1, rhs( 1 ), b( 1, js ), 1,
876  $ f( is, 1 ), ldf )
877  CALL daxpy( js-1, rhs( 2 ), b( 1, jsp1 ), 1,
878  $ f( is, 1 ), ldf )
879  CALL daxpy( js-1, rhs( 3 ), e( 1, js ), 1,
880  $ f( is, 1 ), ldf )
881  CALL daxpy( js-1, rhs( 4 ), e( 1, jsp1 ), 1,
882  $ f( is, 1 ), ldf )
883  END IF
884  IF( i.LT.p ) THEN
885  CALL dger( m-ie, nb, -one, a( is, ie+1 ), lda,
886  $ rhs( 1 ), 1, c( ie+1, js ), ldc )
887  CALL dger( 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 dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
925  IF( ierr.GT.0 )
926  $ info = ierr
927 *
928  CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
929  IF( scaloc.NE.one ) THEN
930  DO 150 k = 1, n
931  CALL dscal( m, scaloc, c( 1, k ), 1 )
932  CALL dscal( 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 dger( mb, js-1, one, rhs( 1 ), 1, b( 1, js ),
949  $ 1, f( is, 1 ), ldf )
950  CALL dger( 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 dgemv( 'T', mb, m-ie, -one, a( is, ie+1 ),
955  $ lda, rhs( 1 ), 1, one, c( ie+1, js ),
956  $ 1 )
957  CALL dgemv( '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 dlaset( '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 dcopy( mb, c( is, js+jj ), 1, rhs( k ), 1 )
1010  CALL dcopy( 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 dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
1019  IF( ierr.GT.0 )
1020  $ info = ierr
1021 *
1022  CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
1023  IF( scaloc.NE.one ) THEN
1024  DO 170 k = 1, n
1025  CALL dscal( m, scaloc, c( 1, k ), 1 )
1026  CALL dscal( 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 dcopy( mb, rhs( k ), 1, c( is, js+jj ), 1 )
1037  CALL dcopy( 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 dgemm( 'N', 'T', mb, js-1, nb, one,
1047  $ c( is, js ), ldc, b( 1, js ), ldb, one,
1048  $ f( is, 1 ), ldf )
1049  CALL dgemm( '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 dgemm( '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 dgemm( '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 DTGSY2
1071 *
logical function lde(RI, RJ, LR)
Definition: dblat2.f:2942
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:110
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:89
subroutine dger(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
DGER
Definition: dger.f:130
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:156
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187
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:111
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:114
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:171
Here is the call graph for this function:
Here is the caller graph for this function: