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

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. More...
 
subroutine zgetc2 (N, A, LDA, IPIV, JPIV, INFO)
 ZGETC2 computes the LU factorization with complete pivoting of the general n-by-n matrix. More...
 
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. More...
 
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. More...
 
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. More...
 

Detailed Description

This is the group of complex16 auxiliary functions for GE matrices

Function Documentation

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 owerflow in the solution.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012
Contributors:
Bo Kagstrom and Peter Poromaa, Department of Computing Science, Umea University, S-901 87 Umea, Sweden.

Definition at line 117 of file zgesc2.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 113 of file zgetc2.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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.
Date
September 2012

Definition at line 117 of file zlange.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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.
Date
September 2012

Definition at line 145 of file zlaqge.f.

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

Here is the caller graph for this function:

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 arrays, 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 arrays, 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 (LDZ,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.
Date
September 2012
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 192 of file ztgex2.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: