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

Functions

real function slansy (NORM, UPLO, N, A, LDA, WORK)
 SLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real symmetric matrix. More...
 
subroutine slaqsy (UPLO, N, A, LDA, S, SCOND, AMAX, EQUED)
 SLAQSY scales a symmetric/Hermitian matrix, using scaling factors computed by spoequ. More...
 
subroutine slasy2 (LTRANL, LTRANR, ISGN, N1, N2, TL, LDTL, TR, LDTR, B, LDB, SCALE, X, LDX, XNORM, INFO)
 SLASY2 solves the Sylvester matrix equation where the matrices are of order 1 or 2. More...
 
subroutine ssyswapr (UPLO, N, A, LDA, I1, I2)
 SSYSWAPR applies an elementary permutation on the rows and columns of a symmetric matrix. More...
 
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). More...
 

Detailed Description

This is the group of real auxiliary functions for SY matrices

Function Documentation

real function slansy ( character  NORM,
character  UPLO,
integer  N,
real, dimension( lda, * )  A,
integer  LDA,
real, dimension( * )  WORK 
)

SLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real symmetric matrix.

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

Purpose:
 SLANSY  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 symmetric matrix A.
Returns
SLANSY
    SLANSY = ( 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 SLANSY as described
          above.
[in]UPLO
          UPLO is CHARACTER*1
          Specifies whether the upper or lower triangular part of the
          symmetric matrix A is to be referenced.
          = 'U':  Upper triangular part of A is referenced
          = 'L':  Lower triangular part of A is referenced
[in]N
          N is INTEGER
          The order of the matrix A.  N >= 0.  When N = 0, SLANSY is
          set to zero.
[in]A
          A is REAL array, dimension (LDA,N)
          The symmetric matrix A.  If UPLO = 'U', the leading n by n
          upper triangular part of A contains the upper triangular part
          of the matrix A, and the strictly lower triangular part of A
          is not referenced.  If UPLO = 'L', the leading n by n lower
          triangular part of A contains the lower triangular part of
          the matrix A, and the strictly upper triangular part of A is
          not referenced.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(N,1).
[out]WORK
          WORK is REAL array, dimension (MAX(1,LWORK)),
          where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise,
          WORK is not referenced.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2015

Definition at line 124 of file slansy.f.

124 *
125 * -- LAPACK auxiliary routine (version 3.6.0) --
126 * -- LAPACK is a software package provided by Univ. of Tennessee, --
127 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
128 * November 2015
129 *
130 * .. Scalar Arguments ..
131  CHARACTER norm, uplo
132  INTEGER lda, n
133 * ..
134 * .. Array Arguments ..
135  REAL a( lda, * ), work( * )
136 * ..
137 *
138 * =====================================================================
139 *
140 * .. Parameters ..
141  REAL one, zero
142  parameter( one = 1.0e+0, zero = 0.0e+0 )
143 * ..
144 * .. Local Scalars ..
145  INTEGER i, j
146  REAL absa, scale, sum, value
147 * ..
148 * .. External Subroutines ..
149  EXTERNAL slassq
150 * ..
151 * .. External Functions ..
152  LOGICAL lsame, sisnan
153  EXTERNAL lsame, sisnan
154 * ..
155 * .. Intrinsic Functions ..
156  INTRINSIC abs, sqrt
157 * ..
158 * .. Executable Statements ..
159 *
160  IF( n.EQ.0 ) THEN
161  VALUE = zero
162  ELSE IF( lsame( norm, 'M' ) ) THEN
163 *
164 * Find max(abs(A(i,j))).
165 *
166  VALUE = zero
167  IF( lsame( uplo, 'U' ) ) THEN
168  DO 20 j = 1, n
169  DO 10 i = 1, j
170  sum = abs( a( i, j ) )
171  IF( VALUE .LT. sum .OR. sisnan( sum ) ) VALUE = sum
172  10 CONTINUE
173  20 CONTINUE
174  ELSE
175  DO 40 j = 1, n
176  DO 30 i = j, n
177  sum = abs( a( i, j ) )
178  IF( VALUE .LT. sum .OR. sisnan( sum ) ) VALUE = sum
179  30 CONTINUE
180  40 CONTINUE
181  END IF
182  ELSE IF( ( lsame( norm, 'I' ) ) .OR. ( lsame( norm, 'O' ) ) .OR.
183  $ ( norm.EQ.'1' ) ) THEN
184 *
185 * Find normI(A) ( = norm1(A), since A is symmetric).
186 *
187  VALUE = zero
188  IF( lsame( uplo, 'U' ) ) THEN
189  DO 60 j = 1, n
190  sum = zero
191  DO 50 i = 1, j - 1
192  absa = abs( a( i, j ) )
193  sum = sum + absa
194  work( i ) = work( i ) + absa
195  50 CONTINUE
196  work( j ) = sum + abs( a( j, j ) )
197  60 CONTINUE
198  DO 70 i = 1, n
199  sum = work( i )
200  IF( VALUE .LT. sum .OR. sisnan( sum ) ) VALUE = sum
201  70 CONTINUE
202  ELSE
203  DO 80 i = 1, n
204  work( i ) = zero
205  80 CONTINUE
206  DO 100 j = 1, n
207  sum = work( j ) + abs( a( j, j ) )
208  DO 90 i = j + 1, n
209  absa = abs( a( i, j ) )
210  sum = sum + absa
211  work( i ) = work( i ) + absa
212  90 CONTINUE
213  IF( VALUE .LT. sum .OR. sisnan( sum ) ) VALUE = sum
214  100 CONTINUE
215  END IF
216  ELSE IF( ( lsame( norm, 'F' ) ) .OR. ( lsame( norm, 'E' ) ) ) THEN
217 *
218 * Find normF(A).
219 *
220  scale = zero
221  sum = one
222  IF( lsame( uplo, 'U' ) ) THEN
223  DO 110 j = 2, n
224  CALL slassq( j-1, a( 1, j ), 1, scale, sum )
225  110 CONTINUE
226  ELSE
227  DO 120 j = 1, n - 1
228  CALL slassq( n-j, a( j+1, j ), 1, scale, sum )
229  120 CONTINUE
230  END IF
231  sum = 2*sum
232  CALL slassq( n, a, lda+1, scale, sum )
233  VALUE = scale*sqrt( sum )
234  END IF
235 *
236  slansy = VALUE
237  RETURN
238 *
239 * End of SLANSY
240 *
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 slansy(NORM, UPLO, N, A, LDA, WORK)
SLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real symmetric matrix.
Definition: slansy.f:124

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine slaqsy ( character  UPLO,
integer  N,
real, dimension( lda, * )  A,
integer  LDA,
real, dimension( * )  S,
real  SCOND,
real  AMAX,
character  EQUED 
)

SLAQSY scales a symmetric/Hermitian matrix, using scaling factors computed by spoequ.

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

Purpose:
 SLAQSY equilibrates a symmetric matrix A using the scaling factors
 in the vector S.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          Specifies whether the upper or lower triangular part of the
          symmetric matrix A is stored.
          = 'U':  Upper triangular
          = 'L':  Lower triangular
[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 symmetric matrix A.  If UPLO = 'U', the leading
          n by n upper triangular part of A contains the upper
          triangular part of the matrix A, and the strictly lower
          triangular part of A is not referenced.  If UPLO = 'L', the
          leading n by n lower triangular part of A contains the lower
          triangular part of the matrix A, and the strictly upper
          triangular part of A is not referenced.

          On exit, if EQUED = 'Y', the equilibrated matrix:
          diag(S) * A * diag(S).
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(N,1).
[in]S
          S is REAL array, dimension (N)
          The scale factors for A.
[in]SCOND
          SCOND is REAL
          Ratio of the smallest S(i) to the largest S(i).
[in]AMAX
          AMAX is REAL
          Absolute value of largest matrix entry.
[out]EQUED
          EQUED is CHARACTER*1
          Specifies whether or not equilibration was done.
          = 'N':  No equilibration.
          = 'Y':  Equilibration was done, i.e., A has been replaced by
                  diag(S) * A * diag(S).
Internal Parameters:
  THRESH is a threshold value used to decide if scaling should be done
  based on the ratio of the scaling factors.  If SCOND < THRESH,
  scaling is done.

  LARGE and SMALL are threshold values used to decide if scaling should
  be done based on the absolute size of the largest matrix element.
  If AMAX > LARGE or AMAX < SMALL, scaling is done.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 135 of file slaqsy.f.

135 *
136 * -- LAPACK auxiliary routine (version 3.4.2) --
137 * -- LAPACK is a software package provided by Univ. of Tennessee, --
138 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
139 * September 2012
140 *
141 * .. Scalar Arguments ..
142  CHARACTER equed, uplo
143  INTEGER lda, n
144  REAL amax, scond
145 * ..
146 * .. Array Arguments ..
147  REAL a( lda, * ), s( * )
148 * ..
149 *
150 * =====================================================================
151 *
152 * .. Parameters ..
153  REAL one, thresh
154  parameter( one = 1.0e+0, thresh = 0.1e+0 )
155 * ..
156 * .. Local Scalars ..
157  INTEGER i, j
158  REAL cj, large, small
159 * ..
160 * .. External Functions ..
161  LOGICAL lsame
162  REAL slamch
163  EXTERNAL lsame, slamch
164 * ..
165 * .. Executable Statements ..
166 *
167 * Quick return if possible
168 *
169  IF( n.LE.0 ) THEN
170  equed = 'N'
171  RETURN
172  END IF
173 *
174 * Initialize LARGE and SMALL.
175 *
176  small = slamch( 'Safe minimum' ) / slamch( 'Precision' )
177  large = one / small
178 *
179  IF( scond.GE.thresh .AND. amax.GE.small .AND. amax.LE.large ) THEN
180 *
181 * No equilibration
182 *
183  equed = 'N'
184  ELSE
185 *
186 * Replace A by diag(S) * A * diag(S).
187 *
188  IF( lsame( uplo, 'U' ) ) THEN
189 *
190 * Upper triangle of A is stored.
191 *
192  DO 20 j = 1, n
193  cj = s( j )
194  DO 10 i = 1, j
195  a( i, j ) = cj*s( i )*a( i, j )
196  10 CONTINUE
197  20 CONTINUE
198  ELSE
199 *
200 * Lower triangle of A is stored.
201 *
202  DO 40 j = 1, n
203  cj = s( j )
204  DO 30 i = j, n
205  a( i, j ) = cj*s( i )*a( i, j )
206  30 CONTINUE
207  40 CONTINUE
208  END IF
209  equed = 'Y'
210  END IF
211 *
212  RETURN
213 *
214 * End of SLAQSY
215 *
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the caller graph for this function:

subroutine slasy2 ( logical  LTRANL,
logical  LTRANR,
integer  ISGN,
integer  N1,
integer  N2,
real, dimension( ldtl, * )  TL,
integer  LDTL,
real, dimension( ldtr, * )  TR,
integer  LDTR,
real, dimension( ldb, * )  B,
integer  LDB,
real  SCALE,
real, dimension( ldx, * )  X,
integer  LDX,
real  XNORM,
integer  INFO 
)

SLASY2 solves the Sylvester matrix equation where the matrices are of order 1 or 2.

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

Purpose:
 SLASY2 solves for the N1 by N2 matrix X, 1 <= N1,N2 <= 2, in

        op(TL)*X + ISGN*X*op(TR) = SCALE*B,

 where TL is N1 by N1, TR is N2 by N2, B is N1 by N2, and ISGN = 1 or
 -1.  op(T) = T or T**T, where T**T denotes the transpose of T.
Parameters
[in]LTRANL
          LTRANL is LOGICAL
          On entry, LTRANL specifies the op(TL):
             = .FALSE., op(TL) = TL,
             = .TRUE., op(TL) = TL**T.
[in]LTRANR
          LTRANR is LOGICAL
          On entry, LTRANR specifies the op(TR):
            = .FALSE., op(TR) = TR,
            = .TRUE., op(TR) = TR**T.
[in]ISGN
          ISGN is INTEGER
          On entry, ISGN specifies the sign of the equation
          as described before. ISGN may only be 1 or -1.
[in]N1
          N1 is INTEGER
          On entry, N1 specifies the order of matrix TL.
          N1 may only be 0, 1 or 2.
[in]N2
          N2 is INTEGER
          On entry, N2 specifies the order of matrix TR.
          N2 may only be 0, 1 or 2.
[in]TL
          TL is REAL array, dimension (LDTL,2)
          On entry, TL contains an N1 by N1 matrix.
[in]LDTL
          LDTL is INTEGER
          The leading dimension of the matrix TL. LDTL >= max(1,N1).
[in]TR
          TR is REAL array, dimension (LDTR,2)
          On entry, TR contains an N2 by N2 matrix.
[in]LDTR
          LDTR is INTEGER
          The leading dimension of the matrix TR. LDTR >= max(1,N2).
[in]B
          B is REAL array, dimension (LDB,2)
          On entry, the N1 by N2 matrix B contains the right-hand
          side of the equation.
[in]LDB
          LDB is INTEGER
          The leading dimension of the matrix B. LDB >= max(1,N1).
[out]SCALE
          SCALE is REAL
          On exit, SCALE contains the scale factor. SCALE is chosen
          less than or equal to 1 to prevent the solution overflowing.
[out]X
          X is REAL array, dimension (LDX,2)
          On exit, X contains the N1 by N2 solution.
[in]LDX
          LDX is INTEGER
          The leading dimension of the matrix X. LDX >= max(1,N1).
[out]XNORM
          XNORM is REAL
          On exit, XNORM is the infinity-norm of the solution.
[out]INFO
          INFO is INTEGER
          On exit, INFO is set to
             0: successful exit.
             1: TL and TR have too close eigenvalues, so TL or
                TR is perturbed to get a nonsingular equation.
          NOTE: In the interests of speed, this routine does not
                check the inputs for errors.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 176 of file slasy2.f.

176 *
177 * -- LAPACK auxiliary routine (version 3.4.2) --
178 * -- LAPACK is a software package provided by Univ. of Tennessee, --
179 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
180 * September 2012
181 *
182 * .. Scalar Arguments ..
183  LOGICAL ltranl, ltranr
184  INTEGER info, isgn, ldb, ldtl, ldtr, ldx, n1, n2
185  REAL scale, xnorm
186 * ..
187 * .. Array Arguments ..
188  REAL b( ldb, * ), tl( ldtl, * ), tr( ldtr, * ),
189  $ x( ldx, * )
190 * ..
191 *
192 * =====================================================================
193 *
194 * .. Parameters ..
195  REAL zero, one
196  parameter( zero = 0.0e+0, one = 1.0e+0 )
197  REAL two, half, eight
198  parameter( two = 2.0e+0, half = 0.5e+0, eight = 8.0e+0 )
199 * ..
200 * .. Local Scalars ..
201  LOGICAL bswap, xswap
202  INTEGER i, ip, ipiv, ipsv, j, jp, jpsv, k
203  REAL bet, eps, gam, l21, sgn, smin, smlnum, tau1,
204  $ temp, u11, u12, u22, xmax
205 * ..
206 * .. Local Arrays ..
207  LOGICAL bswpiv( 4 ), xswpiv( 4 )
208  INTEGER jpiv( 4 ), locl21( 4 ), locu12( 4 ),
209  $ locu22( 4 )
210  REAL btmp( 4 ), t16( 4, 4 ), tmp( 4 ), x2( 2 )
211 * ..
212 * .. External Functions ..
213  INTEGER isamax
214  REAL slamch
215  EXTERNAL isamax, slamch
216 * ..
217 * .. External Subroutines ..
218  EXTERNAL scopy, sswap
219 * ..
220 * .. Intrinsic Functions ..
221  INTRINSIC abs, max
222 * ..
223 * .. Data statements ..
224  DATA locu12 / 3, 4, 1, 2 / , locl21 / 2, 1, 4, 3 / ,
225  $ locu22 / 4, 3, 2, 1 /
226  DATA xswpiv / .false., .false., .true., .true. /
227  DATA bswpiv / .false., .true., .false., .true. /
228 * ..
229 * .. Executable Statements ..
230 *
231 * Do not check the input parameters for errors
232 *
233  info = 0
234 *
235 * Quick return if possible
236 *
237  IF( n1.EQ.0 .OR. n2.EQ.0 )
238  $ RETURN
239 *
240 * Set constants to control overflow
241 *
242  eps = slamch( 'P' )
243  smlnum = slamch( 'S' ) / eps
244  sgn = isgn
245 *
246  k = n1 + n1 + n2 - 2
247  GO TO ( 10, 20, 30, 50 )k
248 *
249 * 1 by 1: TL11*X + SGN*X*TR11 = B11
250 *
251  10 CONTINUE
252  tau1 = tl( 1, 1 ) + sgn*tr( 1, 1 )
253  bet = abs( tau1 )
254  IF( bet.LE.smlnum ) THEN
255  tau1 = smlnum
256  bet = smlnum
257  info = 1
258  END IF
259 *
260  scale = one
261  gam = abs( b( 1, 1 ) )
262  IF( smlnum*gam.GT.bet )
263  $ scale = one / gam
264 *
265  x( 1, 1 ) = ( b( 1, 1 )*scale ) / tau1
266  xnorm = abs( x( 1, 1 ) )
267  RETURN
268 *
269 * 1 by 2:
270 * TL11*[X11 X12] + ISGN*[X11 X12]*op[TR11 TR12] = [B11 B12]
271 * [TR21 TR22]
272 *
273  20 CONTINUE
274 *
275  smin = max( eps*max( abs( tl( 1, 1 ) ), abs( tr( 1, 1 ) ),
276  $ abs( tr( 1, 2 ) ), abs( tr( 2, 1 ) ), abs( tr( 2, 2 ) ) ),
277  $ smlnum )
278  tmp( 1 ) = tl( 1, 1 ) + sgn*tr( 1, 1 )
279  tmp( 4 ) = tl( 1, 1 ) + sgn*tr( 2, 2 )
280  IF( ltranr ) THEN
281  tmp( 2 ) = sgn*tr( 2, 1 )
282  tmp( 3 ) = sgn*tr( 1, 2 )
283  ELSE
284  tmp( 2 ) = sgn*tr( 1, 2 )
285  tmp( 3 ) = sgn*tr( 2, 1 )
286  END IF
287  btmp( 1 ) = b( 1, 1 )
288  btmp( 2 ) = b( 1, 2 )
289  GO TO 40
290 *
291 * 2 by 1:
292 * op[TL11 TL12]*[X11] + ISGN* [X11]*TR11 = [B11]
293 * [TL21 TL22] [X21] [X21] [B21]
294 *
295  30 CONTINUE
296  smin = max( eps*max( abs( tr( 1, 1 ) ), abs( tl( 1, 1 ) ),
297  $ abs( tl( 1, 2 ) ), abs( tl( 2, 1 ) ), abs( tl( 2, 2 ) ) ),
298  $ smlnum )
299  tmp( 1 ) = tl( 1, 1 ) + sgn*tr( 1, 1 )
300  tmp( 4 ) = tl( 2, 2 ) + sgn*tr( 1, 1 )
301  IF( ltranl ) THEN
302  tmp( 2 ) = tl( 1, 2 )
303  tmp( 3 ) = tl( 2, 1 )
304  ELSE
305  tmp( 2 ) = tl( 2, 1 )
306  tmp( 3 ) = tl( 1, 2 )
307  END IF
308  btmp( 1 ) = b( 1, 1 )
309  btmp( 2 ) = b( 2, 1 )
310  40 CONTINUE
311 *
312 * Solve 2 by 2 system using complete pivoting.
313 * Set pivots less than SMIN to SMIN.
314 *
315  ipiv = isamax( 4, tmp, 1 )
316  u11 = tmp( ipiv )
317  IF( abs( u11 ).LE.smin ) THEN
318  info = 1
319  u11 = smin
320  END IF
321  u12 = tmp( locu12( ipiv ) )
322  l21 = tmp( locl21( ipiv ) ) / u11
323  u22 = tmp( locu22( ipiv ) ) - u12*l21
324  xswap = xswpiv( ipiv )
325  bswap = bswpiv( ipiv )
326  IF( abs( u22 ).LE.smin ) THEN
327  info = 1
328  u22 = smin
329  END IF
330  IF( bswap ) THEN
331  temp = btmp( 2 )
332  btmp( 2 ) = btmp( 1 ) - l21*temp
333  btmp( 1 ) = temp
334  ELSE
335  btmp( 2 ) = btmp( 2 ) - l21*btmp( 1 )
336  END IF
337  scale = one
338  IF( ( two*smlnum )*abs( btmp( 2 ) ).GT.abs( u22 ) .OR.
339  $ ( two*smlnum )*abs( btmp( 1 ) ).GT.abs( u11 ) ) THEN
340  scale = half / max( abs( btmp( 1 ) ), abs( btmp( 2 ) ) )
341  btmp( 1 ) = btmp( 1 )*scale
342  btmp( 2 ) = btmp( 2 )*scale
343  END IF
344  x2( 2 ) = btmp( 2 ) / u22
345  x2( 1 ) = btmp( 1 ) / u11 - ( u12 / u11 )*x2( 2 )
346  IF( xswap ) THEN
347  temp = x2( 2 )
348  x2( 2 ) = x2( 1 )
349  x2( 1 ) = temp
350  END IF
351  x( 1, 1 ) = x2( 1 )
352  IF( n1.EQ.1 ) THEN
353  x( 1, 2 ) = x2( 2 )
354  xnorm = abs( x( 1, 1 ) ) + abs( x( 1, 2 ) )
355  ELSE
356  x( 2, 1 ) = x2( 2 )
357  xnorm = max( abs( x( 1, 1 ) ), abs( x( 2, 1 ) ) )
358  END IF
359  RETURN
360 *
361 * 2 by 2:
362 * op[TL11 TL12]*[X11 X12] +ISGN* [X11 X12]*op[TR11 TR12] = [B11 B12]
363 * [TL21 TL22] [X21 X22] [X21 X22] [TR21 TR22] [B21 B22]
364 *
365 * Solve equivalent 4 by 4 system using complete pivoting.
366 * Set pivots less than SMIN to SMIN.
367 *
368  50 CONTINUE
369  smin = max( abs( tr( 1, 1 ) ), abs( tr( 1, 2 ) ),
370  $ abs( tr( 2, 1 ) ), abs( tr( 2, 2 ) ) )
371  smin = max( smin, abs( tl( 1, 1 ) ), abs( tl( 1, 2 ) ),
372  $ abs( tl( 2, 1 ) ), abs( tl( 2, 2 ) ) )
373  smin = max( eps*smin, smlnum )
374  btmp( 1 ) = zero
375  CALL scopy( 16, btmp, 0, t16, 1 )
376  t16( 1, 1 ) = tl( 1, 1 ) + sgn*tr( 1, 1 )
377  t16( 2, 2 ) = tl( 2, 2 ) + sgn*tr( 1, 1 )
378  t16( 3, 3 ) = tl( 1, 1 ) + sgn*tr( 2, 2 )
379  t16( 4, 4 ) = tl( 2, 2 ) + sgn*tr( 2, 2 )
380  IF( ltranl ) THEN
381  t16( 1, 2 ) = tl( 2, 1 )
382  t16( 2, 1 ) = tl( 1, 2 )
383  t16( 3, 4 ) = tl( 2, 1 )
384  t16( 4, 3 ) = tl( 1, 2 )
385  ELSE
386  t16( 1, 2 ) = tl( 1, 2 )
387  t16( 2, 1 ) = tl( 2, 1 )
388  t16( 3, 4 ) = tl( 1, 2 )
389  t16( 4, 3 ) = tl( 2, 1 )
390  END IF
391  IF( ltranr ) THEN
392  t16( 1, 3 ) = sgn*tr( 1, 2 )
393  t16( 2, 4 ) = sgn*tr( 1, 2 )
394  t16( 3, 1 ) = sgn*tr( 2, 1 )
395  t16( 4, 2 ) = sgn*tr( 2, 1 )
396  ELSE
397  t16( 1, 3 ) = sgn*tr( 2, 1 )
398  t16( 2, 4 ) = sgn*tr( 2, 1 )
399  t16( 3, 1 ) = sgn*tr( 1, 2 )
400  t16( 4, 2 ) = sgn*tr( 1, 2 )
401  END IF
402  btmp( 1 ) = b( 1, 1 )
403  btmp( 2 ) = b( 2, 1 )
404  btmp( 3 ) = b( 1, 2 )
405  btmp( 4 ) = b( 2, 2 )
406 *
407 * Perform elimination
408 *
409  DO 100 i = 1, 3
410  xmax = zero
411  DO 70 ip = i, 4
412  DO 60 jp = i, 4
413  IF( abs( t16( ip, jp ) ).GE.xmax ) THEN
414  xmax = abs( t16( ip, jp ) )
415  ipsv = ip
416  jpsv = jp
417  END IF
418  60 CONTINUE
419  70 CONTINUE
420  IF( ipsv.NE.i ) THEN
421  CALL sswap( 4, t16( ipsv, 1 ), 4, t16( i, 1 ), 4 )
422  temp = btmp( i )
423  btmp( i ) = btmp( ipsv )
424  btmp( ipsv ) = temp
425  END IF
426  IF( jpsv.NE.i )
427  $ CALL sswap( 4, t16( 1, jpsv ), 1, t16( 1, i ), 1 )
428  jpiv( i ) = jpsv
429  IF( abs( t16( i, i ) ).LT.smin ) THEN
430  info = 1
431  t16( i, i ) = smin
432  END IF
433  DO 90 j = i + 1, 4
434  t16( j, i ) = t16( j, i ) / t16( i, i )
435  btmp( j ) = btmp( j ) - t16( j, i )*btmp( i )
436  DO 80 k = i + 1, 4
437  t16( j, k ) = t16( j, k ) - t16( j, i )*t16( i, k )
438  80 CONTINUE
439  90 CONTINUE
440  100 CONTINUE
441  IF( abs( t16( 4, 4 ) ).LT.smin )
442  $ t16( 4, 4 ) = smin
443  scale = one
444  IF( ( eight*smlnum )*abs( btmp( 1 ) ).GT.abs( t16( 1, 1 ) ) .OR.
445  $ ( eight*smlnum )*abs( btmp( 2 ) ).GT.abs( t16( 2, 2 ) ) .OR.
446  $ ( eight*smlnum )*abs( btmp( 3 ) ).GT.abs( t16( 3, 3 ) ) .OR.
447  $ ( eight*smlnum )*abs( btmp( 4 ) ).GT.abs( t16( 4, 4 ) ) ) THEN
448  scale = ( one / eight ) / max( abs( btmp( 1 ) ),
449  $ abs( btmp( 2 ) ), abs( btmp( 3 ) ), abs( btmp( 4 ) ) )
450  btmp( 1 ) = btmp( 1 )*scale
451  btmp( 2 ) = btmp( 2 )*scale
452  btmp( 3 ) = btmp( 3 )*scale
453  btmp( 4 ) = btmp( 4 )*scale
454  END IF
455  DO 120 i = 1, 4
456  k = 5 - i
457  temp = one / t16( k, k )
458  tmp( k ) = btmp( k )*temp
459  DO 110 j = k + 1, 4
460  tmp( k ) = tmp( k ) - ( temp*t16( k, j ) )*tmp( j )
461  110 CONTINUE
462  120 CONTINUE
463  DO 130 i = 1, 3
464  IF( jpiv( 4-i ).NE.4-i ) THEN
465  temp = tmp( 4-i )
466  tmp( 4-i ) = tmp( jpiv( 4-i ) )
467  tmp( jpiv( 4-i ) ) = temp
468  END IF
469  130 CONTINUE
470  x( 1, 1 ) = tmp( 1 )
471  x( 2, 1 ) = tmp( 2 )
472  x( 1, 2 ) = tmp( 3 )
473  x( 2, 2 ) = tmp( 4 )
474  xnorm = max( abs( tmp( 1 ) )+abs( tmp( 3 ) ),
475  $ abs( tmp( 2 ) )+abs( tmp( 4 ) ) )
476  RETURN
477 *
478 * End of SLASY2
479 *
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:53
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53
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:

subroutine ssyswapr ( character  UPLO,
integer  N,
real, dimension( lda, n )  A,
integer  LDA,
integer  I1,
integer  I2 
)

SSYSWAPR applies an elementary permutation on the rows and columns of a symmetric matrix.

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

Purpose:
 SSYSWAPR applies an elementary permutation on the rows and the columns of
 a symmetric matrix.
Parameters
[in]UPLO
          UPLO is CHARACTER*1
          Specifies whether the details of the factorization are stored
          as an upper or lower triangular matrix.
          = 'U':  Upper triangular, form is A = U*D*U**T;
          = 'L':  Lower triangular, form is A = L*D*L**T.
[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 NB diagonal matrix D and the multipliers
          used to obtain the factor U or L as computed by SSYTRF.

          On exit, if INFO = 0, the (symmetric) inverse of the original
          matrix.  If UPLO = 'U', the upper triangular part of the
          inverse is formed and the part of A below the diagonal is not
          referenced; if UPLO = 'L' the lower triangular part of the
          inverse is formed and the part of A above the diagonal is
          not referenced.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[in]I1
          I1 is INTEGER
          Index of the first row to swap
[in]I2
          I2 is INTEGER
          Index of the second row to swap
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
September 2012

Definition at line 104 of file ssyswapr.f.

104 *
105 * -- LAPACK auxiliary routine (version 3.4.2) --
106 * -- LAPACK is a software package provided by Univ. of Tennessee, --
107 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
108 * September 2012
109 *
110 * .. Scalar Arguments ..
111  CHARACTER uplo
112  INTEGER i1, i2, lda, n
113 * ..
114 * .. Array Arguments ..
115  REAL a( lda, n )
116 *
117 * =====================================================================
118 *
119 * ..
120 * .. Local Scalars ..
121  LOGICAL upper
122  INTEGER i
123  REAL tmp
124 *
125 * .. External Functions ..
126  LOGICAL lsame
127  EXTERNAL lsame
128 * ..
129 * .. External Subroutines ..
130  EXTERNAL sswap
131 * ..
132 * .. Executable Statements ..
133 *
134  upper = lsame( uplo, 'U' )
135  IF (upper) THEN
136 *
137 * UPPER
138 * first swap
139 * - swap column I1 and I2 from I1 to I1-1
140  CALL sswap( i1-1, a(1,i1), 1, a(1,i2), 1 )
141 *
142 * second swap :
143 * - swap A(I1,I1) and A(I2,I2)
144 * - swap row I1 from I1+1 to I2-1 with col I2 from I1+1 to I2-1
145  tmp=a(i1,i1)
146  a(i1,i1)=a(i2,i2)
147  a(i2,i2)=tmp
148 *
149  DO i=1,i2-i1-1
150  tmp=a(i1,i1+i)
151  a(i1,i1+i)=a(i1+i,i2)
152  a(i1+i,i2)=tmp
153  END DO
154 *
155 * third swap
156 * - swap row I1 and I2 from I2+1 to N
157  DO i=i2+1,n
158  tmp=a(i1,i)
159  a(i1,i)=a(i2,i)
160  a(i2,i)=tmp
161  END DO
162 *
163  ELSE
164 *
165 * LOWER
166 * first swap
167 * - swap row I1 and I2 from I1 to I1-1
168  CALL sswap( i1-1, a(i1,1), lda, a(i2,1), lda )
169 *
170 * second swap :
171 * - swap A(I1,I1) and A(I2,I2)
172 * - swap col I1 from I1+1 to I2-1 with row I2 from I1+1 to I2-1
173  tmp=a(i1,i1)
174  a(i1,i1)=a(i2,i2)
175  a(i2,i2)=tmp
176 *
177  DO i=1,i2-i1-1
178  tmp=a(i1+i,i1)
179  a(i1+i,i1)=a(i2,i1+i)
180  a(i2,i1+i)=tmp
181  END DO
182 *
183 * third swap
184 * - swap col I1 and I2 from I2+1 to N
185  DO i=i2+1,n
186  tmp=a(i,i1)
187  a(i,i1)=a(i,i2)
188  a(i,i2)=tmp
189  END DO
190 *
191  ENDIF
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
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:

subroutine stgsy2 ( character  TRANS,
integer  IJOB,
integer  M,
integer  N,
real, dimension( lda, * )  A,
integer  LDA,
real, dimension( ldb, * )  B,
integer  LDB,
real, dimension( ldc, * )  C,
integer  LDC,
real, dimension( ldd, * )  D,
integer  LDD,
real, dimension( lde, * )  E,
integer  LDE,
real, dimension( ldf, * )  F,
integer  LDF,
real  SCALE,
real  RDSUM,
real  RDSCAL,
integer, dimension( * )  IWORK,
integer  PQ,
integer  INFO 
)

STGSY2 solves the generalized Sylvester equation (unblocked algorithm).

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

Purpose:
 STGSY2 solves the generalized Sylvester equation:

             A * R - L * B = scale * C                (1)
             D * R - L * E = scale * F,

 using Level 1 and 2 BLAS. where R and L are unknown M-by-N matrices,
 (A, D), (B, E) and (C, F) are given matrix pairs of size M-by-M,
 N-by-N and M-by-N, respectively, with real entries. (A, D) and (B, E)
 must be in generalized Schur canonical form, i.e. A, B are upper
 quasi triangular and D, E are upper triangular. The solution (R, L)
 overwrites (C, F). 0 <= SCALE <= 1 is an output scaling factor
 chosen to avoid overflow.

 In matrix notation solving equation (1) corresponds to solve
 Z*x = scale*b, where Z is defined as

        Z = [ kron(In, A)  -kron(B**T, Im) ]             (2)
            [ kron(In, D)  -kron(E**T, Im) ],

 Ik is the identity matrix of size k and X**T is the transpose of X.
 kron(X, Y) is the Kronecker product between the matrices X and Y.
 In the process of solving (1), we solve a number of such systems
 where Dim(In), Dim(In) = 1 or 2.

 If TRANS = 'T', solve the transposed system Z**T*y = scale*b for y,
 which is equivalent to solve for R and L in

             A**T * R  + D**T * L   = scale * C           (3)
             R  * B**T + L  * E**T  = scale * -F

 This case is used to compute an estimate of Dif[(A, D), (B, E)] =
 sigma_min(Z) using reverse communicaton with SLACON.

 STGSY2 also (IJOB >= 1) contributes to the computation in STGSYL
 of an upper bound on the separation between to matrix pairs. Then
 the input (A, D), (B, E) are sub-pencils of the matrix pair in
 STGSYL. See STGSYL for details.
Parameters
[in]TRANS
          TRANS is CHARACTER*1
          = 'N', solve the generalized Sylvester equation (1).
          = 'T': solve the 'transposed' system (3).
[in]IJOB
          IJOB is INTEGER
          Specifies what kind of functionality to be performed.
          = 0: solve (1) only.
          = 1: A contribution from this subsystem to a Frobenius
               norm-based estimate of the separation between two matrix
               pairs is computed. (look ahead strategy is used).
          = 2: A contribution from this subsystem to a Frobenius
               norm-based estimate of the separation between two matrix
               pairs is computed. (SGECON on sub-systems is used.)
          Not referenced if TRANS = 'T'.
[in]M
          M is INTEGER
          On entry, M specifies the order of A and D, and the row
          dimension of C, F, R and L.
[in]N
          N is INTEGER
          On entry, N specifies the order of B and E, and the column
          dimension of C, F, R and L.
[in]A
          A is REAL array, dimension (LDA, M)
          On entry, A contains an upper quasi triangular matrix.
[in]LDA
          LDA is INTEGER
          The leading dimension of the matrix A. LDA >= max(1, M).
[in]B
          B is REAL array, dimension (LDB, N)
          On entry, B contains an upper quasi triangular matrix.
[in]LDB
          LDB is INTEGER
          The leading dimension of the matrix B. LDB >= max(1, N).
[in,out]C
          C is REAL array, dimension (LDC, N)
          On entry, C contains the right-hand-side of the first matrix
          equation in (1).
          On exit, if IJOB = 0, C has been overwritten by the
          solution R.
[in]LDC
          LDC is INTEGER
          The leading dimension of the matrix C. LDC >= max(1, M).
[in]D
          D is REAL array, dimension (LDD, M)
          On entry, D contains an upper triangular matrix.
[in]LDD
          LDD is INTEGER
          The leading dimension of the matrix D. LDD >= max(1, M).
[in]E
          E is REAL array, dimension (LDE, N)
          On entry, E contains an upper triangular matrix.
[in]LDE
          LDE is INTEGER
          The leading dimension of the matrix E. LDE >= max(1, N).
[in,out]F
          F is REAL array, dimension (LDF, N)
          On entry, F contains the right-hand-side of the second matrix
          equation in (1).
          On exit, if IJOB = 0, F has been overwritten by the
          solution L.
[in]LDF
          LDF is INTEGER
          The leading dimension of the matrix F. LDF >= max(1, M).
[out]SCALE
          SCALE is REAL
          On exit, 0 <= SCALE <= 1. If 0 < SCALE < 1, the solutions
          R and L (C and F on entry) will hold the solutions to a
          slightly perturbed system but the input matrices A, B, D and
          E have not been changed. If SCALE = 0, R and L will hold the
          solutions to the homogeneous system with C = F = 0. Normally,
          SCALE = 1.
[in,out]RDSUM
          RDSUM is REAL
          On entry, the sum of squares of computed contributions to
          the Dif-estimate under computation by STGSYL, where the
          scaling factor RDSCAL (see below) has been factored out.
          On exit, the corresponding sum of squares updated with the
          contributions from the current sub-system.
          If TRANS = 'T' RDSUM is not touched.
          NOTE: RDSUM only makes sense when STGSY2 is called by STGSYL.
[in,out]RDSCAL
          RDSCAL is REAL
          On entry, scaling factor used to prevent overflow in RDSUM.
          On exit, RDSCAL is updated w.r.t. the current contributions
          in RDSUM.
          If TRANS = 'T', RDSCAL is not touched.
          NOTE: RDSCAL only makes sense when STGSY2 is called by
                STGSYL.
[out]IWORK
          IWORK is INTEGER array, dimension (M+N+2)
[out]PQ
          PQ is INTEGER
          On exit, the number of subsystems (of size 2-by-2, 4-by-4 and
          8-by-8) solved by this routine.
[out]INFO
          INFO is INTEGER
          On exit, if INFO is set to
            =0: Successful exit
            <0: If INFO = -i, the i-th argument had an illegal value.
            >0: The matrix pairs (A, D) and (B, E) have common or very
                close eigenvalues.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2015
Contributors:
Bo Kagstrom and Peter Poromaa, Department of Computing Science, Umea University, S-901 87 Umea, Sweden.

Definition at line 276 of file stgsy2.f.

276 *
277 * -- LAPACK auxiliary routine (version 3.6.0) --
278 * -- LAPACK is a software package provided by Univ. of Tennessee, --
279 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
280 * November 2015
281 *
282 * .. Scalar Arguments ..
283  CHARACTER trans
284  INTEGER ijob, info, lda, ldb, ldc, ldd, lde, ldf, m, n,
285  $ pq
286  REAL rdscal, rdsum, scale
287 * ..
288 * .. Array Arguments ..
289  INTEGER iwork( * )
290  REAL a( lda, * ), b( ldb, * ), c( ldc, * ),
291  $ d( ldd, * ), e( lde, * ), f( ldf, * )
292 * ..
293 *
294 * =====================================================================
295 * Replaced various illegal calls to SCOPY by calls to SLASET.
296 * Sven Hammarling, 27/5/02.
297 *
298 * .. Parameters ..
299  INTEGER ldz
300  parameter( ldz = 8 )
301  REAL zero, one
302  parameter( zero = 0.0e+0, one = 1.0e+0 )
303 * ..
304 * .. Local Scalars ..
305  LOGICAL notran
306  INTEGER i, ie, ierr, ii, is, isp1, j, je, jj, js, jsp1,
307  $ k, mb, nb, p, q, zdim
308  REAL alpha, scaloc
309 * ..
310 * .. Local Arrays ..
311  INTEGER ipiv( ldz ), jpiv( ldz )
312  REAL rhs( ldz ), z( ldz, ldz )
313 * ..
314 * .. External Functions ..
315  LOGICAL lsame
316  EXTERNAL lsame
317 * ..
318 * .. External Subroutines ..
319  EXTERNAL saxpy, scopy, sgemm, sgemv, sger, sgesc2,
321 * ..
322 * .. Intrinsic Functions ..
323  INTRINSIC max
324 * ..
325 * .. Executable Statements ..
326 *
327 * Decode and test input parameters
328 *
329  info = 0
330  ierr = 0
331  notran = lsame( trans, 'N' )
332  IF( .NOT.notran .AND. .NOT.lsame( trans, 'T' ) ) THEN
333  info = -1
334  ELSE IF( notran ) THEN
335  IF( ( ijob.LT.0 ) .OR. ( ijob.GT.2 ) ) THEN
336  info = -2
337  END IF
338  END IF
339  IF( info.EQ.0 ) THEN
340  IF( m.LE.0 ) THEN
341  info = -3
342  ELSE IF( n.LE.0 ) THEN
343  info = -4
344  ELSE IF( lda.LT.max( 1, m ) ) THEN
345  info = -6
346  ELSE IF( ldb.LT.max( 1, n ) ) THEN
347  info = -8
348  ELSE IF( ldc.LT.max( 1, m ) ) THEN
349  info = -10
350  ELSE IF( ldd.LT.max( 1, m ) ) THEN
351  info = -12
352  ELSE IF( lde.LT.max( 1, n ) ) THEN
353  info = -14
354  ELSE IF( ldf.LT.max( 1, m ) ) THEN
355  info = -16
356  END IF
357  END IF
358  IF( info.NE.0 ) THEN
359  CALL xerbla( 'STGSY2', -info )
360  RETURN
361  END IF
362 *
363 * Determine block structure of A
364 *
365  pq = 0
366  p = 0
367  i = 1
368  10 CONTINUE
369  IF( i.GT.m )
370  $ GO TO 20
371  p = p + 1
372  iwork( p ) = i
373  IF( i.EQ.m )
374  $ GO TO 20
375  IF( a( i+1, i ).NE.zero ) THEN
376  i = i + 2
377  ELSE
378  i = i + 1
379  END IF
380  GO TO 10
381  20 CONTINUE
382  iwork( p+1 ) = m + 1
383 *
384 * Determine block structure of B
385 *
386  q = p + 1
387  j = 1
388  30 CONTINUE
389  IF( j.GT.n )
390  $ GO TO 40
391  q = q + 1
392  iwork( q ) = j
393  IF( j.EQ.n )
394  $ GO TO 40
395  IF( b( j+1, j ).NE.zero ) THEN
396  j = j + 2
397  ELSE
398  j = j + 1
399  END IF
400  GO TO 30
401  40 CONTINUE
402  iwork( q+1 ) = n + 1
403  pq = p*( q-p-1 )
404 *
405  IF( notran ) THEN
406 *
407 * Solve (I, J) - subsystem
408 * A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
409 * D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
410 * for I = P, P - 1, ..., 1; J = 1, 2, ..., Q
411 *
412  scale = one
413  scaloc = one
414  DO 120 j = p + 2, q
415  js = iwork( j )
416  jsp1 = js + 1
417  je = iwork( j+1 ) - 1
418  nb = je - js + 1
419  DO 110 i = p, 1, -1
420 *
421  is = iwork( i )
422  isp1 = is + 1
423  ie = iwork( i+1 ) - 1
424  mb = ie - is + 1
425  zdim = mb*nb*2
426 *
427  IF( ( mb.EQ.1 ) .AND. ( nb.EQ.1 ) ) THEN
428 *
429 * Build a 2-by-2 system Z * x = RHS
430 *
431  z( 1, 1 ) = a( is, is )
432  z( 2, 1 ) = d( is, is )
433  z( 1, 2 ) = -b( js, js )
434  z( 2, 2 ) = -e( js, js )
435 *
436 * Set up right hand side(s)
437 *
438  rhs( 1 ) = c( is, js )
439  rhs( 2 ) = f( is, js )
440 *
441 * Solve Z * x = RHS
442 *
443  CALL sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
444  IF( ierr.GT.0 )
445  $ info = ierr
446 *
447  IF( ijob.EQ.0 ) THEN
448  CALL sgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
449  $ scaloc )
450  IF( scaloc.NE.one ) THEN
451  DO 50 k = 1, n
452  CALL sscal( m, scaloc, c( 1, k ), 1 )
453  CALL sscal( m, scaloc, f( 1, k ), 1 )
454  50 CONTINUE
455  scale = scale*scaloc
456  END IF
457  ELSE
458  CALL slatdf( ijob, zdim, z, ldz, rhs, rdsum,
459  $ rdscal, ipiv, jpiv )
460  END IF
461 *
462 * Unpack solution vector(s)
463 *
464  c( is, js ) = rhs( 1 )
465  f( is, js ) = rhs( 2 )
466 *
467 * Substitute R(I, J) and L(I, J) into remaining
468 * equation.
469 *
470  IF( i.GT.1 ) THEN
471  alpha = -rhs( 1 )
472  CALL saxpy( is-1, alpha, a( 1, is ), 1, c( 1, js ),
473  $ 1 )
474  CALL saxpy( is-1, alpha, d( 1, is ), 1, f( 1, js ),
475  $ 1 )
476  END IF
477  IF( j.LT.q ) THEN
478  CALL saxpy( n-je, rhs( 2 ), b( js, je+1 ), ldb,
479  $ c( is, je+1 ), ldc )
480  CALL saxpy( n-je, rhs( 2 ), e( js, je+1 ), lde,
481  $ f( is, je+1 ), ldf )
482  END IF
483 *
484  ELSE IF( ( mb.EQ.1 ) .AND. ( nb.EQ.2 ) ) THEN
485 *
486 * Build a 4-by-4 system Z * x = RHS
487 *
488  z( 1, 1 ) = a( is, is )
489  z( 2, 1 ) = zero
490  z( 3, 1 ) = d( is, is )
491  z( 4, 1 ) = zero
492 *
493  z( 1, 2 ) = zero
494  z( 2, 2 ) = a( is, is )
495  z( 3, 2 ) = zero
496  z( 4, 2 ) = d( is, is )
497 *
498  z( 1, 3 ) = -b( js, js )
499  z( 2, 3 ) = -b( js, jsp1 )
500  z( 3, 3 ) = -e( js, js )
501  z( 4, 3 ) = -e( js, jsp1 )
502 *
503  z( 1, 4 ) = -b( jsp1, js )
504  z( 2, 4 ) = -b( jsp1, jsp1 )
505  z( 3, 4 ) = zero
506  z( 4, 4 ) = -e( jsp1, jsp1 )
507 *
508 * Set up right hand side(s)
509 *
510  rhs( 1 ) = c( is, js )
511  rhs( 2 ) = c( is, jsp1 )
512  rhs( 3 ) = f( is, js )
513  rhs( 4 ) = f( is, jsp1 )
514 *
515 * Solve Z * x = RHS
516 *
517  CALL sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
518  IF( ierr.GT.0 )
519  $ info = ierr
520 *
521  IF( ijob.EQ.0 ) THEN
522  CALL sgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
523  $ scaloc )
524  IF( scaloc.NE.one ) THEN
525  DO 60 k = 1, n
526  CALL sscal( m, scaloc, c( 1, k ), 1 )
527  CALL sscal( m, scaloc, f( 1, k ), 1 )
528  60 CONTINUE
529  scale = scale*scaloc
530  END IF
531  ELSE
532  CALL slatdf( ijob, zdim, z, ldz, rhs, rdsum,
533  $ rdscal, ipiv, jpiv )
534  END IF
535 *
536 * Unpack solution vector(s)
537 *
538  c( is, js ) = rhs( 1 )
539  c( is, jsp1 ) = rhs( 2 )
540  f( is, js ) = rhs( 3 )
541  f( is, jsp1 ) = rhs( 4 )
542 *
543 * Substitute R(I, J) and L(I, J) into remaining
544 * equation.
545 *
546  IF( i.GT.1 ) THEN
547  CALL sger( is-1, nb, -one, a( 1, is ), 1, rhs( 1 ),
548  $ 1, c( 1, js ), ldc )
549  CALL sger( is-1, nb, -one, d( 1, is ), 1, rhs( 1 ),
550  $ 1, f( 1, js ), ldf )
551  END IF
552  IF( j.LT.q ) THEN
553  CALL saxpy( n-je, rhs( 3 ), b( js, je+1 ), ldb,
554  $ c( is, je+1 ), ldc )
555  CALL saxpy( n-je, rhs( 3 ), e( js, je+1 ), lde,
556  $ f( is, je+1 ), ldf )
557  CALL saxpy( n-je, rhs( 4 ), b( jsp1, je+1 ), ldb,
558  $ c( is, je+1 ), ldc )
559  CALL saxpy( n-je, rhs( 4 ), e( jsp1, je+1 ), lde,
560  $ f( is, je+1 ), ldf )
561  END IF
562 *
563  ELSE IF( ( mb.EQ.2 ) .AND. ( nb.EQ.1 ) ) THEN
564 *
565 * Build a 4-by-4 system Z * x = RHS
566 *
567  z( 1, 1 ) = a( is, is )
568  z( 2, 1 ) = a( isp1, is )
569  z( 3, 1 ) = d( is, is )
570  z( 4, 1 ) = zero
571 *
572  z( 1, 2 ) = a( is, isp1 )
573  z( 2, 2 ) = a( isp1, isp1 )
574  z( 3, 2 ) = d( is, isp1 )
575  z( 4, 2 ) = d( isp1, isp1 )
576 *
577  z( 1, 3 ) = -b( js, js )
578  z( 2, 3 ) = zero
579  z( 3, 3 ) = -e( js, js )
580  z( 4, 3 ) = zero
581 *
582  z( 1, 4 ) = zero
583  z( 2, 4 ) = -b( js, js )
584  z( 3, 4 ) = zero
585  z( 4, 4 ) = -e( js, js )
586 *
587 * Set up right hand side(s)
588 *
589  rhs( 1 ) = c( is, js )
590  rhs( 2 ) = c( isp1, js )
591  rhs( 3 ) = f( is, js )
592  rhs( 4 ) = f( isp1, js )
593 *
594 * Solve Z * x = RHS
595 *
596  CALL sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
597  IF( ierr.GT.0 )
598  $ info = ierr
599  IF( ijob.EQ.0 ) THEN
600  CALL sgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
601  $ scaloc )
602  IF( scaloc.NE.one ) THEN
603  DO 70 k = 1, n
604  CALL sscal( m, scaloc, c( 1, k ), 1 )
605  CALL sscal( m, scaloc, f( 1, k ), 1 )
606  70 CONTINUE
607  scale = scale*scaloc
608  END IF
609  ELSE
610  CALL slatdf( ijob, zdim, z, ldz, rhs, rdsum,
611  $ rdscal, ipiv, jpiv )
612  END IF
613 *
614 * Unpack solution vector(s)
615 *
616  c( is, js ) = rhs( 1 )
617  c( isp1, js ) = rhs( 2 )
618  f( is, js ) = rhs( 3 )
619  f( isp1, js ) = rhs( 4 )
620 *
621 * Substitute R(I, J) and L(I, J) into remaining
622 * equation.
623 *
624  IF( i.GT.1 ) THEN
625  CALL sgemv( 'N', is-1, mb, -one, a( 1, is ), lda,
626  $ rhs( 1 ), 1, one, c( 1, js ), 1 )
627  CALL sgemv( 'N', is-1, mb, -one, d( 1, is ), ldd,
628  $ rhs( 1 ), 1, one, f( 1, js ), 1 )
629  END IF
630  IF( j.LT.q ) THEN
631  CALL sger( mb, n-je, one, rhs( 3 ), 1,
632  $ b( js, je+1 ), ldb, c( is, je+1 ), ldc )
633  CALL sger( mb, n-je, one, rhs( 3 ), 1,
634  $ e( js, je+1 ), lde, f( is, je+1 ), ldf )
635  END IF
636 *
637  ELSE IF( ( mb.EQ.2 ) .AND. ( nb.EQ.2 ) ) THEN
638 *
639 * Build an 8-by-8 system Z * x = RHS
640 *
641  CALL slaset( 'F', ldz, ldz, zero, zero, z, ldz )
642 *
643  z( 1, 1 ) = a( is, is )
644  z( 2, 1 ) = a( isp1, is )
645  z( 5, 1 ) = d( is, is )
646 *
647  z( 1, 2 ) = a( is, isp1 )
648  z( 2, 2 ) = a( isp1, isp1 )
649  z( 5, 2 ) = d( is, isp1 )
650  z( 6, 2 ) = d( isp1, isp1 )
651 *
652  z( 3, 3 ) = a( is, is )
653  z( 4, 3 ) = a( isp1, is )
654  z( 7, 3 ) = d( is, is )
655 *
656  z( 3, 4 ) = a( is, isp1 )
657  z( 4, 4 ) = a( isp1, isp1 )
658  z( 7, 4 ) = d( is, isp1 )
659  z( 8, 4 ) = d( isp1, isp1 )
660 *
661  z( 1, 5 ) = -b( js, js )
662  z( 3, 5 ) = -b( js, jsp1 )
663  z( 5, 5 ) = -e( js, js )
664  z( 7, 5 ) = -e( js, jsp1 )
665 *
666  z( 2, 6 ) = -b( js, js )
667  z( 4, 6 ) = -b( js, jsp1 )
668  z( 6, 6 ) = -e( js, js )
669  z( 8, 6 ) = -e( js, jsp1 )
670 *
671  z( 1, 7 ) = -b( jsp1, js )
672  z( 3, 7 ) = -b( jsp1, jsp1 )
673  z( 7, 7 ) = -e( jsp1, jsp1 )
674 *
675  z( 2, 8 ) = -b( jsp1, js )
676  z( 4, 8 ) = -b( jsp1, jsp1 )
677  z( 8, 8 ) = -e( jsp1, jsp1 )
678 *
679 * Set up right hand side(s)
680 *
681  k = 1
682  ii = mb*nb + 1
683  DO 80 jj = 0, nb - 1
684  CALL scopy( mb, c( is, js+jj ), 1, rhs( k ), 1 )
685  CALL scopy( mb, f( is, js+jj ), 1, rhs( ii ), 1 )
686  k = k + mb
687  ii = ii + mb
688  80 CONTINUE
689 *
690 * Solve Z * x = RHS
691 *
692  CALL sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
693  IF( ierr.GT.0 )
694  $ info = ierr
695  IF( ijob.EQ.0 ) THEN
696  CALL sgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
697  $ scaloc )
698  IF( scaloc.NE.one ) THEN
699  DO 90 k = 1, n
700  CALL sscal( m, scaloc, c( 1, k ), 1 )
701  CALL sscal( m, scaloc, f( 1, k ), 1 )
702  90 CONTINUE
703  scale = scale*scaloc
704  END IF
705  ELSE
706  CALL slatdf( ijob, zdim, z, ldz, rhs, rdsum,
707  $ rdscal, ipiv, jpiv )
708  END IF
709 *
710 * Unpack solution vector(s)
711 *
712  k = 1
713  ii = mb*nb + 1
714  DO 100 jj = 0, nb - 1
715  CALL scopy( mb, rhs( k ), 1, c( is, js+jj ), 1 )
716  CALL scopy( mb, rhs( ii ), 1, f( is, js+jj ), 1 )
717  k = k + mb
718  ii = ii + mb
719  100 CONTINUE
720 *
721 * Substitute R(I, J) and L(I, J) into remaining
722 * equation.
723 *
724  IF( i.GT.1 ) THEN
725  CALL sgemm( 'N', 'N', is-1, nb, mb, -one,
726  $ a( 1, is ), lda, rhs( 1 ), mb, one,
727  $ c( 1, js ), ldc )
728  CALL sgemm( 'N', 'N', is-1, nb, mb, -one,
729  $ d( 1, is ), ldd, rhs( 1 ), mb, one,
730  $ f( 1, js ), ldf )
731  END IF
732  IF( j.LT.q ) THEN
733  k = mb*nb + 1
734  CALL sgemm( 'N', 'N', mb, n-je, nb, one, rhs( k ),
735  $ mb, b( js, je+1 ), ldb, one,
736  $ c( is, je+1 ), ldc )
737  CALL sgemm( 'N', 'N', mb, n-je, nb, one, rhs( k ),
738  $ mb, e( js, je+1 ), lde, one,
739  $ f( is, je+1 ), ldf )
740  END IF
741 *
742  END IF
743 *
744  110 CONTINUE
745  120 CONTINUE
746  ELSE
747 *
748 * Solve (I, J) - subsystem
749 * A(I, I)**T * R(I, J) + D(I, I)**T * L(J, J) = C(I, J)
750 * R(I, I) * B(J, J) + L(I, J) * E(J, J) = -F(I, J)
751 * for I = 1, 2, ..., P, J = Q, Q - 1, ..., 1
752 *
753  scale = one
754  scaloc = one
755  DO 200 i = 1, p
756 *
757  is = iwork( i )
758  isp1 = is + 1
759  ie = iwork( i+1 ) - 1
760  mb = ie - is + 1
761  DO 190 j = q, p + 2, -1
762 *
763  js = iwork( j )
764  jsp1 = js + 1
765  je = iwork( j+1 ) - 1
766  nb = je - js + 1
767  zdim = mb*nb*2
768  IF( ( mb.EQ.1 ) .AND. ( nb.EQ.1 ) ) THEN
769 *
770 * Build a 2-by-2 system Z**T * x = RHS
771 *
772  z( 1, 1 ) = a( is, is )
773  z( 2, 1 ) = -b( js, js )
774  z( 1, 2 ) = d( is, is )
775  z( 2, 2 ) = -e( js, js )
776 *
777 * Set up right hand side(s)
778 *
779  rhs( 1 ) = c( is, js )
780  rhs( 2 ) = f( is, js )
781 *
782 * Solve Z**T * x = RHS
783 *
784  CALL sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
785  IF( ierr.GT.0 )
786  $ info = ierr
787 *
788  CALL sgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
789  IF( scaloc.NE.one ) THEN
790  DO 130 k = 1, n
791  CALL sscal( m, scaloc, c( 1, k ), 1 )
792  CALL sscal( m, scaloc, f( 1, k ), 1 )
793  130 CONTINUE
794  scale = scale*scaloc
795  END IF
796 *
797 * Unpack solution vector(s)
798 *
799  c( is, js ) = rhs( 1 )
800  f( is, js ) = rhs( 2 )
801 *
802 * Substitute R(I, J) and L(I, J) into remaining
803 * equation.
804 *
805  IF( j.GT.p+2 ) THEN
806  alpha = rhs( 1 )
807  CALL saxpy( js-1, alpha, b( 1, js ), 1, f( is, 1 ),
808  $ ldf )
809  alpha = rhs( 2 )
810  CALL saxpy( js-1, alpha, e( 1, js ), 1, f( is, 1 ),
811  $ ldf )
812  END IF
813  IF( i.LT.p ) THEN
814  alpha = -rhs( 1 )
815  CALL saxpy( m-ie, alpha, a( is, ie+1 ), lda,
816  $ c( ie+1, js ), 1 )
817  alpha = -rhs( 2 )
818  CALL saxpy( m-ie, alpha, d( is, ie+1 ), ldd,
819  $ c( ie+1, js ), 1 )
820  END IF
821 *
822  ELSE IF( ( mb.EQ.1 ) .AND. ( nb.EQ.2 ) ) THEN
823 *
824 * Build a 4-by-4 system Z**T * x = RHS
825 *
826  z( 1, 1 ) = a( is, is )
827  z( 2, 1 ) = zero
828  z( 3, 1 ) = -b( js, js )
829  z( 4, 1 ) = -b( jsp1, js )
830 *
831  z( 1, 2 ) = zero
832  z( 2, 2 ) = a( is, is )
833  z( 3, 2 ) = -b( js, jsp1 )
834  z( 4, 2 ) = -b( jsp1, jsp1 )
835 *
836  z( 1, 3 ) = d( is, is )
837  z( 2, 3 ) = zero
838  z( 3, 3 ) = -e( js, js )
839  z( 4, 3 ) = zero
840 *
841  z( 1, 4 ) = zero
842  z( 2, 4 ) = d( is, is )
843  z( 3, 4 ) = -e( js, jsp1 )
844  z( 4, 4 ) = -e( jsp1, jsp1 )
845 *
846 * Set up right hand side(s)
847 *
848  rhs( 1 ) = c( is, js )
849  rhs( 2 ) = c( is, jsp1 )
850  rhs( 3 ) = f( is, js )
851  rhs( 4 ) = f( is, jsp1 )
852 *
853 * Solve Z**T * x = RHS
854 *
855  CALL sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
856  IF( ierr.GT.0 )
857  $ info = ierr
858  CALL sgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
859  IF( scaloc.NE.one ) THEN
860  DO 140 k = 1, n
861  CALL sscal( m, scaloc, c( 1, k ), 1 )
862  CALL sscal( m, scaloc, f( 1, k ), 1 )
863  140 CONTINUE
864  scale = scale*scaloc
865  END IF
866 *
867 * Unpack solution vector(s)
868 *
869  c( is, js ) = rhs( 1 )
870  c( is, jsp1 ) = rhs( 2 )
871  f( is, js ) = rhs( 3 )
872  f( is, jsp1 ) = rhs( 4 )
873 *
874 * Substitute R(I, J) and L(I, J) into remaining
875 * equation.
876 *
877  IF( j.GT.p+2 ) THEN
878  CALL saxpy( js-1, rhs( 1 ), b( 1, js ), 1,
879  $ f( is, 1 ), ldf )
880  CALL saxpy( js-1, rhs( 2 ), b( 1, jsp1 ), 1,
881  $ f( is, 1 ), ldf )
882  CALL saxpy( js-1, rhs( 3 ), e( 1, js ), 1,
883  $ f( is, 1 ), ldf )
884  CALL saxpy( js-1, rhs( 4 ), e( 1, jsp1 ), 1,
885  $ f( is, 1 ), ldf )
886  END IF
887  IF( i.LT.p ) THEN
888  CALL sger( m-ie, nb, -one, a( is, ie+1 ), lda,
889  $ rhs( 1 ), 1, c( ie+1, js ), ldc )
890  CALL sger( m-ie, nb, -one, d( is, ie+1 ), ldd,
891  $ rhs( 3 ), 1, c( ie+1, js ), ldc )
892  END IF
893 *
894  ELSE IF( ( mb.EQ.2 ) .AND. ( nb.EQ.1 ) ) THEN
895 *
896 * Build a 4-by-4 system Z**T * x = RHS
897 *
898  z( 1, 1 ) = a( is, is )
899  z( 2, 1 ) = a( is, isp1 )
900  z( 3, 1 ) = -b( js, js )
901  z( 4, 1 ) = zero
902 *
903  z( 1, 2 ) = a( isp1, is )
904  z( 2, 2 ) = a( isp1, isp1 )
905  z( 3, 2 ) = zero
906  z( 4, 2 ) = -b( js, js )
907 *
908  z( 1, 3 ) = d( is, is )
909  z( 2, 3 ) = d( is, isp1 )
910  z( 3, 3 ) = -e( js, js )
911  z( 4, 3 ) = zero
912 *
913  z( 1, 4 ) = zero
914  z( 2, 4 ) = d( isp1, isp1 )
915  z( 3, 4 ) = zero
916  z( 4, 4 ) = -e( js, js )
917 *
918 * Set up right hand side(s)
919 *
920  rhs( 1 ) = c( is, js )
921  rhs( 2 ) = c( isp1, js )
922  rhs( 3 ) = f( is, js )
923  rhs( 4 ) = f( isp1, js )
924 *
925 * Solve Z**T * x = RHS
926 *
927  CALL sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
928  IF( ierr.GT.0 )
929  $ info = ierr
930 *
931  CALL sgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
932  IF( scaloc.NE.one ) THEN
933  DO 150 k = 1, n
934  CALL sscal( m, scaloc, c( 1, k ), 1 )
935  CALL sscal( m, scaloc, f( 1, k ), 1 )
936  150 CONTINUE
937  scale = scale*scaloc
938  END IF
939 *
940 * Unpack solution vector(s)
941 *
942  c( is, js ) = rhs( 1 )
943  c( isp1, js ) = rhs( 2 )
944  f( is, js ) = rhs( 3 )
945  f( isp1, js ) = rhs( 4 )
946 *
947 * Substitute R(I, J) and L(I, J) into remaining
948 * equation.
949 *
950  IF( j.GT.p+2 ) THEN
951  CALL sger( mb, js-1, one, rhs( 1 ), 1, b( 1, js ),
952  $ 1, f( is, 1 ), ldf )
953  CALL sger( mb, js-1, one, rhs( 3 ), 1, e( 1, js ),
954  $ 1, f( is, 1 ), ldf )
955  END IF
956  IF( i.LT.p ) THEN
957  CALL sgemv( 'T', mb, m-ie, -one, a( is, ie+1 ),
958  $ lda, rhs( 1 ), 1, one, c( ie+1, js ),
959  $ 1 )
960  CALL sgemv( 'T', mb, m-ie, -one, d( is, ie+1 ),
961  $ ldd, rhs( 3 ), 1, one, c( ie+1, js ),
962  $ 1 )
963  END IF
964 *
965  ELSE IF( ( mb.EQ.2 ) .AND. ( nb.EQ.2 ) ) THEN
966 *
967 * Build an 8-by-8 system Z**T * x = RHS
968 *
969  CALL slaset( 'F', ldz, ldz, zero, zero, z, ldz )
970 *
971  z( 1, 1 ) = a( is, is )
972  z( 2, 1 ) = a( is, isp1 )
973  z( 5, 1 ) = -b( js, js )
974  z( 7, 1 ) = -b( jsp1, js )
975 *
976  z( 1, 2 ) = a( isp1, is )
977  z( 2, 2 ) = a( isp1, isp1 )
978  z( 6, 2 ) = -b( js, js )
979  z( 8, 2 ) = -b( jsp1, js )
980 *
981  z( 3, 3 ) = a( is, is )
982  z( 4, 3 ) = a( is, isp1 )
983  z( 5, 3 ) = -b( js, jsp1 )
984  z( 7, 3 ) = -b( jsp1, jsp1 )
985 *
986  z( 3, 4 ) = a( isp1, is )
987  z( 4, 4 ) = a( isp1, isp1 )
988  z( 6, 4 ) = -b( js, jsp1 )
989  z( 8, 4 ) = -b( jsp1, jsp1 )
990 *
991  z( 1, 5 ) = d( is, is )
992  z( 2, 5 ) = d( is, isp1 )
993  z( 5, 5 ) = -e( js, js )
994 *
995  z( 2, 6 ) = d( isp1, isp1 )
996  z( 6, 6 ) = -e( js, js )
997 *
998  z( 3, 7 ) = d( is, is )
999  z( 4, 7 ) = d( is, isp1 )
1000  z( 5, 7 ) = -e( js, jsp1 )
1001  z( 7, 7 ) = -e( jsp1, jsp1 )
1002 *
1003  z( 4, 8 ) = d( isp1, isp1 )
1004  z( 6, 8 ) = -e( js, jsp1 )
1005  z( 8, 8 ) = -e( jsp1, jsp1 )
1006 *
1007 * Set up right hand side(s)
1008 *
1009  k = 1
1010  ii = mb*nb + 1
1011  DO 160 jj = 0, nb - 1
1012  CALL scopy( mb, c( is, js+jj ), 1, rhs( k ), 1 )
1013  CALL scopy( mb, f( is, js+jj ), 1, rhs( ii ), 1 )
1014  k = k + mb
1015  ii = ii + mb
1016  160 CONTINUE
1017 *
1018 *
1019 * Solve Z**T * x = RHS
1020 *
1021  CALL sgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
1022  IF( ierr.GT.0 )
1023  $ info = ierr
1024 *
1025  CALL sgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
1026  IF( scaloc.NE.one ) THEN
1027  DO 170 k = 1, n
1028  CALL sscal( m, scaloc, c( 1, k ), 1 )
1029  CALL sscal( m, scaloc, f( 1, k ), 1 )
1030  170 CONTINUE
1031  scale = scale*scaloc
1032  END IF
1033 *
1034 * Unpack solution vector(s)
1035 *
1036  k = 1
1037  ii = mb*nb + 1
1038  DO 180 jj = 0, nb - 1
1039  CALL scopy( mb, rhs( k ), 1, c( is, js+jj ), 1 )
1040  CALL scopy( mb, rhs( ii ), 1, f( is, js+jj ), 1 )
1041  k = k + mb
1042  ii = ii + mb
1043  180 CONTINUE
1044 *
1045 * Substitute R(I, J) and L(I, J) into remaining
1046 * equation.
1047 *
1048  IF( j.GT.p+2 ) THEN
1049  CALL sgemm( 'N', 'T', mb, js-1, nb, one,
1050  $ c( is, js ), ldc, b( 1, js ), ldb, one,
1051  $ f( is, 1 ), ldf )
1052  CALL sgemm( 'N', 'T', mb, js-1, nb, one,
1053  $ f( is, js ), ldf, e( 1, js ), lde, one,
1054  $ f( is, 1 ), ldf )
1055  END IF
1056  IF( i.LT.p ) THEN
1057  CALL sgemm( 'T', 'N', m-ie, nb, mb, -one,
1058  $ a( is, ie+1 ), lda, c( is, js ), ldc,
1059  $ one, c( ie+1, js ), ldc )
1060  CALL sgemm( 'T', 'N', m-ie, nb, mb, -one,
1061  $ d( is, ie+1 ), ldd, f( is, js ), ldf,
1062  $ one, c( ie+1, js ), ldc )
1063  END IF
1064 *
1065  END IF
1066 *
1067  190 CONTINUE
1068  200 CONTINUE
1069 *
1070  END IF
1071  RETURN
1072 *
1073 * End of STGSY2
1074 *
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:158
subroutine sger(M, N, ALPHA, X, INCX, Y, INCY, A, LDA)
SGER
Definition: sger.f:132
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine sgesc2(N, A, LDA, RHS, IPIV, JPIV, SCALE)
SGESC2 solves a system of linear equations using the LU factorization with complete pivoting computed...
Definition: sgesc2.f:116
subroutine sgetc2(N, A, LDA, IPIV, JPIV, INFO)
SGETC2 computes the LU factorization with complete pivoting of the general n-by-n matrix...
Definition: sgetc2.f:113
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
Definition: saxpy.f:54
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:189
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:55
logical function lde(RI, RJ, LR)
Definition: dblat2.f:2945
subroutine slatdf(IJOB, N, Z, LDZ, RHS, RDSUM, RDSCAL, IPIV, JPIV)
SLATDF uses the LU factorization of the n-by-n matrix computed by sgetc2 and computes a contribution ...
Definition: slatdf.f:173
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: