SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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 ..
149c DATA ( ILNGTH(I), I = 1, 3 ) / 1, 2, 4 /
150c 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
bdlaexc
subroutine bdlaexc(n, t, ldt, j1, n1, n2, itraf, dtraf, work, info)
Definition bdlaexc.f:3
bdtrexc
subroutine bdtrexc(n, t, ldt, ifst, ilst, nitraf, itraf, ndtraf, dtraf, work, info)
Definition bdtrexc.f:3
max
#define max(A, B)
Definition pcgemr.c:180