ergo
template_lapack_tgevc.h
Go to the documentation of this file.
1 /* Ergo, version 3.2, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2012 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  *
18  * Primary academic reference:
19  * Kohn−Sham Density Functional Theory Electronic Structure Calculations
20  * with Linearly Scaling Computational Time and Memory Usage,
21  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
22  * J. Chem. Theory Comput. 7, 340 (2011),
23  * <http://dx.doi.org/10.1021/ct100611z>
24  *
25  * For further information about Ergo, see <http://www.ergoscf.org>.
26  */
27 
28  /* This file belongs to the template_lapack part of the Ergo source
29  * code. The source files in the template_lapack directory are modified
30  * versions of files originally distributed as CLAPACK, see the
31  * Copyright/license notice in the file template_lapack/COPYING.
32  */
33 
34 
35 #ifndef TEMPLATE_LAPACK_TGEVC_HEADER
36 #define TEMPLATE_LAPACK_TGEVC_HEADER
37 
38 
39 template<class Treal>
40 int template_lapack_tgevc(const char *side, const char *howmny, const logical *select,
41  const integer *n, const Treal *a, const integer *lda, const Treal *b, const integer *ldb,
42  Treal *vl, const integer *ldvl, Treal *vr, const integer *ldvr, const integer
43  *mm, integer *m, Treal *work, integer *info)
44 {
45 /* -- LAPACK routine (version 3.0) --
46  Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
47  Courant Institute, Argonne National Lab, and Rice University
48  June 30, 1999
49 
50 
51 
52  Purpose
53  =======
54 
55  DTGEVC computes some or all of the right and/or left generalized
56  eigenvectors of a pair of real upper triangular matrices (A,B).
57 
58  The right generalized eigenvector x and the left generalized
59  eigenvector y of (A,B) corresponding to a generalized eigenvalue
60  w are defined by:
61 
62  (A - wB) * x = 0 and y**H * (A - wB) = 0
63 
64  where y**H denotes the conjugate tranpose of y.
65 
66  If an eigenvalue w is determined by zero diagonal elements of both A
67  and B, a unit vector is returned as the corresponding eigenvector.
68 
69  If all eigenvectors are requested, the routine may either return
70  the matrices X and/or Y of right or left eigenvectors of (A,B), or
71  the products Z*X and/or Q*Y, where Z and Q are input orthogonal
72  matrices. If (A,B) was obtained from the generalized real-Schur
73  factorization of an original pair of matrices
74  (A0,B0) = (Q*A*Z**H,Q*B*Z**H),
75  then Z*X and Q*Y are the matrices of right or left eigenvectors of
76  A.
77 
78  A must be block upper triangular, with 1-by-1 and 2-by-2 diagonal
79  blocks. Corresponding to each 2-by-2 diagonal block is a complex
80  conjugate pair of eigenvalues and eigenvectors; only one
81  eigenvector of the pair is computed, namely the one corresponding
82  to the eigenvalue with positive imaginary part.
83 
84  Arguments
85  =========
86 
87  SIDE (input) CHARACTER*1
88  = 'R': compute right eigenvectors only;
89  = 'L': compute left eigenvectors only;
90  = 'B': compute both right and left eigenvectors.
91 
92  HOWMNY (input) CHARACTER*1
93  = 'A': compute all right and/or left eigenvectors;
94  = 'B': compute all right and/or left eigenvectors, and
95  backtransform them using the input matrices supplied
96  in VR and/or VL;
97  = 'S': compute selected right and/or left eigenvectors,
98  specified by the logical array SELECT.
99 
100  SELECT (input) LOGICAL array, dimension (N)
101  If HOWMNY='S', SELECT specifies the eigenvectors to be
102  computed.
103  If HOWMNY='A' or 'B', SELECT is not referenced.
104  To select the real eigenvector corresponding to the real
105  eigenvalue w(j), SELECT(j) must be set to .TRUE. To select
106  the complex eigenvector corresponding to a complex conjugate
107  pair w(j) and w(j+1), either SELECT(j) or SELECT(j+1) must
108  be set to .TRUE..
109 
110  N (input) INTEGER
111  The order of the matrices A and B. N >= 0.
112 
113  A (input) DOUBLE PRECISION array, dimension (LDA,N)
114  The upper quasi-triangular matrix A.
115 
116  LDA (input) INTEGER
117  The leading dimension of array A. LDA >= max(1, N).
118 
119  B (input) DOUBLE PRECISION array, dimension (LDB,N)
120  The upper triangular matrix B. If A has a 2-by-2 diagonal
121  block, then the corresponding 2-by-2 block of B must be
122  diagonal with positive elements.
123 
124  LDB (input) INTEGER
125  The leading dimension of array B. LDB >= max(1,N).
126 
127  VL (input/output) DOUBLE PRECISION array, dimension (LDVL,MM)
128  On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
129  contain an N-by-N matrix Q (usually the orthogonal matrix Q
130  of left Schur vectors returned by DHGEQZ).
131  On exit, if SIDE = 'L' or 'B', VL contains:
132  if HOWMNY = 'A', the matrix Y of left eigenvectors of (A,B);
133  if HOWMNY = 'B', the matrix Q*Y;
134  if HOWMNY = 'S', the left eigenvectors of (A,B) specified by
135  SELECT, stored consecutively in the columns of
136  VL, in the same order as their eigenvalues.
137  If SIDE = 'R', VL is not referenced.
138 
139  A complex eigenvector corresponding to a complex eigenvalue
140  is stored in two consecutive columns, the first holding the
141  real part, and the second the imaginary part.
142 
143  LDVL (input) INTEGER
144  The leading dimension of array VL.
145  LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise.
146 
147  VR (input/output) DOUBLE PRECISION array, dimension (LDVR,MM)
148  On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
149  contain an N-by-N matrix Q (usually the orthogonal matrix Z
150  of right Schur vectors returned by DHGEQZ).
151  On exit, if SIDE = 'R' or 'B', VR contains:
152  if HOWMNY = 'A', the matrix X of right eigenvectors of (A,B);
153  if HOWMNY = 'B', the matrix Z*X;
154  if HOWMNY = 'S', the right eigenvectors of (A,B) specified by
155  SELECT, stored consecutively in the columns of
156  VR, in the same order as their eigenvalues.
157  If SIDE = 'L', VR is not referenced.
158 
159  A complex eigenvector corresponding to a complex eigenvalue
160  is stored in two consecutive columns, the first holding the
161  real part and the second the imaginary part.
162 
163  LDVR (input) INTEGER
164  The leading dimension of the array VR.
165  LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise.
166 
167  MM (input) INTEGER
168  The number of columns in the arrays VL and/or VR. MM >= M.
169 
170  M (output) INTEGER
171  The number of columns in the arrays VL and/or VR actually
172  used to store the eigenvectors. If HOWMNY = 'A' or 'B', M
173  is set to N. Each selected real eigenvector occupies one
174  column and each selected complex eigenvector occupies two
175  columns.
176 
177  WORK (workspace) DOUBLE PRECISION array, dimension (6*N)
178 
179  INFO (output) INTEGER
180  = 0: successful exit.
181  < 0: if INFO = -i, the i-th argument had an illegal value.
182  > 0: the 2-by-2 block (INFO:INFO+1) does not have a complex
183  eigenvalue.
184 
185  Further Details
186  ===============
187 
188  Allocation of workspace:
189  ---------- -- ---------
190 
191  WORK( j ) = 1-norm of j-th column of A, above the diagonal
192  WORK( N+j ) = 1-norm of j-th column of B, above the diagonal
193  WORK( 2*N+1:3*N ) = real part of eigenvector
194  WORK( 3*N+1:4*N ) = imaginary part of eigenvector
195  WORK( 4*N+1:5*N ) = real part of back-transformed eigenvector
196  WORK( 5*N+1:6*N ) = imaginary part of back-transformed eigenvector
197 
198  Rowwise vs. columnwise solution methods:
199  ------- -- ---------- -------- -------
200 
201  Finding a generalized eigenvector consists basically of solving the
202  singular triangular system
203 
204  (A - w B) x = 0 (for right) or: (A - w B)**H y = 0 (for left)
205 
206  Consider finding the i-th right eigenvector (assume all eigenvalues
207  are real). The equation to be solved is:
208  n i
209  0 = sum C(j,k) v(k) = sum C(j,k) v(k) for j = i,. . .,1
210  k=j k=j
211 
212  where C = (A - w B) (The components v(i+1:n) are 0.)
213 
214  The "rowwise" method is:
215 
216  (1) v(i) := 1
217  for j = i-1,. . .,1:
218  i
219  (2) compute s = - sum C(j,k) v(k) and
220  k=j+1
221 
222  (3) v(j) := s / C(j,j)
223 
224  Step 2 is sometimes called the "dot product" step, since it is an
225  inner product between the j-th row and the portion of the eigenvector
226  that has been computed so far.
227 
228  The "columnwise" method consists basically in doing the sums
229  for all the rows in parallel. As each v(j) is computed, the
230  contribution of v(j) times the j-th column of C is added to the
231  partial sums. Since FORTRAN arrays are stored columnwise, this has
232  the advantage that at each step, the elements of C that are accessed
233  are adjacent to one another, whereas with the rowwise method, the
234  elements accessed at a step are spaced LDA (and LDB) words apart.
235 
236  When finding left eigenvectors, the matrix in question is the
237  transpose of the one in storage, so the rowwise method then
238  actually accesses columns of A and B at each step, and so is the
239  preferred method.
240 
241  =====================================================================
242 
243 
244  Decode and Test the input parameters
245 
246  Parameter adjustments */
247  /* Table of constant values */
248  logical c_true = TRUE_;
249  integer c__2 = 2;
250  Treal c_b35 = 1.;
251  integer c__1 = 1;
252  Treal c_b37 = 0.;
253  logical c_false = FALSE_;
254 
255  /* System generated locals */
256  integer a_dim1, a_offset, b_dim1, b_offset, vl_dim1, vl_offset, vr_dim1,
257  vr_offset, i__1, i__2, i__3, i__4, i__5;
258  Treal d__1, d__2, d__3, d__4, d__5, d__6;
259  /* Local variables */
260  integer ibeg, ieig, iend;
261  Treal dmin__, temp, suma[4] /* was [2][2] */, sumb[4]
262  /* was [2][2] */, xmax;
263  Treal cim2a, cim2b, cre2a, cre2b, temp2, bdiag[2];
264  integer i__, j;
265  Treal acoef, scale;
266  logical ilall;
267  integer iside;
268  Treal sbeta;
269  logical il2by2;
270  integer iinfo;
271  Treal small;
272  logical compl_AAAA;
273  Treal anorm, bnorm;
274  logical compr;
275  Treal temp2i;
276  Treal temp2r;
277  integer ja;
278  logical ilabad, ilbbad;
279  integer jc, je, na;
280  Treal acoefa, bcoefa, cimaga, cimagb;
281  logical ilback;
282  integer im;
283  Treal bcoefi, ascale, bscale, creala;
284  integer jr;
285  Treal crealb;
286  Treal bcoefr;
287  integer jw, nw;
288  Treal salfar, safmin;
289  Treal xscale, bignum;
290  logical ilcomp, ilcplx;
291  integer ihwmny;
292  Treal big;
293  logical lsa, lsb;
294  Treal ulp, sum[4] /* was [2][2] */;
295 #define suma_ref(a_1,a_2) suma[(a_2)*2 + a_1 - 3]
296 #define sumb_ref(a_1,a_2) sumb[(a_2)*2 + a_1 - 3]
297 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
298 #define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]
299 #define vl_ref(a_1,a_2) vl[(a_2)*vl_dim1 + a_1]
300 #define vr_ref(a_1,a_2) vr[(a_2)*vr_dim1 + a_1]
301 #define sum_ref(a_1,a_2) sum[(a_2)*2 + a_1 - 3]
302 
303 
304  --select;
305  a_dim1 = *lda;
306  a_offset = 1 + a_dim1 * 1;
307  a -= a_offset;
308  b_dim1 = *ldb;
309  b_offset = 1 + b_dim1 * 1;
310  b -= b_offset;
311  vl_dim1 = *ldvl;
312  vl_offset = 1 + vl_dim1 * 1;
313  vl -= vl_offset;
314  vr_dim1 = *ldvr;
315  vr_offset = 1 + vr_dim1 * 1;
316  vr -= vr_offset;
317  --work;
318 
319  /* Initialization added by Elias to get rid of compiler warnings. */
320  ilback = 0;
321  /* Function Body */
322  if (template_blas_lsame(howmny, "A")) {
323  ihwmny = 1;
324  ilall = TRUE_;
325  ilback = FALSE_;
326  } else if (template_blas_lsame(howmny, "S")) {
327  ihwmny = 2;
328  ilall = FALSE_;
329  ilback = FALSE_;
330  } else if (template_blas_lsame(howmny, "B") || template_blas_lsame(howmny,
331  "T")) {
332  ihwmny = 3;
333  ilall = TRUE_;
334  ilback = TRUE_;
335  } else {
336  ihwmny = -1;
337  ilall = TRUE_;
338  }
339 
340  if (template_blas_lsame(side, "R")) {
341  iside = 1;
342  compl_AAAA = FALSE_;
343  compr = TRUE_;
344  } else if (template_blas_lsame(side, "L")) {
345  iside = 2;
346  compl_AAAA = TRUE_;
347  compr = FALSE_;
348  } else if (template_blas_lsame(side, "B")) {
349  iside = 3;
350  compl_AAAA = TRUE_;
351  compr = TRUE_;
352  } else {
353  iside = -1;
354  }
355 
356  *info = 0;
357  if (iside < 0) {
358  *info = -1;
359  } else if (ihwmny < 0) {
360  *info = -2;
361  } else if (*n < 0) {
362  *info = -4;
363  } else if (*lda < maxMACRO(1,*n)) {
364  *info = -6;
365  } else if (*ldb < maxMACRO(1,*n)) {
366  *info = -8;
367  }
368  if (*info != 0) {
369  i__1 = -(*info);
370  template_blas_erbla("TGEVC ", &i__1);
371  return 0;
372  }
373 
374 /* Count the number of eigenvectors to be computed */
375 
376  if (! ilall) {
377  im = 0;
378  ilcplx = FALSE_;
379  i__1 = *n;
380  for (j = 1; j <= i__1; ++j) {
381  if (ilcplx) {
382  ilcplx = FALSE_;
383  goto L10;
384  }
385  if (j < *n) {
386  if (a_ref(j + 1, j) != 0.) {
387  ilcplx = TRUE_;
388  }
389  }
390  if (ilcplx) {
391  if (select[j] || select[j + 1]) {
392  im += 2;
393  }
394  } else {
395  if (select[j]) {
396  ++im;
397  }
398  }
399 L10:
400  ;
401  }
402  } else {
403  im = *n;
404  }
405 
406 /* Check 2-by-2 diagonal blocks of A, B */
407 
408  ilabad = FALSE_;
409  ilbbad = FALSE_;
410  i__1 = *n - 1;
411  for (j = 1; j <= i__1; ++j) {
412  if (a_ref(j + 1, j) != 0.) {
413  if (b_ref(j, j) == 0. || b_ref(j + 1, j + 1) == 0. || b_ref(j, j
414  + 1) != 0.) {
415  ilbbad = TRUE_;
416  }
417  if (j < *n - 1) {
418  if (a_ref(j + 2, j + 1) != 0.) {
419  ilabad = TRUE_;
420  }
421  }
422  }
423 /* L20: */
424  }
425 
426  if (ilabad) {
427  *info = -5;
428  } else if (ilbbad) {
429  *info = -7;
430  } else if ( ( compl_AAAA && *ldvl < *n ) || *ldvl < 1) {
431  *info = -10;
432  } else if ( ( compr && *ldvr < *n ) || *ldvr < 1) {
433  *info = -12;
434  } else if (*mm < im) {
435  *info = -13;
436  }
437  if (*info != 0) {
438  i__1 = -(*info);
439  template_blas_erbla("TGEVC ", &i__1);
440  return 0;
441  }
442 
443 /* Quick return if possible */
444 
445  *m = im;
446  if (*n == 0) {
447  return 0;
448  }
449 
450 /* Machine Constants */
451 
452  safmin = template_lapack_lamch("Safe minimum", (Treal)0);
453  big = 1. / safmin;
454  template_lapack_labad(&safmin, &big);
455  ulp = template_lapack_lamch("Epsilon", (Treal)0) * template_lapack_lamch("Base", (Treal)0);
456  small = safmin * *n / ulp;
457  big = 1. / small;
458  bignum = 1. / (safmin * *n);
459 
460 /* Compute the 1-norm of each column of the strictly upper triangular
461  part (i.e., excluding all elements belonging to the diagonal
462  blocks) of A and B to check for possible overflow in the
463  triangular solver. */
464 
465  anorm = (d__1 = a_ref(1, 1), absMACRO(d__1));
466  if (*n > 1) {
467  anorm += (d__1 = a_ref(2, 1), absMACRO(d__1));
468  }
469  bnorm = (d__1 = b_ref(1, 1), absMACRO(d__1));
470  work[1] = 0.;
471  work[*n + 1] = 0.;
472 
473  i__1 = *n;
474  for (j = 2; j <= i__1; ++j) {
475  temp = 0.;
476  temp2 = 0.;
477  if (a_ref(j, j - 1) == 0.) {
478  iend = j - 1;
479  } else {
480  iend = j - 2;
481  }
482  i__2 = iend;
483  for (i__ = 1; i__ <= i__2; ++i__) {
484  temp += (d__1 = a_ref(i__, j), absMACRO(d__1));
485  temp2 += (d__1 = b_ref(i__, j), absMACRO(d__1));
486 /* L30: */
487  }
488  work[j] = temp;
489  work[*n + j] = temp2;
490 /* Computing MIN */
491  i__3 = j + 1;
492  i__2 = minMACRO(i__3,*n);
493  for (i__ = iend + 1; i__ <= i__2; ++i__) {
494  temp += (d__1 = a_ref(i__, j), absMACRO(d__1));
495  temp2 += (d__1 = b_ref(i__, j), absMACRO(d__1));
496 /* L40: */
497  }
498  anorm = maxMACRO(anorm,temp);
499  bnorm = maxMACRO(bnorm,temp2);
500 /* L50: */
501  }
502 
503  ascale = 1. / maxMACRO(anorm,safmin);
504  bscale = 1. / maxMACRO(bnorm,safmin);
505 
506 /* Left eigenvectors */
507 
508  if (compl_AAAA) {
509  ieig = 0;
510 
511 /* Main loop over eigenvalues */
512 
513  ilcplx = FALSE_;
514  i__1 = *n;
515  for (je = 1; je <= i__1; ++je) {
516 
517 /* Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or
518  (b) this would be the second of a complex pair.
519  Check for complex eigenvalue, so as to be sure of which
520  entry(-ies) of SELECT to look at. */
521 
522  if (ilcplx) {
523  ilcplx = FALSE_;
524  goto L220;
525  }
526  nw = 1;
527  if (je < *n) {
528  if (a_ref(je + 1, je) != 0.) {
529  ilcplx = TRUE_;
530  nw = 2;
531  }
532  }
533  if (ilall) {
534  ilcomp = TRUE_;
535  } else if (ilcplx) {
536  ilcomp = select[je] || select[je + 1];
537  } else {
538  ilcomp = select[je];
539  }
540  if (! ilcomp) {
541  goto L220;
542  }
543 
544 /* Decide if (a) singular pencil, (b) real eigenvalue, or
545  (c) complex eigenvalue. */
546 
547  if (! ilcplx) {
548  if ((d__1 = a_ref(je, je), absMACRO(d__1)) <= safmin && (d__2 =
549  b_ref(je, je), absMACRO(d__2)) <= safmin) {
550 
551 /* Singular matrix pencil -- return unit eigenvector */
552 
553  ++ieig;
554  i__2 = *n;
555  for (jr = 1; jr <= i__2; ++jr) {
556  vl_ref(jr, ieig) = 0.;
557 /* L60: */
558  }
559  vl_ref(ieig, ieig) = 1.;
560  goto L220;
561  }
562  }
563 
564 /* Clear vector */
565 
566  i__2 = nw * *n;
567  for (jr = 1; jr <= i__2; ++jr) {
568  work[(*n << 1) + jr] = 0.;
569 /* L70: */
570  }
571 /* T
572  Compute coefficients in ( a A - b B ) y = 0
573  a is ACOEF
574  b is BCOEFR + i*BCOEFI */
575 
576  if (! ilcplx) {
577 
578 /* Real eigenvalue
579 
580  Computing MAX */
581  d__3 = (d__1 = a_ref(je, je), absMACRO(d__1)) * ascale, d__4 = (
582  d__2 = b_ref(je, je), absMACRO(d__2)) * bscale, d__3 = maxMACRO(
583  d__3,d__4);
584  temp = 1. / maxMACRO(d__3,safmin);
585  salfar = temp * a_ref(je, je) * ascale;
586  sbeta = temp * b_ref(je, je) * bscale;
587  acoef = sbeta * ascale;
588  bcoefr = salfar * bscale;
589  bcoefi = 0.;
590 
591 /* Scale to avoid underflow */
592 
593  scale = 1.;
594  lsa = absMACRO(sbeta) >= safmin && absMACRO(acoef) < small;
595  lsb = absMACRO(salfar) >= safmin && absMACRO(bcoefr) < small;
596  if (lsa) {
597  scale = small / absMACRO(sbeta) * minMACRO(anorm,big);
598  }
599  if (lsb) {
600 /* Computing MAX */
601  d__1 = scale, d__2 = small / absMACRO(salfar) * minMACRO(bnorm,big);
602  scale = maxMACRO(d__1,d__2);
603  }
604  if (lsa || lsb) {
605 /* Computing MIN
606  Computing MAX */
607  d__3 = 1., d__4 = absMACRO(acoef), d__3 = maxMACRO(d__3,d__4), d__4
608  = absMACRO(bcoefr);
609  d__1 = scale, d__2 = 1. / (safmin * maxMACRO(d__3,d__4));
610  scale = minMACRO(d__1,d__2);
611  if (lsa) {
612  acoef = ascale * (scale * sbeta);
613  } else {
614  acoef = scale * acoef;
615  }
616  if (lsb) {
617  bcoefr = bscale * (scale * salfar);
618  } else {
619  bcoefr = scale * bcoefr;
620  }
621  }
622  acoefa = absMACRO(acoef);
623  bcoefa = absMACRO(bcoefr);
624 
625 /* First component is 1 */
626 
627  work[(*n << 1) + je] = 1.;
628  xmax = 1.;
629  } else {
630 
631 /* Complex eigenvalue */
632 
633  d__1 = safmin * 100.;
634  template_lapack_lag2(&a_ref(je, je), lda, &b_ref(je, je), ldb, &d__1, &
635  acoef, &temp, &bcoefr, &temp2, &bcoefi);
636  bcoefi = -bcoefi;
637  if (bcoefi == 0.) {
638  *info = je;
639  return 0;
640  }
641 
642 /* Scale to avoid over/underflow */
643 
644  acoefa = absMACRO(acoef);
645  bcoefa = absMACRO(bcoefr) + absMACRO(bcoefi);
646  scale = 1.;
647  if (acoefa * ulp < safmin && acoefa >= safmin) {
648  scale = safmin / ulp / acoefa;
649  }
650  if (bcoefa * ulp < safmin && bcoefa >= safmin) {
651 /* Computing MAX */
652  d__1 = scale, d__2 = safmin / ulp / bcoefa;
653  scale = maxMACRO(d__1,d__2);
654  }
655  if (safmin * acoefa > ascale) {
656  scale = ascale / (safmin * acoefa);
657  }
658  if (safmin * bcoefa > bscale) {
659 /* Computing MIN */
660  d__1 = scale, d__2 = bscale / (safmin * bcoefa);
661  scale = minMACRO(d__1,d__2);
662  }
663  if (scale != 1.) {
664  acoef = scale * acoef;
665  acoefa = absMACRO(acoef);
666  bcoefr = scale * bcoefr;
667  bcoefi = scale * bcoefi;
668  bcoefa = absMACRO(bcoefr) + absMACRO(bcoefi);
669  }
670 
671 /* Compute first two components of eigenvector */
672 
673  temp = acoef * a_ref(je + 1, je);
674  temp2r = acoef * a_ref(je, je) - bcoefr * b_ref(je, je);
675  temp2i = -bcoefi * b_ref(je, je);
676  if (absMACRO(temp) > absMACRO(temp2r) + absMACRO(temp2i)) {
677  work[(*n << 1) + je] = 1.;
678  work[*n * 3 + je] = 0.;
679  work[(*n << 1) + je + 1] = -temp2r / temp;
680  work[*n * 3 + je + 1] = -temp2i / temp;
681  } else {
682  work[(*n << 1) + je + 1] = 1.;
683  work[*n * 3 + je + 1] = 0.;
684  temp = acoef * a_ref(je, je + 1);
685  work[(*n << 1) + je] = (bcoefr * b_ref(je + 1, je + 1) -
686  acoef * a_ref(je + 1, je + 1)) / temp;
687  work[*n * 3 + je] = bcoefi * b_ref(je + 1, je + 1) / temp;
688  }
689 /* Computing MAX */
690  d__5 = (d__1 = work[(*n << 1) + je], absMACRO(d__1)) + (d__2 =
691  work[*n * 3 + je], absMACRO(d__2)), d__6 = (d__3 = work[(*
692  n << 1) + je + 1], absMACRO(d__3)) + (d__4 = work[*n * 3 +
693  je + 1], absMACRO(d__4));
694  xmax = maxMACRO(d__5,d__6);
695  }
696 
697 /* Computing MAX */
698  d__1 = ulp * acoefa * anorm, d__2 = ulp * bcoefa * bnorm, d__1 =
699  maxMACRO(d__1,d__2);
700  dmin__ = maxMACRO(d__1,safmin);
701 
702 /* T
703  Triangular solve of (a A - b B) y = 0
704 
705  T
706  (rowwise in (a A - b B) , or columnwise in (a A - b B) ) */
707 
708  il2by2 = FALSE_;
709 
710  i__2 = *n;
711  for (j = je + nw; j <= i__2; ++j) {
712  if (il2by2) {
713  il2by2 = FALSE_;
714  goto L160;
715  }
716 
717  na = 1;
718  bdiag[0] = b_ref(j, j);
719  if (j < *n) {
720  if (a_ref(j + 1, j) != 0.) {
721  il2by2 = TRUE_;
722  bdiag[1] = b_ref(j + 1, j + 1);
723  na = 2;
724  }
725  }
726 
727 /* Check whether scaling is necessary for dot products */
728 
729  xscale = 1. / maxMACRO(1.,xmax);
730 /* Computing MAX */
731  d__1 = work[j], d__2 = work[*n + j], d__1 = maxMACRO(d__1,d__2),
732  d__2 = acoefa * work[j] + bcoefa * work[*n + j];
733  temp = maxMACRO(d__1,d__2);
734  if (il2by2) {
735 /* Computing MAX */
736  d__1 = temp, d__2 = work[j + 1], d__1 = maxMACRO(d__1,d__2),
737  d__2 = work[*n + j + 1], d__1 = maxMACRO(d__1,d__2),
738  d__2 = acoefa * work[j + 1] + bcoefa * work[*n +
739  j + 1];
740  temp = maxMACRO(d__1,d__2);
741  }
742  if (temp > bignum * xscale) {
743  i__3 = nw - 1;
744  for (jw = 0; jw <= i__3; ++jw) {
745  i__4 = j - 1;
746  for (jr = je; jr <= i__4; ++jr) {
747  work[(jw + 2) * *n + jr] = xscale * work[(jw + 2)
748  * *n + jr];
749 /* L80: */
750  }
751 /* L90: */
752  }
753  xmax *= xscale;
754  }
755 
756 /* Compute dot products
757 
758  j-1
759  SUM = sum conjg( a*A(k,j) - b*B(k,j) )*x(k)
760  k=je
761 
762  To reduce the op count, this is done as
763 
764  _ j-1 _ j-1
765  a*conjg( sum A(k,j)*x(k) ) - b*conjg( sum B(k,j)*x(k) )
766  k=je k=je
767 
768  which may cause underflow problems if A or B are close
769  to underflow. (E.g., less than SMALL.)
770 
771 
772  A series of compiler directives to defeat vectorization
773  for the next loop
774 
775  $PL$ CMCHAR=' '
776  DIR$ NEXTSCALAR
777  $DIR SCALAR
778  DIR$ NEXT SCALAR
779  VD$L NOVECTOR
780  DEC$ NOVECTOR
781  VD$ NOVECTOR
782  VDIR NOVECTOR
783  VOCL LOOP,SCALAR
784  IBM PREFER SCALAR
785  $PL$ CMCHAR='*' */
786 
787  i__3 = nw;
788  for (jw = 1; jw <= i__3; ++jw) {
789 
790 /* $PL$ CMCHAR=' '
791  DIR$ NEXTSCALAR
792  $DIR SCALAR
793  DIR$ NEXT SCALAR
794  VD$L NOVECTOR
795  DEC$ NOVECTOR
796  VD$ NOVECTOR
797  VDIR NOVECTOR
798  VOCL LOOP,SCALAR
799  IBM PREFER SCALAR
800  $PL$ CMCHAR='*' */
801 
802  i__4 = na;
803  for (ja = 1; ja <= i__4; ++ja) {
804  suma_ref(ja, jw) = 0.;
805  sumb_ref(ja, jw) = 0.;
806 
807  i__5 = j - 1;
808  for (jr = je; jr <= i__5; ++jr) {
809  suma_ref(ja, jw) = suma_ref(ja, jw) + a_ref(jr, j
810  + ja - 1) * work[(jw + 1) * *n + jr];
811  sumb_ref(ja, jw) = sumb_ref(ja, jw) + b_ref(jr, j
812  + ja - 1) * work[(jw + 1) * *n + jr];
813 /* L100: */
814  }
815 /* L110: */
816  }
817 /* L120: */
818  }
819 
820 /* $PL$ CMCHAR=' '
821  DIR$ NEXTSCALAR
822  $DIR SCALAR
823  DIR$ NEXT SCALAR
824  VD$L NOVECTOR
825  DEC$ NOVECTOR
826  VD$ NOVECTOR
827  VDIR NOVECTOR
828  VOCL LOOP,SCALAR
829  IBM PREFER SCALAR
830  $PL$ CMCHAR='*' */
831 
832  i__3 = na;
833  for (ja = 1; ja <= i__3; ++ja) {
834  if (ilcplx) {
835  sum_ref(ja, 1) = -acoef * suma_ref(ja, 1) + bcoefr *
836  sumb_ref(ja, 1) - bcoefi * sumb_ref(ja, 2);
837  sum_ref(ja, 2) = -acoef * suma_ref(ja, 2) + bcoefr *
838  sumb_ref(ja, 2) + bcoefi * sumb_ref(ja, 1);
839  } else {
840  sum_ref(ja, 1) = -acoef * suma_ref(ja, 1) + bcoefr *
841  sumb_ref(ja, 1);
842  }
843 /* L130: */
844  }
845 
846 /* T
847  Solve ( a A - b B ) y = SUM(,)
848  with scaling and perturbation of the denominator */
849 
850  template_lapack_laln2(&c_true, &na, &nw, &dmin__, &acoef, &a_ref(j, j), lda,
851  bdiag, &bdiag[1], sum, &c__2, &bcoefr, &bcoefi, &
852  work[(*n << 1) + j], n, &scale, &temp, &iinfo);
853  if (scale < 1.) {
854  i__3 = nw - 1;
855  for (jw = 0; jw <= i__3; ++jw) {
856  i__4 = j - 1;
857  for (jr = je; jr <= i__4; ++jr) {
858  work[(jw + 2) * *n + jr] = scale * work[(jw + 2) *
859  *n + jr];
860 /* L140: */
861  }
862 /* L150: */
863  }
864  xmax = scale * xmax;
865  }
866  xmax = maxMACRO(xmax,temp);
867 L160:
868  ;
869  }
870 
871 /* Copy eigenvector to VL, back transforming if
872  HOWMNY='B'. */
873 
874  ++ieig;
875  if (ilback) {
876  i__2 = nw - 1;
877  for (jw = 0; jw <= i__2; ++jw) {
878  i__3 = *n + 1 - je;
879  template_blas_gemv("N", n, &i__3, &c_b35, &vl_ref(1, je), ldvl, &work[
880  (jw + 2) * *n + je], &c__1, &c_b37, &work[(jw + 4)
881  * *n + 1], &c__1);
882 /* L170: */
883  }
884  template_lapack_lacpy(" ", n, &nw, &work[(*n << 2) + 1], n, &vl_ref(1, je),
885  ldvl);
886  ibeg = 1;
887  } else {
888  template_lapack_lacpy(" ", n, &nw, &work[(*n << 1) + 1], n, &vl_ref(1, ieig)
889  , ldvl);
890  ibeg = je;
891  }
892 
893 /* Scale eigenvector */
894 
895  xmax = 0.;
896  if (ilcplx) {
897  i__2 = *n;
898  for (j = ibeg; j <= i__2; ++j) {
899 /* Computing MAX */
900  d__3 = xmax, d__4 = (d__1 = vl_ref(j, ieig), absMACRO(d__1)) +
901  (d__2 = vl_ref(j, ieig + 1), absMACRO(d__2));
902  xmax = maxMACRO(d__3,d__4);
903 /* L180: */
904  }
905  } else {
906  i__2 = *n;
907  for (j = ibeg; j <= i__2; ++j) {
908 /* Computing MAX */
909  d__2 = xmax, d__3 = (d__1 = vl_ref(j, ieig), absMACRO(d__1));
910  xmax = maxMACRO(d__2,d__3);
911 /* L190: */
912  }
913  }
914 
915  if (xmax > safmin) {
916  xscale = 1. / xmax;
917 
918  i__2 = nw - 1;
919  for (jw = 0; jw <= i__2; ++jw) {
920  i__3 = *n;
921  for (jr = ibeg; jr <= i__3; ++jr) {
922  vl_ref(jr, ieig + jw) = xscale * vl_ref(jr, ieig + jw)
923  ;
924 /* L200: */
925  }
926 /* L210: */
927  }
928  }
929  ieig = ieig + nw - 1;
930 
931 L220:
932  ;
933  }
934  }
935 
936 /* Right eigenvectors */
937 
938  if (compr) {
939  ieig = im + 1;
940 
941 /* Main loop over eigenvalues */
942 
943  ilcplx = FALSE_;
944  for (je = *n; je >= 1; --je) {
945 
946 /* Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or
947  (b) this would be the second of a complex pair.
948  Check for complex eigenvalue, so as to be sure of which
949  entry(-ies) of SELECT to look at -- if complex, SELECT(JE)
950  or SELECT(JE-1).
951  If this is a complex pair, the 2-by-2 diagonal block
952  corresponding to the eigenvalue is in rows/columns JE-1:JE */
953 
954  if (ilcplx) {
955  ilcplx = FALSE_;
956  goto L500;
957  }
958  nw = 1;
959  if (je > 1) {
960  if (a_ref(je, je - 1) != 0.) {
961  ilcplx = TRUE_;
962  nw = 2;
963  }
964  }
965  if (ilall) {
966  ilcomp = TRUE_;
967  } else if (ilcplx) {
968  ilcomp = select[je] || select[je - 1];
969  } else {
970  ilcomp = select[je];
971  }
972  if (! ilcomp) {
973  goto L500;
974  }
975 
976 /* Decide if (a) singular pencil, (b) real eigenvalue, or
977  (c) complex eigenvalue. */
978 
979  if (! ilcplx) {
980  if ((d__1 = a_ref(je, je), absMACRO(d__1)) <= safmin && (d__2 =
981  b_ref(je, je), absMACRO(d__2)) <= safmin) {
982 
983 /* Singular matrix pencil -- unit eigenvector */
984 
985  --ieig;
986  i__1 = *n;
987  for (jr = 1; jr <= i__1; ++jr) {
988  vr_ref(jr, ieig) = 0.;
989 /* L230: */
990  }
991  vr_ref(ieig, ieig) = 1.;
992  goto L500;
993  }
994  }
995 
996 /* Clear vector */
997 
998  i__1 = nw - 1;
999  for (jw = 0; jw <= i__1; ++jw) {
1000  i__2 = *n;
1001  for (jr = 1; jr <= i__2; ++jr) {
1002  work[(jw + 2) * *n + jr] = 0.;
1003 /* L240: */
1004  }
1005 /* L250: */
1006  }
1007 
1008 /* Compute coefficients in ( a A - b B ) x = 0
1009  a is ACOEF
1010  b is BCOEFR + i*BCOEFI */
1011 
1012  if (! ilcplx) {
1013 
1014 /* Real eigenvalue
1015 
1016  Computing MAX */
1017  d__3 = (d__1 = a_ref(je, je), absMACRO(d__1)) * ascale, d__4 = (
1018  d__2 = b_ref(je, je), absMACRO(d__2)) * bscale, d__3 = maxMACRO(
1019  d__3,d__4);
1020  temp = 1. / maxMACRO(d__3,safmin);
1021  salfar = temp * a_ref(je, je) * ascale;
1022  sbeta = temp * b_ref(je, je) * bscale;
1023  acoef = sbeta * ascale;
1024  bcoefr = salfar * bscale;
1025  bcoefi = 0.;
1026 
1027 /* Scale to avoid underflow */
1028 
1029  scale = 1.;
1030  lsa = absMACRO(sbeta) >= safmin && absMACRO(acoef) < small;
1031  lsb = absMACRO(salfar) >= safmin && absMACRO(bcoefr) < small;
1032  if (lsa) {
1033  scale = small / absMACRO(sbeta) * minMACRO(anorm,big);
1034  }
1035  if (lsb) {
1036 /* Computing MAX */
1037  d__1 = scale, d__2 = small / absMACRO(salfar) * minMACRO(bnorm,big);
1038  scale = maxMACRO(d__1,d__2);
1039  }
1040  if (lsa || lsb) {
1041 /* Computing MIN
1042  Computing MAX */
1043  d__3 = 1., d__4 = absMACRO(acoef), d__3 = maxMACRO(d__3,d__4), d__4
1044  = absMACRO(bcoefr);
1045  d__1 = scale, d__2 = 1. / (safmin * maxMACRO(d__3,d__4));
1046  scale = minMACRO(d__1,d__2);
1047  if (lsa) {
1048  acoef = ascale * (scale * sbeta);
1049  } else {
1050  acoef = scale * acoef;
1051  }
1052  if (lsb) {
1053  bcoefr = bscale * (scale * salfar);
1054  } else {
1055  bcoefr = scale * bcoefr;
1056  }
1057  }
1058  acoefa = absMACRO(acoef);
1059  bcoefa = absMACRO(bcoefr);
1060 
1061 /* First component is 1 */
1062 
1063  work[(*n << 1) + je] = 1.;
1064  xmax = 1.;
1065 
1066 /* Compute contribution from column JE of A and B to sum
1067  (See "Further Details", above.) */
1068 
1069  i__1 = je - 1;
1070  for (jr = 1; jr <= i__1; ++jr) {
1071  work[(*n << 1) + jr] = bcoefr * b_ref(jr, je) - acoef *
1072  a_ref(jr, je);
1073 /* L260: */
1074  }
1075  } else {
1076 
1077 /* Complex eigenvalue */
1078 
1079  d__1 = safmin * 100.;
1080  template_lapack_lag2(&a_ref(je - 1, je - 1), lda, &b_ref(je - 1, je - 1),
1081  ldb, &d__1, &acoef, &temp, &bcoefr, &temp2, &bcoefi);
1082  if (bcoefi == 0.) {
1083  *info = je - 1;
1084  return 0;
1085  }
1086 
1087 /* Scale to avoid over/underflow */
1088 
1089  acoefa = absMACRO(acoef);
1090  bcoefa = absMACRO(bcoefr) + absMACRO(bcoefi);
1091  scale = 1.;
1092  if (acoefa * ulp < safmin && acoefa >= safmin) {
1093  scale = safmin / ulp / acoefa;
1094  }
1095  if (bcoefa * ulp < safmin && bcoefa >= safmin) {
1096 /* Computing MAX */
1097  d__1 = scale, d__2 = safmin / ulp / bcoefa;
1098  scale = maxMACRO(d__1,d__2);
1099  }
1100  if (safmin * acoefa > ascale) {
1101  scale = ascale / (safmin * acoefa);
1102  }
1103  if (safmin * bcoefa > bscale) {
1104 /* Computing MIN */
1105  d__1 = scale, d__2 = bscale / (safmin * bcoefa);
1106  scale = minMACRO(d__1,d__2);
1107  }
1108  if (scale != 1.) {
1109  acoef = scale * acoef;
1110  acoefa = absMACRO(acoef);
1111  bcoefr = scale * bcoefr;
1112  bcoefi = scale * bcoefi;
1113  bcoefa = absMACRO(bcoefr) + absMACRO(bcoefi);
1114  }
1115 
1116 /* Compute first two components of eigenvector
1117  and contribution to sums */
1118 
1119  temp = acoef * a_ref(je, je - 1);
1120  temp2r = acoef * a_ref(je, je) - bcoefr * b_ref(je, je);
1121  temp2i = -bcoefi * b_ref(je, je);
1122  if (absMACRO(temp) >= absMACRO(temp2r) + absMACRO(temp2i)) {
1123  work[(*n << 1) + je] = 1.;
1124  work[*n * 3 + je] = 0.;
1125  work[(*n << 1) + je - 1] = -temp2r / temp;
1126  work[*n * 3 + je - 1] = -temp2i / temp;
1127  } else {
1128  work[(*n << 1) + je - 1] = 1.;
1129  work[*n * 3 + je - 1] = 0.;
1130  temp = acoef * a_ref(je - 1, je);
1131  work[(*n << 1) + je] = (bcoefr * b_ref(je - 1, je - 1) -
1132  acoef * a_ref(je - 1, je - 1)) / temp;
1133  work[*n * 3 + je] = bcoefi * b_ref(je - 1, je - 1) / temp;
1134  }
1135 
1136 /* Computing MAX */
1137  d__5 = (d__1 = work[(*n << 1) + je], absMACRO(d__1)) + (d__2 =
1138  work[*n * 3 + je], absMACRO(d__2)), d__6 = (d__3 = work[(*
1139  n << 1) + je - 1], absMACRO(d__3)) + (d__4 = work[*n * 3 +
1140  je - 1], absMACRO(d__4));
1141  xmax = maxMACRO(d__5,d__6);
1142 
1143 /* Compute contribution from columns JE and JE-1
1144  of A and B to the sums. */
1145 
1146  creala = acoef * work[(*n << 1) + je - 1];
1147  cimaga = acoef * work[*n * 3 + je - 1];
1148  crealb = bcoefr * work[(*n << 1) + je - 1] - bcoefi * work[*n
1149  * 3 + je - 1];
1150  cimagb = bcoefi * work[(*n << 1) + je - 1] + bcoefr * work[*n
1151  * 3 + je - 1];
1152  cre2a = acoef * work[(*n << 1) + je];
1153  cim2a = acoef * work[*n * 3 + je];
1154  cre2b = bcoefr * work[(*n << 1) + je] - bcoefi * work[*n * 3
1155  + je];
1156  cim2b = bcoefi * work[(*n << 1) + je] + bcoefr * work[*n * 3
1157  + je];
1158  i__1 = je - 2;
1159  for (jr = 1; jr <= i__1; ++jr) {
1160  work[(*n << 1) + jr] = -creala * a_ref(jr, je - 1) +
1161  crealb * b_ref(jr, je - 1) - cre2a * a_ref(jr, je)
1162  + cre2b * b_ref(jr, je);
1163  work[*n * 3 + jr] = -cimaga * a_ref(jr, je - 1) + cimagb *
1164  b_ref(jr, je - 1) - cim2a * a_ref(jr, je) +
1165  cim2b * b_ref(jr, je);
1166 /* L270: */
1167  }
1168  }
1169 
1170 /* Computing MAX */
1171  d__1 = ulp * acoefa * anorm, d__2 = ulp * bcoefa * bnorm, d__1 =
1172  maxMACRO(d__1,d__2);
1173  dmin__ = maxMACRO(d__1,safmin);
1174 
1175 /* Columnwise triangular solve of (a A - b B) x = 0 */
1176 
1177  il2by2 = FALSE_;
1178  for (j = je - nw; j >= 1; --j) {
1179 
1180 /* If a 2-by-2 block, is in position j-1:j, wait until
1181  next iteration to process it (when it will be j:j+1) */
1182 
1183  if (! il2by2 && j > 1) {
1184  if (a_ref(j, j - 1) != 0.) {
1185  il2by2 = TRUE_;
1186  goto L370;
1187  }
1188  }
1189  bdiag[0] = b_ref(j, j);
1190  if (il2by2) {
1191  na = 2;
1192  bdiag[1] = b_ref(j + 1, j + 1);
1193  } else {
1194  na = 1;
1195  }
1196 
1197 /* Compute x(j) (and x(j+1), if 2-by-2 block) */
1198 
1199  template_lapack_laln2(&c_false, &na, &nw, &dmin__, &acoef, &a_ref(j, j),
1200  lda, bdiag, &bdiag[1], &work[(*n << 1) + j], n, &
1201  bcoefr, &bcoefi, sum, &c__2, &scale, &temp, &iinfo);
1202  if (scale < 1.) {
1203 
1204  i__1 = nw - 1;
1205  for (jw = 0; jw <= i__1; ++jw) {
1206  i__2 = je;
1207  for (jr = 1; jr <= i__2; ++jr) {
1208  work[(jw + 2) * *n + jr] = scale * work[(jw + 2) *
1209  *n + jr];
1210 /* L280: */
1211  }
1212 /* L290: */
1213  }
1214  }
1215 /* Computing MAX */
1216  d__1 = scale * xmax;
1217  xmax = maxMACRO(d__1,temp);
1218 
1219  i__1 = nw;
1220  for (jw = 1; jw <= i__1; ++jw) {
1221  i__2 = na;
1222  for (ja = 1; ja <= i__2; ++ja) {
1223  work[(jw + 1) * *n + j + ja - 1] = sum_ref(ja, jw);
1224 /* L300: */
1225  }
1226 /* L310: */
1227  }
1228 
1229 /* w = w + x(j)*(a A(*,j) - b B(*,j) ) with scaling */
1230 
1231  if (j > 1) {
1232 
1233 /* Check whether scaling is necessary for sum. */
1234 
1235  xscale = 1. / maxMACRO(1.,xmax);
1236  temp = acoefa * work[j] + bcoefa * work[*n + j];
1237  if (il2by2) {
1238 /* Computing MAX */
1239  d__1 = temp, d__2 = acoefa * work[j + 1] + bcoefa *
1240  work[*n + j + 1];
1241  temp = maxMACRO(d__1,d__2);
1242  }
1243 /* Computing MAX */
1244  d__1 = maxMACRO(temp,acoefa);
1245  temp = maxMACRO(d__1,bcoefa);
1246  if (temp > bignum * xscale) {
1247 
1248  i__1 = nw - 1;
1249  for (jw = 0; jw <= i__1; ++jw) {
1250  i__2 = je;
1251  for (jr = 1; jr <= i__2; ++jr) {
1252  work[(jw + 2) * *n + jr] = xscale * work[(jw
1253  + 2) * *n + jr];
1254 /* L320: */
1255  }
1256 /* L330: */
1257  }
1258  xmax *= xscale;
1259  }
1260 
1261 /* Compute the contributions of the off-diagonals of
1262  column j (and j+1, if 2-by-2 block) of A and B to the
1263  sums. */
1264 
1265 
1266  i__1 = na;
1267  for (ja = 1; ja <= i__1; ++ja) {
1268  if (ilcplx) {
1269  creala = acoef * work[(*n << 1) + j + ja - 1];
1270  cimaga = acoef * work[*n * 3 + j + ja - 1];
1271  crealb = bcoefr * work[(*n << 1) + j + ja - 1] -
1272  bcoefi * work[*n * 3 + j + ja - 1];
1273  cimagb = bcoefi * work[(*n << 1) + j + ja - 1] +
1274  bcoefr * work[*n * 3 + j + ja - 1];
1275  i__2 = j - 1;
1276  for (jr = 1; jr <= i__2; ++jr) {
1277  work[(*n << 1) + jr] = work[(*n << 1) + jr] -
1278  creala * a_ref(jr, j + ja - 1) +
1279  crealb * b_ref(jr, j + ja - 1);
1280  work[*n * 3 + jr] = work[*n * 3 + jr] -
1281  cimaga * a_ref(jr, j + ja - 1) +
1282  cimagb * b_ref(jr, j + ja - 1);
1283 /* L340: */
1284  }
1285  } else {
1286  creala = acoef * work[(*n << 1) + j + ja - 1];
1287  crealb = bcoefr * work[(*n << 1) + j + ja - 1];
1288  i__2 = j - 1;
1289  for (jr = 1; jr <= i__2; ++jr) {
1290  work[(*n << 1) + jr] = work[(*n << 1) + jr] -
1291  creala * a_ref(jr, j + ja - 1) +
1292  crealb * b_ref(jr, j + ja - 1);
1293 /* L350: */
1294  }
1295  }
1296 /* L360: */
1297  }
1298  }
1299 
1300  il2by2 = FALSE_;
1301 L370:
1302  ;
1303  }
1304 
1305 /* Copy eigenvector to VR, back transforming if
1306  HOWMNY='B'. */
1307 
1308  ieig -= nw;
1309  if (ilback) {
1310 
1311  i__1 = nw - 1;
1312  for (jw = 0; jw <= i__1; ++jw) {
1313  i__2 = *n;
1314  for (jr = 1; jr <= i__2; ++jr) {
1315  work[(jw + 4) * *n + jr] = work[(jw + 2) * *n + 1] *
1316  vr_ref(jr, 1);
1317 /* L380: */
1318  }
1319 
1320 /* A series of compiler directives to defeat
1321  vectorization for the next loop */
1322 
1323 
1324  i__2 = je;
1325  for (jc = 2; jc <= i__2; ++jc) {
1326  i__3 = *n;
1327  for (jr = 1; jr <= i__3; ++jr) {
1328  work[(jw + 4) * *n + jr] += work[(jw + 2) * *n +
1329  jc] * vr_ref(jr, jc);
1330 /* L390: */
1331  }
1332 /* L400: */
1333  }
1334 /* L410: */
1335  }
1336 
1337  i__1 = nw - 1;
1338  for (jw = 0; jw <= i__1; ++jw) {
1339  i__2 = *n;
1340  for (jr = 1; jr <= i__2; ++jr) {
1341  vr_ref(jr, ieig + jw) = work[(jw + 4) * *n + jr];
1342 /* L420: */
1343  }
1344 /* L430: */
1345  }
1346 
1347  iend = *n;
1348  } else {
1349  i__1 = nw - 1;
1350  for (jw = 0; jw <= i__1; ++jw) {
1351  i__2 = *n;
1352  for (jr = 1; jr <= i__2; ++jr) {
1353  vr_ref(jr, ieig + jw) = work[(jw + 2) * *n + jr];
1354 /* L440: */
1355  }
1356 /* L450: */
1357  }
1358 
1359  iend = je;
1360  }
1361 
1362 /* Scale eigenvector */
1363 
1364  xmax = 0.;
1365  if (ilcplx) {
1366  i__1 = iend;
1367  for (j = 1; j <= i__1; ++j) {
1368 /* Computing MAX */
1369  d__3 = xmax, d__4 = (d__1 = vr_ref(j, ieig), absMACRO(d__1)) +
1370  (d__2 = vr_ref(j, ieig + 1), absMACRO(d__2));
1371  xmax = maxMACRO(d__3,d__4);
1372 /* L460: */
1373  }
1374  } else {
1375  i__1 = iend;
1376  for (j = 1; j <= i__1; ++j) {
1377 /* Computing MAX */
1378  d__2 = xmax, d__3 = (d__1 = vr_ref(j, ieig), absMACRO(d__1));
1379  xmax = maxMACRO(d__2,d__3);
1380 /* L470: */
1381  }
1382  }
1383 
1384  if (xmax > safmin) {
1385  xscale = 1. / xmax;
1386  i__1 = nw - 1;
1387  for (jw = 0; jw <= i__1; ++jw) {
1388  i__2 = iend;
1389  for (jr = 1; jr <= i__2; ++jr) {
1390  vr_ref(jr, ieig + jw) = xscale * vr_ref(jr, ieig + jw)
1391  ;
1392 /* L480: */
1393  }
1394 /* L490: */
1395  }
1396  }
1397 L500:
1398  ;
1399  }
1400  }
1401 
1402  return 0;
1403 
1404 /* End of DTGEVC */
1405 
1406 } /* dtgevc_ */
1407 
1408 #undef sum_ref
1409 #undef vr_ref
1410 #undef vl_ref
1411 #undef b_ref
1412 #undef a_ref
1413 #undef sumb_ref
1414 #undef suma_ref
1415 
1416 
1417 #endif