ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
bslaexc.f
Go to the documentation of this file.
1  SUBROUTINE bslaexc( N, T, LDT, J1, N1, N2, ITRAF, DTRAF, WORK,
2  $ INFO )
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  $ 1 )-t33 ) ).GT.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  $ abs( d( 4, 2 ) ) ).GT.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  IF( n2.EQ.2 ) 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  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 slanv2( 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 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 *
367  END
max
#define max(A, B)
Definition: pcgemr.c:180
bslaexc
subroutine bslaexc(N, T, LDT, J1, N1, N2, ITRAF, DTRAF, WORK, INFO)
Definition: bslaexc.f:3