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

Functions

subroutine sgesc2 (n, a, lda, rhs, ipiv, jpiv, scale)
 SGESC2 solves a system of linear equations using the LU factorization with complete pivoting computed by sgetc2.
subroutine sgetc2 (n, a, lda, ipiv, jpiv, info)
 SGETC2 computes the LU factorization with complete pivoting of the general n-by-n matrix.
real function slange (norm, m, n, a, lda, work)
 SLANGE 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 slaqge (m, n, a, lda, r, c, rowcnd, colcnd, amax, equed)
 SLAQGE scales a general rectangular matrix, using row and column scaling factors computed by sgeequ.
subroutine stgex2 (wantq, wantz, n, a, lda, b, ldb, q, ldq, z, ldz, j1, n1, n2, work, lwork, info)
 STGEX2 swaps adjacent diagonal blocks in an upper (quasi) triangular matrix pair by an orthogonal equivalence transformation.

Detailed Description

This is the group of real auxiliary functions for GE matrices

Function Documentation

◆ sgesc2()

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

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

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

Purpose:
!>
!> SGESC2 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 SGETC2.
!>
Parameters
[in]N
!>          N is INTEGER
!>          The order of the matrix A.
!>
[in]A
!>          A is REAL array, dimension (LDA,N)
!>          On entry, the  LU part of the factorization of the n-by-n
!>          matrix A computed by SGETC2:  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 REAL 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 113 of file sgesc2.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 REAL SCALE
122* ..
123* .. Array Arguments ..
124 INTEGER IPIV( * ), JPIV( * )
125 REAL A( LDA, * ), RHS( * )
126* ..
127*
128* =====================================================================
129*
130* .. Parameters ..
131 REAL ONE, TWO
132 parameter( one = 1.0e+0, two = 2.0e+0 )
133* ..
134* .. Local Scalars ..
135 INTEGER I, J
136 REAL BIGNUM, EPS, SMLNUM, TEMP
137* ..
138* .. External Subroutines ..
139 EXTERNAL slabad, slaswp, sscal
140* ..
141* .. External Functions ..
142 INTEGER ISAMAX
143 REAL SLAMCH
144 EXTERNAL isamax, slamch
145* ..
146* .. Intrinsic Functions ..
147 INTRINSIC abs
148* ..
149* .. Executable Statements ..
150*
151* Set constant to control overflow
152*
153 eps = slamch( 'P' )
154 smlnum = slamch( 'S' ) / eps
155 bignum = one / smlnum
156 CALL slabad( smlnum, bignum )
157*
158* Apply permutations IPIV to RHS
159*
160 CALL slaswp( 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 = isamax( 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 sscal( 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 slaswp( 1, rhs, lda, 1, n-1, jpiv, -1 )
194 RETURN
195*
196* End of SGESC2
197*
slabad
subroutine slabad(small, large)
SLABAD
Definition slabad.f:74
isamax
integer function isamax(n, sx, incx)
ISAMAX
Definition isamax.f:71
slaswp
subroutine slaswp(n, a, lda, k1, k2, ipiv, incx)
SLASWP performs a series of row interchanges on a general rectangular matrix.
Definition slaswp.f:115
sscal
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68

◆ sgetc2()

subroutine sgetc2 ( integer n,
real, dimension( lda, * ) a,
integer lda,
integer, dimension( * ) ipiv,
integer, dimension( * ) jpiv,
integer info )

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

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

Purpose:
!>
!> SGETC2 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 REAL 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 sgetc2.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 REAL 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 sger, slabad, sswap
136* ..
137* .. External Functions ..
138 REAL SLAMCH
139 EXTERNAL slamch
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 = 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 ) = 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 sswap( 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 sswap( 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 sger( 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 SGETC2
230*
sswap
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82
sger
subroutine sger(m, n, alpha, x, incx, y, incy, a, lda)
SGER
Definition sger.f:130
max
#define max(a, b)
Definition macros.h:21

◆ slange()

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

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

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

◆ slaqge()

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

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

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

Purpose:
!>
!> SLAQGE 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 REAL 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 140 of file slaqge.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 REAL AMAX, COLCND, ROWCND
151* ..
152* .. Array Arguments ..
153 REAL A( LDA, * ), C( * ), R( * )
154* ..
155*
156* =====================================================================
157*
158* .. Parameters ..
159 REAL ONE, THRESH
160 parameter( one = 1.0e+0, thresh = 0.1e+0 )
161* ..
162* .. Local Scalars ..
163 INTEGER I, J
164 REAL CJ, LARGE, SMALL
165* ..
166* .. External Functions ..
167 REAL SLAMCH
168 EXTERNAL slamch
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 = slamch( 'Safe minimum' ) / slamch( '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 SLAQGE
232*

◆ stgex2()

subroutine stgex2 ( logical wantq,
logical wantz,
integer n,
real, dimension( lda, * ) a,
integer lda,
real, dimension( ldb, * ) b,
integer ldb,
real, dimension( ldq, * ) q,
integer ldq,
real, dimension( ldz, * ) z,
integer ldz,
integer j1,
integer n1,
integer n2,
real, dimension( * ) work,
integer lwork,
integer info )

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

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

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