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

Functions

subroutine cgesc2 (n, a, lda, rhs, ipiv, jpiv, scale)
 CGESC2 solves a system of linear equations using the LU factorization with complete pivoting computed by sgetc2.
subroutine cgetc2 (n, a, lda, ipiv, jpiv, info)
 CGETC2 computes the LU factorization with complete pivoting of the general n-by-n matrix.
real function clange (norm, m, n, a, lda, work)
 CLANGE 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 claqge (m, n, a, lda, r, c, rowcnd, colcnd, amax, equed)
 CLAQGE scales a general rectangular matrix, using row and column scaling factors computed by sgeequ.
subroutine ctgex2 (wantq, wantz, n, a, lda, b, ldb, q, ldq, z, ldz, j1, info)
 CTGEX2 swaps adjacent diagonal blocks in an upper (quasi) triangular matrix pair by an unitary equivalence transformation.

Detailed Description

This is the group of complex auxiliary functions for GE matrices

Function Documentation

◆ cgesc2()

subroutine cgesc2 ( integer n,
complex, dimension( lda, * ) a,
integer lda,
complex, dimension( * ) rhs,
integer, dimension( * ) ipiv,
integer, dimension( * ) jpiv,
real scale )

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

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

Purpose:
!>
!> CGESC2 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 CGETC2.
!>
!>
Parameters
[in]N
!>          N is INTEGER
!>          The number of columns of the matrix A.
!>
[in]A
!>          A is COMPLEX array, dimension (LDA, N)
!>          On entry, the  LU part of the factorization of the n-by-n
!>          matrix A computed by CGETC2:  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 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 REAL
!>           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 cgesc2.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 REAL SCALE
123* ..
124* .. Array Arguments ..
125 INTEGER IPIV( * ), JPIV( * )
126 COMPLEX A( LDA, * ), RHS( * )
127* ..
128*
129* =====================================================================
130*
131* .. Parameters ..
132 REAL ZERO, ONE, TWO
133 parameter( zero = 0.0e+0, one = 1.0e+0, two = 2.0e+0 )
134* ..
135* .. Local Scalars ..
136 INTEGER I, J
137 REAL BIGNUM, EPS, SMLNUM
138 COMPLEX TEMP
139* ..
140* .. External Subroutines ..
141 EXTERNAL claswp, cscal, slabad
142* ..
143* .. External Functions ..
144 INTEGER ICAMAX
145 REAL SLAMCH
146 EXTERNAL icamax, slamch
147* ..
148* .. Intrinsic Functions ..
149 INTRINSIC abs, cmplx, real
150* ..
151* .. Executable Statements ..
152*
153* Set constant to control overflow
154*
155 eps = slamch( 'P' )
156 smlnum = slamch( 'S' ) / eps
157 bignum = one / smlnum
158 CALL slabad( smlnum, bignum )
159*
160* Apply permutations IPIV to RHS
161*
162 CALL claswp( 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 = icamax( n, rhs, 1 )
179 IF( two*smlnum*abs( rhs( i ) ).GT.abs( a( n, n ) ) ) THEN
180 temp = cmplx( one / two, zero ) / abs( rhs( i ) )
181 CALL cscal( n, temp, rhs( 1 ), 1 )
182 scale = scale*real( temp )
183 END IF
184 DO 40 i = n, 1, -1
185 temp = cmplx( 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 claswp( 1, rhs, lda, 1, n-1, jpiv, -1 )
195 RETURN
196*
197* End of CGESC2
198*
cmplx
float cmplx[2]
Definition pblas.h:136
slabad
subroutine slabad(small, large)
SLABAD
Definition slabad.f:74
icamax
integer function icamax(n, cx, incx)
ICAMAX
Definition icamax.f:71
claswp
subroutine claswp(n, a, lda, k1, k2, ipiv, incx)
CLASWP performs a series of row interchanges on a general rectangular matrix.
Definition claswp.f:115
cscal
subroutine cscal(n, ca, cx, incx)
CSCAL
Definition cscal.f:78
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68

◆ cgetc2()

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

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

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

Purpose:
!>
!> CGETC2 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 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 cgetc2.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 A( LDA, * )
122* ..
123*
124* =====================================================================
125*
126* .. Parameters ..
127 REAL ZERO, ONE
128 parameter( zero = 0.0e+0, one = 1.0e+0 )
129* ..
130* .. Local Scalars ..
131 INTEGER I, IP, IPV, J, JP, JPV
132 REAL BIGNUM, EPS, SMIN, SMLNUM, XMAX
133* ..
134* .. External Subroutines ..
135 EXTERNAL cgeru, cswap, slabad
136* ..
137* .. External Functions ..
138 REAL SLAMCH
139 EXTERNAL slamch
140* ..
141* .. Intrinsic Functions ..
142 INTRINSIC abs, cmplx, 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 = slamch( 'P' )
156 smlnum = slamch( 'S' ) / eps
157 bignum = one / smlnum
158 CALL slabad( 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 ) = cmplx( 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 cswap( 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 cswap( 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 ) = cmplx( 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 cgeru( n-i, n-i, -cmplx( 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 ) = cmplx( 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 CGETC2
230*
cswap
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81
cgeru
subroutine cgeru(m, n, alpha, x, incx, y, incy, a, lda)
CGERU
Definition cgeru.f:130
max
#define max(a, b)
Definition macros.h:21

◆ clange()

real function clange ( character norm,
integer m,
integer n,
complex, dimension( lda, * ) a,
integer lda,
real, dimension( * ) work )

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

Purpose:
!>
!> CLANGE  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
CLANGE 
!>
!>    CLANGE = ( 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 CLANGE as described
!>          above.
!>
[in]M
!>          M is INTEGER
!>          The number of rows of the matrix A.  M >= 0.  When M = 0,
!>          CLANGE is set to zero.
!>
[in]N
!>          N is INTEGER
!>          The number of columns of the matrix A.  N >= 0.  When N = 0,
!>          CLANGE is set to zero.
!>
[in]A
!>          A is COMPLEX 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 REAL 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 clange.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 REAL WORK( * )
126 COMPLEX A( LDA, * )
127* ..
128*
129* =====================================================================
130*
131* .. Parameters ..
132 REAL ONE, ZERO
133 parameter( one = 1.0e+0, zero = 0.0e+0 )
134* ..
135* .. Local Scalars ..
136 INTEGER I, J
137 REAL SCALE, SUM, VALUE, TEMP
138* ..
139* .. External Functions ..
140 LOGICAL LSAME, SISNAN
141 EXTERNAL lsame, sisnan
142* ..
143* .. External Subroutines ..
144 EXTERNAL classq
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. sisnan( 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. sisnan( 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. sisnan( 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 classq( m, a( 1, j ), 1, scale, sum )
201 90 CONTINUE
202 VALUE = scale*sqrt( sum )
203 END IF
204*
205 clange = VALUE
206 RETURN
207*
208* End of CLANGE
209*
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
sisnan
logical function sisnan(sin)
SISNAN tests input for NaN.
Definition sisnan.f:59
classq
subroutine classq(n, x, incx, scl, sumsq)
CLASSQ updates a sum of squares represented in scaled form.
Definition classq.f90:137
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
clange
real function clange(norm, m, n, a, lda, work)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition clange.f:115
min
#define min(a, b)
Definition macros.h:20

◆ claqge()

subroutine claqge ( integer m,
integer n,
complex, dimension( lda, * ) a,
integer lda,
real, dimension( * ) r,
real, dimension( * ) c,
real rowcnd,
real colcnd,
real amax,
character equed )

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

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

Purpose:
!>
!> CLAQGE 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 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 REAL array, dimension (M)
!>          The row scale factors for A.
!>
[in]C
!>          C is REAL array, dimension (N)
!>          The column scale factors for A.
!>
[in]ROWCND
!>          ROWCND is REAL
!>          Ratio of the smallest R(i) to the largest R(i).
!>
[in]COLCND
!>          COLCND is REAL
!>          Ratio of the smallest C(i) to the largest C(i).
!>
[in]AMAX
!>          AMAX is REAL
!>          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 claqge.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 REAL AMAX, COLCND, ROWCND
152* ..
153* .. Array Arguments ..
154 REAL C( * ), R( * )
155 COMPLEX A( LDA, * )
156* ..
157*
158* =====================================================================
159*
160* .. Parameters ..
161 REAL ONE, THRESH
162 parameter( one = 1.0e+0, thresh = 0.1e+0 )
163* ..
164* .. Local Scalars ..
165 INTEGER I, J
166 REAL CJ, LARGE, SMALL
167* ..
168* .. External Functions ..
169 REAL SLAMCH
170 EXTERNAL slamch
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 = slamch( 'Safe minimum' ) / slamch( '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 CLAQGE
234*

◆ ctgex2()

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

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

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

Purpose:
!>
!> CTGEX2 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 array, dimension (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 array, dimension (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 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 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 ctgex2.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 A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
201 $ Z( LDZ, * )
202* ..
203*
204* =====================================================================
205*
206* .. Parameters ..
207 COMPLEX CZERO, CONE
208 parameter( czero = ( 0.0e+0, 0.0e+0 ),
209 $ cone = ( 1.0e+0, 0.0e+0 ) )
210 REAL TWENTY
211 parameter( twenty = 2.0e+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 REAL CQ, CZ, EPS, SA, SB, SCALE, SMLNUM, SUM,
221 $ THRESHA, THRESHB
222 COMPLEX CDUM, F, G, SQ, SZ
223* ..
224* .. Local Arrays ..
225 COMPLEX S( LDST, LDST ), T( LDST, LDST ), WORK( 8 )
226* ..
227* .. External Functions ..
228 REAL SLAMCH
229 EXTERNAL slamch
230* ..
231* .. External Subroutines ..
232 EXTERNAL clacpy, clartg, classq, crot
233* ..
234* .. Intrinsic Functions ..
235 INTRINSIC abs, conjg, max, real, 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 clacpy( 'Full', m, m, a( j1, j1 ), lda, s, ldst )
253 CALL clacpy( 'Full', m, m, b( j1, j1 ), ldb, t, ldst )
254*
255* Compute the threshold for testing the acceptance of swapping.
256*
257 eps = slamch( 'P' )
258 smlnum = slamch( 'S' ) / eps
259 scale = real( czero )
260 sum = real( cone )
261 CALL clacpy( 'Full', m, m, s, ldst, work, m )
262 CALL clacpy( 'Full', m, m, t, ldst, work( m*m+1 ), m )
263 CALL classq( m*m, work, 1, scale, sum )
264 sa = scale*sqrt( sum )
265 scale = dble( czero )
266 sum = dble( cone )
267 CALL classq( 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 clartg( g, f, cz, sz, cdum )
289 sz = -sz
290 CALL crot( 2, s( 1, 1 ), 1, s( 1, 2 ), 1, cz, conjg( sz ) )
291 CALL crot( 2, t( 1, 1 ), 1, t( 1, 2 ), 1, cz, conjg( sz ) )
292 IF( sa.GE.sb ) THEN
293 CALL clartg( s( 1, 1 ), s( 2, 1 ), cq, sq, cdum )
294 ELSE
295 CALL clartg( t( 1, 1 ), t( 2, 1 ), cq, sq, cdum )
296 END IF
297 CALL crot( 2, s( 1, 1 ), ldst, s( 2, 1 ), ldst, cq, sq )
298 CALL crot( 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, B-QL**H*T*QR)) <= O(EPS*F-norm((A, B)))
312*
313 CALL clacpy( 'Full', m, m, s, ldst, work, m )
314 CALL clacpy( 'Full', m, m, t, ldst, work( m*m+1 ), m )
315 CALL crot( 2, work, 1, work( 3 ), 1, cz, -conjg( sz ) )
316 CALL crot( 2, work( 5 ), 1, work( 7 ), 1, cz, -conjg( sz ) )
317 CALL crot( 2, work, 2, work( 2 ), 2, cq, -sq )
318 CALL crot( 2, work( 5 ), 2, work( 6 ), 2, cq, -sq )
319 DO 10 i = 1, 2
320 work( i ) = work( i ) - a( j1+i-1, j1 )
321 work( i+2 ) = work( i+2 ) - a( j1+i-1, j1+1 )
322 work( i+4 ) = work( i+4 ) - b( j1+i-1, j1 )
323 work( i+6 ) = work( i+6 ) - b( j1+i-1, j1+1 )
324 10 CONTINUE
325 scale = dble( czero )
326 sum = dble( cone )
327 CALL classq( m*m, work, 1, scale, sum )
328 sa = scale*sqrt( sum )
329 scale = dble( czero )
330 sum = dble( cone )
331 CALL classq( m*m, work(m*m+1), 1, scale, sum )
332 sb = scale*sqrt( sum )
333 strong = sa.LE.thresha .AND. sb.LE.threshb
334 IF( .NOT.strong )
335 $ GO TO 20
336 END IF
337*
338* If the swap is accepted ("weakly" and "strongly"), apply the
339* equivalence transformations to the original matrix pair (A,B)
340*
341 CALL crot( j1+1, a( 1, j1 ), 1, a( 1, j1+1 ), 1, cz, conjg( sz ) )
342 CALL crot( j1+1, b( 1, j1 ), 1, b( 1, j1+1 ), 1, cz, conjg( sz ) )
343 CALL crot( n-j1+1, a( j1, j1 ), lda, a( j1+1, j1 ), lda, cq, sq )
344 CALL crot( n-j1+1, b( j1, j1 ), ldb, b( j1+1, j1 ), ldb, cq, sq )
345*
346* Set N1 by N2 (2,1) blocks to 0
347*
348 a( j1+1, j1 ) = czero
349 b( j1+1, j1 ) = czero
350*
351* Accumulate transformations into Q and Z if requested.
352*
353 IF( wantz )
354 $ CALL crot( n, z( 1, j1 ), 1, z( 1, j1+1 ), 1, cz, conjg( sz ) )
355 IF( wantq )
356 $ CALL crot( n, q( 1, j1 ), 1, q( 1, j1+1 ), 1, cq, conjg( sq ) )
357*
358* Exit with INFO = 0 if swap was successfully performed.
359*
360 RETURN
361*
362* Exit with INFO = 1 if swap was rejected.
363*
364 20 CONTINUE
365 info = 1
366 RETURN
367*
368* End of CTGEX2
369*
clartg
subroutine clartg(f, g, c, s, r)
CLARTG generates a plane rotation with real cosine and complex sine.
Definition clartg.f90:118
crot
subroutine crot(n, cx, incx, cy, incy, c, s)
CROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors.
Definition crot.f:103
clacpy
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:103