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

Go to the source code of this file.

Functions/Subroutines

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

Function/Subroutine Documentation

◆ bslaexc()

subroutine bslaexc ( integer n,
real, dimension( ldt, * ) t,
integer ldt,
integer j1,
integer n1,
integer n2,
integer, dimension( * ) itraf,
real, dimension( * ) dtraf,
real, dimension( * ) work,
integer info )

Definition at line 1 of file bslaexc.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 REAL DTRAF( * ), T( LDT, * ), WORK( * )
15* ..
16*
17* Purpose
18* =======
19*
20* BSLAEXC 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 BSTREXC 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) REAL 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 BSTREXC.
63*
64* DTRAF (output) REAL 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 BSTREXC.
70*
71* WORK (workspace) REAL 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 REAL ZERO, ONE
83 parameter( zero = 0.0e+0, one = 1.0e+0 )
84 REAL TEN
85 parameter( ten = 10.0 )
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 REAL 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 REAL D( LDD, 4 ), X( LDX, 2 )
97* ..
98* .. External Functions ..
99 REAL SLAMCH, SLANGE
100 EXTERNAL slamch, slange
101* ..
102* .. External Subroutines ..
103 EXTERNAL slamov, slanv2, slarfg, slarfx, slartg, slasy2,
104 $ srot
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 slartg( t( j1, j2 ), t22-t11, cs, sn, temp )
134*
135* Apply transformation to the matrix T.
136*
137 IF( j3.LE.n )
138 $ CALL srot( n-j1-1, t( j1, j3 ), ldt, t( j2, j3 ), ldt, cs,
139 $ sn )
140 CALL srot( 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 slamov( 'Full', nd, nd, t( j1, j1 ), ldt, d, ldd )
158 dnorm = slange( 'Max', nd, nd, d, ldd, work )
159*
160* Compute machine-dependent threshold for test for accepting
161* swap.
162*
163 eps = slamch( 'P' )
164 smlnum = slamch( 'S' ) / eps
165 thresh = max( ten*eps*dnorm, smlnum )
166*
167* Solve T11*X - X*T22 = scale*T12 for X.
168*
169 CALL slasy2( .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 slarfg( 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 slarfx( 'Left', 3, 3, dtraf, tau, d, ldd, work )
194 CALL slarfx( '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 slarfx( 'Left', 3, n-j1+1, dtraf, tau, t( j1, j1 ), ldt,
204 $ work )
205 CALL slarfx( '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 slarfg( 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 slarfx( 'Left', 3, 3, dtraf, tau, d, ldd, work )
236 CALL slarfx( '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.GT. $ 1 )-T33 ) )THRESH )GO TO 50
242*
243* Accept swap: apply transformation to the entire matrix T.
244*
245 CALL SLARFX( 'right', J3, 3, DTRAF, TAU, T( 1, J1 ), LDT,
246 $ WORK )
247 CALL SLARFX( '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 SLARFG( 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 SLARFG( 3, DTRAF( 4 ), DTRAF( 5 ), 1, TAU2 )
281 DTRAF( 4 ) = ONE
282*
283* Perform swap provisionally on diagonal block in D.
284*
285 CALL SLARFX( 'left', 3, 4, DTRAF, TAU1, D, LDD, WORK )
286 CALL SLARFX( 'right', 4, 3, DTRAF, TAU1, D, LDD, WORK )
287 CALL SLARFX( 'left', 3, 4, DTRAF( 4 ), TAU2, D( 2, 1 ), LDD,
288 $ WORK )
289 CALL SLARFX( '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.GT. $ ABS( D( 4, 2 ) ) )THRESH )GO TO 50
296*
297* Accept swap: apply transformation to the entire matrix T.
298*
299 CALL SLARFX( 'left', 3, N-J1+1, DTRAF, TAU1, T( J1, J1 ), LDT,
300 $ WORK )
301 CALL SLARFX( 'right', J4, 3, DTRAF, TAU1, T( 1, J1 ), LDT,
302 $ WORK )
303 CALL SLARFX( 'left', 3, N-J1+1, DTRAF( 4 ), TAU2, T( J2, J1 ),
304 $ LDT, WORK )
305 CALL SLARFX( '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.EQ. IF( N22 ) THEN
324*
325* Standardize new 2-by-2 block T11
326*
327 CALL SLANV2( T( J1, J1 ), T( J1, J2 ), T( J2, J1 ),
328 $ T( J2, J2 ), WR1, WI1, WR2, WI2, CS, SN )
329 CALL SROT( N-J1-1, T( J1, J1+2 ), LDT, T( J2, J1+2 ), LDT,
330 $ CS, SN )
331 CALL SROT( 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.EQ. IF( N12 ) THEN
340*
341* Standardize new 2-by-2 block T22
342*
343 J3 = J1 + N2
344 J4 = J3 + 1
345 CALL SLANV2( T( J3, J3 ), T( J3, J4 ), T( J4, J3 ),
346 $ T( J4, J4 ), WR1, WI1, WR2, WI2, CS, SN )
347.LE. IF( J3+2N )
348 $ CALL SROT( N-J3-1, T( J3, J3+2 ), LDT, T( J4, J3+2 ),
349 $ LDT, CS, SN )
350 CALL SROT( 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 BSLAEXC
366*
slartg
subroutine slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
Definition slartg.f90:113
slange
real function slange(norm, m, n, a, lda, work)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition slange.f:114
slanv2
subroutine slanv2(a, b, c, d, rt1r, rt1i, rt2r, rt2i, cs, sn)
SLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric matrix in standard form.
Definition slanv2.f:127
slarfg
subroutine slarfg(n, alpha, x, incx, tau)
SLARFG generates an elementary reflector (Householder matrix).
Definition slarfg.f:106
slarfx
subroutine slarfx(side, m, n, v, tau, c, ldc, work)
SLARFX applies an elementary reflector to a general rectangular matrix, with loop unrolling when the ...
Definition slarfx.f:120
slasy2
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.
Definition slasy2.f:174
srot
subroutine srot(n, sx, incx, sy, incy, c, s)
SROT
Definition srot.f:92
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
max
#define max(a, b)
Definition macros.h:21