Annotation of rpl/lapack/lapack/dlaexc.f, revision 1.9
1.9 ! bertrand 1: *> \brief \b DLAEXC
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DLAEXC + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaexc.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaexc.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaexc.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DLAEXC( WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK,
! 22: * INFO )
! 23: *
! 24: * .. Scalar Arguments ..
! 25: * LOGICAL WANTQ
! 26: * INTEGER INFO, J1, LDQ, LDT, N, N1, N2
! 27: * ..
! 28: * .. Array Arguments ..
! 29: * DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WORK( * )
! 30: * ..
! 31: *
! 32: *
! 33: *> \par Purpose:
! 34: * =============
! 35: *>
! 36: *> \verbatim
! 37: *>
! 38: *> DLAEXC swaps adjacent diagonal blocks T11 and T22 of order 1 or 2 in
! 39: *> an upper quasi-triangular matrix T by an orthogonal similarity
! 40: *> transformation.
! 41: *>
! 42: *> T must be in Schur canonical form, that is, block upper triangular
! 43: *> with 1-by-1 and 2-by-2 diagonal blocks; each 2-by-2 diagonal block
! 44: *> has its diagonal elemnts equal and its off-diagonal elements of
! 45: *> opposite sign.
! 46: *> \endverbatim
! 47: *
! 48: * Arguments:
! 49: * ==========
! 50: *
! 51: *> \param[in] WANTQ
! 52: *> \verbatim
! 53: *> WANTQ is LOGICAL
! 54: *> = .TRUE. : accumulate the transformation in the matrix Q;
! 55: *> = .FALSE.: do not accumulate the transformation.
! 56: *> \endverbatim
! 57: *>
! 58: *> \param[in] N
! 59: *> \verbatim
! 60: *> N is INTEGER
! 61: *> The order of the matrix T. N >= 0.
! 62: *> \endverbatim
! 63: *>
! 64: *> \param[in,out] T
! 65: *> \verbatim
! 66: *> T is DOUBLE PRECISION array, dimension (LDT,N)
! 67: *> On entry, the upper quasi-triangular matrix T, in Schur
! 68: *> canonical form.
! 69: *> On exit, the updated matrix T, again in Schur canonical form.
! 70: *> \endverbatim
! 71: *>
! 72: *> \param[in] LDT
! 73: *> \verbatim
! 74: *> LDT is INTEGER
! 75: *> The leading dimension of the array T. LDT >= max(1,N).
! 76: *> \endverbatim
! 77: *>
! 78: *> \param[in,out] Q
! 79: *> \verbatim
! 80: *> Q is DOUBLE PRECISION array, dimension (LDQ,N)
! 81: *> On entry, if WANTQ is .TRUE., the orthogonal matrix Q.
! 82: *> On exit, if WANTQ is .TRUE., the updated matrix Q.
! 83: *> If WANTQ is .FALSE., Q is not referenced.
! 84: *> \endverbatim
! 85: *>
! 86: *> \param[in] LDQ
! 87: *> \verbatim
! 88: *> LDQ is INTEGER
! 89: *> The leading dimension of the array Q.
! 90: *> LDQ >= 1; and if WANTQ is .TRUE., LDQ >= N.
! 91: *> \endverbatim
! 92: *>
! 93: *> \param[in] J1
! 94: *> \verbatim
! 95: *> J1 is INTEGER
! 96: *> The index of the first row of the first block T11.
! 97: *> \endverbatim
! 98: *>
! 99: *> \param[in] N1
! 100: *> \verbatim
! 101: *> N1 is INTEGER
! 102: *> The order of the first block T11. N1 = 0, 1 or 2.
! 103: *> \endverbatim
! 104: *>
! 105: *> \param[in] N2
! 106: *> \verbatim
! 107: *> N2 is INTEGER
! 108: *> The order of the second block T22. N2 = 0, 1 or 2.
! 109: *> \endverbatim
! 110: *>
! 111: *> \param[out] WORK
! 112: *> \verbatim
! 113: *> WORK is DOUBLE PRECISION array, dimension (N)
! 114: *> \endverbatim
! 115: *>
! 116: *> \param[out] INFO
! 117: *> \verbatim
! 118: *> INFO is INTEGER
! 119: *> = 0: successful exit
! 120: *> = 1: the transformed matrix T would be too far from Schur
! 121: *> form; the blocks are not swapped and T and Q are
! 122: *> unchanged.
! 123: *> \endverbatim
! 124: *
! 125: * Authors:
! 126: * ========
! 127: *
! 128: *> \author Univ. of Tennessee
! 129: *> \author Univ. of California Berkeley
! 130: *> \author Univ. of Colorado Denver
! 131: *> \author NAG Ltd.
! 132: *
! 133: *> \date November 2011
! 134: *
! 135: *> \ingroup doubleOTHERauxiliary
! 136: *
! 137: * =====================================================================
1.1 bertrand 138: SUBROUTINE DLAEXC( WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK,
139: $ INFO )
140: *
1.9 ! bertrand 141: * -- LAPACK auxiliary routine (version 3.4.0) --
1.1 bertrand 142: * -- LAPACK is a software package provided by Univ. of Tennessee, --
143: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 144: * November 2011
1.1 bertrand 145: *
146: * .. Scalar Arguments ..
147: LOGICAL WANTQ
148: INTEGER INFO, J1, LDQ, LDT, N, N1, N2
149: * ..
150: * .. Array Arguments ..
151: DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WORK( * )
152: * ..
153: *
154: * =====================================================================
155: *
156: * .. Parameters ..
157: DOUBLE PRECISION ZERO, ONE
158: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
159: DOUBLE PRECISION TEN
160: PARAMETER ( TEN = 1.0D+1 )
161: INTEGER LDD, LDX
162: PARAMETER ( LDD = 4, LDX = 2 )
163: * ..
164: * .. Local Scalars ..
165: INTEGER IERR, J2, J3, J4, K, ND
166: DOUBLE PRECISION CS, DNORM, EPS, SCALE, SMLNUM, SN, T11, T22,
167: $ T33, TAU, TAU1, TAU2, TEMP, THRESH, WI1, WI2,
168: $ WR1, WR2, XNORM
169: * ..
170: * .. Local Arrays ..
171: DOUBLE PRECISION D( LDD, 4 ), U( 3 ), U1( 3 ), U2( 3 ),
172: $ X( LDX, 2 )
173: * ..
174: * .. External Functions ..
175: DOUBLE PRECISION DLAMCH, DLANGE
176: EXTERNAL DLAMCH, DLANGE
177: * ..
178: * .. External Subroutines ..
179: EXTERNAL DLACPY, DLANV2, DLARFG, DLARFX, DLARTG, DLASY2,
180: $ DROT
181: * ..
182: * .. Intrinsic Functions ..
183: INTRINSIC ABS, MAX
184: * ..
185: * .. Executable Statements ..
186: *
187: INFO = 0
188: *
189: * Quick return if possible
190: *
191: IF( N.EQ.0 .OR. N1.EQ.0 .OR. N2.EQ.0 )
192: $ RETURN
193: IF( J1+N1.GT.N )
194: $ RETURN
195: *
196: J2 = J1 + 1
197: J3 = J1 + 2
198: J4 = J1 + 3
199: *
200: IF( N1.EQ.1 .AND. N2.EQ.1 ) THEN
201: *
202: * Swap two 1-by-1 blocks.
203: *
204: T11 = T( J1, J1 )
205: T22 = T( J2, J2 )
206: *
207: * Determine the transformation to perform the interchange.
208: *
209: CALL DLARTG( T( J1, J2 ), T22-T11, CS, SN, TEMP )
210: *
211: * Apply transformation to the matrix T.
212: *
213: IF( J3.LE.N )
214: $ CALL DROT( N-J1-1, T( J1, J3 ), LDT, T( J2, J3 ), LDT, CS,
215: $ SN )
216: CALL DROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN )
217: *
218: T( J1, J1 ) = T22
219: T( J2, J2 ) = T11
220: *
221: IF( WANTQ ) THEN
222: *
223: * Accumulate transformation in the matrix Q.
224: *
225: CALL DROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN )
226: END IF
227: *
228: ELSE
229: *
230: * Swapping involves at least one 2-by-2 block.
231: *
232: * Copy the diagonal block of order N1+N2 to the local array D
233: * and compute its norm.
234: *
235: ND = N1 + N2
236: CALL DLACPY( 'Full', ND, ND, T( J1, J1 ), LDT, D, LDD )
237: DNORM = DLANGE( 'Max', ND, ND, D, LDD, WORK )
238: *
239: * Compute machine-dependent threshold for test for accepting
240: * swap.
241: *
242: EPS = DLAMCH( 'P' )
243: SMLNUM = DLAMCH( 'S' ) / EPS
244: THRESH = MAX( TEN*EPS*DNORM, SMLNUM )
245: *
246: * Solve T11*X - X*T22 = scale*T12 for X.
247: *
248: CALL DLASY2( .FALSE., .FALSE., -1, N1, N2, D, LDD,
249: $ D( N1+1, N1+1 ), LDD, D( 1, N1+1 ), LDD, SCALE, X,
250: $ LDX, XNORM, IERR )
251: *
252: * Swap the adjacent diagonal blocks.
253: *
254: K = N1 + N1 + N2 - 3
255: GO TO ( 10, 20, 30 )K
256: *
257: 10 CONTINUE
258: *
259: * N1 = 1, N2 = 2: generate elementary reflector H so that:
260: *
261: * ( scale, X11, X12 ) H = ( 0, 0, * )
262: *
263: U( 1 ) = SCALE
264: U( 2 ) = X( 1, 1 )
265: U( 3 ) = X( 1, 2 )
266: CALL DLARFG( 3, U( 3 ), U, 1, TAU )
267: U( 3 ) = ONE
268: T11 = T( J1, J1 )
269: *
270: * Perform swap provisionally on diagonal block in D.
271: *
272: CALL DLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK )
273: CALL DLARFX( 'R', 3, 3, U, TAU, D, LDD, WORK )
274: *
275: * Test whether to reject swap.
276: *
277: IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 3,
278: $ 3 )-T11 ) ).GT.THRESH )GO TO 50
279: *
280: * Accept swap: apply transformation to the entire matrix T.
281: *
282: CALL DLARFX( 'L', 3, N-J1+1, U, TAU, T( J1, J1 ), LDT, WORK )
283: CALL DLARFX( 'R', J2, 3, U, TAU, T( 1, J1 ), LDT, WORK )
284: *
285: T( J3, J1 ) = ZERO
286: T( J3, J2 ) = ZERO
287: T( J3, J3 ) = T11
288: *
289: IF( WANTQ ) THEN
290: *
291: * Accumulate transformation in the matrix Q.
292: *
293: CALL DLARFX( 'R', N, 3, U, TAU, Q( 1, J1 ), LDQ, WORK )
294: END IF
295: GO TO 40
296: *
297: 20 CONTINUE
298: *
299: * N1 = 2, N2 = 1: generate elementary reflector H so that:
300: *
301: * H ( -X11 ) = ( * )
302: * ( -X21 ) = ( 0 )
303: * ( scale ) = ( 0 )
304: *
305: U( 1 ) = -X( 1, 1 )
306: U( 2 ) = -X( 2, 1 )
307: U( 3 ) = SCALE
308: CALL DLARFG( 3, U( 1 ), U( 2 ), 1, TAU )
309: U( 1 ) = ONE
310: T33 = T( J3, J3 )
311: *
312: * Perform swap provisionally on diagonal block in D.
313: *
314: CALL DLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK )
315: CALL DLARFX( 'R', 3, 3, U, TAU, D, LDD, WORK )
316: *
317: * Test whether to reject swap.
318: *
319: IF( MAX( ABS( D( 2, 1 ) ), ABS( D( 3, 1 ) ), ABS( D( 1,
320: $ 1 )-T33 ) ).GT.THRESH )GO TO 50
321: *
322: * Accept swap: apply transformation to the entire matrix T.
323: *
324: CALL DLARFX( 'R', J3, 3, U, TAU, T( 1, J1 ), LDT, WORK )
325: CALL DLARFX( 'L', 3, N-J1, U, TAU, T( J1, J2 ), LDT, WORK )
326: *
327: T( J1, J1 ) = T33
328: T( J2, J1 ) = ZERO
329: T( J3, J1 ) = ZERO
330: *
331: IF( WANTQ ) THEN
332: *
333: * Accumulate transformation in the matrix Q.
334: *
335: CALL DLARFX( 'R', N, 3, U, TAU, Q( 1, J1 ), LDQ, WORK )
336: END IF
337: GO TO 40
338: *
339: 30 CONTINUE
340: *
341: * N1 = 2, N2 = 2: generate elementary reflectors H(1) and H(2) so
342: * that:
343: *
344: * H(2) H(1) ( -X11 -X12 ) = ( * * )
345: * ( -X21 -X22 ) ( 0 * )
346: * ( scale 0 ) ( 0 0 )
347: * ( 0 scale ) ( 0 0 )
348: *
349: U1( 1 ) = -X( 1, 1 )
350: U1( 2 ) = -X( 2, 1 )
351: U1( 3 ) = SCALE
352: CALL DLARFG( 3, U1( 1 ), U1( 2 ), 1, TAU1 )
353: U1( 1 ) = ONE
354: *
355: TEMP = -TAU1*( X( 1, 2 )+U1( 2 )*X( 2, 2 ) )
356: U2( 1 ) = -TEMP*U1( 2 ) - X( 2, 2 )
357: U2( 2 ) = -TEMP*U1( 3 )
358: U2( 3 ) = SCALE
359: CALL DLARFG( 3, U2( 1 ), U2( 2 ), 1, TAU2 )
360: U2( 1 ) = ONE
361: *
362: * Perform swap provisionally on diagonal block in D.
363: *
364: CALL DLARFX( 'L', 3, 4, U1, TAU1, D, LDD, WORK )
365: CALL DLARFX( 'R', 4, 3, U1, TAU1, D, LDD, WORK )
366: CALL DLARFX( 'L', 3, 4, U2, TAU2, D( 2, 1 ), LDD, WORK )
367: CALL DLARFX( 'R', 4, 3, U2, TAU2, D( 1, 2 ), LDD, WORK )
368: *
369: * Test whether to reject swap.
370: *
371: IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 4, 1 ) ),
372: $ ABS( D( 4, 2 ) ) ).GT.THRESH )GO TO 50
373: *
374: * Accept swap: apply transformation to the entire matrix T.
375: *
376: CALL DLARFX( 'L', 3, N-J1+1, U1, TAU1, T( J1, J1 ), LDT, WORK )
377: CALL DLARFX( 'R', J4, 3, U1, TAU1, T( 1, J1 ), LDT, WORK )
378: CALL DLARFX( 'L', 3, N-J1+1, U2, TAU2, T( J2, J1 ), LDT, WORK )
379: CALL DLARFX( 'R', J4, 3, U2, TAU2, T( 1, J2 ), LDT, WORK )
380: *
381: T( J3, J1 ) = ZERO
382: T( J3, J2 ) = ZERO
383: T( J4, J1 ) = ZERO
384: T( J4, J2 ) = ZERO
385: *
386: IF( WANTQ ) THEN
387: *
388: * Accumulate transformation in the matrix Q.
389: *
390: CALL DLARFX( 'R', N, 3, U1, TAU1, Q( 1, J1 ), LDQ, WORK )
391: CALL DLARFX( 'R', N, 3, U2, TAU2, Q( 1, J2 ), LDQ, WORK )
392: END IF
393: *
394: 40 CONTINUE
395: *
396: IF( N2.EQ.2 ) THEN
397: *
398: * Standardize new 2-by-2 block T11
399: *
400: CALL DLANV2( T( J1, J1 ), T( J1, J2 ), T( J2, J1 ),
401: $ T( J2, J2 ), WR1, WI1, WR2, WI2, CS, SN )
402: CALL DROT( N-J1-1, T( J1, J1+2 ), LDT, T( J2, J1+2 ), LDT,
403: $ CS, SN )
404: CALL DROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN )
405: IF( WANTQ )
406: $ CALL DROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN )
407: END IF
408: *
409: IF( N1.EQ.2 ) THEN
410: *
411: * Standardize new 2-by-2 block T22
412: *
413: J3 = J1 + N2
414: J4 = J3 + 1
415: CALL DLANV2( T( J3, J3 ), T( J3, J4 ), T( J4, J3 ),
416: $ T( J4, J4 ), WR1, WI1, WR2, WI2, CS, SN )
417: IF( J3+2.LE.N )
418: $ CALL DROT( N-J3-1, T( J3, J3+2 ), LDT, T( J4, J3+2 ),
419: $ LDT, CS, SN )
420: CALL DROT( J3-1, T( 1, J3 ), 1, T( 1, J4 ), 1, CS, SN )
421: IF( WANTQ )
422: $ CALL DROT( N, Q( 1, J3 ), 1, Q( 1, J4 ), 1, CS, SN )
423: END IF
424: *
425: END IF
426: RETURN
427: *
428: * Exit with INFO = 1 if swap was rejected.
429: *
430: 50 CONTINUE
431: INFO = 1
432: RETURN
433: *
434: * End of DLAEXC
435: *
436: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaexc.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack