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

Go to the source code of this file.

Functions/Subroutines

subroutine bdtrexc (n, t, ldt, ifst, ilst, nitraf, itraf, ndtraf, dtraf, work, info)

Function/Subroutine Documentation

◆ bdtrexc()

subroutine bdtrexc ( integer n,
double precision, dimension( ldt, * ) t,
integer ldt,
integer ifst,
integer ilst,
integer nitraf,
integer, dimension( * ) itraf,
integer ndtraf,
double precision, dimension( * ) dtraf,
double precision, dimension( * ) work,
integer info )

Definition at line 1 of file bdtrexc.f.

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.LT. IF( HEREILST )
387 $ GO TO 10
388*
389 ELSE
390*
391 HERE = IFST
392 20 CONTINUE
393*
394* Swap block with next one above
395*
396.EQ..OR..EQ. IF( NBF1 NBF2 ) THEN
397*
398* Current block either 1 by 1 or 2 by 2
399*
400 NBNEXT = 1
401.GE. IF( HERE3 ) THEN
402.NE. IF( T( HERE-1, HERE-2 )ZERO )
403 $ NBNEXT = 2
404 END IF
405*
406 LITRAF = ILNGTH(NBF+NBNEXT-1)
407 LDTRAF = DLNGTH(NBF+NBNEXT-1)
408.GT. IF( CITRAF+LITRAF-1NITRAF ) THEN
409 INFO = -6
410 CALL XERBLA( 'bdtrexc', -INFO )
411 RETURN
412 END IF
413.GT. IF( CDTRAF+LDTRAF-1NDTRAF ) 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.NE. IF( INFO0 ) 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.EQ. IF( NBF2 ) THEN
433.EQ. IF( T( HERE+1, HERE )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.GE. IF( HERE3 ) THEN
444.NE. IF( T( HERE-1, HERE-2 )ZERO )
445 $ NBNEXT = 2
446 END IF
447 LITRAF = ILNGTH(NBNEXT)
448 LDTRAF = DLNGTH(NBNEXT)
449.GT. IF( CITRAF+LITRAF-1NITRAF ) THEN
450 INFO = -6
451 CALL XERBLA( 'bdtrexc', -INFO )
452 RETURN
453 END IF
454.GT. IF( CDTRAF+LDTRAF-1NDTRAF ) 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.NE. IF( INFO0 ) 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.EQ. IF( NBNEXT1 ) THEN
471*
472* Swap two 1 by 1 blocks, no problems possible
473*
474 LITRAF = ILNGTH(1)
475 LDTRAF = DLNGTH(1)
476.GT. IF( CITRAF+LITRAF-1NITRAF ) THEN
477 INFO = -6
478 CALL XERBLA( 'bdtrexc', -INFO )
479 RETURN
480 END IF
481.GT. IF( CDTRAF+LDTRAF-1NDTRAF ) 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.EQ. IF( T( HERE, HERE-1 )ZERO )
496 $ NBNEXT = 1
497.EQ. IF( NBNEXT2 ) THEN
498*
499* 2 by 2 Block did not split
500*
501 LITRAF = ILNGTH(2)
502 LDTRAF = DLNGTH(2)
503.GT. IF( CITRAF+LITRAF-1NITRAF ) THEN
504 INFO = -6
505 CALL XERBLA( 'bdtrexc', -INFO )
506 RETURN
507 END IF
508.GT. IF( CDTRAF+LDTRAF-1NDTRAF ) 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.NE. IF( INFO0 ) 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.GT. IF( CITRAF+2*LITRAF-1NITRAF ) THEN
532 INFO = -6
533 CALL XERBLA( 'bdtrexc', -INFO )
534 RETURN
535 END IF
536.GT. IF( CDTRAF+2*LDTRAF-1NDTRAF ) 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.GT. IF( HEREILST )
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*
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
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
max
#define max(a, b)
Definition macros.h:21