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

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

Detailed Description

This is the group of real auxiliary functions for GE matrices

Function Documentation

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 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 116 of file sgesc2.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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 owerflow 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.
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 sgetc2.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  REAL a( lda, * )
125 * ..
126 *
127 * =====================================================================
128 *
129 * .. Parameters ..
130  REAL zero, one
131  parameter( zero = 0.0e+0, one = 1.0e+0 )
132 * ..
133 * .. Local Scalars ..
134  INTEGER i, ip, ipv, j, jp, jpv
135  REAL bignum, eps, smin, smlnum, xmax
136 * ..
137 * .. External Subroutines ..
138  EXTERNAL sger, slabad, sswap
139 * ..
140 * .. External Functions ..
141  REAL slamch
142  EXTERNAL slamch
143 * ..
144 * .. Intrinsic Functions ..
145  INTRINSIC abs, max
146 * ..
147 * .. Executable Statements ..
148 *
149 * Set constants to control overflow
150 *
151  info = 0
152  eps = slamch( 'P' )
153  smlnum = slamch( 'S' ) / eps
154  bignum = one / smlnum
155  CALL slabad( 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 sswap( 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 sswap( 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 ) = smin
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 sger( n-i, n-i, -one, a( i+1, i ), 1, a( i, i+1 ), lda,
199  $ 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 ) = smin
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 SGETC2
215 *
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine sger(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
SGER
Definition: sger.f:132
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:53

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 116 of file slange.f.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 144 of file slaqge.f.

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

Here is the caller graph for this function:

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 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 REAL 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 REAL array, dimension (LDZ,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.
Date
November 2015
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 223 of file stgex2.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: