OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
real
LAPACK » Symmetric Matrix » Auxiliary routines

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.
subroutine slaqsy (uplo, n, a, lda, s, scond, amax, equed)
 SLAQSY scales a symmetric/Hermitian matrix, using scaling factors computed by spoequ.
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.
subroutine ssyswapr (uplo, n, a, lda, i1, i2)
 SSYSWAPR applies an elementary permutation on the rows and columns of a symmetric matrix.
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).

Detailed Description

This is the group of real auxiliary functions for SY matrices

Function Documentation

◆ slansy()

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.

Definition at line 121 of file slansy.f.

122*
123* -- LAPACK auxiliary routine --
124* -- LAPACK is a software package provided by Univ. of Tennessee, --
125* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
126*
127* .. Scalar Arguments ..
128 CHARACTER NORM, UPLO
129 INTEGER LDA, N
130* ..
131* .. Array Arguments ..
132 REAL A( LDA, * ), WORK( * )
133* ..
134*
135* =====================================================================
136*
137* .. Parameters ..
138 REAL ONE, ZERO
139 parameter( one = 1.0e+0, zero = 0.0e+0 )
140* ..
141* .. Local Scalars ..
142 INTEGER I, J
143 REAL ABSA, SCALE, SUM, VALUE
144* ..
145* .. External Subroutines ..
146 EXTERNAL slassq
147* ..
148* .. External Functions ..
149 LOGICAL LSAME, SISNAN
150 EXTERNAL lsame, sisnan
151* ..
152* .. Intrinsic Functions ..
153 INTRINSIC abs, sqrt
154* ..
155* .. Executable Statements ..
156*
157 IF( n.EQ.0 ) THEN
158 VALUE = zero
159 ELSE IF( lsame( norm, 'M' ) ) THEN
160*
161* Find max(abs(A(i,j))).
162*
163 VALUE = zero
164 IF( lsame( uplo, 'U' ) ) THEN
165 DO 20 j = 1, n
166 DO 10 i = 1, j
167 sum = abs( a( i, j ) )
168 IF( VALUE .LT. sum .OR. sisnan( sum ) ) VALUE = sum
169 10 CONTINUE
170 20 CONTINUE
171 ELSE
172 DO 40 j = 1, n
173 DO 30 i = j, n
174 sum = abs( a( i, j ) )
175 IF( VALUE .LT. sum .OR. sisnan( sum ) ) VALUE = sum
176 30 CONTINUE
177 40 CONTINUE
178 END IF
179 ELSE IF( ( lsame( norm, 'I' ) ) .OR. ( lsame( norm, 'O' ) ) .OR.
180 $ ( norm.EQ.'1' ) ) THEN
181*
182* Find normI(A) ( = norm1(A), since A is symmetric).
183*
184 VALUE = zero
185 IF( lsame( uplo, 'U' ) ) THEN
186 DO 60 j = 1, n
187 sum = zero
188 DO 50 i = 1, j - 1
189 absa = abs( a( i, j ) )
190 sum = sum + absa
191 work( i ) = work( i ) + absa
192 50 CONTINUE
193 work( j ) = sum + abs( a( j, j ) )
194 60 CONTINUE
195 DO 70 i = 1, n
196 sum = work( i )
197 IF( VALUE .LT. sum .OR. sisnan( sum ) ) VALUE = sum
198 70 CONTINUE
199 ELSE
200 DO 80 i = 1, n
201 work( i ) = zero
202 80 CONTINUE
203 DO 100 j = 1, n
204 sum = work( j ) + abs( a( j, j ) )
205 DO 90 i = j + 1, n
206 absa = abs( a( i, j ) )
207 sum = sum + absa
208 work( i ) = work( i ) + absa
209 90 CONTINUE
210 IF( VALUE .LT. sum .OR. sisnan( sum ) ) VALUE = sum
211 100 CONTINUE
212 END IF
213 ELSE IF( ( lsame( norm, 'F' ) ) .OR. ( lsame( norm, 'E' ) ) ) THEN
214*
215* Find normF(A).
216*
217 scale = zero
218 sum = one
219 IF( lsame( uplo, 'U' ) ) THEN
220 DO 110 j = 2, n
221 CALL slassq( j-1, a( 1, j ), 1, scale, sum )
222 110 CONTINUE
223 ELSE
224 DO 120 j = 1, n - 1
225 CALL slassq( n-j, a( j+1, j ), 1, scale, sum )
226 120 CONTINUE
227 END IF
228 sum = 2*sum
229 CALL slassq( n, a, lda+1, scale, sum )
230 VALUE = scale*sqrt( sum )
231 END IF
232*
233 slansy = VALUE
234 RETURN
235*
236* End of SLANSY
237*
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
slassq
subroutine slassq(n, x, incx, scl, sumsq)
SLASSQ updates a sum of squares represented in scaled form.
Definition slassq.f90:137
sisnan
logical function sisnan(sin)
SISNAN tests input for NaN.
Definition sisnan.f:59
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
slansy
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,...
Definition slansy.f:122

◆ slaqsy()

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.

Definition at line 132 of file slaqsy.f.

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

◆ slasy2()

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.

Definition at line 172 of file slasy2.f.

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

◆ ssyswapr()

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.

Definition at line 101 of file ssyswapr.f.

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

◆ stgsy2()

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

Definition at line 271 of file stgsy2.f.

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