ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
bdtrexc.f
Go to the documentation of this file.
1  SUBROUTINE bdtrexc( N, T, LDT, IFST, ILST, NITRAF, ITRAF,
2  $ NDTRAF, DTRAF, WORK, INFO )
3  IMPLICIT NONE
4 *
5 *
6 * -- LAPACK routine (version 3.0) --
7 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
8 * Courant Institute, Argonne National Lab, and Rice University
9 * March 31, 1993
10 *
11 * .. Scalar Arguments ..
12  INTEGER IFST, ILST, INFO, LDT, N, NDTRAF, NITRAF
13 * ..
14 * .. Array Arguments ..
15  INTEGER ITRAF( * )
16  DOUBLE PRECISION DTRAF( * ), T( LDT, * ), WORK( * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * BDTREXC reorders the real Schur factorization of a real matrix
23 * A = Q*T*Q**T, so that the diagonal block of T with row index IFST is
24 * moved to row ILST.
25 *
26 * The real Schur form T is reordered by an orthogonal similarity
27 * transformation Z**T*T*Z. In contrast to the LAPACK routine DTREXC,
28 * the orthogonal matrix Z is not explicitly constructed but
29 * represented by paramaters contained in the arrays ITRAF and DTRAF,
30 * see further details.
31 *
32 * T must be in Schur canonical form (as returned by DHSEQR), that is,
33 * block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
34 * 2-by-2 diagonal block has its diagonal elements equal and its
35 * off-diagonal elements of opposite sign.
36 *
37 * Arguments
38 * =========
39 *
40 * N (input) INTEGER
41 * The order of the matrix T. N >= 0.
42 *
43 * T (input/output) DOUBLE PRECISION array, dimension (LDT,N)
44 * On entry, the upper quasi-triangular matrix T, in Schur
45 * Schur canonical form.
46 * On exit, the reordered upper quasi-triangular matrix, again
47 * in Schur canonical form.
48 *
49 * LDT (input) INTEGER
50 * The leading dimension of the array T. LDT >= max(1,N).
51 *
52 * IFST (input/output) INTEGER
53 * ILST (input/output) INTEGER
54 * Specify the reordering of the diagonal blocks of T.
55 * The block with row index IFST is moved to row ILST, by a
56 * sequence of transpositions between adjacent blocks.
57 * On exit, if IFST pointed on entry to the second row of a
58 * 2-by-2 block, it is changed to point to the first row; ILST
59 * always points to the first row of the block in its final
60 * position (which may differ from its input value by +1 or -1).
61 * 1 <= IFST <= N; 1 <= ILST <= N.
62 *
63 * NITRAF (input/output) INTEGER
64 * On entry, length of the array ITRAF.
65 * As a minimum requirement, NITRAF >= max(1,|ILST-IFST|).
66 * If there are 2-by-2 blocks in T then NITRAF must be larger;
67 * a safe choice is NITRAF >= max(1,2*|ILST-IFST|).
68 * On exit, actual length of the array ITRAF.
69 *
70 * ITRAF (output) INTEGER array, length NITRAF
71 * List of parameters for representing the transformation
72 * matrix Z, see further details.
73 *
74 * NDTRAF (input/output) INTEGER
75 * On entry, length of the array DTRAF.
76 * As a minimum requirement, NDTRAF >= max(1,2*|ILST-IFST|).
77 * If there are 2-by-2 blocks in T then NDTRAF must be larger;
78 * a safe choice is NDTRAF >= max(1,5*|ILST-IFST|).
79 * On exit, actual length of the array DTRAF.
80 *
81 * DTRAF (output) DOUBLE PRECISION array, length NDTRAF
82 * List of parameters for representing the transformation
83 * matrix Z, see further details.
84 *
85 * WORK (workspace) DOUBLE PRECISION array, dimension (N)
86 *
87 * INFO (output) INTEGER
88 * = 0: successful exit
89 * < 0: if INFO = -i, the i-th argument had an illegal value
90 * = 1: two adjacent blocks were too close to swap (the problem
91 * is very ill-conditioned); T may have been partially
92 * reordered, and ILST points to the first row of the
93 * current position of the block being moved.
94 * = 2: the 2 by 2 block to be reordered split into two 1 by 1
95 * blocks and the second block failed to swap with an
96 * adjacent block. ILST points to the first row of the
97 * current position of the whole block being moved.
98 *
99 * Further Details
100 * ===============
101 *
102 * The orthogonal transformation matrix Z is a product of NITRAF
103 * elementary orthogonal transformations. The parameters defining these
104 * transformations are stored in the arrays ITRAF and DTRAF as follows:
105 *
106 * Consider the i-th transformation acting on rows/columns POS,
107 * POS+1, ... If this transformation is
108 *
109 * (1) a Givens rotation with cosine C and sine S then
110 *
111 * ITRAF(i) = POS,
112 * DTRAF(j) = C, DTRAF(j+1) = S;
113 *
114 * (2) a Householder reflector H = I - tau * v * v' with
115 * v = [ 1; v2; v3 ] then
116 *
117 * ITRAF(i) = N + POS,
118 * DTRAF(j) = tau, DTRAF(j+1) = v2, DTRAF(j+2) = v3;
119 *
120 * (3) a Householder reflector H = I - tau * v * v' with
121 * v = [ v1; v2; 1 ] then
122 *
123 * ITRAF(i) = 2*N + POS,
124 * DTRAF(j) = v1, DTRAF(j+1) = v2, DTRAF(j+2) = tau;
125 *
126 * Note that the parameters in DTRAF are stored consecutively.
127 *
128 * =====================================================================
129 *
130 * .. Parameters ..
131  DOUBLE PRECISION ZERO
132  parameter( zero = 0.0d+0 )
133  INTEGER DLNGTH(3), ILNGTH(3)
134 * ..
135 * .. Local Scalars ..
136  INTEGER CDTRAF, CITRAF, LDTRAF, LITRAF, HERE, I, NBF,
137  $ nbl, nbnext
138 * ..
139 * .. External Functions ..
140  LOGICAL LSAME
141  EXTERNAL lsame
142 * ..
143 * .. External Subroutines ..
144  EXTERNAL bdlaexc, xerbla
145 * ..
146 * .. Intrinsic Functions ..
147  INTRINSIC abs, max
148 * .. Data Statements ..
149 c DATA ( ILNGTH(I), I = 1, 3 ) / 1, 2, 4 /
150 c DATA ( DLNGTH(I), I = 1, 3 ) / 2, 5, 10 /
151  DATA ilngth(1)/1/, ilngth(2)/2/, ilngth(3)/4/
152  DATA dlngth(1)/2/, dlngth(2)/5/, dlngth(3)/10/
153 * ..
154 * .. Executable Statements ..
155 *
156 * Decode and test the input arguments.
157 *
158  info = 0
159  IF( n.LT.0 ) THEN
160  info = -1
161  ELSE IF( ldt.LT.max( 1, n ) ) THEN
162  info = -3
163  ELSE IF( ifst.LT.1 .OR. ifst.GT.n ) THEN
164  info = -4
165  ELSE IF( ilst.LT.1 .OR. ilst.GT.n ) THEN
166  info = -5
167  ELSE IF ( nitraf.LT.max( 1, abs( ilst-ifst ) ) ) THEN
168  info = -6
169  ELSE IF ( ndtraf.LT.max( 1, 2*abs( ilst-ifst ) ) ) THEN
170  info = -8
171  END IF
172  IF( info.NE.0 ) THEN
173  CALL xerbla( 'DTREXC', -info )
174  RETURN
175  END IF
176 *
177 * Quick return if possible
178 *
179  IF( n.LE.1 )
180  $ RETURN
181  citraf = 1
182  cdtraf = 1
183 *
184 * Determine the first row of specified block
185 * and find out it is 1 by 1 or 2 by 2.
186 *
187  IF( ifst.GT.1 ) THEN
188  IF( t( ifst, ifst-1 ).NE.zero )
189  $ ifst = ifst - 1
190  END IF
191  nbf = 1
192  IF( ifst.LT.n ) THEN
193  IF( t( ifst+1, ifst ).NE.zero )
194  $ nbf = 2
195  END IF
196 *
197 * Determine the first row of the final block
198 * and find out it is 1 by 1 or 2 by 2.
199 *
200  IF( ilst.GT.1 ) THEN
201  IF( t( ilst, ilst-1 ).NE.zero )
202  $ ilst = ilst - 1
203  END IF
204  nbl = 1
205  IF( ilst.LT.n ) THEN
206  IF( t( ilst+1, ilst ).NE.zero )
207  $ nbl = 2
208  END IF
209 *
210  IF( ifst.EQ.ilst )
211  $ RETURN
212 *
213  IF( ifst.LT.ilst ) THEN
214 *
215 * Update ILST
216 *
217  IF( nbf.EQ.2 .AND. nbl.EQ.1 )
218  $ ilst = ilst - 1
219  IF( nbf.EQ.1 .AND. nbl.EQ.2 )
220  $ ilst = ilst + 1
221 *
222  here = ifst
223 *
224  10 CONTINUE
225 *
226 * Swap block with next one below
227 *
228  IF( nbf.EQ.1 .OR. nbf.EQ.2 ) THEN
229 *
230 * Current block either 1 by 1 or 2 by 2
231 *
232  nbnext = 1
233  IF( here+nbf+1.LE.n ) THEN
234  IF( t( here+nbf+1, here+nbf ).NE.zero )
235  $ nbnext = 2
236  END IF
237 *
238  litraf = ilngth(nbf+nbnext-1)
239  ldtraf = dlngth(nbf+nbnext-1)
240  IF( citraf+litraf-1.GT.nitraf ) THEN
241  info = -6
242  CALL xerbla( 'BDTREXC', -info )
243  RETURN
244  END IF
245  IF( cdtraf+ldtraf-1.GT.ndtraf ) THEN
246  info = -8
247  CALL xerbla( 'BDTREXC', -info )
248  RETURN
249  END IF
250  CALL bdlaexc( n, t, ldt, here, nbf, nbnext, itraf(citraf),
251  $ dtraf(cdtraf), work, info )
252  IF( info.NE.0 ) THEN
253  ilst = here
254  nitraf = citraf - 1
255  ndtraf = cdtraf - 1
256  RETURN
257  END IF
258  citraf = citraf + litraf
259  cdtraf = cdtraf + ldtraf
260  here = here + nbnext
261 *
262 * Test if 2 by 2 block breaks into two 1 by 1 blocks
263 *
264  IF( nbf.EQ.2 ) THEN
265  IF( t( here+1, here ).EQ.zero )
266  $ nbf = 3
267  END IF
268 *
269  ELSE
270 *
271 * Current block consists of two 1 by 1 blocks each of which
272 * must be swapped individually
273 *
274  nbnext = 1
275  IF( here+3.LE.n ) THEN
276  IF( t( here+3, here+2 ).NE.zero )
277  $ nbnext = 2
278  END IF
279  litraf = ilngth(nbnext)
280  ldtraf = dlngth(nbnext)
281  IF( citraf+litraf-1.GT.nitraf ) THEN
282  info = -6
283  CALL xerbla( 'BDTREXC', -info )
284  RETURN
285  END IF
286  IF( cdtraf+ldtraf-1.GT.ndtraf ) THEN
287  info = -8
288  CALL xerbla( 'BDTREXC', -info )
289  RETURN
290  END IF
291  CALL bdlaexc( n, t, ldt, here+1, 1, nbnext, itraf(citraf),
292  $ dtraf(cdtraf), work, info )
293  IF( info.NE.0 ) THEN
294  ilst = here
295  nitraf = citraf - 1
296  ndtraf = cdtraf - 1
297  RETURN
298  END IF
299  citraf = citraf + litraf
300  cdtraf = cdtraf + ldtraf
301 *
302  IF( nbnext.EQ.1 ) THEN
303 *
304 * Swap two 1 by 1 blocks, no problems possible
305 *
306  litraf = ilngth(1)
307  ldtraf = dlngth(1)
308  IF( citraf+litraf-1.GT.nitraf ) THEN
309  info = -6
310  CALL xerbla( 'BDTREXC', -info )
311  RETURN
312  END IF
313  IF( cdtraf+ldtraf-1.GT.ndtraf ) THEN
314  info = -8
315  CALL xerbla( 'BDTREXC', -info )
316  RETURN
317  END IF
318  CALL bdlaexc( n, t, ldt, here, 1, nbnext, itraf(citraf),
319  $ dtraf(cdtraf), work, info )
320  citraf = citraf + litraf
321  cdtraf = cdtraf + ldtraf
322  here = here + 1
323  ELSE
324 *
325 * Recompute NBNEXT in case 2 by 2 split
326 *
327  IF( t( here+2, here+1 ).EQ.zero )
328  $ nbnext = 1
329  IF( nbnext.EQ.2 ) THEN
330 *
331 * 2 by 2 Block did not split
332 *
333  litraf = ilngth(2)
334  ldtraf = dlngth(2)
335  IF( citraf+litraf-1.GT.nitraf ) THEN
336  info = -6
337  CALL xerbla( 'BDTREXC', -info )
338  RETURN
339  END IF
340  IF( cdtraf+ldtraf-1.GT.ndtraf ) THEN
341  info = -8
342  CALL xerbla( 'BDTREXC', -info )
343  RETURN
344  END IF
345  CALL bdlaexc( n, t, ldt, here, 1, nbnext,
346  $ itraf(citraf), dtraf(cdtraf), work,
347  $ info )
348  IF( info.NE.0 ) THEN
349  info = 2
350  ilst = here
351  nitraf = citraf - 1
352  ndtraf = cdtraf - 1
353  RETURN
354  END IF
355  citraf = citraf + litraf
356  cdtraf = cdtraf + ldtraf
357  here = here + 2
358  ELSE
359 *
360 * 2 by 2 Block did split
361 *
362  litraf = ilngth(1)
363  ldtraf = dlngth(1)
364  IF( citraf+2*litraf-1.GT.nitraf ) THEN
365  info = -6
366  CALL xerbla( 'BDTREXC', -info )
367  RETURN
368  END IF
369  IF( cdtraf+2*ldtraf-1.GT.ndtraf ) THEN
370  info = -8
371  CALL xerbla( 'BDTREXC', -info )
372  RETURN
373  END IF
374  CALL bdlaexc( n, t, ldt, here, 1, 1, itraf(citraf),
375  $ dtraf(cdtraf), work, info )
376  citraf = citraf + litraf
377  cdtraf = cdtraf + ldtraf
378  CALL bdlaexc( n, t, ldt, here+1, 1, 1, itraf(citraf),
379  $ dtraf(cdtraf), work, info )
380  citraf = citraf + litraf
381  cdtraf = cdtraf + ldtraf
382  here = here + 2
383  END IF
384  END IF
385  END IF
386  IF( here.LT.ilst )
387  $ GO TO 10
388 *
389  ELSE
390 *
391  here = ifst
392  20 CONTINUE
393 *
394 * Swap block with next one above
395 *
396  IF( nbf.EQ.1 .OR. nbf.EQ.2 ) THEN
397 *
398 * Current block either 1 by 1 or 2 by 2
399 *
400  nbnext = 1
401  IF( here.GE.3 ) THEN
402  IF( t( here-1, here-2 ).NE.zero )
403  $ nbnext = 2
404  END IF
405 *
406  litraf = ilngth(nbf+nbnext-1)
407  ldtraf = dlngth(nbf+nbnext-1)
408  IF( citraf+litraf-1.GT.nitraf ) THEN
409  info = -6
410  CALL xerbla( 'BDTREXC', -info )
411  RETURN
412  END IF
413  IF( cdtraf+ldtraf-1.GT.ndtraf ) THEN
414  info = -8
415  CALL xerbla( 'BDTREXC', -info )
416  RETURN
417  END IF
418  CALL bdlaexc( n, t, ldt, here-nbnext, nbnext, nbf,
419  $ itraf(citraf), dtraf(cdtraf), work, info )
420  IF( info.NE.0 ) THEN
421  ilst = here
422  nitraf = citraf - 1
423  ndtraf = cdtraf - 1
424  RETURN
425  END IF
426  citraf = citraf + litraf
427  cdtraf = cdtraf + ldtraf
428  here = here - nbnext
429 *
430 * Test if 2 by 2 block breaks into two 1 by 1 blocks
431 *
432  IF( nbf.EQ.2 ) THEN
433  IF( t( here+1, here ).EQ.zero )
434  $ nbf = 3
435  END IF
436 *
437  ELSE
438 *
439 * Current block consists of two 1 by 1 blocks each of which
440 * must be swapped individually
441 *
442  nbnext = 1
443  IF( here.GE.3 ) THEN
444  IF( t( here-1, here-2 ).NE.zero )
445  $ nbnext = 2
446  END IF
447  litraf = ilngth(nbnext)
448  ldtraf = dlngth(nbnext)
449  IF( citraf+litraf-1.GT.nitraf ) THEN
450  info = -6
451  CALL xerbla( 'BDTREXC', -info )
452  RETURN
453  END IF
454  IF( cdtraf+ldtraf-1.GT.ndtraf ) THEN
455  info = -8
456  CALL xerbla( 'BDTREXC', -info )
457  RETURN
458  END IF
459  CALL bdlaexc( n, t, ldt, here-nbnext, nbnext, 1,
460  $ itraf(citraf), dtraf(cdtraf), work, info )
461  IF( info.NE.0 ) THEN
462  ilst = here
463  nitraf = citraf - 1
464  ndtraf = cdtraf - 1
465  RETURN
466  END IF
467  citraf = citraf + litraf
468  cdtraf = cdtraf + ldtraf
469 *
470  IF( nbnext.EQ.1 ) THEN
471 *
472 * Swap two 1 by 1 blocks, no problems possible
473 *
474  litraf = ilngth(1)
475  ldtraf = dlngth(1)
476  IF( citraf+litraf-1.GT.nitraf ) THEN
477  info = -6
478  CALL xerbla( 'BDTREXC', -info )
479  RETURN
480  END IF
481  IF( cdtraf+ldtraf-1.GT.ndtraf ) THEN
482  info = -8
483  CALL xerbla( 'BDTREXC', -info )
484  RETURN
485  END IF
486  CALL bdlaexc( n, t, ldt, here, nbnext, 1, itraf(citraf),
487  $ dtraf(cdtraf), work, info )
488  citraf = citraf + litraf
489  cdtraf = cdtraf + ldtraf
490  here = here - 1
491  ELSE
492 *
493 * Recompute NBNEXT in case 2 by 2 split
494 *
495  IF( t( here, here-1 ).EQ.zero )
496  $ nbnext = 1
497  IF( nbnext.EQ.2 ) THEN
498 *
499 * 2 by 2 Block did not split
500 *
501  litraf = ilngth(2)
502  ldtraf = dlngth(2)
503  IF( citraf+litraf-1.GT.nitraf ) THEN
504  info = -6
505  CALL xerbla( 'BDTREXC', -info )
506  RETURN
507  END IF
508  IF( cdtraf+ldtraf-1.GT.ndtraf ) THEN
509  info = -8
510  CALL xerbla( 'BDTREXC', -info )
511  RETURN
512  END IF
513  CALL bdlaexc( n, t, ldt, here-1, 2, 1, itraf(citraf),
514  $ dtraf(cdtraf), work, info )
515  IF( info.NE.0 ) THEN
516  info = 2
517  ilst = here
518  nitraf = citraf - 1
519  ndtraf = cdtraf - 1
520  RETURN
521  END IF
522  citraf = citraf + litraf
523  cdtraf = cdtraf + ldtraf
524  here = here - 2
525  ELSE
526 *
527 * 2 by 2 Block did split
528 *
529  litraf = ilngth(1)
530  ldtraf = dlngth(1)
531  IF( citraf+2*litraf-1.GT.nitraf ) THEN
532  info = -6
533  CALL xerbla( 'BDTREXC', -info )
534  RETURN
535  END IF
536  IF( cdtraf+2*ldtraf-1.GT.ndtraf ) THEN
537  info = -8
538  CALL xerbla( 'BDTREXC', -info )
539  RETURN
540  END IF
541  CALL bdlaexc( n, t, ldt, here, 1, 1, itraf(citraf),
542  $ dtraf(cdtraf), work, info )
543  citraf = citraf + litraf
544  cdtraf = cdtraf + ldtraf
545  CALL bdlaexc( n, t, ldt, here-1, 1, 1, itraf(citraf),
546  $ dtraf(cdtraf), work, info )
547  citraf = citraf + litraf
548  cdtraf = cdtraf + ldtraf
549  here = here - 2
550  END IF
551  END IF
552  END IF
553  IF( here.GT.ilst )
554  $ GO TO 20
555  END IF
556  ilst = here
557  nitraf = citraf - 1
558  ndtraf = cdtraf - 1
559 *
560  RETURN
561 *
562 * End of BDTREXC
563 *
564  END
max
#define max(A, B)
Definition: pcgemr.c:180
bdtrexc
subroutine bdtrexc(N, T, LDT, IFST, ILST, NITRAF, ITRAF, NDTRAF, DTRAF, WORK, INFO)
Definition: bdtrexc.f:3
bdlaexc
subroutine bdlaexc(N, T, LDT, J1, N1, N2, ITRAF, DTRAF, WORK, INFO)
Definition: bdlaexc.f:3