LAPACK  3.6.0
LAPACK: Linear Algebra PACKage
Collaboration diagram for double:

Functions

subroutine dgesc2 (N, A, LDA, RHS, IPIV, JPIV, SCALE)
 DGESC2 solves a system of linear equations using the LU factorization with complete pivoting computed by sgetc2. More...
 
subroutine dgetc2 (N, A, LDA, IPIV, JPIV, INFO)
 DGETC2 computes the LU factorization with complete pivoting of the general n-by-n matrix. More...
 
double precision function dlange (NORM, M, N, A, LDA, WORK)
 DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value of any element of a general rectangular matrix. More...
 
subroutine dlaqge (M, N, A, LDA, R, C, ROWCND, COLCND, AMAX, EQUED)
 DLAQGE scales a general rectangular matrix, using row and column scaling factors computed by sgeequ. More...
 
subroutine dtgex2 (WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, J1, N1, N2, WORK, LWORK, INFO)
 DTGEX2 swaps adjacent diagonal blocks in an upper (quasi) triangular matrix pair by an orthogonal equivalence transformation. More...
 

Detailed Description

This is the group of double auxiliary functions for GE matrices

Function Documentation

subroutine dgesc2 ( integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  RHS,
integer, dimension( * )  IPIV,
integer, dimension( * )  JPIV,
double precision  SCALE 
)

DGESC2 solves a system of linear equations using the LU factorization with complete pivoting computed by sgetc2.

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

Purpose:
 DGESC2 solves a system of linear equations

           A * X = scale* RHS

 with a general N-by-N matrix A using the LU factorization with
 complete pivoting computed by DGETC2.
Parameters
[in]N
          N is INTEGER
          The order of the matrix A.
[in]A
          A is DOUBLE PRECISION array, dimension (LDA,N)
          On entry, the  LU part of the factorization of the n-by-n
          matrix A computed by DGETC2:  A = P * L * U * Q
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1, N).
[in,out]RHS
          RHS is DOUBLE PRECISION array, dimension (N).
          On entry, the right hand side vector b.
          On exit, the solution vector X.
[in]IPIV
          IPIV is INTEGER array, dimension (N).
          The pivot indices; for 1 <= i <= N, row i of the
          matrix has been interchanged with row IPIV(i).
[in]JPIV
          JPIV is INTEGER array, dimension (N).
          The pivot indices; for 1 <= j <= N, column j of the
          matrix has been interchanged with column JPIV(j).
[out]SCALE
          SCALE is DOUBLE PRECISION
          On exit, SCALE contains the scale factor. SCALE is chosen
          0 <= SCALE <= 1 to prevent owerflow in the solution.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012
Contributors:
Bo Kagstrom and Peter Poromaa, Department of Computing Science, Umea University, S-901 87 Umea, Sweden.

Definition at line 116 of file dgesc2.f.

116 *
117 * -- LAPACK auxiliary routine (version 3.4.2) --
118 * -- LAPACK is a software package provided by Univ. of Tennessee, --
119 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
120 * September 2012
121 *
122 * .. Scalar Arguments ..
123  INTEGER lda, n
124  DOUBLE PRECISION scale
125 * ..
126 * .. Array Arguments ..
127  INTEGER ipiv( * ), jpiv( * )
128  DOUBLE PRECISION a( lda, * ), rhs( * )
129 * ..
130 *
131 * =====================================================================
132 *
133 * .. Parameters ..
134  DOUBLE PRECISION one, two
135  parameter( one = 1.0d+0, two = 2.0d+0 )
136 * ..
137 * .. Local Scalars ..
138  INTEGER i, j
139  DOUBLE PRECISION bignum, eps, smlnum, temp
140 * ..
141 * .. External Subroutines ..
142  EXTERNAL dlaswp, dscal
143 * ..
144 * .. External Functions ..
145  INTEGER idamax
146  DOUBLE PRECISION dlamch
147  EXTERNAL idamax, dlamch
148 * ..
149 * .. Intrinsic Functions ..
150  INTRINSIC abs
151 * ..
152 * .. Executable Statements ..
153 *
154 * Set constant to control owerflow
155 *
156  eps = dlamch( 'P' )
157  smlnum = dlamch( 'S' ) / eps
158  bignum = one / smlnum
159  CALL dlabad( smlnum, bignum )
160 *
161 * Apply permutations IPIV to RHS
162 *
163  CALL dlaswp( 1, rhs, lda, 1, n-1, ipiv, 1 )
164 *
165 * Solve for L part
166 *
167  DO 20 i = 1, n - 1
168  DO 10 j = i + 1, n
169  rhs( j ) = rhs( j ) - a( j, i )*rhs( i )
170  10 CONTINUE
171  20 CONTINUE
172 *
173 * Solve for U part
174 *
175  scale = one
176 *
177 * Check for scaling
178 *
179  i = idamax( n, rhs, 1 )
180  IF( two*smlnum*abs( rhs( i ) ).GT.abs( a( n, n ) ) ) THEN
181  temp = ( one / two ) / abs( rhs( i ) )
182  CALL dscal( n, temp, rhs( 1 ), 1 )
183  scale = scale*temp
184  END IF
185 *
186  DO 40 i = n, 1, -1
187  temp = one / a( i, i )
188  rhs( i ) = rhs( i )*temp
189  DO 30 j = i + 1, n
190  rhs( i ) = rhs( i ) - rhs( j )*( a( i, j )*temp )
191  30 CONTINUE
192  40 CONTINUE
193 *
194 * Apply permutations JPIV to the solution (RHS)
195 *
196  CALL dlaswp( 1, rhs, lda, 1, n-1, jpiv, -1 )
197  RETURN
198 *
199 * End of DGESC2
200 *
subroutine dlaswp(N, A, LDA, K1, K2, IPIV, INCX)
DLASWP performs a series of row interchanges on a general rectangular matrix.
Definition: dlaswp.f:116
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
integer function idamax(N, DX, INCX)
IDAMAX
Definition: idamax.f:53
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dgetc2 ( integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
integer, dimension( * )  IPIV,
integer, dimension( * )  JPIV,
integer  INFO 
)

DGETC2 computes the LU factorization with complete pivoting of the general n-by-n matrix.

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

Purpose:
 DGETC2 computes an LU factorization with complete pivoting of the
 n-by-n matrix A. The factorization has the form A = P * L * U * Q,
 where P and Q are permutation matrices, L is lower triangular with
 unit diagonal elements and U is upper triangular.

 This is the Level 2 BLAS algorithm.
Parameters
[in]N
          N is INTEGER
          The order of the matrix A. N >= 0.
[in,out]A
          A is DOUBLE PRECISION array, dimension (LDA, N)
          On entry, the n-by-n matrix A to be factored.
          On exit, the factors L and U from the factorization
          A = P*L*U*Q; the unit diagonal elements of L are not stored.
          If U(k, k) appears to be less than SMIN, U(k, k) is given the
          value of SMIN, i.e., giving a nonsingular perturbed system.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[out]IPIV
          IPIV is INTEGER array, dimension(N).
          The pivot indices; for 1 <= i <= N, row i of the
          matrix has been interchanged with row IPIV(i).
[out]JPIV
          JPIV is INTEGER array, dimension(N).
          The pivot indices; for 1 <= j <= N, column j of the
          matrix has been interchanged with column JPIV(j).
[out]INFO
          INFO is INTEGER
           = 0: successful exit
           > 0: if INFO = k, U(k, k) is likely to produce owerflow if
                we try to solve for x in Ax = b. So U is perturbed to
                avoid the overflow.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2013
Contributors:
Bo Kagstrom and Peter Poromaa, Department of Computing Science, Umea University, S-901 87 Umea, Sweden.

Definition at line 113 of file dgetc2.f.

113 *
114 * -- LAPACK auxiliary routine (version 3.5.0) --
115 * -- LAPACK is a software package provided by Univ. of Tennessee, --
116 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
117 * November 2013
118 *
119 * .. Scalar Arguments ..
120  INTEGER info, lda, n
121 * ..
122 * .. Array Arguments ..
123  INTEGER ipiv( * ), jpiv( * )
124  DOUBLE PRECISION a( lda, * )
125 * ..
126 *
127 * =====================================================================
128 *
129 * .. Parameters ..
130  DOUBLE PRECISION zero, one
131  parameter( zero = 0.0d+0, one = 1.0d+0 )
132 * ..
133 * .. Local Scalars ..
134  INTEGER i, ip, ipv, j, jp, jpv
135  DOUBLE PRECISION bignum, eps, smin, smlnum, xmax
136 * ..
137 * .. External Subroutines ..
138  EXTERNAL dger, dswap
139 * ..
140 * .. External Functions ..
141  DOUBLE PRECISION dlamch
142  EXTERNAL dlamch
143 * ..
144 * .. Intrinsic Functions ..
145  INTRINSIC abs, max
146 * ..
147 * .. Executable Statements ..
148 *
149 * Set constants to control overflow
150 *
151  info = 0
152  eps = dlamch( 'P' )
153  smlnum = dlamch( 'S' ) / eps
154  bignum = one / smlnum
155  CALL dlabad( smlnum, bignum )
156 *
157 * Factorize A using complete pivoting.
158 * Set pivots less than SMIN to SMIN.
159 *
160  DO 40 i = 1, n - 1
161 *
162 * Find max element in matrix A
163 *
164  xmax = zero
165  DO 20 ip = i, n
166  DO 10 jp = i, n
167  IF( abs( a( ip, jp ) ).GE.xmax ) THEN
168  xmax = abs( a( ip, jp ) )
169  ipv = ip
170  jpv = jp
171  END IF
172  10 CONTINUE
173  20 CONTINUE
174  IF( i.EQ.1 )
175  $ smin = max( eps*xmax, smlnum )
176 *
177 * Swap rows
178 *
179  IF( ipv.NE.i )
180  $ CALL dswap( n, a( ipv, 1 ), lda, a( i, 1 ), lda )
181  ipiv( i ) = ipv
182 *
183 * Swap columns
184 *
185  IF( jpv.NE.i )
186  $ CALL dswap( n, a( 1, jpv ), 1, a( 1, i ), 1 )
187  jpiv( i ) = jpv
188 *
189 * Check for singularity
190 *
191  IF( abs( a( i, i ) ).LT.smin ) THEN
192  info = i
193  a( i, i ) = smin
194  END IF
195  DO 30 j = i + 1, n
196  a( j, i ) = a( j, i ) / a( i, i )
197  30 CONTINUE
198  CALL dger( n-i, n-i, -one, a( i+1, i ), 1, a( i, i+1 ), lda,
199  $ a( i+1, i+1 ), lda )
200  40 CONTINUE
201 *
202  IF( abs( a( n, n ) ).LT.smin ) THEN
203  info = n
204  a( n, n ) = smin
205  END IF
206 *
207 * Set last pivots to N
208 *
209  ipiv( n ) = n
210  jpiv( n ) = n
211 *
212  RETURN
213 *
214 * End of DGETC2
215 *
subroutine dger(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
DGER
Definition: dger.f:132
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:53
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65

Here is the call graph for this function:

Here is the caller graph for this function:

double precision function dlange ( character  NORM,
integer  M,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  WORK 
)

DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value of any element of a general rectangular matrix.

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

Purpose:
 DLANGE  returns the value of the one norm,  or the Frobenius norm, or
 the  infinity norm,  or the  element of  largest absolute value  of a
 real matrix A.
Returns
DLANGE
    DLANGE = ( max(abs(A(i,j))), NORM = 'M' or 'm'
             (
             ( norm1(A),         NORM = '1', 'O' or 'o'
             (
             ( normI(A),         NORM = 'I' or 'i'
             (
             ( normF(A),         NORM = 'F', 'f', 'E' or 'e'

 where  norm1  denotes the  one norm of a matrix (maximum column sum),
 normI  denotes the  infinity norm  of a matrix  (maximum row sum) and
 normF  denotes the  Frobenius norm of a matrix (square root of sum of
 squares).  Note that  max(abs(A(i,j)))  is not a consistent matrix norm.
Parameters
[in]NORM
          NORM is CHARACTER*1
          Specifies the value to be returned in DLANGE as described
          above.
[in]M
          M is INTEGER
          The number of rows of the matrix A.  M >= 0.  When M = 0,
          DLANGE is set to zero.
[in]N
          N is INTEGER
          The number of columns of the matrix A.  N >= 0.  When N = 0,
          DLANGE is set to zero.
[in]A
          A is DOUBLE PRECISION array, dimension (LDA,N)
          The m by n matrix A.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(M,1).
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)),
          where LWORK >= M when NORM = 'I'; otherwise, WORK is not
          referenced.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 116 of file dlange.f.

116 *
117 * -- LAPACK auxiliary routine (version 3.4.2) --
118 * -- LAPACK is a software package provided by Univ. of Tennessee, --
119 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
120 * September 2012
121 *
122 * .. Scalar Arguments ..
123  CHARACTER norm
124  INTEGER lda, m, n
125 * ..
126 * .. Array Arguments ..
127  DOUBLE PRECISION a( lda, * ), work( * )
128 * ..
129 *
130 * =====================================================================
131 *
132 * .. Parameters ..
133  DOUBLE PRECISION one, zero
134  parameter( one = 1.0d+0, zero = 0.0d+0 )
135 * ..
136 * .. Local Scalars ..
137  INTEGER i, j
138  DOUBLE PRECISION scale, sum, VALUE, temp
139 * ..
140 * .. External Subroutines ..
141  EXTERNAL dlassq
142 * ..
143 * .. External Functions ..
144  LOGICAL lsame, disnan
145  EXTERNAL lsame, disnan
146 * ..
147 * .. Intrinsic Functions ..
148  INTRINSIC abs, min, sqrt
149 * ..
150 * .. Executable Statements ..
151 *
152  IF( min( m, n ).EQ.0 ) THEN
153  VALUE = zero
154  ELSE IF( lsame( norm, 'M' ) ) THEN
155 *
156 * Find max(abs(A(i,j))).
157 *
158  VALUE = zero
159  DO 20 j = 1, n
160  DO 10 i = 1, m
161  temp = abs( a( i, j ) )
162  IF( VALUE.LT.temp .OR. disnan( temp ) ) VALUE = temp
163  10 CONTINUE
164  20 CONTINUE
165  ELSE IF( ( lsame( norm, 'O' ) ) .OR. ( norm.EQ.'1' ) ) THEN
166 *
167 * Find norm1(A).
168 *
169  VALUE = zero
170  DO 40 j = 1, n
171  sum = zero
172  DO 30 i = 1, m
173  sum = sum + abs( a( i, j ) )
174  30 CONTINUE
175  IF( VALUE.LT.sum .OR. disnan( sum ) ) VALUE = sum
176  40 CONTINUE
177  ELSE IF( lsame( norm, 'I' ) ) THEN
178 *
179 * Find normI(A).
180 *
181  DO 50 i = 1, m
182  work( i ) = zero
183  50 CONTINUE
184  DO 70 j = 1, n
185  DO 60 i = 1, m
186  work( i ) = work( i ) + abs( a( i, j ) )
187  60 CONTINUE
188  70 CONTINUE
189  VALUE = zero
190  DO 80 i = 1, m
191  temp = work( i )
192  IF( VALUE.LT.temp .OR. disnan( temp ) ) VALUE = temp
193  80 CONTINUE
194  ELSE IF( ( lsame( norm, 'F' ) ) .OR. ( lsame( norm, 'E' ) ) ) THEN
195 *
196 * Find normF(A).
197 *
198  scale = zero
199  sum = one
200  DO 90 j = 1, n
201  CALL dlassq( m, a( 1, j ), 1, scale, sum )
202  90 CONTINUE
203  VALUE = scale*sqrt( sum )
204  END IF
205 *
206  dlange = VALUE
207  RETURN
208 *
209 * End of DLANGE
210 *
logical function disnan(DIN)
DISNAN tests input for NaN.
Definition: disnan.f:61
subroutine dlassq(N, X, INCX, SCALE, SUMSQ)
DLASSQ updates a sum of squares represented in scaled form.
Definition: dlassq.f:105
double precision function dlange(NORM, M, N, A, LDA, WORK)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: dlange.f:116
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine dlaqge ( integer  M,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  R,
double precision, dimension( * )  C,
double precision  ROWCND,
double precision  COLCND,
double precision  AMAX,
character  EQUED 
)

DLAQGE scales a general rectangular matrix, using row and column scaling factors computed by sgeequ.

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

Purpose:
 DLAQGE equilibrates a general M by N matrix A using the row and
 column scaling factors in the vectors R and C.
Parameters
[in]M
          M is INTEGER
          The number of rows of the matrix A.  M >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrix A.  N >= 0.
[in,out]A
          A is DOUBLE PRECISION array, dimension (LDA,N)
          On entry, the M by N matrix A.
          On exit, the equilibrated matrix.  See EQUED for the form of
          the equilibrated matrix.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(M,1).
[in]R
          R is DOUBLE PRECISION array, dimension (M)
          The row scale factors for A.
[in]C
          C is DOUBLE PRECISION array, dimension (N)
          The column scale factors for A.
[in]ROWCND
          ROWCND is DOUBLE PRECISION
          Ratio of the smallest R(i) to the largest R(i).
[in]COLCND
          COLCND is DOUBLE PRECISION
          Ratio of the smallest C(i) to the largest C(i).
[in]AMAX
          AMAX is DOUBLE PRECISION
          Absolute value of largest matrix entry.
[out]EQUED
          EQUED is CHARACTER*1
          Specifies the form of equilibration that was done.
          = 'N':  No equilibration
          = 'R':  Row equilibration, i.e., A has been premultiplied by
                  diag(R).
          = 'C':  Column equilibration, i.e., A has been postmultiplied
                  by diag(C).
          = 'B':  Both row and column equilibration, i.e., A has been
                  replaced by diag(R) * A * diag(C).
Internal Parameters:
  THRESH is a threshold value used to decide if row or column scaling
  should be done based on the ratio of the row or column scaling
  factors.  If ROWCND < THRESH, row scaling is done, and if
  COLCND < THRESH, column scaling is done.

  LARGE and SMALL are threshold values used to decide if row scaling
  should be done based on the absolute size of the largest matrix
  element.  If AMAX > LARGE or AMAX < SMALL, row scaling is done.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 144 of file dlaqge.f.

144 *
145 * -- LAPACK auxiliary routine (version 3.4.2) --
146 * -- LAPACK is a software package provided by Univ. of Tennessee, --
147 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
148 * September 2012
149 *
150 * .. Scalar Arguments ..
151  CHARACTER equed
152  INTEGER lda, m, n
153  DOUBLE PRECISION amax, colcnd, rowcnd
154 * ..
155 * .. Array Arguments ..
156  DOUBLE PRECISION a( lda, * ), c( * ), r( * )
157 * ..
158 *
159 * =====================================================================
160 *
161 * .. Parameters ..
162  DOUBLE PRECISION one, thresh
163  parameter( one = 1.0d+0, thresh = 0.1d+0 )
164 * ..
165 * .. Local Scalars ..
166  INTEGER i, j
167  DOUBLE PRECISION cj, large, small
168 * ..
169 * .. External Functions ..
170  DOUBLE PRECISION dlamch
171  EXTERNAL dlamch
172 * ..
173 * .. Executable Statements ..
174 *
175 * Quick return if possible
176 *
177  IF( m.LE.0 .OR. n.LE.0 ) THEN
178  equed = 'N'
179  RETURN
180  END IF
181 *
182 * Initialize LARGE and SMALL.
183 *
184  small = dlamch( 'Safe minimum' ) / dlamch( 'Precision' )
185  large = one / small
186 *
187  IF( rowcnd.GE.thresh .AND. amax.GE.small .AND. amax.LE.large )
188  $ THEN
189 *
190 * No row scaling
191 *
192  IF( colcnd.GE.thresh ) THEN
193 *
194 * No column scaling
195 *
196  equed = 'N'
197  ELSE
198 *
199 * Column scaling
200 *
201  DO 20 j = 1, n
202  cj = c( j )
203  DO 10 i = 1, m
204  a( i, j ) = cj*a( i, j )
205  10 CONTINUE
206  20 CONTINUE
207  equed = 'C'
208  END IF
209  ELSE IF( colcnd.GE.thresh ) THEN
210 *
211 * Row scaling, no column scaling
212 *
213  DO 40 j = 1, n
214  DO 30 i = 1, m
215  a( i, j ) = r( i )*a( i, j )
216  30 CONTINUE
217  40 CONTINUE
218  equed = 'R'
219  ELSE
220 *
221 * Row and column scaling
222 *
223  DO 60 j = 1, n
224  cj = c( j )
225  DO 50 i = 1, m
226  a( i, j ) = cj*r( i )*a( i, j )
227  50 CONTINUE
228  60 CONTINUE
229  equed = 'B'
230  END IF
231 *
232  RETURN
233 *
234 * End of DLAQGE
235 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65

Here is the caller graph for this function:

subroutine dtgex2 ( logical  WANTQ,
logical  WANTZ,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( ldb, * )  B,
integer  LDB,
double precision, dimension( ldq, * )  Q,
integer  LDQ,
double precision, dimension( ldz, * )  Z,
integer  LDZ,
integer  J1,
integer  N1,
integer  N2,
double precision, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

DTGEX2 swaps adjacent diagonal blocks in an upper (quasi) triangular matrix pair by an orthogonal equivalence transformation.

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

Purpose:
 DTGEX2 swaps adjacent diagonal blocks (A11, B11) and (A22, B22)
 of size 1-by-1 or 2-by-2 in an upper (quasi) triangular matrix pair
 (A, B) by an orthogonal equivalence transformation.

 (A, B) must be in generalized real Schur canonical form (as returned
 by DGGES), i.e. A is block upper triangular with 1-by-1 and 2-by-2
 diagonal blocks. B is upper triangular.

 Optionally, the matrices Q and Z of generalized Schur vectors are
 updated.

        Q(in) * A(in) * Z(in)**T = Q(out) * A(out) * Z(out)**T
        Q(in) * B(in) * Z(in)**T = Q(out) * B(out) * Z(out)**T
Parameters
[in]WANTQ
          WANTQ is LOGICAL
          .TRUE. : update the left transformation matrix Q;
          .FALSE.: do not update Q.
[in]WANTZ
          WANTZ is LOGICAL
          .TRUE. : update the right transformation matrix Z;
          .FALSE.: do not update Z.
[in]N
          N is INTEGER
          The order of the matrices A and B. N >= 0.
[in,out]A
          A is DOUBLE PRECISION array, dimensions (LDA,N)
          On entry, the matrix A in the pair (A, B).
          On exit, the updated matrix A.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A. LDA >= max(1,N).
[in,out]B
          B is DOUBLE PRECISION array, dimensions (LDB,N)
          On entry, the matrix B in the pair (A, B).
          On exit, the updated matrix B.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B. LDB >= max(1,N).
[in,out]Q
          Q is DOUBLE PRECISION array, dimension (LDQ,N)
          On entry, if WANTQ = .TRUE., the orthogonal matrix Q.
          On exit, the updated matrix Q.
          Not referenced if WANTQ = .FALSE..
[in]LDQ
          LDQ is INTEGER
          The leading dimension of the array Q. LDQ >= 1.
          If WANTQ = .TRUE., LDQ >= N.
[in,out]Z
          Z is DOUBLE PRECISION array, dimension (LDZ,N)
          On entry, if WANTZ =.TRUE., the orthogonal matrix Z.
          On exit, the updated matrix Z.
          Not referenced if WANTZ = .FALSE..
[in]LDZ
          LDZ is INTEGER
          The leading dimension of the array Z. LDZ >= 1.
          If WANTZ = .TRUE., LDZ >= N.
[in]J1
          J1 is INTEGER
          The index to the first block (A11, B11). 1 <= J1 <= N.
[in]N1
          N1 is INTEGER
          The order of the first block (A11, B11). N1 = 0, 1 or 2.
[in]N2
          N2 is INTEGER
          The order of the second block (A22, B22). N2 = 0, 1 or 2.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)).
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.
          LWORK >=  MAX( 1, N*(N2+N1), (N2+N1)*(N2+N1)*2 )
[out]INFO
          INFO is INTEGER
            =0: Successful exit
            >0: If INFO = 1, the transformed matrix (A, B) would be
                too far from generalized Schur form; the blocks are
                not swapped and (A, B) and (Q, Z) are unchanged.
                The problem of swapping is too ill-conditioned.
            <0: If INFO = -16: LWORK is too small. Appropriate value
                for LWORK is returned in WORK(1).
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2015
Further Details:
In the current code both weak and strong stability tests are performed. The user can omit the strong stability test by changing the internal logical parameter WANDS to .FALSE.. See ref. [2] for details.
Contributors:
Bo Kagstrom and Peter Poromaa, Department of Computing Science, Umea University, S-901 87 Umea, Sweden.
References:
  [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
      Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
      M.S. Moonen et al (eds), Linear Algebra for Large Scale and
      Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.

  [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified
      Eigenvalues of a Regular Matrix Pair (A, B) and Condition
      Estimation: Theory, Algorithms and Software,
      Report UMINF - 94.04, Department of Computing Science, Umea
      University, S-901 87 Umea, Sweden, 1994. Also as LAPACK Working
      Note 87. To appear in Numerical Algorithms, 1996.

Definition at line 223 of file dtgex2.f.

223 *
224 * -- LAPACK auxiliary routine (version 3.6.0) --
225 * -- LAPACK is a software package provided by Univ. of Tennessee, --
226 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
227 * November 2015
228 *
229 * .. Scalar Arguments ..
230  LOGICAL wantq, wantz
231  INTEGER info, j1, lda, ldb, ldq, ldz, lwork, n, n1, n2
232 * ..
233 * .. Array Arguments ..
234  DOUBLE PRECISION a( lda, * ), b( ldb, * ), q( ldq, * ),
235  $ work( * ), z( ldz, * )
236 * ..
237 *
238 * =====================================================================
239 * Replaced various illegal calls to DCOPY by calls to DLASET, or by DO
240 * loops. Sven Hammarling, 1/5/02.
241 *
242 * .. Parameters ..
243  DOUBLE PRECISION zero, one
244  parameter( zero = 0.0d+0, one = 1.0d+0 )
245  DOUBLE PRECISION twenty
246  parameter( twenty = 2.0d+01 )
247  INTEGER ldst
248  parameter( ldst = 4 )
249  LOGICAL wands
250  parameter( wands = .true. )
251 * ..
252 * .. Local Scalars ..
253  LOGICAL dtrong, weak
254  INTEGER i, idum, linfo, m
255  DOUBLE PRECISION bqra21, brqa21, ddum, dnorm, dscale, dsum, eps,
256  $ f, g, sa, sb, scale, smlnum, ss, thresh, ws
257 * ..
258 * .. Local Arrays ..
259  INTEGER iwork( ldst )
260  DOUBLE PRECISION ai( 2 ), ar( 2 ), be( 2 ), ir( ldst, ldst ),
261  $ ircop( ldst, ldst ), li( ldst, ldst ),
262  $ licop( ldst, ldst ), s( ldst, ldst ),
263  $ scpy( ldst, ldst ), t( ldst, ldst ),
264  $ taul( ldst ), taur( ldst ), tcpy( ldst, ldst )
265 * ..
266 * .. External Functions ..
267  DOUBLE PRECISION dlamch
268  EXTERNAL dlamch
269 * ..
270 * .. External Subroutines ..
271  EXTERNAL dgemm, dgeqr2, dgerq2, dlacpy, dlagv2, dlartg,
273  $ drot, dscal, dtgsy2
274 * ..
275 * .. Intrinsic Functions ..
276  INTRINSIC abs, max, sqrt
277 * ..
278 * .. Executable Statements ..
279 *
280  info = 0
281 *
282 * Quick return if possible
283 *
284  IF( n.LE.1 .OR. n1.LE.0 .OR. n2.LE.0 )
285  $ RETURN
286  IF( n1.GT.n .OR. ( j1+n1 ).GT.n )
287  $ RETURN
288  m = n1 + n2
289  IF( lwork.LT.max( 1, n*m, m*m*2 ) ) THEN
290  info = -16
291  work( 1 ) = max( 1, n*m, m*m*2 )
292  RETURN
293  END IF
294 *
295  weak = .false.
296  dtrong = .false.
297 *
298 * Make a local copy of selected block
299 *
300  CALL dlaset( 'Full', ldst, ldst, zero, zero, li, ldst )
301  CALL dlaset( 'Full', ldst, ldst, zero, zero, ir, ldst )
302  CALL dlacpy( 'Full', m, m, a( j1, j1 ), lda, s, ldst )
303  CALL dlacpy( 'Full', m, m, b( j1, j1 ), ldb, t, ldst )
304 *
305 * Compute threshold for testing acceptance of swapping.
306 *
307  eps = dlamch( 'P' )
308  smlnum = dlamch( 'S' ) / eps
309  dscale = zero
310  dsum = one
311  CALL dlacpy( 'Full', m, m, s, ldst, work, m )
312  CALL dlassq( m*m, work, 1, dscale, dsum )
313  CALL dlacpy( 'Full', m, m, t, ldst, work, m )
314  CALL dlassq( m*m, work, 1, dscale, dsum )
315  dnorm = dscale*sqrt( dsum )
316 *
317 * THRES has been changed from
318 * THRESH = MAX( TEN*EPS*SA, SMLNUM )
319 * to
320 * THRESH = MAX( TWENTY*EPS*SA, SMLNUM )
321 * on 04/01/10.
322 * "Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by
323 * Jim Demmel and Guillaume Revy. See forum post 1783.
324 *
325  thresh = max( twenty*eps*dnorm, smlnum )
326 *
327  IF( m.EQ.2 ) THEN
328 *
329 * CASE 1: Swap 1-by-1 and 1-by-1 blocks.
330 *
331 * Compute orthogonal QL and RQ that swap 1-by-1 and 1-by-1 blocks
332 * using Givens rotations and perform the swap tentatively.
333 *
334  f = s( 2, 2 )*t( 1, 1 ) - t( 2, 2 )*s( 1, 1 )
335  g = s( 2, 2 )*t( 1, 2 ) - t( 2, 2 )*s( 1, 2 )
336  sb = abs( t( 2, 2 ) )
337  sa = abs( s( 2, 2 ) )
338  CALL dlartg( f, g, ir( 1, 2 ), ir( 1, 1 ), ddum )
339  ir( 2, 1 ) = -ir( 1, 2 )
340  ir( 2, 2 ) = ir( 1, 1 )
341  CALL drot( 2, s( 1, 1 ), 1, s( 1, 2 ), 1, ir( 1, 1 ),
342  $ ir( 2, 1 ) )
343  CALL drot( 2, t( 1, 1 ), 1, t( 1, 2 ), 1, ir( 1, 1 ),
344  $ ir( 2, 1 ) )
345  IF( sa.GE.sb ) THEN
346  CALL dlartg( s( 1, 1 ), s( 2, 1 ), li( 1, 1 ), li( 2, 1 ),
347  $ ddum )
348  ELSE
349  CALL dlartg( t( 1, 1 ), t( 2, 1 ), li( 1, 1 ), li( 2, 1 ),
350  $ ddum )
351  END IF
352  CALL drot( 2, s( 1, 1 ), ldst, s( 2, 1 ), ldst, li( 1, 1 ),
353  $ li( 2, 1 ) )
354  CALL drot( 2, t( 1, 1 ), ldst, t( 2, 1 ), ldst, li( 1, 1 ),
355  $ li( 2, 1 ) )
356  li( 2, 2 ) = li( 1, 1 )
357  li( 1, 2 ) = -li( 2, 1 )
358 *
359 * Weak stability test:
360 * |S21| + |T21| <= O(EPS * F-norm((S, T)))
361 *
362  ws = abs( s( 2, 1 ) ) + abs( t( 2, 1 ) )
363  weak = ws.LE.thresh
364  IF( .NOT.weak )
365  $ GO TO 70
366 *
367  IF( wands ) THEN
368 *
369 * Strong stability test:
370 * F-norm((A-QL**T*S*QR, B-QL**T*T*QR)) <= O(EPS*F-norm((A,B)))
371 *
372  CALL dlacpy( 'Full', m, m, a( j1, j1 ), lda, work( m*m+1 ),
373  $ m )
374  CALL dgemm( 'N', 'N', m, m, m, one, li, ldst, s, ldst, zero,
375  $ work, m )
376  CALL dgemm( 'N', 'T', m, m, m, -one, work, m, ir, ldst, one,
377  $ work( m*m+1 ), m )
378  dscale = zero
379  dsum = one
380  CALL dlassq( m*m, work( m*m+1 ), 1, dscale, dsum )
381 *
382  CALL dlacpy( 'Full', m, m, b( j1, j1 ), ldb, work( m*m+1 ),
383  $ m )
384  CALL dgemm( 'N', 'N', m, m, m, one, li, ldst, t, ldst, zero,
385  $ work, m )
386  CALL dgemm( 'N', 'T', m, m, m, -one, work, m, ir, ldst, one,
387  $ work( m*m+1 ), m )
388  CALL dlassq( m*m, work( m*m+1 ), 1, dscale, dsum )
389  ss = dscale*sqrt( dsum )
390  dtrong = ss.LE.thresh
391  IF( .NOT.dtrong )
392  $ GO TO 70
393  END IF
394 *
395 * Update (A(J1:J1+M-1, M+J1:N), B(J1:J1+M-1, M+J1:N)) and
396 * (A(1:J1-1, J1:J1+M), B(1:J1-1, J1:J1+M)).
397 *
398  CALL drot( j1+1, a( 1, j1 ), 1, a( 1, j1+1 ), 1, ir( 1, 1 ),
399  $ ir( 2, 1 ) )
400  CALL drot( j1+1, b( 1, j1 ), 1, b( 1, j1+1 ), 1, ir( 1, 1 ),
401  $ ir( 2, 1 ) )
402  CALL drot( n-j1+1, a( j1, j1 ), lda, a( j1+1, j1 ), lda,
403  $ li( 1, 1 ), li( 2, 1 ) )
404  CALL drot( n-j1+1, b( j1, j1 ), ldb, b( j1+1, j1 ), ldb,
405  $ li( 1, 1 ), li( 2, 1 ) )
406 *
407 * Set N1-by-N2 (2,1) - blocks to ZERO.
408 *
409  a( j1+1, j1 ) = zero
410  b( j1+1, j1 ) = zero
411 *
412 * Accumulate transformations into Q and Z if requested.
413 *
414  IF( wantz )
415  $ CALL drot( n, z( 1, j1 ), 1, z( 1, j1+1 ), 1, ir( 1, 1 ),
416  $ ir( 2, 1 ) )
417  IF( wantq )
418  $ CALL drot( n, q( 1, j1 ), 1, q( 1, j1+1 ), 1, li( 1, 1 ),
419  $ li( 2, 1 ) )
420 *
421 * Exit with INFO = 0 if swap was successfully performed.
422 *
423  RETURN
424 *
425  ELSE
426 *
427 * CASE 2: Swap 1-by-1 and 2-by-2 blocks, or 2-by-2
428 * and 2-by-2 blocks.
429 *
430 * Solve the generalized Sylvester equation
431 * S11 * R - L * S22 = SCALE * S12
432 * T11 * R - L * T22 = SCALE * T12
433 * for R and L. Solutions in LI and IR.
434 *
435  CALL dlacpy( 'Full', n1, n2, t( 1, n1+1 ), ldst, li, ldst )
436  CALL dlacpy( 'Full', n1, n2, s( 1, n1+1 ), ldst,
437  $ ir( n2+1, n1+1 ), ldst )
438  CALL dtgsy2( 'N', 0, n1, n2, s, ldst, s( n1+1, n1+1 ), ldst,
439  $ ir( n2+1, n1+1 ), ldst, t, ldst, t( n1+1, n1+1 ),
440  $ ldst, li, ldst, scale, dsum, dscale, iwork, idum,
441  $ linfo )
442 *
443 * Compute orthogonal matrix QL:
444 *
445 * QL**T * LI = [ TL ]
446 * [ 0 ]
447 * where
448 * LI = [ -L ]
449 * [ SCALE * identity(N2) ]
450 *
451  DO 10 i = 1, n2
452  CALL dscal( n1, -one, li( 1, i ), 1 )
453  li( n1+i, i ) = scale
454  10 CONTINUE
455  CALL dgeqr2( m, n2, li, ldst, taul, work, linfo )
456  IF( linfo.NE.0 )
457  $ GO TO 70
458  CALL dorg2r( m, m, n2, li, ldst, taul, work, linfo )
459  IF( linfo.NE.0 )
460  $ GO TO 70
461 *
462 * Compute orthogonal matrix RQ:
463 *
464 * IR * RQ**T = [ 0 TR],
465 *
466 * where IR = [ SCALE * identity(N1), R ]
467 *
468  DO 20 i = 1, n1
469  ir( n2+i, i ) = scale
470  20 CONTINUE
471  CALL dgerq2( n1, m, ir( n2+1, 1 ), ldst, taur, work, linfo )
472  IF( linfo.NE.0 )
473  $ GO TO 70
474  CALL dorgr2( m, m, n1, ir, ldst, taur, work, linfo )
475  IF( linfo.NE.0 )
476  $ GO TO 70
477 *
478 * Perform the swapping tentatively:
479 *
480  CALL dgemm( 'T', 'N', m, m, m, one, li, ldst, s, ldst, zero,
481  $ work, m )
482  CALL dgemm( 'N', 'T', m, m, m, one, work, m, ir, ldst, zero, s,
483  $ ldst )
484  CALL dgemm( 'T', 'N', m, m, m, one, li, ldst, t, ldst, zero,
485  $ work, m )
486  CALL dgemm( 'N', 'T', m, m, m, one, work, m, ir, ldst, zero, t,
487  $ ldst )
488  CALL dlacpy( 'F', m, m, s, ldst, scpy, ldst )
489  CALL dlacpy( 'F', m, m, t, ldst, tcpy, ldst )
490  CALL dlacpy( 'F', m, m, ir, ldst, ircop, ldst )
491  CALL dlacpy( 'F', m, m, li, ldst, licop, ldst )
492 *
493 * Triangularize the B-part by an RQ factorization.
494 * Apply transformation (from left) to A-part, giving S.
495 *
496  CALL dgerq2( m, m, t, ldst, taur, work, linfo )
497  IF( linfo.NE.0 )
498  $ GO TO 70
499  CALL dormr2( 'R', 'T', m, m, m, t, ldst, taur, s, ldst, work,
500  $ linfo )
501  IF( linfo.NE.0 )
502  $ GO TO 70
503  CALL dormr2( 'L', 'N', m, m, m, t, ldst, taur, ir, ldst, work,
504  $ linfo )
505  IF( linfo.NE.0 )
506  $ GO TO 70
507 *
508 * Compute F-norm(S21) in BRQA21. (T21 is 0.)
509 *
510  dscale = zero
511  dsum = one
512  DO 30 i = 1, n2
513  CALL dlassq( n1, s( n2+1, i ), 1, dscale, dsum )
514  30 CONTINUE
515  brqa21 = dscale*sqrt( dsum )
516 *
517 * Triangularize the B-part by a QR factorization.
518 * Apply transformation (from right) to A-part, giving S.
519 *
520  CALL dgeqr2( m, m, tcpy, ldst, taul, work, linfo )
521  IF( linfo.NE.0 )
522  $ GO TO 70
523  CALL dorm2r( 'L', 'T', m, m, m, tcpy, ldst, taul, scpy, ldst,
524  $ work, info )
525  CALL dorm2r( 'R', 'N', m, m, m, tcpy, ldst, taul, licop, ldst,
526  $ work, info )
527  IF( linfo.NE.0 )
528  $ GO TO 70
529 *
530 * Compute F-norm(S21) in BQRA21. (T21 is 0.)
531 *
532  dscale = zero
533  dsum = one
534  DO 40 i = 1, n2
535  CALL dlassq( n1, scpy( n2+1, i ), 1, dscale, dsum )
536  40 CONTINUE
537  bqra21 = dscale*sqrt( dsum )
538 *
539 * Decide which method to use.
540 * Weak stability test:
541 * F-norm(S21) <= O(EPS * F-norm((S, T)))
542 *
543  IF( bqra21.LE.brqa21 .AND. bqra21.LE.thresh ) THEN
544  CALL dlacpy( 'F', m, m, scpy, ldst, s, ldst )
545  CALL dlacpy( 'F', m, m, tcpy, ldst, t, ldst )
546  CALL dlacpy( 'F', m, m, ircop, ldst, ir, ldst )
547  CALL dlacpy( 'F', m, m, licop, ldst, li, ldst )
548  ELSE IF( brqa21.GE.thresh ) THEN
549  GO TO 70
550  END IF
551 *
552 * Set lower triangle of B-part to zero
553 *
554  CALL dlaset( 'Lower', m-1, m-1, zero, zero, t(2,1), ldst )
555 *
556  IF( wands ) THEN
557 *
558 * Strong stability test:
559 * F-norm((A-QL*S*QR**T, B-QL*T*QR**T)) <= O(EPS*F-norm((A,B)))
560 *
561  CALL dlacpy( 'Full', m, m, a( j1, j1 ), lda, work( m*m+1 ),
562  $ m )
563  CALL dgemm( 'N', 'N', m, m, m, one, li, ldst, s, ldst, zero,
564  $ work, m )
565  CALL dgemm( 'N', 'N', m, m, m, -one, work, m, ir, ldst, one,
566  $ work( m*m+1 ), m )
567  dscale = zero
568  dsum = one
569  CALL dlassq( m*m, work( m*m+1 ), 1, dscale, dsum )
570 *
571  CALL dlacpy( 'Full', m, m, b( j1, j1 ), ldb, work( m*m+1 ),
572  $ m )
573  CALL dgemm( 'N', 'N', m, m, m, one, li, ldst, t, ldst, zero,
574  $ work, m )
575  CALL dgemm( 'N', 'N', m, m, m, -one, work, m, ir, ldst, one,
576  $ work( m*m+1 ), m )
577  CALL dlassq( m*m, work( m*m+1 ), 1, dscale, dsum )
578  ss = dscale*sqrt( dsum )
579  dtrong = ( ss.LE.thresh )
580  IF( .NOT.dtrong )
581  $ GO TO 70
582 *
583  END IF
584 *
585 * If the swap is accepted ("weakly" and "strongly"), apply the
586 * transformations and set N1-by-N2 (2,1)-block to zero.
587 *
588  CALL dlaset( 'Full', n1, n2, zero, zero, s(n2+1,1), ldst )
589 *
590 * copy back M-by-M diagonal block starting at index J1 of (A, B)
591 *
592  CALL dlacpy( 'F', m, m, s, ldst, a( j1, j1 ), lda )
593  CALL dlacpy( 'F', m, m, t, ldst, b( j1, j1 ), ldb )
594  CALL dlaset( 'Full', ldst, ldst, zero, zero, t, ldst )
595 *
596 * Standardize existing 2-by-2 blocks.
597 *
598  CALL dlaset( 'Full', m, m, zero, zero, work, m )
599  work( 1 ) = one
600  t( 1, 1 ) = one
601  idum = lwork - m*m - 2
602  IF( n2.GT.1 ) THEN
603  CALL dlagv2( a( j1, j1 ), lda, b( j1, j1 ), ldb, ar, ai, be,
604  $ work( 1 ), work( 2 ), t( 1, 1 ), t( 2, 1 ) )
605  work( m+1 ) = -work( 2 )
606  work( m+2 ) = work( 1 )
607  t( n2, n2 ) = t( 1, 1 )
608  t( 1, 2 ) = -t( 2, 1 )
609  END IF
610  work( m*m ) = one
611  t( m, m ) = one
612 *
613  IF( n1.GT.1 ) THEN
614  CALL dlagv2( a( j1+n2, j1+n2 ), lda, b( j1+n2, j1+n2 ), ldb,
615  $ taur, taul, work( m*m+1 ), work( n2*m+n2+1 ),
616  $ work( n2*m+n2+2 ), t( n2+1, n2+1 ),
617  $ t( m, m-1 ) )
618  work( m*m ) = work( n2*m+n2+1 )
619  work( m*m-1 ) = -work( n2*m+n2+2 )
620  t( m, m ) = t( n2+1, n2+1 )
621  t( m-1, m ) = -t( m, m-1 )
622  END IF
623  CALL dgemm( 'T', 'N', n2, n1, n2, one, work, m, a( j1, j1+n2 ),
624  $ lda, zero, work( m*m+1 ), n2 )
625  CALL dlacpy( 'Full', n2, n1, work( m*m+1 ), n2, a( j1, j1+n2 ),
626  $ lda )
627  CALL dgemm( 'T', 'N', n2, n1, n2, one, work, m, b( j1, j1+n2 ),
628  $ ldb, zero, work( m*m+1 ), n2 )
629  CALL dlacpy( 'Full', n2, n1, work( m*m+1 ), n2, b( j1, j1+n2 ),
630  $ ldb )
631  CALL dgemm( 'N', 'N', m, m, m, one, li, ldst, work, m, zero,
632  $ work( m*m+1 ), m )
633  CALL dlacpy( 'Full', m, m, work( m*m+1 ), m, li, ldst )
634  CALL dgemm( 'N', 'N', n2, n1, n1, one, a( j1, j1+n2 ), lda,
635  $ t( n2+1, n2+1 ), ldst, zero, work, n2 )
636  CALL dlacpy( 'Full', n2, n1, work, n2, a( j1, j1+n2 ), lda )
637  CALL dgemm( 'N', 'N', n2, n1, n1, one, b( j1, j1+n2 ), ldb,
638  $ t( n2+1, n2+1 ), ldst, zero, work, n2 )
639  CALL dlacpy( 'Full', n2, n1, work, n2, b( j1, j1+n2 ), ldb )
640  CALL dgemm( 'T', 'N', m, m, m, one, ir, ldst, t, ldst, zero,
641  $ work, m )
642  CALL dlacpy( 'Full', m, m, work, m, ir, ldst )
643 *
644 * Accumulate transformations into Q and Z if requested.
645 *
646  IF( wantq ) THEN
647  CALL dgemm( 'N', 'N', n, m, m, one, q( 1, j1 ), ldq, li,
648  $ ldst, zero, work, n )
649  CALL dlacpy( 'Full', n, m, work, n, q( 1, j1 ), ldq )
650 *
651  END IF
652 *
653  IF( wantz ) THEN
654  CALL dgemm( 'N', 'N', n, m, m, one, z( 1, j1 ), ldz, ir,
655  $ ldst, zero, work, n )
656  CALL dlacpy( 'Full', n, m, work, n, z( 1, j1 ), ldz )
657 *
658  END IF
659 *
660 * Update (A(J1:J1+M-1, M+J1:N), B(J1:J1+M-1, M+J1:N)) and
661 * (A(1:J1-1, J1:J1+M), B(1:J1-1, J1:J1+M)).
662 *
663  i = j1 + m
664  IF( i.LE.n ) THEN
665  CALL dgemm( 'T', 'N', m, n-i+1, m, one, li, ldst,
666  $ a( j1, i ), lda, zero, work, m )
667  CALL dlacpy( 'Full', m, n-i+1, work, m, a( j1, i ), lda )
668  CALL dgemm( 'T', 'N', m, n-i+1, m, one, li, ldst,
669  $ b( j1, i ), ldb, zero, work, m )
670  CALL dlacpy( 'Full', m, n-i+1, work, m, b( j1, i ), ldb )
671  END IF
672  i = j1 - 1
673  IF( i.GT.0 ) THEN
674  CALL dgemm( 'N', 'N', i, m, m, one, a( 1, j1 ), lda, ir,
675  $ ldst, zero, work, i )
676  CALL dlacpy( 'Full', i, m, work, i, a( 1, j1 ), lda )
677  CALL dgemm( 'N', 'N', i, m, m, one, b( 1, j1 ), ldb, ir,
678  $ ldst, zero, work, i )
679  CALL dlacpy( 'Full', i, m, work, i, b( 1, j1 ), ldb )
680  END IF
681 *
682 * Exit with INFO = 0 if swap was successfully performed.
683 *
684  RETURN
685 *
686  END IF
687 *
688 * Exit with INFO = 1 if swap was rejected.
689 *
690  70 CONTINUE
691 *
692  info = 1
693  RETURN
694 *
695 * End of DTGEX2
696 *
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 dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:105
subroutine dlartg(F, G, CS, SN, R)
DLARTG generates a plane rotation with real cosine and real sine.
Definition: dlartg.f:99
subroutine dlassq(N, X, INCX, SCALE, SUMSQ)
DLASSQ updates a sum of squares represented in scaled form.
Definition: dlassq.f:105
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:189
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
Definition: drot.f:53
subroutine dorg2r(M, N, K, A, LDA, TAU, WORK, INFO)
DORG2R generates all or part of the orthogonal matrix Q from a QR factorization determined by sgeqrf ...
Definition: dorg2r.f:116
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine dtgsy2(TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL, IWORK, PQ, INFO)
DTGSY2 solves the generalized Sylvester equation (unblocked algorithm).
Definition: dtgsy2.f:276
subroutine dgeqr2(M, N, A, LDA, TAU, WORK, INFO)
DGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm...
Definition: dgeqr2.f:123
subroutine dlagv2(A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, CSL, SNL, CSR, SNR)
DLAGV2 computes the Generalized Schur factorization of a real 2-by-2 matrix pencil (A...
Definition: dlagv2.f:159
subroutine dorm2r(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
DORM2R multiplies a general matrix by the orthogonal matrix from a QR factorization determined by sge...
Definition: dorm2r.f:161
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55
subroutine dgerq2(M, N, A, LDA, TAU, WORK, INFO)
DGERQ2 computes the RQ factorization of a general rectangular matrix using an unblocked algorithm...
Definition: dgerq2.f:125
subroutine dormr2(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
DORMR2 multiplies a general matrix by the orthogonal matrix from a RQ factorization determined by sge...
Definition: dormr2.f:161
subroutine dorgr2(M, N, K, A, LDA, TAU, WORK, INFO)
DORGR2 generates all or part of the orthogonal matrix Q from an RQ factorization determined by sgerqf...
Definition: dorgr2.f:116

Here is the call graph for this function:

Here is the caller graph for this function: