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

Functions

subroutine zgesc2 (n, a, lda, rhs, ipiv, jpiv, scale)
 ZGESC2 solves a system of linear equations using the LU factorization with complete pivoting computed by sgetc2.
subroutine zgetc2 (n, a, lda, ipiv, jpiv, info)
 ZGETC2 computes the LU factorization with complete pivoting of the general n-by-n matrix.
double precision function zlange (norm, m, n, a, lda, work)
 ZLANGE 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 zlaqge (m, n, a, lda, r, c, rowcnd, colcnd, amax, equed)
 ZLAQGE scales a general rectangular matrix, using row and column scaling factors computed by sgeequ.
subroutine ztgex2 (wantq, wantz, n, a, lda, b, ldb, q, ldq, z, ldz, j1, info)
 ZTGEX2 swaps adjacent diagonal blocks in an upper (quasi) triangular matrix pair by an unitary equivalence transformation.

Detailed Description

This is the group of complex16 auxiliary functions for GE matrices

Function Documentation

◆ zgesc2()

subroutine zgesc2 ( integer n,
complex*16, dimension( lda, * ) a,
integer lda,
complex*16, dimension( * ) rhs,
integer, dimension( * ) ipiv,
integer, dimension( * ) jpiv,
double precision scale )

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

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

Purpose:
!>
!> ZGESC2 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 ZGETC2.
!>
!>
Parameters
[in]N
!>          N is INTEGER
!>          The number of columns of the matrix A.
!>
[in]A
!>          A is COMPLEX*16 array, dimension (LDA, N)
!>          On entry, the  LU part of the factorization of the n-by-n
!>          matrix A computed by ZGETC2:  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 COMPLEX*16 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 114 of file zgesc2.f.

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

◆ zgetc2()

subroutine zgetc2 ( integer n,
complex*16, dimension( lda, * ) a,
integer lda,
integer, dimension( * ) ipiv,
integer, dimension( * ) jpiv,
integer info )

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

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

Purpose:
!>
!> ZGETC2 computes an LU factorization, using 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 a level 1 BLAS version of the algorithm.
!>
Parameters
[in]N
!>          N is INTEGER
!>          The order of the matrix A. N >= 0.
!>
[in,out]A
!>          A is COMPLEX*16 array, dimension (LDA, N)
!>          On entry, the n-by-n matrix 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, 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
!>                one tries 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 zgetc2.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 COMPLEX*16 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 zgeru, zswap, dlabad
136* ..
137* .. External Functions ..
138 DOUBLE PRECISION DLAMCH
139 EXTERNAL dlamch
140* ..
141* .. Intrinsic Functions ..
142 INTRINSIC abs, dcmplx, 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 ) = dcmplx( smlnum, zero )
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 zswap( 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 zswap( 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 ) = dcmplx( smin, zero )
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 zgeru( n-i, n-i, -dcmplx( one ), a( i+1, i ), 1,
214 $ a( i, i+1 ), lda, 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 ) = dcmplx( smin, zero )
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 ZGETC2
230*
zswap
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP
Definition zswap.f:81
zgeru
subroutine zgeru(m, n, alpha, x, incx, y, incy, a, lda)
ZGERU
Definition zgeru.f:130
max
#define max(a, b)
Definition macros.h:21

◆ zlange()

double precision function zlange ( character norm,
integer m,
integer n,
complex*16, dimension( lda, * ) a,
integer lda,
double precision, dimension( * ) work )

ZLANGE 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 ZLANGE + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> ZLANGE  returns the value of the one norm,  or the Frobenius norm, or
!> the  infinity norm,  or the  element of  largest absolute value  of a
!> complex matrix A.
!>
Returns
ZLANGE 
!>
!>    ZLANGE = ( 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 ZLANGE as described
!>          above.
!>
[in]M
!>          M is INTEGER
!>          The number of rows of the matrix A.  M >= 0.  When M = 0,
!>          ZLANGE is set to zero.
!>
[in]N
!>          N is INTEGER
!>          The number of columns of the matrix A.  N >= 0.  When N = 0,
!>          ZLANGE is set to zero.
!>
[in]A
!>          A is COMPLEX*16 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 114 of file zlange.f.

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

◆ zlaqge()

subroutine zlaqge ( integer m,
integer n,
complex*16, dimension( lda, * ) a,
integer lda,
double precision, dimension( * ) r,
double precision, dimension( * ) c,
double precision rowcnd,
double precision colcnd,
double precision amax,
character equed )

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

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

Purpose:
!>
!> ZLAQGE 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 COMPLEX*16 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 141 of file zlaqge.f.

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

◆ ztgex2()

subroutine ztgex2 ( logical wantq,
logical wantz,
integer n,
complex*16, dimension( lda, * ) a,
integer lda,
complex*16, dimension( ldb, * ) b,
integer ldb,
complex*16, dimension( ldq, * ) q,
integer ldq,
complex*16, dimension( ldz, * ) z,
integer ldz,
integer j1,
integer info )

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

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

Purpose:
!>
!> ZTGEX2 swaps adjacent diagonal 1 by 1 blocks (A11,B11) and (A22,B22)
!> in an upper triangular matrix pair (A, B) by an unitary equivalence
!> transformation.
!>
!> (A, B) must be in generalized Schur canonical form, that is, A and
!> B are both upper triangular.
!>
!> Optionally, the matrices Q and Z of generalized Schur vectors are
!> updated.
!>
!>        Q(in) * A(in) * Z(in)**H = Q(out) * A(out) * Z(out)**H
!>        Q(in) * B(in) * Z(in)**H = Q(out) * B(out) * Z(out)**H
!>
!>
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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 array, dimension (LDQ,N)
!>          If WANTQ = .TRUE, on entry, the unitary 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 COMPLEX*16 array, dimension (LDZ,N)
!>          If WANTZ = .TRUE, on entry, the unitary 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).
!>
[out]INFO
!>          INFO is INTEGER
!>           =0:  Successful exit.
!>           =1:  The transformed matrix pair (A, B) would be too far
!>                from generalized Schur form; the problem is ill-
!>                conditioned.
!>
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 188 of file ztgex2.f.

190*
191* -- LAPACK auxiliary routine --
192* -- LAPACK is a software package provided by Univ. of Tennessee, --
193* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
194*
195* .. Scalar Arguments ..
196 LOGICAL WANTQ, WANTZ
197 INTEGER INFO, J1, LDA, LDB, LDQ, LDZ, N
198* ..
199* .. Array Arguments ..
200 COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
201 $ Z( LDZ, * )
202* ..
203*
204* =====================================================================
205*
206* .. Parameters ..
207 COMPLEX*16 CZERO, CONE
208 parameter( czero = ( 0.0d+0, 0.0d+0 ),
209 $ cone = ( 1.0d+0, 0.0d+0 ) )
210 DOUBLE PRECISION TWENTY
211 parameter( twenty = 2.0d+1 )
212 INTEGER LDST
213 parameter( ldst = 2 )
214 LOGICAL WANDS
215 parameter( wands = .true. )
216* ..
217* .. Local Scalars ..
218 LOGICAL STRONG, WEAK
219 INTEGER I, M
220 DOUBLE PRECISION CQ, CZ, EPS, SA, SB, SCALE, SMLNUM, SUM,
221 $ THRESHA, THRESHB
222 COMPLEX*16 CDUM, F, G, SQ, SZ
223* ..
224* .. Local Arrays ..
225 COMPLEX*16 S( LDST, LDST ), T( LDST, LDST ), WORK( 8 )
226* ..
227* .. External Functions ..
228 DOUBLE PRECISION DLAMCH
229 EXTERNAL dlamch
230* ..
231* .. External Subroutines ..
232 EXTERNAL zlacpy, zlartg, zlassq, zrot
233* ..
234* .. Intrinsic Functions ..
235 INTRINSIC abs, dble, dconjg, max, sqrt
236* ..
237* .. Executable Statements ..
238*
239 info = 0
240*
241* Quick return if possible
242*
243 IF( n.LE.1 )
244 $ RETURN
245*
246 m = ldst
247 weak = .false.
248 strong = .false.
249*
250* Make a local copy of selected block in (A, B)
251*
252 CALL zlacpy( 'Full', m, m, a( j1, j1 ), lda, s, ldst )
253 CALL zlacpy( 'Full', m, m, b( j1, j1 ), ldb, t, ldst )
254*
255* Compute the threshold for testing the acceptance of swapping.
256*
257 eps = dlamch( 'P' )
258 smlnum = dlamch( 'S' ) / eps
259 scale = dble( czero )
260 sum = dble( cone )
261 CALL zlacpy( 'Full', m, m, s, ldst, work, m )
262 CALL zlacpy( 'Full', m, m, t, ldst, work( m*m+1 ), m )
263 CALL zlassq( m*m, work, 1, scale, sum )
264 sa = scale*sqrt( sum )
265 scale = dble( czero )
266 sum = dble( cone )
267 CALL zlassq( m*m, work(m*m+1), 1, scale, sum )
268 sb = scale*sqrt( sum )
269*
270* THRES has been changed from
271* THRESH = MAX( TEN*EPS*SA, SMLNUM )
272* to
273* THRESH = MAX( TWENTY*EPS*SA, SMLNUM )
274* on 04/01/10.
275* "Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by
276* Jim Demmel and Guillaume Revy. See forum post 1783.
277*
278 thresha = max( twenty*eps*sa, smlnum )
279 threshb = max( twenty*eps*sb, smlnum )
280*
281* Compute unitary QL and RQ that swap 1-by-1 and 1-by-1 blocks
282* using Givens rotations and perform the swap tentatively.
283*
284 f = s( 2, 2 )*t( 1, 1 ) - t( 2, 2 )*s( 1, 1 )
285 g = s( 2, 2 )*t( 1, 2 ) - t( 2, 2 )*s( 1, 2 )
286 sa = abs( s( 2, 2 ) ) * abs( t( 1, 1 ) )
287 sb = abs( s( 1, 1 ) ) * abs( t( 2, 2 ) )
288 CALL zlartg( g, f, cz, sz, cdum )
289 sz = -sz
290 CALL zrot( 2, s( 1, 1 ), 1, s( 1, 2 ), 1, cz, dconjg( sz ) )
291 CALL zrot( 2, t( 1, 1 ), 1, t( 1, 2 ), 1, cz, dconjg( sz ) )
292 IF( sa.GE.sb ) THEN
293 CALL zlartg( s( 1, 1 ), s( 2, 1 ), cq, sq, cdum )
294 ELSE
295 CALL zlartg( t( 1, 1 ), t( 2, 1 ), cq, sq, cdum )
296 END IF
297 CALL zrot( 2, s( 1, 1 ), ldst, s( 2, 1 ), ldst, cq, sq )
298 CALL zrot( 2, t( 1, 1 ), ldst, t( 2, 1 ), ldst, cq, sq )
299*
300* Weak stability test: |S21| <= O(EPS F-norm((A)))
301* and |T21| <= O(EPS F-norm((B)))
302*
303 weak = abs( s( 2, 1 ) ).LE.thresha .AND.
304 $ abs( t( 2, 1 ) ).LE.threshb
305 IF( .NOT.weak )
306 $ GO TO 20
307*
308 IF( wands ) THEN
309*
310* Strong stability test:
311* F-norm((A-QL**H*S*QR)) <= O(EPS*F-norm((A)))
312* and
313* F-norm((B-QL**H*T*QR)) <= O(EPS*F-norm((B)))
314*
315 CALL zlacpy( 'Full', m, m, s, ldst, work, m )
316 CALL zlacpy( 'Full', m, m, t, ldst, work( m*m+1 ), m )
317 CALL zrot( 2, work, 1, work( 3 ), 1, cz, -dconjg( sz ) )
318 CALL zrot( 2, work( 5 ), 1, work( 7 ), 1, cz, -dconjg( sz ) )
319 CALL zrot( 2, work, 2, work( 2 ), 2, cq, -sq )
320 CALL zrot( 2, work( 5 ), 2, work( 6 ), 2, cq, -sq )
321 DO 10 i = 1, 2
322 work( i ) = work( i ) - a( j1+i-1, j1 )
323 work( i+2 ) = work( i+2 ) - a( j1+i-1, j1+1 )
324 work( i+4 ) = work( i+4 ) - b( j1+i-1, j1 )
325 work( i+6 ) = work( i+6 ) - b( j1+i-1, j1+1 )
326 10 CONTINUE
327 scale = dble( czero )
328 sum = dble( cone )
329 CALL zlassq( m*m, work, 1, scale, sum )
330 sa = scale*sqrt( sum )
331 scale = dble( czero )
332 sum = dble( cone )
333 CALL zlassq( m*m, work(m*m+1), 1, scale, sum )
334 sb = scale*sqrt( sum )
335 strong = sa.LE.thresha .AND. sb.LE.threshb
336 IF( .NOT.strong )
337 $ GO TO 20
338 END IF
339*
340* If the swap is accepted ("weakly" and "strongly"), apply the
341* equivalence transformations to the original matrix pair (A,B)
342*
343 CALL zrot( j1+1, a( 1, j1 ), 1, a( 1, j1+1 ), 1, cz,
344 $ dconjg( sz ) )
345 CALL zrot( j1+1, b( 1, j1 ), 1, b( 1, j1+1 ), 1, cz,
346 $ dconjg( sz ) )
347 CALL zrot( n-j1+1, a( j1, j1 ), lda, a( j1+1, j1 ), lda, cq, sq )
348 CALL zrot( n-j1+1, b( j1, j1 ), ldb, b( j1+1, j1 ), ldb, cq, sq )
349*
350* Set N1 by N2 (2,1) blocks to 0
351*
352 a( j1+1, j1 ) = czero
353 b( j1+1, j1 ) = czero
354*
355* Accumulate transformations into Q and Z if requested.
356*
357 IF( wantz )
358 $ CALL zrot( n, z( 1, j1 ), 1, z( 1, j1+1 ), 1, cz,
359 $ dconjg( sz ) )
360 IF( wantq )
361 $ CALL zrot( n, q( 1, j1 ), 1, q( 1, j1+1 ), 1, cq,
362 $ dconjg( sq ) )
363*
364* Exit with INFO = 0 if swap was successfully performed.
365*
366 RETURN
367*
368* Exit with INFO = 1 if swap was rejected.
369*
370 20 CONTINUE
371 info = 1
372 RETURN
373*
374* End of ZTGEX2
375*
zlartg
subroutine zlartg(f, g, c, s, r)
ZLARTG generates a plane rotation with real cosine and complex sine.
Definition zlartg.f90:118
zlacpy
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:103
zrot
subroutine zrot(n, cx, incx, cy, incy, c, s)
ZROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors.
Definition zrot.f:103