OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
double
LAPACK » General Matrices » Auxiliary routines

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.
subroutine dgetc2 (n, a, lda, ipiv, jpiv, info)
 DGETC2 computes the LU factorization with complete pivoting of the general n-by-n matrix.
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.
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.
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.

Detailed Description

This is the group of double auxiliary functions for GE matrices

Function Documentation

◆ dgesc2()

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 overflow in the solution.
!>
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 113 of file dgesc2.f.

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

◆ dgetc2()

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

Definition at line 110 of file dgetc2.f.

111*
112* -- LAPACK auxiliary routine --
113* -- LAPACK is a software package provided by Univ. of Tennessee, --
114* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
115*
116* .. Scalar Arguments ..
117 INTEGER INFO, LDA, N
118* ..
119* .. Array Arguments ..
120 INTEGER IPIV( * ), JPIV( * )
121 DOUBLE PRECISION A( LDA, * )
122* ..
123*
124* =====================================================================
125*
126* .. Parameters ..
127 DOUBLE PRECISION ZERO, ONE
128 parameter( zero = 0.0d+0, one = 1.0d+0 )
129* ..
130* .. Local Scalars ..
131 INTEGER I, IP, IPV, J, JP, JPV
132 DOUBLE PRECISION BIGNUM, EPS, SMIN, SMLNUM, XMAX
133* ..
134* .. External Subroutines ..
135 EXTERNAL dger, dswap, dlabad
136* ..
137* .. External Functions ..
138 DOUBLE PRECISION DLAMCH
139 EXTERNAL dlamch
140* ..
141* .. Intrinsic Functions ..
142 INTRINSIC abs, max
143* ..
144* .. Executable Statements ..
145*
146 info = 0
147*
148* Quick return if possible
149*
150 IF( n.EQ.0 )
151 $ RETURN
152*
153* Set constants to control overflow
154*
155 eps = dlamch( 'P' )
156 smlnum = dlamch( 'S' ) / eps
157 bignum = one / smlnum
158 CALL dlabad( smlnum, bignum )
159*
160* Handle the case N=1 by itself
161*
162 IF( n.EQ.1 ) THEN
163 ipiv( 1 ) = 1
164 jpiv( 1 ) = 1
165 IF( abs( a( 1, 1 ) ).LT.smlnum ) THEN
166 info = 1
167 a( 1, 1 ) = smlnum
168 END IF
169 RETURN
170 END IF
171*
172* Factorize A using complete pivoting.
173* Set pivots less than SMIN to SMIN.
174*
175 DO 40 i = 1, n - 1
176*
177* Find max element in matrix A
178*
179 xmax = zero
180 DO 20 ip = i, n
181 DO 10 jp = i, n
182 IF( abs( a( ip, jp ) ).GE.xmax ) THEN
183 xmax = abs( a( ip, jp ) )
184 ipv = ip
185 jpv = jp
186 END IF
187 10 CONTINUE
188 20 CONTINUE
189 IF( i.EQ.1 )
190 $ smin = max( eps*xmax, smlnum )
191*
192* Swap rows
193*
194 IF( ipv.NE.i )
195 $ CALL dswap( n, a( ipv, 1 ), lda, a( i, 1 ), lda )
196 ipiv( i ) = ipv
197*
198* Swap columns
199*
200 IF( jpv.NE.i )
201 $ CALL dswap( n, a( 1, jpv ), 1, a( 1, i ), 1 )
202 jpiv( i ) = jpv
203*
204* Check for singularity
205*
206 IF( abs( a( i, i ) ).LT.smin ) THEN
207 info = i
208 a( i, i ) = smin
209 END IF
210 DO 30 j = i + 1, n
211 a( j, i ) = a( j, i ) / a( i, i )
212 30 CONTINUE
213 CALL dger( n-i, n-i, -one, a( i+1, i ), 1, a( i, i+1 ), lda,
214 $ a( i+1, i+1 ), lda )
215 40 CONTINUE
216*
217 IF( abs( a( n, n ) ).LT.smin ) THEN
218 info = n
219 a( n, n ) = smin
220 END IF
221*
222* Set last pivots to N
223*
224 ipiv( n ) = n
225 jpiv( n ) = n
226*
227 RETURN
228*
229* End of DGETC2
230*
dswap
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82
dger
subroutine dger(m, n, alpha, x, incx, y, incy, a, lda)
DGER
Definition dger.f:130
max
#define max(a, b)
Definition macros.h:21

◆ dlange()

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.

Definition at line 113 of file dlange.f.

114*
115* -- LAPACK auxiliary routine --
116* -- LAPACK is a software package provided by Univ. of Tennessee, --
117* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
118*
119* .. Scalar Arguments ..
120 CHARACTER NORM
121 INTEGER LDA, M, N
122* ..
123* .. Array Arguments ..
124 DOUBLE PRECISION A( LDA, * ), WORK( * )
125* ..
126*
127* =====================================================================
128*
129* .. Parameters ..
130 DOUBLE PRECISION ONE, ZERO
131 parameter( one = 1.0d+0, zero = 0.0d+0 )
132* ..
133* .. Local Scalars ..
134 INTEGER I, J
135 DOUBLE PRECISION SCALE, SUM, VALUE, TEMP
136* ..
137* .. External Subroutines ..
138 EXTERNAL dlassq
139* ..
140* .. External Functions ..
141 LOGICAL LSAME, DISNAN
142 EXTERNAL lsame, disnan
143* ..
144* .. Intrinsic Functions ..
145 INTRINSIC abs, min, sqrt
146* ..
147* .. Executable Statements ..
148*
149 IF( min( m, n ).EQ.0 ) THEN
150 VALUE = zero
151 ELSE IF( lsame( norm, 'M' ) ) THEN
152*
153* Find max(abs(A(i,j))).
154*
155 VALUE = zero
156 DO 20 j = 1, n
157 DO 10 i = 1, m
158 temp = abs( a( i, j ) )
159 IF( VALUE.LT.temp .OR. disnan( temp ) ) VALUE = temp
160 10 CONTINUE
161 20 CONTINUE
162 ELSE IF( ( lsame( norm, 'O' ) ) .OR. ( norm.EQ.'1' ) ) THEN
163*
164* Find norm1(A).
165*
166 VALUE = zero
167 DO 40 j = 1, n
168 sum = zero
169 DO 30 i = 1, m
170 sum = sum + abs( a( i, j ) )
171 30 CONTINUE
172 IF( VALUE.LT.sum .OR. disnan( sum ) ) VALUE = sum
173 40 CONTINUE
174 ELSE IF( lsame( norm, 'I' ) ) THEN
175*
176* Find normI(A).
177*
178 DO 50 i = 1, m
179 work( i ) = zero
180 50 CONTINUE
181 DO 70 j = 1, n
182 DO 60 i = 1, m
183 work( i ) = work( i ) + abs( a( i, j ) )
184 60 CONTINUE
185 70 CONTINUE
186 VALUE = zero
187 DO 80 i = 1, m
188 temp = work( i )
189 IF( VALUE.LT.temp .OR. disnan( temp ) ) VALUE = temp
190 80 CONTINUE
191 ELSE IF( ( lsame( norm, 'F' ) ) .OR. ( lsame( norm, 'E' ) ) ) THEN
192*
193* Find normF(A).
194*
195 scale = zero
196 sum = one
197 DO 90 j = 1, n
198 CALL dlassq( m, a( 1, j ), 1, scale, sum )
199 90 CONTINUE
200 VALUE = scale*sqrt( sum )
201 END IF
202*
203 dlange = VALUE
204 RETURN
205*
206* End of DLANGE
207*
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
dlassq
subroutine dlassq(n, x, incx, scl, sumsq)
DLASSQ updates a sum of squares represented in scaled form.
Definition dlassq.f90:137
disnan
logical function disnan(din)
DISNAN tests input for NaN.
Definition disnan.f:59
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
dlange
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:114
min
#define min(a, b)
Definition macros.h:20

◆ dlaqge()

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.

Definition at line 140 of file dlaqge.f.

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

◆ dtgex2()

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.
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 219 of file dtgex2.f.

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