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

Functions

double precision function dlansy (norm, uplo, n, a, lda, work)
 DLANSY 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 dlaqsy (uplo, n, a, lda, s, scond, amax, equed)
 DLAQSY scales a symmetric/Hermitian matrix, using scaling factors computed by spoequ.
subroutine dlasy2 (ltranl, ltranr, isgn, n1, n2, tl, ldtl, tr, ldtr, b, ldb, scale, x, ldx, xnorm, info)
 DLASY2 solves the Sylvester matrix equation where the matrices are of order 1 or 2.
subroutine dsyswapr (uplo, n, a, lda, i1, i2)
 DSYSWAPR applies an elementary permutation on the rows and columns of a symmetric matrix.
subroutine dtgsy2 (trans, ijob, m, n, a, lda, b, ldb, c, ldc, d, ldd, e, lde, f, ldf, scale, rdsum, rdscal, iwork, pq, info)
 DTGSY2 solves the generalized Sylvester equation (unblocked algorithm).

Detailed Description

This is the group of double auxiliary functions for SY matrices

Function Documentation

◆ dlansy()

double precision function dlansy ( character norm,
character uplo,
integer n,
double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( * ) work )

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

Purpose:
!>
!> DLANSY  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
DLANSY 
!>
!>    DLANSY = ( 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 DLANSY 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, DLANSY is
!>          set to zero.
!>
[in]A
!>          A is DOUBLE PRECISION 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 DOUBLE PRECISION 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 dlansy.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 DOUBLE PRECISION A( LDA, * ), WORK( * )
133* ..
134*
135* =====================================================================
136*
137* .. Parameters ..
138 DOUBLE PRECISION ONE, ZERO
139 parameter( one = 1.0d+0, zero = 0.0d+0 )
140* ..
141* .. Local Scalars ..
142 INTEGER I, J
143 DOUBLE PRECISION ABSA, SCALE, SUM, VALUE
144* ..
145* .. External Subroutines ..
146 EXTERNAL dlassq
147* ..
148* .. External Functions ..
149 LOGICAL LSAME, DISNAN
150 EXTERNAL lsame, disnan
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. disnan( 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. disnan( 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. disnan( 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. disnan( 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 dlassq( j-1, a( 1, j ), 1, scale, sum )
222 110 CONTINUE
223 ELSE
224 DO 120 j = 1, n - 1
225 CALL dlassq( n-j, a( j+1, j ), 1, scale, sum )
226 120 CONTINUE
227 END IF
228 sum = 2*sum
229 CALL dlassq( n, a, lda+1, scale, sum )
230 VALUE = scale*sqrt( sum )
231 END IF
232*
233 dlansy = VALUE
234 RETURN
235*
236* End of DLANSY
237*
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
dlassq
subroutine dlassq(n, x, incx, scl, sumsq)
DLASSQ updates a sum of squares represented in scaled form.
Definition dlassq.f90:137
disnan
logical function disnan(din)
DISNAN tests input for NaN.
Definition disnan.f:59
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
dlansy
double precision function dlansy(norm, uplo, n, a, lda, work)
DLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition dlansy.f:122

◆ dlaqsy()

subroutine dlaqsy ( character uplo,
integer n,
double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( * ) s,
double precision scond,
double precision amax,
character equed )

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

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

Purpose:
!>
!> DLAQSY 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N)
!>          The scale factors for A.
!>
[in]SCOND
!>          SCOND is DOUBLE PRECISION
!>          Ratio of the smallest S(i) to the largest S(i).
!>
[in]AMAX
!>          AMAX is DOUBLE PRECISION
!>          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 dlaqsy.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 DOUBLE PRECISION AMAX, SCOND
142* ..
143* .. Array Arguments ..
144 DOUBLE PRECISION A( LDA, * ), S( * )
145* ..
146*
147* =====================================================================
148*
149* .. Parameters ..
150 DOUBLE PRECISION ONE, THRESH
151 parameter( one = 1.0d+0, thresh = 0.1d+0 )
152* ..
153* .. Local Scalars ..
154 INTEGER I, J
155 DOUBLE PRECISION CJ, LARGE, SMALL
156* ..
157* .. External Functions ..
158 LOGICAL LSAME
159 DOUBLE PRECISION DLAMCH
160 EXTERNAL lsame, dlamch
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 = dlamch( 'Safe minimum' ) / dlamch( '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 DLAQSY
212*
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69

◆ dlasy2()

subroutine dlasy2 ( logical ltranl,
logical ltranr,
integer isgn,
integer n1,
integer n2,
double precision, dimension( ldtl, * ) tl,
integer ldtl,
double precision, dimension( ldtr, * ) tr,
integer ldtr,
double precision, dimension( ldb, * ) b,
integer ldb,
double precision scale,
double precision, dimension( ldx, * ) x,
integer ldx,
double precision xnorm,
integer info )

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

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

Purpose:
!>
!> DLASY2 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
!>          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 DOUBLE PRECISION 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 DOUBLE PRECISION
!>          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 dlasy2.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 DOUBLE PRECISION SCALE, XNORM
183* ..
184* .. Array Arguments ..
185 DOUBLE PRECISION B( LDB, * ), TL( LDTL, * ), TR( LDTR, * ),
186 $ X( LDX, * )
187* ..
188*
189* =====================================================================
190*
191* .. Parameters ..
192 DOUBLE PRECISION ZERO, ONE
193 parameter( zero = 0.0d+0, one = 1.0d+0 )
194 DOUBLE PRECISION TWO, HALF, EIGHT
195 parameter( two = 2.0d+0, half = 0.5d+0, eight = 8.0d+0 )
196* ..
197* .. Local Scalars ..
198 LOGICAL BSWAP, XSWAP
199 INTEGER I, IP, IPIV, IPSV, J, JP, JPSV, K
200 DOUBLE PRECISION 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 DOUBLE PRECISION BTMP( 4 ), T16( 4, 4 ), TMP( 4 ), X2( 2 )
208* ..
209* .. External Functions ..
210 INTEGER IDAMAX
211 DOUBLE PRECISION DLAMCH
212 EXTERNAL idamax, dlamch
213* ..
214* .. External Subroutines ..
215 EXTERNAL dcopy, dswap
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 = dlamch( 'P' )
240 smlnum = dlamch( '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 = idamax( 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 dcopy( 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 dswap( 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 dswap( 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 DLASY2
478*
idamax
integer function idamax(n, dx, incx)
IDAMAX
Definition idamax.f:71
dswap
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
max
#define max(a, b)
Definition macros.h:21

◆ dsyswapr()

subroutine dsyswapr ( character uplo,
integer n,
double precision, dimension( lda, n ) a,
integer lda,
integer i1,
integer i2 )

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

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

Purpose:
!>
!> DSYSWAPR 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 DOUBLE PRECISION 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 DSYTRF.
!>
!>          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 dsyswapr.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 DOUBLE PRECISION A( LDA, N )
113*
114* =====================================================================
115*
116* ..
117* .. Local Scalars ..
118 LOGICAL UPPER
119 INTEGER I
120 DOUBLE PRECISION TMP
121*
122* .. External Functions ..
123 LOGICAL LSAME
124 EXTERNAL lsame
125* ..
126* .. External Subroutines ..
127 EXTERNAL dswap
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 dswap( 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 dswap( 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

◆ dtgsy2()

subroutine dtgsy2 ( character trans,
integer ijob,
integer m,
integer n,
double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( ldb, * ) b,
integer ldb,
double precision, dimension( ldc, * ) c,
integer ldc,
double precision, dimension( ldd, * ) d,
integer ldd,
double precision, dimension( lde, * ) e,
integer lde,
double precision, dimension( ldf, * ) f,
integer ldf,
double precision scale,
double precision rdsum,
double precision rdscal,
integer, dimension( * ) iwork,
integer pq,
integer info )

DTGSY2 solves the generalized Sylvester equation (unblocked algorithm).

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

Purpose:
!>
!> DTGSY2 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 DLACON.
!>
!> DTGSY2 also (IJOB >= 1) contributes to the computation in DTGSYL
!> 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
!> DTGSYL. See DTGSYL 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. (DGECON 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION
!>          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 DOUBLE PRECISION
!>          On entry, the sum of squares of computed contributions to
!>          the Dif-estimate under computation by DTGSYL, 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 DTGSY2 is called by DTGSYL.
!>
[in,out]RDSCAL
!>          RDSCAL is DOUBLE PRECISION
!>          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 DTGSY2 is called by
!>                DTGSYL.
!>
[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 dtgsy2.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 DOUBLE PRECISION RDSCAL, RDSUM, SCALE
284* ..
285* .. Array Arguments ..
286 INTEGER IWORK( * )
287 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ),
288 $ D( LDD, * ), E( LDE, * ), F( LDF, * )
289* ..
290*
291* =====================================================================
292* Replaced various illegal calls to DCOPY by calls to DLASET.
293* Sven Hammarling, 27/5/02.
294*
295* .. Parameters ..
296 INTEGER LDZ
297 parameter( ldz = 8 )
298 DOUBLE PRECISION ZERO, ONE
299 parameter( zero = 0.0d+0, one = 1.0d+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 DOUBLE PRECISION ALPHA, SCALOC
306* ..
307* .. Local Arrays ..
308 INTEGER IPIV( LDZ ), JPIV( LDZ )
309 DOUBLE PRECISION RHS( LDZ ), Z( LDZ, LDZ )
310* ..
311* .. External Functions ..
312 LOGICAL LSAME
313 EXTERNAL lsame
314* ..
315* .. External Subroutines ..
316 EXTERNAL daxpy, dcopy, dgemm, dgemv, dger, dgesc2,
317 $ dgetc2, dlaset, dlatdf, dscal, 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( 'DTGSY2', -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 dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
441 IF( ierr.GT.0 )
442 $ info = ierr
443*
444 IF( ijob.EQ.0 ) THEN
445 CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
446 $ scaloc )
447 IF( scaloc.NE.one ) THEN
448 DO 50 k = 1, n
449 CALL dscal( m, scaloc, c( 1, k ), 1 )
450 CALL dscal( m, scaloc, f( 1, k ), 1 )
451 50 CONTINUE
452 scale = scale*scaloc
453 END IF
454 ELSE
455 CALL dlatdf( 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 daxpy( is-1, alpha, a( 1, is ), 1, c( 1, js ),
470 $ 1 )
471 CALL daxpy( is-1, alpha, d( 1, is ), 1, f( 1, js ),
472 $ 1 )
473 END IF
474 IF( j.LT.q ) THEN
475 CALL daxpy( n-je, rhs( 2 ), b( js, je+1 ), ldb,
476 $ c( is, je+1 ), ldc )
477 CALL daxpy( 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 dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
515 IF( ierr.GT.0 )
516 $ info = ierr
517*
518 IF( ijob.EQ.0 ) THEN
519 CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
520 $ scaloc )
521 IF( scaloc.NE.one ) THEN
522 DO 60 k = 1, n
523 CALL dscal( m, scaloc, c( 1, k ), 1 )
524 CALL dscal( m, scaloc, f( 1, k ), 1 )
525 60 CONTINUE
526 scale = scale*scaloc
527 END IF
528 ELSE
529 CALL dlatdf( 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 dger( is-1, nb, -one, a( 1, is ), 1, rhs( 1 ),
545 $ 1, c( 1, js ), ldc )
546 CALL dger( 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 daxpy( n-je, rhs( 3 ), b( js, je+1 ), ldb,
551 $ c( is, je+1 ), ldc )
552 CALL daxpy( n-je, rhs( 3 ), e( js, je+1 ), lde,
553 $ f( is, je+1 ), ldf )
554 CALL daxpy( n-je, rhs( 4 ), b( jsp1, je+1 ), ldb,
555 $ c( is, je+1 ), ldc )
556 CALL daxpy( 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 dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
594 IF( ierr.GT.0 )
595 $ info = ierr
596 IF( ijob.EQ.0 ) THEN
597 CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
598 $ scaloc )
599 IF( scaloc.NE.one ) THEN
600 DO 70 k = 1, n
601 CALL dscal( m, scaloc, c( 1, k ), 1 )
602 CALL dscal( m, scaloc, f( 1, k ), 1 )
603 70 CONTINUE
604 scale = scale*scaloc
605 END IF
606 ELSE
607 CALL dlatdf( 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 dgemv( 'N', is-1, mb, -one, a( 1, is ), lda,
623 $ rhs( 1 ), 1, one, c( 1, js ), 1 )
624 CALL dgemv( '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 dger( mb, n-je, one, rhs( 3 ), 1,
629 $ b( js, je+1 ), ldb, c( is, je+1 ), ldc )
630 CALL dger( 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 dlaset( '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 dcopy( mb, c( is, js+jj ), 1, rhs( k ), 1 )
682 CALL dcopy( 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 dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
690 IF( ierr.GT.0 )
691 $ info = ierr
692 IF( ijob.EQ.0 ) THEN
693 CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv,
694 $ scaloc )
695 IF( scaloc.NE.one ) THEN
696 DO 90 k = 1, n
697 CALL dscal( m, scaloc, c( 1, k ), 1 )
698 CALL dscal( m, scaloc, f( 1, k ), 1 )
699 90 CONTINUE
700 scale = scale*scaloc
701 END IF
702 ELSE
703 CALL dlatdf( 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 dcopy( mb, rhs( k ), 1, c( is, js+jj ), 1 )
713 CALL dcopy( 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 dgemm( 'N', 'N', is-1, nb, mb, -one,
723 $ a( 1, is ), lda, rhs( 1 ), mb, one,
724 $ c( 1, js ), ldc )
725 CALL dgemm( '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 dgemm( '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 dgemm( '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 dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
782 IF( ierr.GT.0 )
783 $ info = ierr
784*
785 CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
786 IF( scaloc.NE.one ) THEN
787 DO 130 k = 1, n
788 CALL dscal( m, scaloc, c( 1, k ), 1 )
789 CALL dscal( 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 daxpy( js-1, alpha, b( 1, js ), 1, f( is, 1 ),
805 $ ldf )
806 alpha = rhs( 2 )
807 CALL daxpy( 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 daxpy( m-ie, alpha, a( is, ie+1 ), lda,
813 $ c( ie+1, js ), 1 )
814 alpha = -rhs( 2 )
815 CALL daxpy( 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 dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
853 IF( ierr.GT.0 )
854 $ info = ierr
855 CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
856 IF( scaloc.NE.one ) THEN
857 DO 140 k = 1, n
858 CALL dscal( m, scaloc, c( 1, k ), 1 )
859 CALL dscal( 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 daxpy( js-1, rhs( 1 ), b( 1, js ), 1,
876 $ f( is, 1 ), ldf )
877 CALL daxpy( js-1, rhs( 2 ), b( 1, jsp1 ), 1,
878 $ f( is, 1 ), ldf )
879 CALL daxpy( js-1, rhs( 3 ), e( 1, js ), 1,
880 $ f( is, 1 ), ldf )
881 CALL daxpy( js-1, rhs( 4 ), e( 1, jsp1 ), 1,
882 $ f( is, 1 ), ldf )
883 END IF
884 IF( i.LT.p ) THEN
885 CALL dger( m-ie, nb, -one, a( is, ie+1 ), lda,
886 $ rhs( 1 ), 1, c( ie+1, js ), ldc )
887 CALL dger( 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 dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
925 IF( ierr.GT.0 )
926 $ info = ierr
927*
928 CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
929 IF( scaloc.NE.one ) THEN
930 DO 150 k = 1, n
931 CALL dscal( m, scaloc, c( 1, k ), 1 )
932 CALL dscal( 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 dger( mb, js-1, one, rhs( 1 ), 1, b( 1, js ),
949 $ 1, f( is, 1 ), ldf )
950 CALL dger( 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 dgemv( 'T', mb, m-ie, -one, a( is, ie+1 ),
955 $ lda, rhs( 1 ), 1, one, c( ie+1, js ),
956 $ 1 )
957 CALL dgemv( '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 dlaset( '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 dcopy( mb, c( is, js+jj ), 1, rhs( k ), 1 )
1010 CALL dcopy( 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 dgetc2( zdim, z, ldz, ipiv, jpiv, ierr )
1019 IF( ierr.GT.0 )
1020 $ info = ierr
1021*
1022 CALL dgesc2( zdim, z, ldz, rhs, ipiv, jpiv, scaloc )
1023 IF( scaloc.NE.one ) THEN
1024 DO 170 k = 1, n
1025 CALL dscal( m, scaloc, c( 1, k ), 1 )
1026 CALL dscal( 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 dcopy( mb, rhs( k ), 1, c( is, js+jj ), 1 )
1037 CALL dcopy( 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 dgemm( 'N', 'T', mb, js-1, nb, one,
1047 $ c( is, js ), ldc, b( 1, js ), ldb, one,
1048 $ f( is, 1 ), ldf )
1049 CALL dgemm( '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 dgemm( '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 dgemm( '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 DTGSY2
1071*
lde
logical function lde(ri, rj, lr)
Definition dblat2.f:2942
alpha
#define alpha
Definition eval.h:35
dlaset
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:110
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
dgetc2
subroutine dgetc2(n, a, lda, ipiv, jpiv, info)
DGETC2 computes the LU factorization with complete pivoting of the general n-by-n matrix.
Definition dgetc2.f:111
dgesc2
subroutine dgesc2(n, a, lda, rhs, ipiv, jpiv, scale)
DGESC2 solves a system of linear equations using the LU factorization with complete pivoting computed...
Definition dgesc2.f:114
dlatdf
subroutine dlatdf(ijob, n, z, ldz, rhs, rdsum, rdscal, ipiv, jpiv)
DLATDF uses the LU factorization of the n-by-n matrix computed by sgetc2 and computes a contribution ...
Definition dlatdf.f:171
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
daxpy
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
dgemv
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:156
dger
subroutine dger(m, n, alpha, x, incx, y, incy, a, lda)
DGER
Definition dger.f:130
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:187