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
43
44
45template<class Treal>
46int 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 }
405L10:
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);
873L160:
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
937L220:
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_;
1307L370:
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 }
1403L500:
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
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:146
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:46
int integer
Definition: template_blas_common.h:40
#define absMACRO(x)
Definition: template_blas_common.h:47
#define minMACRO(a, b)
Definition: template_blas_common.h:46
#define maxMACRO(a, b)
Definition: template_blas_common.h:45
bool logical
Definition: template_blas_common.h:41
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
#define TRUE_
Definition: template_lapack_common.h:42
#define FALSE_
Definition: template_lapack_common.h:43
int template_lapack_labad(Treal *small, Treal *large)
Definition: template_lapack_labad.h:42
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
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
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
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition: template_lapack_lamch.h:202
#define vl_ref(a_1, a_2)
#define sum_ref(a_1, a_2)
#define a_ref(a_1, a_2)
#define b_ref(a_1, a_2)
#define sumb_ref(a_1, a_2)
#define suma_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:46
#define vr_ref(a_1, a_2)