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