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