Mercurial > octave-nkf

comparison libcruft/lapack/dtgevc.f @ 7034:68db500cb558

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-10-16 18:54:19 by jwe]
author jwe
date Tue, 16 Oct 2007 18:54:23 +0000
parents 15cddaacbc2d
children
comparison
equal deleted inserted replaced
7033:f0142f2afdc6 7034:68db500cb558
1 SUBROUTINE DTGEVC( SIDE, HOWMNY, SELECT, N, A, LDA, B, LDB, VL, 1 SUBROUTINE DTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
2 $ LDVL, VR, LDVR, MM, M, WORK, INFO ) 2 $ LDVL, VR, LDVR, MM, M, WORK, INFO )
3 * 3 *
4 * -- LAPACK routine (version 3.0) -- 4 * -- LAPACK routine (version 3.1) --
5 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 5 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
6 * Courant Institute, Argonne National Lab, and Rice University 6 * November 2006
7 * June 30, 1999
8 * 7 *
9 * .. Scalar Arguments .. 8 * .. Scalar Arguments ..
10 CHARACTER HOWMNY, SIDE 9 CHARACTER HOWMNY, SIDE
11 INTEGER INFO, LDA, LDB, LDVL, LDVR, M, MM, N 10 INTEGER INFO, LDP, LDS, LDVL, LDVR, M, MM, N
12 * .. 11 * ..
13 * .. Array Arguments .. 12 * .. Array Arguments ..
14 LOGICAL SELECT( * ) 13 LOGICAL SELECT( * )
15 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), VL( LDVL, * ), 14 DOUBLE PRECISION P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
16 $ VR( LDVR, * ), WORK( * ) 15 $ VR( LDVR, * ), WORK( * )
17 * .. 16 * ..
18 * 17 *
19 * 18 *
20 * Purpose 19 * Purpose
21 * ======= 20 * =======
22 * 21 *
23 * DTGEVC computes some or all of the right and/or left generalized 22 * DTGEVC computes some or all of the right and/or left eigenvectors of
24 * eigenvectors of a pair of real upper triangular matrices (A,B). 23 * a pair of real matrices (S,P), where S is a quasi-triangular matrix
25 * 24 * and P is upper triangular. Matrix pairs of this type are produced by
26 * The right generalized eigenvector x and the left generalized 25 * the generalized Schur factorization of a matrix pair (A,B):
27 * eigenvector y of (A,B) corresponding to a generalized eigenvalue 26 *
28 * w are defined by: 27 * A = Q*S*Z**T, B = Q*P*Z**T
29 * 28 *
30 * (A - wB) * x = 0 and y**H * (A - wB) = 0 29 * as computed by DGGHRD + DHGEQZ.
31 * 30 *
31 * The right eigenvector x and the left eigenvector y of (S,P)
32 * corresponding to an eigenvalue w are defined by:
33 *
34 * S*x = w*P*x, (y**H)*S = w*(y**H)*P,
35 *
32 * where y**H denotes the conjugate tranpose of y. 36 * where y**H denotes the conjugate tranpose of y.
33 * 37 * The eigenvalues are not input to this routine, but are computed
34 * If an eigenvalue w is determined by zero diagonal elements of both A 38 * directly from the diagonal blocks of S and P.
35 * and B, a unit vector is returned as the corresponding eigenvector. 39 *
36 * 40 * This routine returns the matrices X and/or Y of right and left
37 * If all eigenvectors are requested, the routine may either return 41 * eigenvectors of (S,P), or the products Z*X and/or Q*Y,
38 * the matrices X and/or Y of right or left eigenvectors of (A,B), or 42 * where Z and Q are input matrices.
39 * the products Z*X and/or Q*Y, where Z and Q are input orthogonal 43 * If Q and Z are the orthogonal factors from the generalized Schur
40 * matrices. If (A,B) was obtained from the generalized real-Schur 44 * factorization of a matrix pair (A,B), then Z*X and Q*Y
41 * factorization of an original pair of matrices 45 * are the matrices of right and left eigenvectors of (A,B).
42 * (A0,B0) = (Q*A*Z**H,Q*B*Z**H), 46 *
43 * then Z*X and Q*Y are the matrices of right or left eigenvectors of
44 * A.
45 *
46 * A must be block upper triangular, with 1-by-1 and 2-by-2 diagonal
47 * blocks. Corresponding to each 2-by-2 diagonal block is a complex
48 * conjugate pair of eigenvalues and eigenvectors; only one
49 * eigenvector of the pair is computed, namely the one corresponding
50 * to the eigenvalue with positive imaginary part.
51 *
52 * Arguments 47 * Arguments
53 * ========= 48 * =========
54 * 49 *
55 * SIDE (input) CHARACTER*1 50 * SIDE (input) CHARACTER*1
56 * = 'R': compute right eigenvectors only; 51 * = 'R': compute right eigenvectors only;
57 * = 'L': compute left eigenvectors only; 52 * = 'L': compute left eigenvectors only;
58 * = 'B': compute both right and left eigenvectors. 53 * = 'B': compute both right and left eigenvectors.
59 * 54 *
60 * HOWMNY (input) CHARACTER*1 55 * HOWMNY (input) CHARACTER*1
61 * = 'A': compute all right and/or left eigenvectors; 56 * = 'A': compute all right and/or left eigenvectors;
62 * = 'B': compute all right and/or left eigenvectors, and 57 * = 'B': compute all right and/or left eigenvectors,
63 * backtransform them using the input matrices supplied 58 * backtransformed by the matrices in VR and/or VL;
64 * in VR and/or VL;
65 * = 'S': compute selected right and/or left eigenvectors, 59 * = 'S': compute selected right and/or left eigenvectors,
66 * specified by the logical array SELECT. 60 * specified by the logical array SELECT.
67 * 61 *
68 * SELECT (input) LOGICAL array, dimension (N) 62 * SELECT (input) LOGICAL array, dimension (N)
69 * If HOWMNY='S', SELECT specifies the eigenvectors to be 63 * If HOWMNY='S', SELECT specifies the eigenvectors to be
70 * computed. 64 * computed. If w(j) is a real eigenvalue, the corresponding
71 * If HOWMNY='A' or 'B', SELECT is not referenced. 65 * real eigenvector is computed if SELECT(j) is .TRUE..
72 * To select the real eigenvector corresponding to the real 66 * If w(j) and w(j+1) are the real and imaginary parts of a
73 * eigenvalue w(j), SELECT(j) must be set to .TRUE. To select 67 * complex eigenvalue, the corresponding complex eigenvector
74 * the complex eigenvector corresponding to a complex conjugate 68 * is computed if either SELECT(j) or SELECT(j+1) is .TRUE.,
75 * pair w(j) and w(j+1), either SELECT(j) or SELECT(j+1) must 69 * and on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is
76 * be set to .TRUE.. 70 * set to .FALSE..
71 * Not referenced if HOWMNY = 'A' or 'B'.
77 * 72 *
78 * N (input) INTEGER 73 * N (input) INTEGER
79 * The order of the matrices A and B. N >= 0. 74 * The order of the matrices S and P. N >= 0.
80 * 75 *
81 * A (input) DOUBLE PRECISION array, dimension (LDA,N) 76 * S (input) DOUBLE PRECISION array, dimension (LDS,N)
82 * The upper quasi-triangular matrix A. 77 * The upper quasi-triangular matrix S from a generalized Schur
83 * 78 * factorization, as computed by DHGEQZ.
84 * LDA (input) INTEGER 79 *
85 * The leading dimension of array A. LDA >= max(1, N). 80 * LDS (input) INTEGER
86 * 81 * The leading dimension of array S. LDS >= max(1,N).
87 * B (input) DOUBLE PRECISION array, dimension (LDB,N) 82 *
88 * The upper triangular matrix B. If A has a 2-by-2 diagonal 83 * P (input) DOUBLE PRECISION array, dimension (LDP,N)
89 * block, then the corresponding 2-by-2 block of B must be 84 * The upper triangular matrix P from a generalized Schur
90 * diagonal with positive elements. 85 * factorization, as computed by DHGEQZ.
91 * 86 * 2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks
92 * LDB (input) INTEGER 87 * of S must be in positive diagonal form.
93 * The leading dimension of array B. LDB >= max(1,N). 88 *
89 * LDP (input) INTEGER
90 * The leading dimension of array P. LDP >= max(1,N).
94 * 91 *
95 * VL (input/output) DOUBLE PRECISION array, dimension (LDVL,MM) 92 * VL (input/output) DOUBLE PRECISION array, dimension (LDVL,MM)
96 * On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must 93 * On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
97 * contain an N-by-N matrix Q (usually the orthogonal matrix Q 94 * contain an N-by-N matrix Q (usually the orthogonal matrix Q
98 * of left Schur vectors returned by DHGEQZ). 95 * of left Schur vectors returned by DHGEQZ).
99 * On exit, if SIDE = 'L' or 'B', VL contains: 96 * On exit, if SIDE = 'L' or 'B', VL contains:
100 * if HOWMNY = 'A', the matrix Y of left eigenvectors of (A,B); 97 * if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P);
101 * if HOWMNY = 'B', the matrix Q*Y; 98 * if HOWMNY = 'B', the matrix Q*Y;
102 * if HOWMNY = 'S', the left eigenvectors of (A,B) specified by 99 * if HOWMNY = 'S', the left eigenvectors of (S,P) specified by
103 * SELECT, stored consecutively in the columns of 100 * SELECT, stored consecutively in the columns of
104 * VL, in the same order as their eigenvalues. 101 * VL, in the same order as their eigenvalues.
105 * If SIDE = 'R', VL is not referenced.
106 * 102 *
107 * A complex eigenvector corresponding to a complex eigenvalue 103 * A complex eigenvector corresponding to a complex eigenvalue
108 * is stored in two consecutive columns, the first holding the 104 * is stored in two consecutive columns, the first holding the
109 * real part, and the second the imaginary part. 105 * real part, and the second the imaginary part.
110 * 106 *
107 * Not referenced if SIDE = 'R'.
108 *
111 * LDVL (input) INTEGER 109 * LDVL (input) INTEGER
112 * The leading dimension of array VL. 110 * The leading dimension of array VL. LDVL >= 1, and if
113 * LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise. 111 * SIDE = 'L' or 'B', LDVL >= N.
114 * 112 *
115 * VR (input/output) DOUBLE PRECISION array, dimension (LDVR,MM) 113 * VR (input/output) DOUBLE PRECISION array, dimension (LDVR,MM)
116 * On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must 114 * On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
117 * contain an N-by-N matrix Q (usually the orthogonal matrix Z 115 * contain an N-by-N matrix Z (usually the orthogonal matrix Z
118 * of right Schur vectors returned by DHGEQZ). 116 * of right Schur vectors returned by DHGEQZ).
117 *
119 * On exit, if SIDE = 'R' or 'B', VR contains: 118 * On exit, if SIDE = 'R' or 'B', VR contains:
120 * if HOWMNY = 'A', the matrix X of right eigenvectors of (A,B); 119 * if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P);
121 * if HOWMNY = 'B', the matrix Z*X; 120 * if HOWMNY = 'B' or 'b', the matrix Z*X;
122 * if HOWMNY = 'S', the right eigenvectors of (A,B) specified by 121 * if HOWMNY = 'S' or 's', the right eigenvectors of (S,P)
123 * SELECT, stored consecutively in the columns of 122 * specified by SELECT, stored consecutively in the
124 * VR, in the same order as their eigenvalues. 123 * columns of VR, in the same order as their
125 * If SIDE = 'L', VR is not referenced. 124 * eigenvalues.
126 * 125 *
127 * A complex eigenvector corresponding to a complex eigenvalue 126 * A complex eigenvector corresponding to a complex eigenvalue
128 * is stored in two consecutive columns, the first holding the 127 * is stored in two consecutive columns, the first holding the
129 * real part and the second the imaginary part. 128 * real part and the second the imaginary part.
129 *
130 * Not referenced if SIDE = 'L'.
130 * 131 *
131 * LDVR (input) INTEGER 132 * LDVR (input) INTEGER
132 * The leading dimension of the array VR. 133 * The leading dimension of the array VR. LDVR >= 1, and if
133 * LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise. 134 * SIDE = 'R' or 'B', LDVR >= N.
134 * 135 *
135 * MM (input) INTEGER 136 * MM (input) INTEGER
136 * The number of columns in the arrays VL and/or VR. MM >= M. 137 * The number of columns in the arrays VL and/or VR. MM >= M.
137 * 138 *
138 * M (output) INTEGER 139 * M (output) INTEGER
197 * for all the rows in parallel. As each v(j) is computed, the 198 * for all the rows in parallel. As each v(j) is computed, the
198 * contribution of v(j) times the j-th column of C is added to the 199 * contribution of v(j) times the j-th column of C is added to the
199 * partial sums. Since FORTRAN arrays are stored columnwise, this has 200 * partial sums. Since FORTRAN arrays are stored columnwise, this has
200 * the advantage that at each step, the elements of C that are accessed 201 * the advantage that at each step, the elements of C that are accessed
201 * are adjacent to one another, whereas with the rowwise method, the 202 * are adjacent to one another, whereas with the rowwise method, the
202 * elements accessed at a step are spaced LDA (and LDB) words apart. 203 * elements accessed at a step are spaced LDS (and LDP) words apart.
203 * 204 *
204 * When finding left eigenvectors, the matrix in question is the 205 * When finding left eigenvectors, the matrix in question is the
205 * transpose of the one in storage, so the rowwise method then 206 * transpose of the one in storage, so the rowwise method then
206 * actually accesses columns of A and B at each step, and so is the 207 * actually accesses columns of A and B at each step, and so is the
207 * preferred method. 208 * preferred method.
224 $ CREALB, DMIN, SAFMIN, SALFAR, SBETA, SCALE, 225 $ CREALB, DMIN, SAFMIN, SALFAR, SBETA, SCALE,
225 $ SMALL, TEMP, TEMP2, TEMP2I, TEMP2R, ULP, XMAX, 226 $ SMALL, TEMP, TEMP2, TEMP2I, TEMP2R, ULP, XMAX,
226 $ XSCALE 227 $ XSCALE
227 * .. 228 * ..
228 * .. Local Arrays .. 229 * .. Local Arrays ..
229 DOUBLE PRECISION BDIAG( 2 ), SUM( 2, 2 ), SUMA( 2, 2 ), 230 DOUBLE PRECISION BDIAG( 2 ), SUM( 2, 2 ), SUMS( 2, 2 ),
230 $ SUMB( 2, 2 ) 231 $ SUMP( 2, 2 )
231 * .. 232 * ..
232 * .. External Functions .. 233 * .. External Functions ..
233 LOGICAL LSAME 234 LOGICAL LSAME
234 DOUBLE PRECISION DLAMCH 235 DOUBLE PRECISION DLAMCH
235 EXTERNAL LSAME, DLAMCH 236 EXTERNAL LSAME, DLAMCH
236 * .. 237 * ..
237 * .. External Subroutines .. 238 * .. External Subroutines ..
238 EXTERNAL DGEMV, DLACPY, DLAG2, DLALN2, XERBLA 239 EXTERNAL DGEMV, DLABAD, DLACPY, DLAG2, DLALN2, XERBLA
239 * .. 240 * ..
240 * .. Intrinsic Functions .. 241 * .. Intrinsic Functions ..
241 INTRINSIC ABS, MAX, MIN 242 INTRINSIC ABS, MAX, MIN
242 * .. 243 * ..
243 * .. Executable Statements .. 244 * .. Executable Statements ..
250 ILBACK = .FALSE. 251 ILBACK = .FALSE.
251 ELSE IF( LSAME( HOWMNY, 'S' ) ) THEN 252 ELSE IF( LSAME( HOWMNY, 'S' ) ) THEN
252 IHWMNY = 2 253 IHWMNY = 2
253 ILALL = .FALSE. 254 ILALL = .FALSE.
254 ILBACK = .FALSE. 255 ILBACK = .FALSE.
255 ELSE IF( LSAME( HOWMNY, 'B' ) .OR. LSAME( HOWMNY, 'T' ) ) THEN 256 ELSE IF( LSAME( HOWMNY, 'B' ) ) THEN
256 IHWMNY = 3 257 IHWMNY = 3
257 ILALL = .TRUE. 258 ILALL = .TRUE.
258 ILBACK = .TRUE. 259 ILBACK = .TRUE.
259 ELSE 260 ELSE
260 IHWMNY = -1 261 IHWMNY = -1
282 INFO = -1 283 INFO = -1
283 ELSE IF( IHWMNY.LT.0 ) THEN 284 ELSE IF( IHWMNY.LT.0 ) THEN
284 INFO = -2 285 INFO = -2
285 ELSE IF( N.LT.0 ) THEN 286 ELSE IF( N.LT.0 ) THEN
286 INFO = -4 287 INFO = -4
287 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 288 ELSE IF( LDS.LT.MAX( 1, N ) ) THEN
288 INFO = -6 289 INFO = -6
289 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 290 ELSE IF( LDP.LT.MAX( 1, N ) ) THEN
290 INFO = -8 291 INFO = -8
291 END IF 292 END IF
292 IF( INFO.NE.0 ) THEN 293 IF( INFO.NE.0 ) THEN
293 CALL XERBLA( 'DTGEVC', -INFO ) 294 CALL XERBLA( 'DTGEVC', -INFO )
294 RETURN 295 RETURN
303 IF( ILCPLX ) THEN 304 IF( ILCPLX ) THEN
304 ILCPLX = .FALSE. 305 ILCPLX = .FALSE.
305 GO TO 10 306 GO TO 10
306 END IF 307 END IF
307 IF( J.LT.N ) THEN 308 IF( J.LT.N ) THEN
308 IF( A( J+1, J ).NE.ZERO ) 309 IF( S( J+1, J ).NE.ZERO )
309 $ ILCPLX = .TRUE. 310 $ ILCPLX = .TRUE.
310 END IF 311 END IF
311 IF( ILCPLX ) THEN 312 IF( ILCPLX ) THEN
312 IF( SELECT( J ) .OR. SELECT( J+1 ) ) 313 IF( SELECT( J ) .OR. SELECT( J+1 ) )
313 $ IM = IM + 2 314 $ IM = IM + 2
323 * Check 2-by-2 diagonal blocks of A, B 324 * Check 2-by-2 diagonal blocks of A, B
324 * 325 *
325 ILABAD = .FALSE. 326 ILABAD = .FALSE.
326 ILBBAD = .FALSE. 327 ILBBAD = .FALSE.
327 DO 20 J = 1, N - 1 328 DO 20 J = 1, N - 1
328 IF( A( J+1, J ).NE.ZERO ) THEN 329 IF( S( J+1, J ).NE.ZERO ) THEN
329 IF( B( J, J ).EQ.ZERO .OR. B( J+1, J+1 ).EQ.ZERO .OR. 330 IF( P( J, J ).EQ.ZERO .OR. P( J+1, J+1 ).EQ.ZERO .OR.
330 $ B( J, J+1 ).NE.ZERO )ILBBAD = .TRUE. 331 $ P( J, J+1 ).NE.ZERO )ILBBAD = .TRUE.
331 IF( J.LT.N-1 ) THEN 332 IF( J.LT.N-1 ) THEN
332 IF( A( J+2, J+1 ).NE.ZERO ) 333 IF( S( J+2, J+1 ).NE.ZERO )
333 $ ILABAD = .TRUE. 334 $ ILABAD = .TRUE.
334 END IF 335 END IF
335 END IF 336 END IF
336 20 CONTINUE 337 20 CONTINUE
337 * 338 *
370 * Compute the 1-norm of each column of the strictly upper triangular 371 * Compute the 1-norm of each column of the strictly upper triangular
371 * part (i.e., excluding all elements belonging to the diagonal 372 * part (i.e., excluding all elements belonging to the diagonal
372 * blocks) of A and B to check for possible overflow in the 373 * blocks) of A and B to check for possible overflow in the
373 * triangular solver. 374 * triangular solver.
374 * 375 *
375 ANORM = ABS( A( 1, 1 ) ) 376 ANORM = ABS( S( 1, 1 ) )
376 IF( N.GT.1 ) 377 IF( N.GT.1 )
377 $ ANORM = ANORM + ABS( A( 2, 1 ) ) 378 $ ANORM = ANORM + ABS( S( 2, 1 ) )
378 BNORM = ABS( B( 1, 1 ) ) 379 BNORM = ABS( P( 1, 1 ) )
379 WORK( 1 ) = ZERO 380 WORK( 1 ) = ZERO
380 WORK( N+1 ) = ZERO 381 WORK( N+1 ) = ZERO
381 * 382 *
382 DO 50 J = 2, N 383 DO 50 J = 2, N
383 TEMP = ZERO 384 TEMP = ZERO
384 TEMP2 = ZERO 385 TEMP2 = ZERO
385 IF( A( J, J-1 ).EQ.ZERO ) THEN 386 IF( S( J, J-1 ).EQ.ZERO ) THEN
386 IEND = J - 1 387 IEND = J - 1
387 ELSE 388 ELSE
388 IEND = J - 2 389 IEND = J - 2
389 END IF 390 END IF
390 DO 30 I = 1, IEND 391 DO 30 I = 1, IEND
391 TEMP = TEMP + ABS( A( I, J ) ) 392 TEMP = TEMP + ABS( S( I, J ) )
392 TEMP2 = TEMP2 + ABS( B( I, J ) ) 393 TEMP2 = TEMP2 + ABS( P( I, J ) )
393 30 CONTINUE 394 30 CONTINUE
394 WORK( J ) = TEMP 395 WORK( J ) = TEMP
395 WORK( N+J ) = TEMP2 396 WORK( N+J ) = TEMP2
396 DO 40 I = IEND + 1, MIN( J+1, N ) 397 DO 40 I = IEND + 1, MIN( J+1, N )
397 TEMP = TEMP + ABS( A( I, J ) ) 398 TEMP = TEMP + ABS( S( I, J ) )
398 TEMP2 = TEMP2 + ABS( B( I, J ) ) 399 TEMP2 = TEMP2 + ABS( P( I, J ) )
399 40 CONTINUE 400 40 CONTINUE
400 ANORM = MAX( ANORM, TEMP ) 401 ANORM = MAX( ANORM, TEMP )
401 BNORM = MAX( BNORM, TEMP2 ) 402 BNORM = MAX( BNORM, TEMP2 )
402 50 CONTINUE 403 50 CONTINUE
403 * 404 *
423 ILCPLX = .FALSE. 424 ILCPLX = .FALSE.
424 GO TO 220 425 GO TO 220
425 END IF 426 END IF
426 NW = 1 427 NW = 1
427 IF( JE.LT.N ) THEN 428 IF( JE.LT.N ) THEN
428 IF( A( JE+1, JE ).NE.ZERO ) THEN 429 IF( S( JE+1, JE ).NE.ZERO ) THEN
429 ILCPLX = .TRUE. 430 ILCPLX = .TRUE.
430 NW = 2 431 NW = 2
431 END IF 432 END IF
432 END IF 433 END IF
433 IF( ILALL ) THEN 434 IF( ILALL ) THEN
442 * 443 *
443 * Decide if (a) singular pencil, (b) real eigenvalue, or 444 * Decide if (a) singular pencil, (b) real eigenvalue, or
444 * (c) complex eigenvalue. 445 * (c) complex eigenvalue.
445 * 446 *
446 IF( .NOT.ILCPLX ) THEN 447 IF( .NOT.ILCPLX ) THEN
447 IF( ABS( A( JE, JE ) ).LE.SAFMIN .AND. 448 IF( ABS( S( JE, JE ) ).LE.SAFMIN .AND.
448 $ ABS( B( JE, JE ) ).LE.SAFMIN ) THEN 449 $ ABS( P( JE, JE ) ).LE.SAFMIN ) THEN
449 * 450 *
450 * Singular matrix pencil -- return unit eigenvector 451 * Singular matrix pencil -- return unit eigenvector
451 * 452 *
452 IEIG = IEIG + 1 453 IEIG = IEIG + 1
453 DO 60 JR = 1, N 454 DO 60 JR = 1, N
470 * 471 *
471 IF( .NOT.ILCPLX ) THEN 472 IF( .NOT.ILCPLX ) THEN
472 * 473 *
473 * Real eigenvalue 474 * Real eigenvalue
474 * 475 *
475 TEMP = ONE / MAX( ABS( A( JE, JE ) )*ASCALE, 476 TEMP = ONE / MAX( ABS( S( JE, JE ) )*ASCALE,
476 $ ABS( B( JE, JE ) )*BSCALE, SAFMIN ) 477 $ ABS( P( JE, JE ) )*BSCALE, SAFMIN )
477 SALFAR = ( TEMP*A( JE, JE ) )*ASCALE 478 SALFAR = ( TEMP*S( JE, JE ) )*ASCALE
478 SBETA = ( TEMP*B( JE, JE ) )*BSCALE 479 SBETA = ( TEMP*P( JE, JE ) )*BSCALE
479 ACOEF = SBETA*ASCALE 480 ACOEF = SBETA*ASCALE
480 BCOEFR = SALFAR*BSCALE 481 BCOEFR = SALFAR*BSCALE
481 BCOEFI = ZERO 482 BCOEFI = ZERO
482 * 483 *
483 * Scale to avoid underflow 484 * Scale to avoid underflow
515 XMAX = ONE 516 XMAX = ONE
516 ELSE 517 ELSE
517 * 518 *
518 * Complex eigenvalue 519 * Complex eigenvalue
519 * 520 *
520 CALL DLAG2( A( JE, JE ), LDA, B( JE, JE ), LDB, 521 CALL DLAG2( S( JE, JE ), LDS, P( JE, JE ), LDP,
521 $ SAFMIN*SAFETY, ACOEF, TEMP, BCOEFR, TEMP2, 522 $ SAFMIN*SAFETY, ACOEF, TEMP, BCOEFR, TEMP2,
522 $ BCOEFI ) 523 $ BCOEFI )
523 BCOEFI = -BCOEFI 524 BCOEFI = -BCOEFI
524 IF( BCOEFI.EQ.ZERO ) THEN 525 IF( BCOEFI.EQ.ZERO ) THEN
525 INFO = JE 526 INFO = JE
547 BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI ) 548 BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI )
548 END IF 549 END IF
549 * 550 *
550 * Compute first two components of eigenvector 551 * Compute first two components of eigenvector
551 * 552 *
552 TEMP = ACOEF*A( JE+1, JE ) 553 TEMP = ACOEF*S( JE+1, JE )
553 TEMP2R = ACOEF*A( JE, JE ) - BCOEFR*B( JE, JE ) 554 TEMP2R = ACOEF*S( JE, JE ) - BCOEFR*P( JE, JE )
554 TEMP2I = -BCOEFI*B( JE, JE ) 555 TEMP2I = -BCOEFI*P( JE, JE )
555 IF( ABS( TEMP ).GT.ABS( TEMP2R )+ABS( TEMP2I ) ) THEN 556 IF( ABS( TEMP ).GT.ABS( TEMP2R )+ABS( TEMP2I ) ) THEN
556 WORK( 2*N+JE ) = ONE 557 WORK( 2*N+JE ) = ONE
557 WORK( 3*N+JE ) = ZERO 558 WORK( 3*N+JE ) = ZERO
558 WORK( 2*N+JE+1 ) = -TEMP2R / TEMP 559 WORK( 2*N+JE+1 ) = -TEMP2R / TEMP
559 WORK( 3*N+JE+1 ) = -TEMP2I / TEMP 560 WORK( 3*N+JE+1 ) = -TEMP2I / TEMP
560 ELSE 561 ELSE
561 WORK( 2*N+JE+1 ) = ONE 562 WORK( 2*N+JE+1 ) = ONE
562 WORK( 3*N+JE+1 ) = ZERO 563 WORK( 3*N+JE+1 ) = ZERO
563 TEMP = ACOEF*A( JE, JE+1 ) 564 TEMP = ACOEF*S( JE, JE+1 )
564 WORK( 2*N+JE ) = ( BCOEFR*B( JE+1, JE+1 )-ACOEF* 565 WORK( 2*N+JE ) = ( BCOEFR*P( JE+1, JE+1 )-ACOEF*
565 $ A( JE+1, JE+1 ) ) / TEMP 566 $ S( JE+1, JE+1 ) ) / TEMP
566 WORK( 3*N+JE ) = BCOEFI*B( JE+1, JE+1 ) / TEMP 567 WORK( 3*N+JE ) = BCOEFI*P( JE+1, JE+1 ) / TEMP
567 END IF 568 END IF
568 XMAX = MAX( ABS( WORK( 2*N+JE ) )+ABS( WORK( 3*N+JE ) ), 569 XMAX = MAX( ABS( WORK( 2*N+JE ) )+ABS( WORK( 3*N+JE ) ),
569 $ ABS( WORK( 2*N+JE+1 ) )+ABS( WORK( 3*N+JE+1 ) ) ) 570 $ ABS( WORK( 2*N+JE+1 ) )+ABS( WORK( 3*N+JE+1 ) ) )
570 END IF 571 END IF
571 * 572 *
584 IL2BY2 = .FALSE. 585 IL2BY2 = .FALSE.
585 GO TO 160 586 GO TO 160
586 END IF 587 END IF
587 * 588 *
588 NA = 1 589 NA = 1
589 BDIAG( 1 ) = B( J, J ) 590 BDIAG( 1 ) = P( J, J )
590 IF( J.LT.N ) THEN 591 IF( J.LT.N ) THEN
591 IF( A( J+1, J ).NE.ZERO ) THEN 592 IF( S( J+1, J ).NE.ZERO ) THEN
592 IL2BY2 = .TRUE. 593 IL2BY2 = .TRUE.
593 BDIAG( 2 ) = B( J+1, J+1 ) 594 BDIAG( 2 ) = P( J+1, J+1 )
594 NA = 2 595 NA = 2
595 END IF 596 END IF
596 END IF 597 END IF
597 * 598 *
598 * Check whether scaling is necessary for dot products 599 * Check whether scaling is necessary for dot products
614 END IF 615 END IF
615 * 616 *
616 * Compute dot products 617 * Compute dot products
617 * 618 *
618 * j-1 619 * j-1
619 * SUM = sum conjg( a*A(k,j) - b*B(k,j) )*x(k) 620 * SUM = sum conjg( a*S(k,j) - b*P(k,j) )*x(k)
620 * k=je 621 * k=je
621 * 622 *
622 * To reduce the op count, this is done as 623 * To reduce the op count, this is done as
623 * 624 *
624 * _ j-1 _ j-1 625 * _ j-1 _ j-1
625 * a*conjg( sum A(k,j)*x(k) ) - b*conjg( sum B(k,j)*x(k) ) 626 * a*conjg( sum S(k,j)*x(k) ) - b*conjg( sum P(k,j)*x(k) )
626 * k=je k=je 627 * k=je k=je
627 * 628 *
628 * which may cause underflow problems if A or B are close 629 * which may cause underflow problems if A or B are close
629 * to underflow. (E.g., less than SMALL.) 630 * to underflow. (E.g., less than SMALL.)
630 * 631 *
657 *VOCL LOOP,SCALAR 658 *VOCL LOOP,SCALAR
658 CIBM PREFER SCALAR 659 CIBM PREFER SCALAR
659 *$PL$ CMCHAR='*' 660 *$PL$ CMCHAR='*'
660 * 661 *
661 DO 110 JA = 1, NA 662 DO 110 JA = 1, NA
662 SUMA( JA, JW ) = ZERO 663 SUMS( JA, JW ) = ZERO
663 SUMB( JA, JW ) = ZERO 664 SUMP( JA, JW ) = ZERO
664 * 665 *
665 DO 100 JR = JE, J - 1 666 DO 100 JR = JE, J - 1
666 SUMA( JA, JW ) = SUMA( JA, JW ) + 667 SUMS( JA, JW ) = SUMS( JA, JW ) +
667 $ A( JR, J+JA-1 )* 668 $ S( JR, J+JA-1 )*
668 $ WORK( ( JW+1 )*N+JR ) 669 $ WORK( ( JW+1 )*N+JR )
669 SUMB( JA, JW ) = SUMB( JA, JW ) + 670 SUMP( JA, JW ) = SUMP( JA, JW ) +
670 $ B( JR, J+JA-1 )* 671 $ P( JR, J+JA-1 )*
671 $ WORK( ( JW+1 )*N+JR ) 672 $ WORK( ( JW+1 )*N+JR )
672 100 CONTINUE 673 100 CONTINUE
673 110 CONTINUE 674 110 CONTINUE
674 120 CONTINUE 675 120 CONTINUE
675 * 676 *
685 CIBM PREFER SCALAR 686 CIBM PREFER SCALAR
686 *$PL$ CMCHAR='*' 687 *$PL$ CMCHAR='*'
687 * 688 *
688 DO 130 JA = 1, NA 689 DO 130 JA = 1, NA
689 IF( ILCPLX ) THEN 690 IF( ILCPLX ) THEN
690 SUM( JA, 1 ) = -ACOEF*SUMA( JA, 1 ) + 691 SUM( JA, 1 ) = -ACOEF*SUMS( JA, 1 ) +
691 $ BCOEFR*SUMB( JA, 1 ) - 692 $ BCOEFR*SUMP( JA, 1 ) -
692 $ BCOEFI*SUMB( JA, 2 ) 693 $ BCOEFI*SUMP( JA, 2 )
693 SUM( JA, 2 ) = -ACOEF*SUMA( JA, 2 ) + 694 SUM( JA, 2 ) = -ACOEF*SUMS( JA, 2 ) +
694 $ BCOEFR*SUMB( JA, 2 ) + 695 $ BCOEFR*SUMP( JA, 2 ) +
695 $ BCOEFI*SUMB( JA, 1 ) 696 $ BCOEFI*SUMP( JA, 1 )
696 ELSE 697 ELSE
697 SUM( JA, 1 ) = -ACOEF*SUMA( JA, 1 ) + 698 SUM( JA, 1 ) = -ACOEF*SUMS( JA, 1 ) +
698 $ BCOEFR*SUMB( JA, 1 ) 699 $ BCOEFR*SUMP( JA, 1 )
699 END IF 700 END IF
700 130 CONTINUE 701 130 CONTINUE
701 * 702 *
702 * T 703 * T
703 * Solve ( a A - b B ) y = SUM(,) 704 * Solve ( a A - b B ) y = SUM(,)
704 * with scaling and perturbation of the denominator 705 * with scaling and perturbation of the denominator
705 * 706 *
706 CALL DLALN2( .TRUE., NA, NW, DMIN, ACOEF, A( J, J ), LDA, 707 CALL DLALN2( .TRUE., NA, NW, DMIN, ACOEF, S( J, J ), LDS,
707 $ BDIAG( 1 ), BDIAG( 2 ), SUM, 2, BCOEFR, 708 $ BDIAG( 1 ), BDIAG( 2 ), SUM, 2, BCOEFR,
708 $ BCOEFI, WORK( 2*N+J ), N, SCALE, TEMP, 709 $ BCOEFI, WORK( 2*N+J ), N, SCALE, TEMP,
709 $ IINFO ) 710 $ IINFO )
710 IF( SCALE.LT.ONE ) THEN 711 IF( SCALE.LT.ONE ) THEN
711 DO 150 JW = 0, NW - 1 712 DO 150 JW = 0, NW - 1
788 ILCPLX = .FALSE. 789 ILCPLX = .FALSE.
789 GO TO 500 790 GO TO 500
790 END IF 791 END IF
791 NW = 1 792 NW = 1
792 IF( JE.GT.1 ) THEN 793 IF( JE.GT.1 ) THEN
793 IF( A( JE, JE-1 ).NE.ZERO ) THEN 794 IF( S( JE, JE-1 ).NE.ZERO ) THEN
794 ILCPLX = .TRUE. 795 ILCPLX = .TRUE.
795 NW = 2 796 NW = 2
796 END IF 797 END IF
797 END IF 798 END IF
798 IF( ILALL ) THEN 799 IF( ILALL ) THEN
807 * 808 *
808 * Decide if (a) singular pencil, (b) real eigenvalue, or 809 * Decide if (a) singular pencil, (b) real eigenvalue, or
809 * (c) complex eigenvalue. 810 * (c) complex eigenvalue.
810 * 811 *
811 IF( .NOT.ILCPLX ) THEN 812 IF( .NOT.ILCPLX ) THEN
812 IF( ABS( A( JE, JE ) ).LE.SAFMIN .AND. 813 IF( ABS( S( JE, JE ) ).LE.SAFMIN .AND.
813 $ ABS( B( JE, JE ) ).LE.SAFMIN ) THEN 814 $ ABS( P( JE, JE ) ).LE.SAFMIN ) THEN
814 * 815 *
815 * Singular matrix pencil -- unit eigenvector 816 * Singular matrix pencil -- unit eigenvector
816 * 817 *
817 IEIG = IEIG - 1 818 IEIG = IEIG - 1
818 DO 230 JR = 1, N 819 DO 230 JR = 1, N
837 * 838 *
838 IF( .NOT.ILCPLX ) THEN 839 IF( .NOT.ILCPLX ) THEN
839 * 840 *
840 * Real eigenvalue 841 * Real eigenvalue
841 * 842 *
842 TEMP = ONE / MAX( ABS( A( JE, JE ) )*ASCALE, 843 TEMP = ONE / MAX( ABS( S( JE, JE ) )*ASCALE,
843 $ ABS( B( JE, JE ) )*BSCALE, SAFMIN ) 844 $ ABS( P( JE, JE ) )*BSCALE, SAFMIN )
844 SALFAR = ( TEMP*A( JE, JE ) )*ASCALE 845 SALFAR = ( TEMP*S( JE, JE ) )*ASCALE
845 SBETA = ( TEMP*B( JE, JE ) )*BSCALE 846 SBETA = ( TEMP*P( JE, JE ) )*BSCALE
846 ACOEF = SBETA*ASCALE 847 ACOEF = SBETA*ASCALE
847 BCOEFR = SALFAR*BSCALE 848 BCOEFR = SALFAR*BSCALE
848 BCOEFI = ZERO 849 BCOEFI = ZERO
849 * 850 *
850 * Scale to avoid underflow 851 * Scale to avoid underflow
883 * 884 *
884 * Compute contribution from column JE of A and B to sum 885 * Compute contribution from column JE of A and B to sum
885 * (See "Further Details", above.) 886 * (See "Further Details", above.)
886 * 887 *
887 DO 260 JR = 1, JE - 1 888 DO 260 JR = 1, JE - 1
888 WORK( 2*N+JR ) = BCOEFR*B( JR, JE ) - 889 WORK( 2*N+JR ) = BCOEFR*P( JR, JE ) -
889 $ ACOEF*A( JR, JE ) 890 $ ACOEF*S( JR, JE )
890 260 CONTINUE 891 260 CONTINUE
891 ELSE 892 ELSE
892 * 893 *
893 * Complex eigenvalue 894 * Complex eigenvalue
894 * 895 *
895 CALL DLAG2( A( JE-1, JE-1 ), LDA, B( JE-1, JE-1 ), LDB, 896 CALL DLAG2( S( JE-1, JE-1 ), LDS, P( JE-1, JE-1 ), LDP,
896 $ SAFMIN*SAFETY, ACOEF, TEMP, BCOEFR, TEMP2, 897 $ SAFMIN*SAFETY, ACOEF, TEMP, BCOEFR, TEMP2,
897 $ BCOEFI ) 898 $ BCOEFI )
898 IF( BCOEFI.EQ.ZERO ) THEN 899 IF( BCOEFI.EQ.ZERO ) THEN
899 INFO = JE - 1 900 INFO = JE - 1
900 RETURN 901 RETURN
922 END IF 923 END IF
923 * 924 *
924 * Compute first two components of eigenvector 925 * Compute first two components of eigenvector
925 * and contribution to sums 926 * and contribution to sums
926 * 927 *
927 TEMP = ACOEF*A( JE, JE-1 ) 928 TEMP = ACOEF*S( JE, JE-1 )
928 TEMP2R = ACOEF*A( JE, JE ) - BCOEFR*B( JE, JE ) 929 TEMP2R = ACOEF*S( JE, JE ) - BCOEFR*P( JE, JE )
929 TEMP2I = -BCOEFI*B( JE, JE ) 930 TEMP2I = -BCOEFI*P( JE, JE )
930 IF( ABS( TEMP ).GE.ABS( TEMP2R )+ABS( TEMP2I ) ) THEN 931 IF( ABS( TEMP ).GE.ABS( TEMP2R )+ABS( TEMP2I ) ) THEN
931 WORK( 2*N+JE ) = ONE 932 WORK( 2*N+JE ) = ONE
932 WORK( 3*N+JE ) = ZERO 933 WORK( 3*N+JE ) = ZERO
933 WORK( 2*N+JE-1 ) = -TEMP2R / TEMP 934 WORK( 2*N+JE-1 ) = -TEMP2R / TEMP
934 WORK( 3*N+JE-1 ) = -TEMP2I / TEMP 935 WORK( 3*N+JE-1 ) = -TEMP2I / TEMP
935 ELSE 936 ELSE
936 WORK( 2*N+JE-1 ) = ONE 937 WORK( 2*N+JE-1 ) = ONE
937 WORK( 3*N+JE-1 ) = ZERO 938 WORK( 3*N+JE-1 ) = ZERO
938 TEMP = ACOEF*A( JE-1, JE ) 939 TEMP = ACOEF*S( JE-1, JE )
939 WORK( 2*N+JE ) = ( BCOEFR*B( JE-1, JE-1 )-ACOEF* 940 WORK( 2*N+JE ) = ( BCOEFR*P( JE-1, JE-1 )-ACOEF*
940 $ A( JE-1, JE-1 ) ) / TEMP 941 $ S( JE-1, JE-1 ) ) / TEMP
941 WORK( 3*N+JE ) = BCOEFI*B( JE-1, JE-1 ) / TEMP 942 WORK( 3*N+JE ) = BCOEFI*P( JE-1, JE-1 ) / TEMP
942 END IF 943 END IF
943 * 944 *
944 XMAX = MAX( ABS( WORK( 2*N+JE ) )+ABS( WORK( 3*N+JE ) ), 945 XMAX = MAX( ABS( WORK( 2*N+JE ) )+ABS( WORK( 3*N+JE ) ),
945 $ ABS( WORK( 2*N+JE-1 ) )+ABS( WORK( 3*N+JE-1 ) ) ) 946 $ ABS( WORK( 2*N+JE-1 ) )+ABS( WORK( 3*N+JE-1 ) ) )
946 * 947 *
956 CRE2A = ACOEF*WORK( 2*N+JE ) 957 CRE2A = ACOEF*WORK( 2*N+JE )
957 CIM2A = ACOEF*WORK( 3*N+JE ) 958 CIM2A = ACOEF*WORK( 3*N+JE )
958 CRE2B = BCOEFR*WORK( 2*N+JE ) - BCOEFI*WORK( 3*N+JE ) 959 CRE2B = BCOEFR*WORK( 2*N+JE ) - BCOEFI*WORK( 3*N+JE )
959 CIM2B = BCOEFI*WORK( 2*N+JE ) + BCOEFR*WORK( 3*N+JE ) 960 CIM2B = BCOEFI*WORK( 2*N+JE ) + BCOEFR*WORK( 3*N+JE )
960 DO 270 JR = 1, JE - 2 961 DO 270 JR = 1, JE - 2
961 WORK( 2*N+JR ) = -CREALA*A( JR, JE-1 ) + 962 WORK( 2*N+JR ) = -CREALA*S( JR, JE-1 ) +
962 $ CREALB*B( JR, JE-1 ) - 963 $ CREALB*P( JR, JE-1 ) -
963 $ CRE2A*A( JR, JE ) + CRE2B*B( JR, JE ) 964 $ CRE2A*S( JR, JE ) + CRE2B*P( JR, JE )
964 WORK( 3*N+JR ) = -CIMAGA*A( JR, JE-1 ) + 965 WORK( 3*N+JR ) = -CIMAGA*S( JR, JE-1 ) +
965 $ CIMAGB*B( JR, JE-1 ) - 966 $ CIMAGB*P( JR, JE-1 ) -
966 $ CIM2A*A( JR, JE ) + CIM2B*B( JR, JE ) 967 $ CIM2A*S( JR, JE ) + CIM2B*P( JR, JE )
967 270 CONTINUE 968 270 CONTINUE
968 END IF 969 END IF
969 * 970 *
970 DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN ) 971 DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN )
971 * 972 *
976 * 977 *
977 * If a 2-by-2 block, is in position j-1:j, wait until 978 * If a 2-by-2 block, is in position j-1:j, wait until
978 * next iteration to process it (when it will be j:j+1) 979 * next iteration to process it (when it will be j:j+1)
979 * 980 *
980 IF( .NOT.IL2BY2 .AND. J.GT.1 ) THEN 981 IF( .NOT.IL2BY2 .AND. J.GT.1 ) THEN
981 IF( A( J, J-1 ).NE.ZERO ) THEN 982 IF( S( J, J-1 ).NE.ZERO ) THEN
982 IL2BY2 = .TRUE. 983 IL2BY2 = .TRUE.
983 GO TO 370 984 GO TO 370
984 END IF 985 END IF
985 END IF 986 END IF
986 BDIAG( 1 ) = B( J, J ) 987 BDIAG( 1 ) = P( J, J )
987 IF( IL2BY2 ) THEN 988 IF( IL2BY2 ) THEN
988 NA = 2 989 NA = 2
989 BDIAG( 2 ) = B( J+1, J+1 ) 990 BDIAG( 2 ) = P( J+1, J+1 )
990 ELSE 991 ELSE
991 NA = 1 992 NA = 1
992 END IF 993 END IF
993 * 994 *
994 * Compute x(j) (and x(j+1), if 2-by-2 block) 995 * Compute x(j) (and x(j+1), if 2-by-2 block)
995 * 996 *
996 CALL DLALN2( .FALSE., NA, NW, DMIN, ACOEF, A( J, J ), 997 CALL DLALN2( .FALSE., NA, NW, DMIN, ACOEF, S( J, J ),
997 $ LDA, BDIAG( 1 ), BDIAG( 2 ), WORK( 2*N+J ), 998 $ LDS, BDIAG( 1 ), BDIAG( 2 ), WORK( 2*N+J ),
998 $ N, BCOEFR, BCOEFI, SUM, 2, SCALE, TEMP, 999 $ N, BCOEFR, BCOEFI, SUM, 2, SCALE, TEMP,
999 $ IINFO ) 1000 $ IINFO )
1000 IF( SCALE.LT.ONE ) THEN 1001 IF( SCALE.LT.ONE ) THEN
1001 * 1002 *
1002 DO 290 JW = 0, NW - 1 1003 DO 290 JW = 0, NW - 1
1012 DO 300 JA = 1, NA 1013 DO 300 JA = 1, NA
1013 WORK( ( JW+1 )*N+J+JA-1 ) = SUM( JA, JW ) 1014 WORK( ( JW+1 )*N+J+JA-1 ) = SUM( JA, JW )
1014 300 CONTINUE 1015 300 CONTINUE
1015 310 CONTINUE 1016 310 CONTINUE
1016 * 1017 *
1017 * w = w + x(j)*(a A(*,j) - b B(*,j) ) with scaling 1018 * w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling
1018 * 1019 *
1019 IF( J.GT.1 ) THEN 1020 IF( J.GT.1 ) THEN
1020 * 1021 *
1021 * Check whether scaling is necessary for sum. 1022 * Check whether scaling is necessary for sum.
1022 * 1023 *
1050 $ BCOEFI*WORK( 3*N+J+JA-1 ) 1051 $ BCOEFI*WORK( 3*N+J+JA-1 )
1051 CIMAGB = BCOEFI*WORK( 2*N+J+JA-1 ) + 1052 CIMAGB = BCOEFI*WORK( 2*N+J+JA-1 ) +
1052 $ BCOEFR*WORK( 3*N+J+JA-1 ) 1053 $ BCOEFR*WORK( 3*N+J+JA-1 )
1053 DO 340 JR = 1, J - 1 1054 DO 340 JR = 1, J - 1
1054 WORK( 2*N+JR ) = WORK( 2*N+JR ) - 1055 WORK( 2*N+JR ) = WORK( 2*N+JR ) -
1055 $ CREALA*A( JR, J+JA-1 ) + 1056 $ CREALA*S( JR, J+JA-1 ) +
1056 $ CREALB*B( JR, J+JA-1 ) 1057 $ CREALB*P( JR, J+JA-1 )
1057 WORK( 3*N+JR ) = WORK( 3*N+JR ) - 1058 WORK( 3*N+JR ) = WORK( 3*N+JR ) -
1058 $ CIMAGA*A( JR, J+JA-1 ) + 1059 $ CIMAGA*S( JR, J+JA-1 ) +
1059 $ CIMAGB*B( JR, J+JA-1 ) 1060 $ CIMAGB*P( JR, J+JA-1 )
1060 340 CONTINUE 1061 340 CONTINUE
1061 ELSE 1062 ELSE
1062 CREALA = ACOEF*WORK( 2*N+J+JA-1 ) 1063 CREALA = ACOEF*WORK( 2*N+J+JA-1 )
1063 CREALB = BCOEFR*WORK( 2*N+J+JA-1 ) 1064 CREALB = BCOEFR*WORK( 2*N+J+JA-1 )
1064 DO 350 JR = 1, J - 1 1065 DO 350 JR = 1, J - 1
1065 WORK( 2*N+JR ) = WORK( 2*N+JR ) - 1066 WORK( 2*N+JR ) = WORK( 2*N+JR ) -
1066 $ CREALA*A( JR, J+JA-1 ) + 1067 $ CREALA*S( JR, J+JA-1 ) +
1067 $ CREALB*B( JR, J+JA-1 ) 1068 $ CREALB*P( JR, J+JA-1 )
1068 350 CONTINUE 1069 350 CONTINUE
1069 END IF 1070 END IF
1070 360 CONTINUE 1071 360 CONTINUE
1071 END IF 1072 END IF
1072 * 1073 *