OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
bdlaexc.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine bdlaexc (n, t, ldt, j1, n1, n2, itraf, dtraf, work, info)

Function/Subroutine Documentation

◆ bdlaexc()

subroutine bdlaexc ( integer n,
double precision, dimension( ldt, * ) t,
integer ldt,
integer j1,
integer n1,
integer n2,
integer, dimension( * ) itraf,
double precision, dimension( * ) dtraf,
double precision, dimension( * ) work,
integer info )

Definition at line 1 of file bdlaexc.f.

3 IMPLICIT NONE
4*
5* -- ScaLAPACK auxiliary routine (version 2.0.2) --
6* Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver
7* May 1 2012
8*
9* .. Scalar Arguments ..
10 INTEGER INFO, J1, LDT, N, N1, N2
11* ..
12* .. Array Arguments ..
13 INTEGER ITRAF( * )
14 DOUBLE PRECISION DTRAF( * ), T( LDT, * ), WORK( * )
15* ..
16*
17* Purpose
18* =======
19*
20* BDLAEXC swaps adjacent diagonal blocks T11 and T22 of order 1 or 2 in
21* an upper quasi-triangular matrix T by an orthogonal similarity
22* transformation.
23*
24* In contrast to the LAPACK routine DLAEXC, the orthogonal
25* transformation matrix Q is not explicitly constructed but
26* represented by paramaters contained in the arrays ITRAF and DTRAF,
27* see the description of BDTREXC for more details.
28*
29* T must be in Schur canonical form, that is, block upper triangular
30* with 1-by-1 and 2-by-2 diagonal blocks; each 2-by-2 diagonal block
31* has its diagonal elemnts equal and its off-diagonal elements of
32* opposite sign.
33*
34* Arguments
35* =========
36*
37* N (input) INTEGER
38* The order of the matrix T. N >= 0.
39*
40* T (input/output) DOUBLE PRECISION array, dimension (LDT,N)
41* On entry, the upper quasi-triangular matrix T, in Schur
42* canonical form.
43* On exit, the updated matrix T, again in Schur canonical form.
44*
45* LDT (input) INTEGER
46* The leading dimension of the array T. LDT >= max(1,N).
47*
48* J1 (input) INTEGER
49* The index of the first row of the first block T11.
50*
51* N1 (input) INTEGER
52* The order of the first block T11. N1 = 0, 1 or 2.
53*
54* N2 (input) INTEGER
55* The order of the second block T22. N2 = 0, 1 or 2.
56*
57* ITRAF (output) INTEGER array, length k, where
58* k = 1, if N1+N2 = 2;
59* k = 2, if N1+N2 = 3;
60* k = 4, if N1+N2 = 4.
61* List of parameters for representing the transformation
62* matrix Q, see BDTREXC.
63*
64* DTRAF (output) DOUBLE PRECISION array, length k, where
65* k = 2, if N1+N2 = 2;
66* k = 5, if N1+N2 = 3;
67* k = 10, if N1+N2 = 4.
68* List of parameters for representing the transformation
69* matrix Q, see BDTREXC.
70*
71* WORK (workspace) DOUBLE PRECISION array, dimension (N)
72*
73* INFO (output) INTEGER
74* = 0: successful exit
75* = 1: the transformed matrix T would be too far from Schur
76* form; the blocks are not swapped and T and Q are
77* unchanged.
78*
79* =====================================================================
80*
81* .. Parameters ..
82 DOUBLE PRECISION ZERO, ONE
83 parameter( zero = 0.0d+0, one = 1.0d+0 )
84 DOUBLE PRECISION TEN
85 parameter( ten = 1.0d+1 )
86 INTEGER LDD, LDX
87 parameter( ldd = 4, ldx = 2 )
88* ..
89* .. Local Scalars ..
90 INTEGER IERR, J2, J3, J4, K, LD, LI, ND
91 DOUBLE PRECISION CS, DNORM, EPS, SCALE, SMLNUM, SN, T11, T22,
92 $ T33, TAU, TAU1, TAU2, TEMP, THRESH, WI1, WI2,
93 $ WR1, WR2, XNORM
94* ..
95* .. Local Arrays ..
96 DOUBLE PRECISION D( LDD, 4 ), X( LDX, 2 )
97* ..
98* .. External Functions ..
99 DOUBLE PRECISION DLAMCH, DLANGE
100 EXTERNAL dlamch, dlange
101* ..
102* .. External Subroutines ..
103 EXTERNAL dlamov, dlanv2, dlarfg, dlarfx, dlartg, dlasy2,
104 $ drot
105* ..
106* .. Intrinsic Functions ..
107 INTRINSIC abs, max
108* ..
109* .. Executable Statements ..
110*
111 info = 0
112*
113* Quick return if possible
114*
115 IF( n.EQ.0 .OR. n1.EQ.0 .OR. n2.EQ.0 )
116 $ RETURN
117 IF( j1+n1.GT.n )
118 $ RETURN
119*
120 j2 = j1 + 1
121 j3 = j1 + 2
122 j4 = j1 + 3
123*
124 IF( n1.EQ.1 .AND. n2.EQ.1 ) THEN
125*
126* Swap two 1-by-1 blocks.
127*
128 t11 = t( j1, j1 )
129 t22 = t( j2, j2 )
130*
131* Determine the transformation to perform the interchange.
132*
133 CALL dlartg( t( j1, j2 ), t22-t11, cs, sn, temp )
134*
135* Apply transformation to the matrix T.
136*
137 IF( j3.LE.n )
138 $ CALL drot( n-j1-1, t( j1, j3 ), ldt, t( j2, j3 ), ldt, cs,
139 $ sn )
140 CALL drot( j1-1, t( 1, j1 ), 1, t( 1, j2 ), 1, cs, sn )
141*
142 t( j1, j1 ) = t22
143 t( j2, j2 ) = t11
144*
145 itraf( 1 ) = j1
146 dtraf( 1 ) = cs
147 dtraf( 2 ) = sn
148*
149 ELSE
150*
151* Swapping involves at least one 2-by-2 block.
152*
153* Copy the diagonal block of order N1+N2 to the local array D
154* and compute its norm.
155*
156 nd = n1 + n2
157 CALL dlamov( 'Full', nd, nd, t( j1, j1 ), ldt, d, ldd )
158 dnorm = dlange( 'Max', nd, nd, d, ldd, work )
159*
160* Compute machine-dependent threshold for test for accepting
161* swap.
162*
163 eps = dlamch( 'P' )
164 smlnum = dlamch( 'S' ) / eps
165 thresh = max( ten*eps*dnorm, smlnum )
166*
167* Solve T11*X - X*T22 = scale*T12 for X.
168*
169 CALL dlasy2( .false., .false., -1, n1, n2, d, ldd,
170 $ d( n1+1, n1+1 ), ldd, d( 1, n1+1 ), ldd, scale, x,
171 $ ldx, xnorm, ierr )
172*
173* Swap the adjacent diagonal blocks.
174*
175 k = n1 + n1 + n2 - 3
176 GO TO ( 10, 20, 30 )k
177*
178 10 CONTINUE
179*
180* N1 = 1, N2 = 2: generate elementary reflector H so that:
181*
182* ( scale, X11, X12 ) H = ( 0, 0, * )
183*
184 dtraf( 1 ) = scale
185 dtraf( 2 ) = x( 1, 1 )
186 dtraf( 3 ) = x( 1, 2 )
187 CALL dlarfg( 3, dtraf( 3 ), dtraf, 1, tau )
188 dtraf( 3 ) = one
189 t11 = t( j1, j1 )
190*
191* Perform swap provisionally on diagonal block in D.
192*
193 CALL dlarfx( 'Left', 3, 3, dtraf, tau, d, ldd, work )
194 CALL dlarfx( 'Right', 3, 3, dtraf, tau, d, ldd, work )
195*
196* Test whether to reject swap.
197*
198 IF( max( abs( d( 3, 1 ) ), abs( d( 3, 2 ) ), abs( d( 3,
199 $ 3 )-t11 ) ).GT.thresh )GO TO 50
200*
201* Accept swap: apply transformation to the entire matrix T.
202*
203 CALL dlarfx( 'Left', 3, n-j1+1, dtraf, tau, t( j1, j1 ), ldt,
204 $ work )
205 CALL dlarfx( 'Right', j2, 3, dtraf, tau, t( 1, j1 ), ldt,
206 $ work )
207*
208 t( j3, j1 ) = zero
209 t( j3, j2 ) = zero
210 t( j3, j3 ) = t11
211*
212 itraf( 1 ) = 2*n + j1
213 li = 2
214 dtraf( 3 ) = tau
215 ld = 4
216 GO TO 40
217*
218 20 CONTINUE
219*
220* N1 = 2, N2 = 1: generate elementary reflector H so that:
221*
222* H ( -X11 ) = ( * )
223* ( -X21 ) = ( 0 )
224* ( scale ) = ( 0 )
225*
226 dtraf( 1 ) = -x( 1, 1 )
227 dtraf( 2 ) = -x( 2, 1 )
228 dtraf( 3 ) = scale
229 CALL dlarfg( 3, dtraf( 1 ), dtraf( 2 ), 1, tau )
230 dtraf( 1 ) = one
231 t33 = t( j3, j3 )
232*
233* Perform swap provisionally on diagonal block in D.
234*
235 CALL dlarfx( 'Left', 3, 3, dtraf, tau, d, ldd, work )
236 CALL dlarfx( 'Right', 3, 3, dtraf, tau, d, ldd, work )
237*
238* Test whether to reject swap.
239*
240 IF( max( abs( d( 2, 1 ) ), abs( d( 3, 1 ) ), abs( d( 1,
241 $ 1 )-t33 ) ).GT.thresh )GO TO 50
242*
243* Accept swap: apply transformation to the entire matrix T.
244*
245 CALL dlarfx( 'Right', j3, 3, dtraf, tau, t( 1, j1 ), ldt,
246 $ work )
247 CALL dlarfx( 'Left', 3, n-j1, dtraf, tau, t( j1, j2 ), ldt,
248 $ work )
249*
250 t( j1, j1 ) = t33
251 t( j2, j1 ) = zero
252 t( j3, j1 ) = zero
253*
254 itraf( 1 ) = n + j1
255 li = 2
256 dtraf( 1 ) = tau
257 ld = 4
258 GO TO 40
259*
260 30 CONTINUE
261*
262* N1 = 2, N2 = 2: generate elementary reflectors H(1) and H(2) so
263* that:
264*
265* H(2) H(1) ( -X11 -X12 ) = ( * * )
266* ( -X21 -X22 ) ( 0 * )
267* ( scale 0 ) ( 0 0 )
268* ( 0 scale ) ( 0 0 )
269*
270 dtraf( 1 ) = -x( 1, 1 )
271 dtraf( 2 ) = -x( 2, 1 )
272 dtraf( 3 ) = scale
273 CALL dlarfg( 3, dtraf( 1 ), dtraf( 2 ), 1, tau1 )
274 dtraf( 1 ) = one
275*
276 temp = -tau1*( x( 1, 2 )+dtraf( 2 )*x( 2, 2 ) )
277 dtraf( 4 ) = -temp*dtraf( 2 ) - x( 2, 2 )
278 dtraf( 5 ) = -temp*dtraf( 3 )
279 dtraf( 6 ) = scale
280 CALL dlarfg( 3, dtraf( 4 ), dtraf( 5 ), 1, tau2 )
281 dtraf( 4 ) = one
282*
283* Perform swap provisionally on diagonal block in D.
284*
285 CALL dlarfx( 'Left', 3, 4, dtraf, tau1, d, ldd, work )
286 CALL dlarfx( 'Right', 4, 3, dtraf, tau1, d, ldd, work )
287 CALL dlarfx( 'Left', 3, 4, dtraf( 4 ), tau2, d( 2, 1 ), ldd,
288 $ work )
289 CALL dlarfx( 'Right', 4, 3, dtraf( 4 ), tau2, d( 1, 2 ), ldd,
290 $ work )
291*
292* Test whether to reject swap.
293*
294 IF( max( abs( d( 3, 1 ) ), abs( d( 3, 2 ) ), abs( d( 4, 1 ) ),
295 $ abs( d( 4, 2 ) ) ).GT.thresh )GO TO 50
296*
297* Accept swap: apply transformation to the entire matrix T.
298*
299 CALL dlarfx( 'Left', 3, n-j1+1, dtraf, tau1, t( j1, j1 ), ldt,
300 $ work )
301 CALL dlarfx( 'Right', j4, 3, dtraf, tau1, t( 1, j1 ), ldt,
302 $ work )
303 CALL dlarfx( 'Left', 3, n-j1+1, dtraf( 4 ), tau2, t( j2, j1 ),
304 $ ldt, work )
305 CALL dlarfx( 'Right', j4, 3, dtraf( 4 ), tau2, t( 1, j2 ), ldt,
306 $ work )
307*
308 t( j3, j1 ) = zero
309 t( j3, j2 ) = zero
310 t( j4, j1 ) = zero
311 t( j4, j2 ) = zero
312*
313 itraf( 1 ) = n + j1
314 itraf( 2 ) = n + j2
315 li = 3
316 dtraf( 1 ) = tau1
317 dtraf( 4 ) = tau2
318 ld = 7
319 GO TO 40
320*
321 40 CONTINUE
322*
323 IF( n2.EQ.2 ) THEN
324*
325* Standardize new 2-by-2 block T11
326*
327 CALL dlanv2( t( j1, j1 ), t( j1, j2 ), t( j2, j1 ),
328 $ t( j2, j2 ), wr1, wi1, wr2, wi2, cs, sn )
329 CALL drot( n-j1-1, t( j1, j1+2 ), ldt, t( j2, j1+2 ), ldt,
330 $ cs, sn )
331 CALL drot( j1-1, t( 1, j1 ), 1, t( 1, j2 ), 1, cs, sn )
332 itraf( li ) = j1
333 li = li + 1
334 dtraf( ld ) = cs
335 dtraf( ld+1 ) = sn
336 ld = ld + 2
337 END IF
338*
339 IF( n1.EQ.2 ) THEN
340*
341* Standardize new 2-by-2 block T22
342*
343 j3 = j1 + n2
344 j4 = j3 + 1
345 CALL dlanv2( t( j3, j3 ), t( j3, j4 ), t( j4, j3 ),
346 $ t( j4, j4 ), wr1, wi1, wr2, wi2, cs, sn )
347 IF( j3+2.LE.n )
348 $ CALL drot( n-j3-1, t( j3, j3+2 ), ldt, t( j4, j3+2 ),
349 $ ldt, cs, sn )
350 CALL drot( j3-1, t( 1, j3 ), 1, t( 1, j4 ), 1, cs, sn )
351 itraf( li ) = j3
352 dtraf( ld ) = cs
353 dtraf( ld+1 ) = sn
354 END IF
355*
356 END IF
357 RETURN
358*
359* Exit with INFO = 1 if swap was rejected.
360*
361 50 CONTINUE
362 info = 1
363 RETURN
364*
365* End of BDLAEXC
366*
dlartg
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition dlartg.f90:113
dlange
double precision function dlange(norm, m, n, a, lda, work)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition dlange.f:114
dlarfx
subroutine dlarfx(side, m, n, v, tau, c, ldc, work)
DLARFX applies an elementary reflector to a general rectangular matrix, with loop unrolling when the ...
Definition dlarfx.f:120
dlanv2
subroutine dlanv2(a, b, c, d, rt1r, rt1i, rt2r, rt2i, cs, sn)
DLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric matrix in standard form.
Definition dlanv2.f:127
dlarfg
subroutine dlarfg(n, alpha, x, incx, tau)
DLARFG generates an elementary reflector (Householder matrix).
Definition dlarfg.f:106
dlasy2
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.
Definition dlasy2.f:174
drot
subroutine drot(n, dx, incx, dy, incy, c, s)
DROT
Definition drot.f:92
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
max
#define max(a, b)
Definition macros.h:21