Annotation of rpl/lapack/lapack/dgges.f, revision 1.4
1.1 bertrand 1: SUBROUTINE DGGES( JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B, LDB,
2: $ SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR,
3: $ LDVSR, WORK, LWORK, BWORK, INFO )
4: *
5: * -- LAPACK driver routine (version 3.2) --
6: * -- LAPACK is a software package provided by Univ. of Tennessee, --
7: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8: * November 2006
9: *
10: * .. Scalar Arguments ..
11: CHARACTER JOBVSL, JOBVSR, SORT
12: INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM
13: * ..
14: * .. Array Arguments ..
15: LOGICAL BWORK( * )
16: DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
17: $ B( LDB, * ), BETA( * ), VSL( LDVSL, * ),
18: $ VSR( LDVSR, * ), WORK( * )
19: * ..
20: * .. Function Arguments ..
21: LOGICAL SELCTG
22: EXTERNAL SELCTG
23: * ..
24: *
25: * Purpose
26: * =======
27: *
28: * DGGES computes for a pair of N-by-N real nonsymmetric matrices (A,B),
29: * the generalized eigenvalues, the generalized real Schur form (S,T),
30: * optionally, the left and/or right matrices of Schur vectors (VSL and
31: * VSR). This gives the generalized Schur factorization
32: *
33: * (A,B) = ( (VSL)*S*(VSR)**T, (VSL)*T*(VSR)**T )
34: *
35: * Optionally, it also orders the eigenvalues so that a selected cluster
36: * of eigenvalues appears in the leading diagonal blocks of the upper
37: * quasi-triangular matrix S and the upper triangular matrix T.The
38: * leading columns of VSL and VSR then form an orthonormal basis for the
39: * corresponding left and right eigenspaces (deflating subspaces).
40: *
41: * (If only the generalized eigenvalues are needed, use the driver
42: * DGGEV instead, which is faster.)
43: *
44: * A generalized eigenvalue for a pair of matrices (A,B) is a scalar w
45: * or a ratio alpha/beta = w, such that A - w*B is singular. It is
46: * usually represented as the pair (alpha,beta), as there is a
47: * reasonable interpretation for beta=0 or both being zero.
48: *
49: * A pair of matrices (S,T) is in generalized real Schur form if T is
50: * upper triangular with non-negative diagonal and S is block upper
51: * triangular with 1-by-1 and 2-by-2 blocks. 1-by-1 blocks correspond
52: * to real generalized eigenvalues, while 2-by-2 blocks of S will be
53: * "standardized" by making the corresponding elements of T have the
54: * form:
55: * [ a 0 ]
56: * [ 0 b ]
57: *
58: * and the pair of corresponding 2-by-2 blocks in S and T will have a
59: * complex conjugate pair of generalized eigenvalues.
60: *
61: *
62: * Arguments
63: * =========
64: *
65: * JOBVSL (input) CHARACTER*1
66: * = 'N': do not compute the left Schur vectors;
67: * = 'V': compute the left Schur vectors.
68: *
69: * JOBVSR (input) CHARACTER*1
70: * = 'N': do not compute the right Schur vectors;
71: * = 'V': compute the right Schur vectors.
72: *
73: * SORT (input) CHARACTER*1
74: * Specifies whether or not to order the eigenvalues on the
75: * diagonal of the generalized Schur form.
76: * = 'N': Eigenvalues are not ordered;
77: * = 'S': Eigenvalues are ordered (see SELCTG);
78: *
79: * SELCTG (external procedure) LOGICAL FUNCTION of three DOUBLE PRECISION arguments
80: * SELCTG must be declared EXTERNAL in the calling subroutine.
81: * If SORT = 'N', SELCTG is not referenced.
82: * If SORT = 'S', SELCTG is used to select eigenvalues to sort
83: * to the top left of the Schur form.
84: * An eigenvalue (ALPHAR(j)+ALPHAI(j))/BETA(j) is selected if
85: * SELCTG(ALPHAR(j),ALPHAI(j),BETA(j)) is true; i.e. if either
86: * one of a complex conjugate pair of eigenvalues is selected,
87: * then both complex eigenvalues are selected.
88: *
89: * Note that in the ill-conditioned case, a selected complex
90: * eigenvalue may no longer satisfy SELCTG(ALPHAR(j),ALPHAI(j),
91: * BETA(j)) = .TRUE. after ordering. INFO is to be set to N+2
92: * in this case.
93: *
94: * N (input) INTEGER
95: * The order of the matrices A, B, VSL, and VSR. N >= 0.
96: *
97: * A (input/output) DOUBLE PRECISION array, dimension (LDA, N)
98: * On entry, the first of the pair of matrices.
99: * On exit, A has been overwritten by its generalized Schur
100: * form S.
101: *
102: * LDA (input) INTEGER
103: * The leading dimension of A. LDA >= max(1,N).
104: *
105: * B (input/output) DOUBLE PRECISION array, dimension (LDB, N)
106: * On entry, the second of the pair of matrices.
107: * On exit, B has been overwritten by its generalized Schur
108: * form T.
109: *
110: * LDB (input) INTEGER
111: * The leading dimension of B. LDB >= max(1,N).
112: *
113: * SDIM (output) INTEGER
114: * If SORT = 'N', SDIM = 0.
115: * If SORT = 'S', SDIM = number of eigenvalues (after sorting)
116: * for which SELCTG is true. (Complex conjugate pairs for which
117: * SELCTG is true for either eigenvalue count as 2.)
118: *
119: * ALPHAR (output) DOUBLE PRECISION array, dimension (N)
120: * ALPHAI (output) DOUBLE PRECISION array, dimension (N)
121: * BETA (output) DOUBLE PRECISION array, dimension (N)
122: * On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
123: * be the generalized eigenvalues. ALPHAR(j) + ALPHAI(j)*i,
124: * and BETA(j),j=1,...,N are the diagonals of the complex Schur
125: * form (S,T) that would result if the 2-by-2 diagonal blocks of
126: * the real Schur form of (A,B) were further reduced to
127: * triangular form using 2-by-2 complex unitary transformations.
128: * If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
129: * positive, then the j-th and (j+1)-st eigenvalues are a
130: * complex conjugate pair, with ALPHAI(j+1) negative.
131: *
132: * Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)
133: * may easily over- or underflow, and BETA(j) may even be zero.
134: * Thus, the user should avoid naively computing the ratio.
135: * However, ALPHAR and ALPHAI will be always less than and
136: * usually comparable with norm(A) in magnitude, and BETA always
137: * less than and usually comparable with norm(B).
138: *
139: * VSL (output) DOUBLE PRECISION array, dimension (LDVSL,N)
140: * If JOBVSL = 'V', VSL will contain the left Schur vectors.
141: * Not referenced if JOBVSL = 'N'.
142: *
143: * LDVSL (input) INTEGER
144: * The leading dimension of the matrix VSL. LDVSL >=1, and
145: * if JOBVSL = 'V', LDVSL >= N.
146: *
147: * VSR (output) DOUBLE PRECISION array, dimension (LDVSR,N)
148: * If JOBVSR = 'V', VSR will contain the right Schur vectors.
149: * Not referenced if JOBVSR = 'N'.
150: *
151: * LDVSR (input) INTEGER
152: * The leading dimension of the matrix VSR. LDVSR >= 1, and
153: * if JOBVSR = 'V', LDVSR >= N.
154: *
155: * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
156: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
157: *
158: * LWORK (input) INTEGER
159: * The dimension of the array WORK.
160: * If N = 0, LWORK >= 1, else LWORK >= 8*N+16.
161: * For good performance , LWORK must generally be larger.
162: *
163: * If LWORK = -1, then a workspace query is assumed; the routine
164: * only calculates the optimal size of the WORK array, returns
165: * this value as the first entry of the WORK array, and no error
166: * message related to LWORK is issued by XERBLA.
167: *
168: * BWORK (workspace) LOGICAL array, dimension (N)
169: * Not referenced if SORT = 'N'.
170: *
171: * INFO (output) INTEGER
172: * = 0: successful exit
173: * < 0: if INFO = -i, the i-th argument had an illegal value.
174: * = 1,...,N:
175: * The QZ iteration failed. (A,B) are not in Schur
176: * form, but ALPHAR(j), ALPHAI(j), and BETA(j) should
177: * be correct for j=INFO+1,...,N.
178: * > N: =N+1: other than QZ iteration failed in DHGEQZ.
179: * =N+2: after reordering, roundoff changed values of
180: * some complex eigenvalues so that leading
181: * eigenvalues in the Generalized Schur form no
182: * longer satisfy SELCTG=.TRUE. This could also
183: * be caused due to scaling.
184: * =N+3: reordering failed in DTGSEN.
185: *
186: * =====================================================================
187: *
188: * .. Parameters ..
189: DOUBLE PRECISION ZERO, ONE
190: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
191: * ..
192: * .. Local Scalars ..
193: LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
194: $ LQUERY, LST2SL, WANTST
195: INTEGER I, ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT,
196: $ ILO, IP, IRIGHT, IROWS, ITAU, IWRK, MAXWRK,
197: $ MINWRK
198: DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PVSL,
199: $ PVSR, SAFMAX, SAFMIN, SMLNUM
200: * ..
201: * .. Local Arrays ..
202: INTEGER IDUM( 1 )
203: DOUBLE PRECISION DIF( 2 )
204: * ..
205: * .. External Subroutines ..
206: EXTERNAL DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ, DLABAD,
207: $ DLACPY, DLASCL, DLASET, DORGQR, DORMQR, DTGSEN,
208: $ XERBLA
209: * ..
210: * .. External Functions ..
211: LOGICAL LSAME
212: INTEGER ILAENV
213: DOUBLE PRECISION DLAMCH, DLANGE
214: EXTERNAL LSAME, ILAENV, DLAMCH, DLANGE
215: * ..
216: * .. Intrinsic Functions ..
217: INTRINSIC ABS, MAX, SQRT
218: * ..
219: * .. Executable Statements ..
220: *
221: * Decode the input arguments
222: *
223: IF( LSAME( JOBVSL, 'N' ) ) THEN
224: IJOBVL = 1
225: ILVSL = .FALSE.
226: ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
227: IJOBVL = 2
228: ILVSL = .TRUE.
229: ELSE
230: IJOBVL = -1
231: ILVSL = .FALSE.
232: END IF
233: *
234: IF( LSAME( JOBVSR, 'N' ) ) THEN
235: IJOBVR = 1
236: ILVSR = .FALSE.
237: ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
238: IJOBVR = 2
239: ILVSR = .TRUE.
240: ELSE
241: IJOBVR = -1
242: ILVSR = .FALSE.
243: END IF
244: *
245: WANTST = LSAME( SORT, 'S' )
246: *
247: * Test the input arguments
248: *
249: INFO = 0
250: LQUERY = ( LWORK.EQ.-1 )
251: IF( IJOBVL.LE.0 ) THEN
252: INFO = -1
253: ELSE IF( IJOBVR.LE.0 ) THEN
254: INFO = -2
255: ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
256: INFO = -3
257: ELSE IF( N.LT.0 ) THEN
258: INFO = -5
259: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
260: INFO = -7
261: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
262: INFO = -9
263: ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
264: INFO = -15
265: ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
266: INFO = -17
267: END IF
268: *
269: * Compute workspace
270: * (Note: Comments in the code beginning "Workspace:" describe the
271: * minimal amount of workspace needed at that point in the code,
272: * as well as the preferred amount for good performance.
273: * NB refers to the optimal block size for the immediately
274: * following subroutine, as returned by ILAENV.)
275: *
276: IF( INFO.EQ.0 ) THEN
277: IF( N.GT.0 )THEN
278: MINWRK = MAX( 8*N, 6*N + 16 )
279: MAXWRK = MINWRK - N +
280: $ N*ILAENV( 1, 'DGEQRF', ' ', N, 1, N, 0 )
281: MAXWRK = MAX( MAXWRK, MINWRK - N +
282: $ N*ILAENV( 1, 'DORMQR', ' ', N, 1, N, -1 ) )
283: IF( ILVSL ) THEN
284: MAXWRK = MAX( MAXWRK, MINWRK - N +
285: $ N*ILAENV( 1, 'DORGQR', ' ', N, 1, N, -1 ) )
286: END IF
287: ELSE
288: MINWRK = 1
289: MAXWRK = 1
290: END IF
291: WORK( 1 ) = MAXWRK
292: *
293: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY )
294: $ INFO = -19
295: END IF
296: *
297: IF( INFO.NE.0 ) THEN
298: CALL XERBLA( 'DGGES ', -INFO )
299: RETURN
300: ELSE IF( LQUERY ) THEN
301: RETURN
302: END IF
303: *
304: * Quick return if possible
305: *
306: IF( N.EQ.0 ) THEN
307: SDIM = 0
308: RETURN
309: END IF
310: *
311: * Get machine constants
312: *
313: EPS = DLAMCH( 'P' )
314: SAFMIN = DLAMCH( 'S' )
315: SAFMAX = ONE / SAFMIN
316: CALL DLABAD( SAFMIN, SAFMAX )
317: SMLNUM = SQRT( SAFMIN ) / EPS
318: BIGNUM = ONE / SMLNUM
319: *
320: * Scale A if max element outside range [SMLNUM,BIGNUM]
321: *
322: ANRM = DLANGE( 'M', N, N, A, LDA, WORK )
323: ILASCL = .FALSE.
324: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
325: ANRMTO = SMLNUM
326: ILASCL = .TRUE.
327: ELSE IF( ANRM.GT.BIGNUM ) THEN
328: ANRMTO = BIGNUM
329: ILASCL = .TRUE.
330: END IF
331: IF( ILASCL )
332: $ CALL DLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
333: *
334: * Scale B if max element outside range [SMLNUM,BIGNUM]
335: *
336: BNRM = DLANGE( 'M', N, N, B, LDB, WORK )
337: ILBSCL = .FALSE.
338: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
339: BNRMTO = SMLNUM
340: ILBSCL = .TRUE.
341: ELSE IF( BNRM.GT.BIGNUM ) THEN
342: BNRMTO = BIGNUM
343: ILBSCL = .TRUE.
344: END IF
345: IF( ILBSCL )
346: $ CALL DLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
347: *
348: * Permute the matrix to make it more nearly triangular
349: * (Workspace: need 6*N + 2*N space for storing balancing factors)
350: *
351: ILEFT = 1
352: IRIGHT = N + 1
353: IWRK = IRIGHT + N
354: CALL DGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ),
355: $ WORK( IRIGHT ), WORK( IWRK ), IERR )
356: *
357: * Reduce B to triangular form (QR decomposition of B)
358: * (Workspace: need N, prefer N*NB)
359: *
360: IROWS = IHI + 1 - ILO
361: ICOLS = N + 1 - ILO
362: ITAU = IWRK
363: IWRK = ITAU + IROWS
364: CALL DGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
365: $ WORK( IWRK ), LWORK+1-IWRK, IERR )
366: *
367: * Apply the orthogonal transformation to matrix A
368: * (Workspace: need N, prefer N*NB)
369: *
370: CALL DORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
371: $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
372: $ LWORK+1-IWRK, IERR )
373: *
374: * Initialize VSL
375: * (Workspace: need N, prefer N*NB)
376: *
377: IF( ILVSL ) THEN
378: CALL DLASET( 'Full', N, N, ZERO, ONE, VSL, LDVSL )
379: IF( IROWS.GT.1 ) THEN
380: CALL DLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
381: $ VSL( ILO+1, ILO ), LDVSL )
382: END IF
383: CALL DORGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
384: $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
385: END IF
386: *
387: * Initialize VSR
388: *
389: IF( ILVSR )
390: $ CALL DLASET( 'Full', N, N, ZERO, ONE, VSR, LDVSR )
391: *
392: * Reduce to generalized Hessenberg form
393: * (Workspace: none needed)
394: *
395: CALL DGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
396: $ LDVSL, VSR, LDVSR, IERR )
397: *
398: * Perform QZ algorithm, computing Schur vectors if desired
399: * (Workspace: need N)
400: *
401: IWRK = ITAU
402: CALL DHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
403: $ ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR,
404: $ WORK( IWRK ), LWORK+1-IWRK, IERR )
405: IF( IERR.NE.0 ) THEN
406: IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
407: INFO = IERR
408: ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
409: INFO = IERR - N
410: ELSE
411: INFO = N + 1
412: END IF
413: GO TO 50
414: END IF
415: *
416: * Sort eigenvalues ALPHA/BETA if desired
417: * (Workspace: need 4*N+16 )
418: *
419: SDIM = 0
420: IF( WANTST ) THEN
421: *
422: * Undo scaling on eigenvalues before SELCTGing
423: *
424: IF( ILASCL ) THEN
425: CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N,
426: $ IERR )
427: CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N,
428: $ IERR )
429: END IF
430: IF( ILBSCL )
431: $ CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
432: *
433: * Select eigenvalues
434: *
435: DO 10 I = 1, N
436: BWORK( I ) = SELCTG( ALPHAR( I ), ALPHAI( I ), BETA( I ) )
437: 10 CONTINUE
438: *
439: CALL DTGSEN( 0, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB, ALPHAR,
440: $ ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, SDIM, PVSL,
441: $ PVSR, DIF, WORK( IWRK ), LWORK-IWRK+1, IDUM, 1,
442: $ IERR )
443: IF( IERR.EQ.1 )
444: $ INFO = N + 3
445: *
446: END IF
447: *
448: * Apply back-permutation to VSL and VSR
449: * (Workspace: none needed)
450: *
451: IF( ILVSL )
452: $ CALL DGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ),
453: $ WORK( IRIGHT ), N, VSL, LDVSL, IERR )
454: *
455: IF( ILVSR )
456: $ CALL DGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ),
457: $ WORK( IRIGHT ), N, VSR, LDVSR, IERR )
458: *
459: * Check if unscaling would cause over/underflow, if so, rescale
460: * (ALPHAR(I),ALPHAI(I),BETA(I)) so BETA(I) is on the order of
461: * B(I,I) and ALPHAR(I) and ALPHAI(I) are on the order of A(I,I)
462: *
463: IF( ILASCL ) THEN
464: DO 20 I = 1, N
465: IF( ALPHAI( I ).NE.ZERO ) THEN
466: IF( ( ALPHAR( I ) / SAFMAX ).GT.( ANRMTO / ANRM ) .OR.
467: $ ( SAFMIN / ALPHAR( I ) ).GT.( ANRM / ANRMTO ) ) THEN
468: WORK( 1 ) = ABS( A( I, I ) / ALPHAR( I ) )
469: BETA( I ) = BETA( I )*WORK( 1 )
470: ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
471: ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
472: ELSE IF( ( ALPHAI( I ) / SAFMAX ).GT.
473: $ ( ANRMTO / ANRM ) .OR.
474: $ ( SAFMIN / ALPHAI( I ) ).GT.( ANRM / ANRMTO ) )
475: $ THEN
476: WORK( 1 ) = ABS( A( I, I+1 ) / ALPHAI( I ) )
477: BETA( I ) = BETA( I )*WORK( 1 )
478: ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
479: ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
480: END IF
481: END IF
482: 20 CONTINUE
483: END IF
484: *
485: IF( ILBSCL ) THEN
486: DO 30 I = 1, N
487: IF( ALPHAI( I ).NE.ZERO ) THEN
488: IF( ( BETA( I ) / SAFMAX ).GT.( BNRMTO / BNRM ) .OR.
489: $ ( SAFMIN / BETA( I ) ).GT.( BNRM / BNRMTO ) ) THEN
490: WORK( 1 ) = ABS( B( I, I ) / BETA( I ) )
491: BETA( I ) = BETA( I )*WORK( 1 )
492: ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
493: ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
494: END IF
495: END IF
496: 30 CONTINUE
497: END IF
498: *
499: * Undo scaling
500: *
501: IF( ILASCL ) THEN
502: CALL DLASCL( 'H', 0, 0, ANRMTO, ANRM, N, N, A, LDA, IERR )
503: CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, IERR )
504: CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N, IERR )
505: END IF
506: *
507: IF( ILBSCL ) THEN
508: CALL DLASCL( 'U', 0, 0, BNRMTO, BNRM, N, N, B, LDB, IERR )
509: CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
510: END IF
511: *
512: IF( WANTST ) THEN
513: *
514: * Check if reordering is correct
515: *
516: LASTSL = .TRUE.
517: LST2SL = .TRUE.
518: SDIM = 0
519: IP = 0
520: DO 40 I = 1, N
521: CURSL = SELCTG( ALPHAR( I ), ALPHAI( I ), BETA( I ) )
522: IF( ALPHAI( I ).EQ.ZERO ) THEN
523: IF( CURSL )
524: $ SDIM = SDIM + 1
525: IP = 0
526: IF( CURSL .AND. .NOT.LASTSL )
527: $ INFO = N + 2
528: ELSE
529: IF( IP.EQ.1 ) THEN
530: *
531: * Last eigenvalue of conjugate pair
532: *
533: CURSL = CURSL .OR. LASTSL
534: LASTSL = CURSL
535: IF( CURSL )
536: $ SDIM = SDIM + 2
537: IP = -1
538: IF( CURSL .AND. .NOT.LST2SL )
539: $ INFO = N + 2
540: ELSE
541: *
542: * First eigenvalue of conjugate pair
543: *
544: IP = 1
545: END IF
546: END IF
547: LST2SL = LASTSL
548: LASTSL = CURSL
549: 40 CONTINUE
550: *
551: END IF
552: *
553: 50 CONTINUE
554: *
555: WORK( 1 ) = MAXWRK
556: *
557: RETURN
558: *
559: * End of DGGES
560: *
561: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgges.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack