Context Navigation

eispack.cxx

Visit:

Last change on this file was 5390, checked in by westram, 16 years ago
TAB-Ex
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 481.8 KB

Line
1	/* eispack.f -- translated by f2c (version 19950110).
2	You must link the resulting object file with the libraries:
3	-lf2c -lm (in that order)
4	*/
5
6	#ifdef __cplusplus
7	extern "C" {
8	#endif
9	#include "f2c.h"
10
11	/* Table of constant values */
12
13	static doublereal c_b141 = 1.;
14	static doublereal c_b550 = 0.;
15
16	/* Subroutine / int cdiv_(doublereal ar, doublereal ai, doublereal br,
17	doublereal bi, doublereal cr, doublereal *ci)
18	{
19	/* System generated locals */
20	doublereal d_1, d_2;
21
22	/* Local variables */
23	static doublereal s, ais, bis, ars, brs;
24
25
26	/* COMPLEX DIVISION, (CR,CI) = (AR,AI)/(BR,BI) */
27
28	s = abs(br) + abs(bi);
29	ars = *ar / s;
30	ais = *ai / s;
31	brs = *br / s;
32	bis = *bi / s;
33	/* Computing 2nd power */
34	d_1 = brs;
35	/* Computing 2nd power */
36	d_2 = bis;
37	s = d_1 * d_1 + d_2 * d_2;
38	cr = (ars brs + ais * bis) / s;
39	ci = (ais brs - ars * bis) / s;
40	return 0;
41	} /* cdiv_ */
42
43	/* Subroutine / int csroot_(doublereal xr, doublereal xi, doublereal yr,
44	doublereal *yi)
45	{
46	/* Builtin functions */
47	double sqrt(doublereal);
48
49	/* Local variables */
50	static doublereal s, ti, tr;
51	extern doublereal pythag_(doublereal , doublereal );
52
53
54	/* (YR,YI) = COMPLEX DSQRT(XR,XI) */
55	/* BRANCH CHOSEN SO THAT YR .GE. 0.0 AND SIGN(YI) .EQ. SIGN(XI) */
56
57	tr = *xr;
58	ti = *xi;
59	s = sqrt((pythag_(&tr, &ti) + abs(tr)) * .5);
60	if (tr >= 0.) {
61	*yr = s;
62	}
63	if (ti < 0.) {
64	s = -s;
65	}
66	if (tr <= 0.) {
67	*yi = s;
68	}
69	if (tr < 0.) {
70	yr = ti / yi * .5;
71	}
72	if (tr > 0.) {
73	yi = ti / yr * .5;
74	}
75	return 0;
76	} /* csroot_ */
77
78	doublereal epslon_(doublereal *x)
79	{
80	/* System generated locals */
81	doublereal ret_val, d_1;
82
83	/* Local variables */
84	static doublereal a, b, c, eps;
85
86
87	/* ESTIMATE UNIT ROUNDOFF IN QUANTITIES OF SIZE X. */
88
89
90	/* THIS PROGRAM SHOULD FUNCTION PROPERLY ON ALL SYSTEMS */
91	/* SATISFYING THE FOLLOWING TWO ASSUMPTIONS, */
92	/* 1. THE BASE USED IN REPRESENTING FLOATING POINT */
93	/* NUMBERS IS NOT A POWER OF THREE. */
94	/* 2. THE QUANTITY A IN STATEMENT 10 IS REPRESENTED TO */
95	/* THE ACCURACY USED IN FLOATING POINT VARIABLES */
96	/* THAT ARE STORED IN MEMORY. */
97	/* THE STATEMENT NUMBER 10 AND THE GO TO 10 ARE INTENDED TO */
98	/* FORCE OPTIMIZING COMPILERS TO GENERATE CODE SATISFYING */
99	/* ASSUMPTION 2. */
100	/* UNDER THESE ASSUMPTIONS, IT SHOULD BE TRUE THAT, */
101	/* A IS NOT EXACTLY EQUAL TO FOUR-THIRDS, */
102	/* B HAS A ZERO FOR ITS LAST BIT OR DIGIT, */
103	/* C IS NOT EXACTLY EQUAL TO ONE, */
104	/* EPS MEASURES THE SEPARATION OF 1.0 FROM */
105	/* THE NEXT LARGER FLOATING POINT NUMBER. */
106	/* THE DEVELOPERS OF EISPACK WOULD APPRECIATE BEING INFORMED */
107	/* ABOUT ANY SYSTEMS WHERE THESE ASSUMPTIONS DO NOT HOLD. */
108
109	/* THIS VERSION DATED 4/6/83. */
110
111	a = 1.3333333333333333;
112	L10:
113	b = a - 1.;
114	c = b + b + b;
115	eps = (d_1 = c - 1., abs(d_1));
116	if (eps == 0.) {
117	goto L10;
118	}
119	ret_val = eps * abs(*x);
120	return ret_val;
121	} /* epslon_ */
122
123	doublereal pythag_(doublereal a, doublereal b)
124	{
125	/* System generated locals */
126	doublereal ret_val, d_1, d_2, d_3;
127
128	/* Local variables */
129	static doublereal p, r, s, t, u;
130
131
132	/* FINDS DSQRT(A2+B2) WITHOUT OVERFLOW OR DESTRUCTIVE UNDERFLOW */
133
134	/* Computing MAX */
135	d_1 = abs(a), d_2 = abs(b);
136	p = max(d_1,d_2);
137	if (p == 0.) {
138	goto L20;
139	}
140	/* Computing MIN */
141	d_2 = abs(a), d_3 = abs(b);
142	/* Computing 2nd power */
143	d_1 = min(d_2,d_3) / p;
144	r = d_1 * d_1;
145	L10:
146	t = r + 4.;
147	if (t == 4.) {
148	goto L20;
149	}
150	s = r / t;
151	u = s * 2. + 1.;
152	p = u * p;
153	/* Computing 2nd power */
154	d_1 = s / u;
155	r = d_1 * d_1 * r;
156	goto L10;
157	L20:
158	ret_val = p;
159	return ret_val;
160	} /* pythag_ */
161
162	/* Subroutine / int bakvec_(integer nm, integer n, doublereal t,
163	doublereal e, integer m, doublereal z, integer ierr)
164	{
165	/* System generated locals */
166	integer t_dim1, t_offset, z_dim1, z_offset, i_1, i_2;
167
168	/* Local variables */
169	static integer i, j;
170
171
172
173	/* THIS SUBROUTINE FORMS THE EIGENVECTORS OF A NONSYMMETRIC */
174	/* TRIDIAGONAL MATRIX BY BACK TRANSFORMING THOSE OF THE */
175	/* CORRESPONDING SYMMETRIC MATRIX DETERMINED BY FIGI. */
176
177	/* ON INPUT */
178
179	/* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
180	/* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
181	/* DIMENSION STATEMENT. */
182
183	/* N IS THE ORDER OF THE MATRIX. */
184
185	/* T CONTAINS THE NONSYMMETRIC MATRIX. ITS SUBDIAGONAL IS */
186	/* STORED IN THE LAST N-1 POSITIONS OF THE FIRST COLUMN, */
187	/* ITS DIAGONAL IN THE N POSITIONS OF THE SECOND COLUMN, */
188	/* AND ITS SUPERDIAGONAL IN THE FIRST N-1 POSITIONS OF */
189	/* THE THIRD COLUMN. T(1,1) AND T(N,3) ARE ARBITRARY. */
190
191	/* E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE SYMMETRIC */
192	/* MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY. */
193
194	/* M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED. */
195
196	/* Z CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED */
197	/* IN ITS FIRST M COLUMNS. */
198
199	/* ON OUTPUT */
200
201	/* T IS UNALTERED. */
202
203	/* E IS DESTROYED. */
204
205	/* Z CONTAINS THE TRANSFORMED EIGENVECTORS */
206	/* IN ITS FIRST M COLUMNS. */
207
208	/* IERR IS SET TO */
209	/* ZERO FOR NORMAL RETURN, */
210	/* 2*N+I IF E(I) IS ZERO WITH T(I,1) OR T(I-1,3) NON-ZERO.
211	*/
212	/* IN THIS CASE, THE SYMMETRIC MATRIX IS NOT SIMILAR
213	*/
214	/* TO THE ORIGINAL MATRIX, AND THE EIGENVECTORS */
215	/* CANNOT BE FOUND BY THIS PROGRAM. */
216
217	/* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
218	/* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
219	*/
220
221	/* THIS VERSION DATED AUGUST 1983. */
222
223	/* ------------------------------------------------------------------
224	*/
225
226	/* Parameter adjustments */
227	t_dim1 = *nm;
228	t_offset = t_dim1 + 1;
229	t -= t_offset;
230	--e;
231	z_dim1 = *nm;
232	z_offset = z_dim1 + 1;
233	z -= z_offset;
234
235	/* Function Body */
236	*ierr = 0;
237	if (*m == 0) {
238	goto L1001;
239	}
240	e[1] = 1.;
241	if (*n == 1) {
242	goto L1001;
243	}
244
245	i_1 = *n;
246	for (i = 2; i <= i_1; ++i) {
247	if (e[i] != 0.) {
248	goto L80;
249	}
250	if (t[i + t_dim1] != 0. \|\| t[i - 1 + t_dim1 * 3] != 0.) {
251	goto L1000;
252	}
253	e[i] = 1.;
254	goto L100;
255	L80:
256	e[i] = e[i - 1] * e[i] / t[i - 1 + t_dim1 * 3];
257	L100:
258	;
259	}
260
261	i_1 = *m;
262	for (j = 1; j <= i_1; ++j) {
263
264	i_2 = *n;
265	for (i = 2; i <= i_2; ++i) {
266	z[i + j * z_dim1] *= e[i];
267	/* L120: */
268	}
269	}
270
271	goto L1001;
272	/* .......... SET ERROR -- EIGENVECTORS CANNOT BE */
273	/* FOUND BY THIS PROGRAM .......... */
274	L1000:
275	ierr = (n << 1) + i;
276	L1001:
277	return 0;
278	} /* bakvec_ */
279
280	/* Subroutine / int balanc_(integer nm, integer n, doublereal a, integer *
281	low, integer igh, doublereal scale)
282	{
283	/* System generated locals */
284	integer a_dim1, a_offset, i_1, i_2;
285	doublereal d_1;
286
287	/* Local variables */
288	static integer iexc;
289	static doublereal c, f, g;
290	static integer i, j, k, l, m;
291	static doublereal r, s, radix, b2;
292	static integer jj;
293	static logical noconv;
294
295
296
297	/* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BALANCE, */
298	/* NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH. */
299	/* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971). */
300
301	/* THIS SUBROUTINE BALANCES A REAL MATRIX AND ISOLATES */
302	/* EIGENVALUES WHENEVER POSSIBLE. */
303
304	/* ON INPUT */
305
306	/* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
307	/* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
308	/* DIMENSION STATEMENT. */
309
310	/* N IS THE ORDER OF THE MATRIX. */
311
312	/* A CONTAINS THE INPUT MATRIX TO BE BALANCED. */
313
314	/* ON OUTPUT */
315
316	/* A CONTAINS THE BALANCED MATRIX. */
317
318	/* LOW AND IGH ARE TWO INTEGERS SUCH THAT A(I,J) */
319	/* IS EQUAL TO ZERO IF */
320	/* (1) I IS GREATER THAN J AND */
321	/* (2) J=1,...,LOW-1 OR I=IGH+1,...,N. */
322
323	/* SCALE CONTAINS INFORMATION DETERMINING THE */
324	/* PERMUTATIONS AND SCALING FACTORS USED. */
325
326	/* SUPPOSE THAT THE PRINCIPAL SUBMATRIX IN ROWS LOW THROUGH IGH */
327	/* HAS BEEN BALANCED, THAT P(J) DENOTES THE INDEX INTERCHANGED */
328	/* WITH J DURING THE PERMUTATION STEP, AND THAT THE ELEMENTS */
329	/* OF THE DIAGONAL MATRIX USED ARE DENOTED BY D(I,J). THEN */
330	/* SCALE(J) = P(J), FOR J = 1,...,LOW-1 */
331	/* = D(J,J), J = LOW,...,IGH */
332	/* = P(J) J = IGH+1,...,N. */
333	/* THE ORDER IN WHICH THE INTERCHANGES ARE MADE IS N TO IGH+1, */
334	/* THEN 1 TO LOW-1. */
335
336	/* NOTE THAT 1 IS RETURNED FOR IGH IF IGH IS ZERO FORMALLY. */
337
338	/* THE ALGOL PROCEDURE EXC CONTAINED IN BALANCE APPEARS IN */
339	/* BALANC IN LINE. (NOTE THAT THE ALGOL ROLES OF IDENTIFIERS */
340	/* K,L HAVE BEEN REVERSED.) */
341
342	/* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
343	/* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
344	*/
345
346	/* THIS VERSION DATED AUGUST 1983. */
347
348	/* ------------------------------------------------------------------
349	*/
350
351	/* Parameter adjustments */
352	--scale;
353	a_dim1 = *nm;
354	a_offset = a_dim1 + 1;
355	a -= a_offset;
356
357	/* Function Body */
358	radix = 16.;
359
360	b2 = radix * radix;
361	k = 1;
362	l = *n;
363	goto L100;
364	/* .......... IN-LINE PROCEDURE FOR ROW AND */
365	/* COLUMN EXCHANGE .......... */
366	L20:
367	scale[m] = (doublereal) j;
368	if (j == m) {
369	goto L50;
370	}
371
372	i_1 = l;
373	for (i = 1; i <= i_1; ++i) {
374	f = a[i + j * a_dim1];
375	a[i + j * a_dim1] = a[i + m * a_dim1];
376	a[i + m * a_dim1] = f;
377	/* L30: */
378	}
379
380	i_1 = *n;
381	for (i = k; i <= i_1; ++i) {
382	f = a[j + i * a_dim1];
383	a[j + i * a_dim1] = a[m + i * a_dim1];
384	a[m + i * a_dim1] = f;
385	/* L40: */
386	}
387
388	L50:
389	switch (iexc) {
390	case 1: goto L80;
391	case 2: goto L130;
392	}
393	/* .......... SEARCH FOR ROWS ISOLATING AN EIGENVALUE */
394	/* AND PUSH THEM DOWN .......... */
395	L80:
396	if (l == 1) {
397	goto L280;
398	}
399	--l;
400	/* .......... FOR J=L STEP -1 UNTIL 1 DO -- .......... */
401	L100:
402	i_1 = l;
403	for (jj = 1; jj <= i_1; ++jj) {
404	j = l + 1 - jj;
405
406	i_2 = l;
407	for (i = 1; i <= i_2; ++i) {
408	if (i == j) {
409	goto L110;
410	}
411	if (a[j + i * a_dim1] != 0.) {
412	goto L120;
413	}
414	L110:
415	;
416	}
417
418	m = l;
419	iexc = 1;
420	goto L20;
421	L120:
422	;
423	}
424
425	goto L140;
426	/* .......... SEARCH FOR COLUMNS ISOLATING AN EIGENVALUE */
427	/* AND PUSH THEM LEFT .......... */
428	L130:
429	++k;
430
431	L140:
432	i_1 = l;
433	for (j = k; j <= i_1; ++j) {
434
435	i_2 = l;
436	for (i = k; i <= i_2; ++i) {
437	if (i == j) {
438	goto L150;
439	}
440	if (a[i + j * a_dim1] != 0.) {
441	goto L170;
442	}
443	L150:
444	;
445	}
446
447	m = k;
448	iexc = 2;
449	goto L20;
450	L170:
451	;
452	}
453	/* .......... NOW BALANCE THE SUBMATRIX IN ROWS K TO L .......... */
454	i_1 = l;
455	for (i = k; i <= i_1; ++i) {
456	/* L180: */
457	scale[i] = 1.;
458	}
459	/* .......... ITERATIVE LOOP FOR NORM REDUCTION .......... */
460	L190:
461	noconv = FALSE_;
462
463	i_1 = l;
464	for (i = k; i <= i_1; ++i) {
465	c = 0.;
466	r = 0.;
467
468	i_2 = l;
469	for (j = k; j <= i_2; ++j) {
470	if (j == i) {
471	goto L200;
472	}
473	c += (d_1 = a[j + i * a_dim1], abs(d_1));
474	r += (d_1 = a[i + j * a_dim1], abs(d_1));
475	L200:
476	;
477	}
478	/* .......... GUARD AGAINST ZERO C OR R DUE TO UNDERFLOW .........
479	. */
480	if (c == 0. \|\| r == 0.) {
481	goto L270;
482	}
483	g = r / radix;
484	f = 1.;
485	s = c + r;
486	L210:
487	if (c >= g) {
488	goto L220;
489	}
490	f *= radix;
491	c *= b2;
492	goto L210;
493	L220:
494	g = r * radix;
495	L230:
496	if (c < g) {
497	goto L240;
498	}
499	f /= radix;
500	c /= b2;
501	goto L230;
502	/* .......... NOW BALANCE .......... */
503	L240:
504	if ((c + r) / f >= s * .95) {
505	goto L270;
506	}
507	g = 1. / f;
508	scale[i] *= f;
509	noconv = TRUE_;
510
511	i_2 = *n;
512	for (j = k; j <= i_2; ++j) {
513	/* L250: */
514	a[i + j * a_dim1] *= g;
515	}
516
517	i_2 = l;
518	for (j = 1; j <= i_2; ++j) {
519	/* L260: */
520	a[j + i * a_dim1] *= f;
521	}
522
523	L270:
524	;
525	}
526
527	if (noconv) {
528	goto L190;
529	}
530
531	L280:
532	*low = k;
533	*igh = l;
534	return 0;
535	} /* balanc_ */
536
537	/* Subroutine / int balbak_(integer nm, integer n, integer low, integer *
538	igh, doublereal scale, integer m, doublereal *z)
539	{
540	/* System generated locals */
541	integer z_dim1, z_offset, i_1, i_2;
542
543	/* Local variables */
544	static integer i, j, k;
545	static doublereal s;
546	static integer ii;
547
548
549
550	/* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BALBAK, */
551	/* NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH. */
552	/* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971). */
553
554	/* THIS SUBROUTINE FORMS THE EIGENVECTORS OF A REAL GENERAL */
555	/* MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING */
556	/* BALANCED MATRIX DETERMINED BY BALANC. */
557
558	/* ON INPUT */
559
560	/* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
561	/* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
562	/* DIMENSION STATEMENT. */
563
564	/* N IS THE ORDER OF THE MATRIX. */
565
566	/* LOW AND IGH ARE INTEGERS DETERMINED BY BALANC. */
567
568	/* SCALE CONTAINS INFORMATION DETERMINING THE PERMUTATIONS */
569	/* AND SCALING FACTORS USED BY BALANC. */
570
571	/* M IS THE NUMBER OF COLUMNS OF Z TO BE BACK TRANSFORMED. */
572
573	/* Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGEN- */
574	/* VECTORS TO BE BACK TRANSFORMED IN ITS FIRST M COLUMNS. */
575
576	/* ON OUTPUT */
577
578	/* Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE */
579	/* TRANSFORMED EIGENVECTORS IN ITS FIRST M COLUMNS. */
580
581	/* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
582	/* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
583	*/
584
585	/* THIS VERSION DATED AUGUST 1983. */
586
587	/* ------------------------------------------------------------------
588	*/
589
590	/* Parameter adjustments */
591	--scale;
592	z_dim1 = *nm;
593	z_offset = z_dim1 + 1;
594	z -= z_offset;
595
596	/* Function Body */
597	if (*m == 0) {
598	goto L200;
599	}
600	if (igh == low) {
601	goto L120;
602	}
603
604	i_1 = *igh;
605	for (i = *low; i <= i_1; ++i) {
606	s = scale[i];
607	/* .......... LEFT HAND EIGENVECTORS ARE BACK TRANSFORMED */
608	/* IF THE FOREGOING STATEMENT IS REPLACED BY */
609	/* S=1.0D0/SCALE(I). .......... */
610	i_2 = *m;
611	for (j = 1; j <= i_2; ++j) {
612	/* L100: */
613	z[i + j * z_dim1] *= s;
614	}
615
616	/* L110: */
617	}
618	/* ......... FOR I=LOW-1 STEP -1 UNTIL 1, */
619	/* IGH+1 STEP 1 UNTIL N DO -- .......... */
620	L120:
621	i_1 = *n;
622	for (ii = 1; ii <= i_1; ++ii) {
623	i = ii;
624	if (i >= low && i <= igh) {
625	goto L140;
626	}
627	if (i < *low) {
628	i = *low - ii;
629	}
630	k = (integer) scale[i];
631	if (k == i) {
632	goto L140;
633	}
634
635	i_2 = *m;
636	for (j = 1; j <= i_2; ++j) {
637	s = z[i + j * z_dim1];
638	z[i + j * z_dim1] = z[k + j * z_dim1];
639	z[k + j * z_dim1] = s;
640	/* L130: */
641	}
642
643	L140:
644	;
645	}
646
647	L200:
648	return 0;
649	} /* balbak_ */
650
651	/* Subroutine / int bandr_(integer nm, integer n, integer mb, doublereal *
652	a, doublereal d, doublereal e, doublereal e2, logical matz,
653	doublereal *z)
654	{
655	/* System generated locals */
656	integer a_dim1, a_offset, z_dim1, z_offset, i_1, i_2, i_3, i_4, i_5,
657	i_6;
658	doublereal d_1;
659
660	/* Builtin functions */
661	double sqrt(doublereal);
662
663	/* Local variables */
664	static doublereal dmin_;
665	static integer maxl, maxr;
666	static doublereal g;
667	static integer j, k, l, r;
668	static doublereal u, b1, b2, c2, f1, f2;
669	static integer i1, i2, j1, j2, m1, n2, r1;
670	static doublereal s2;
671	static integer kr, mr;
672	static doublereal dminrt;
673	static integer ugl;
674
675
676
677	/* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BANDRD, */
678	/* NUM. MATH. 12, 231-241(1968) BY SCHWARZ. */
679	/* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 273-283(1971). */
680
681	/* THIS SUBROUTINE REDUCES A REAL SYMMETRIC BAND MATRIX */
682	/* TO A SYMMETRIC TRIDIAGONAL MATRIX USING AND OPTIONALLY */
683	/* ACCUMULATING ORTHOGONAL SIMILARITY TRANSFORMATIONS. */
684
685	/* ON INPUT */
686
687	/* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
688	/* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
689	/* DIMENSION STATEMENT. */
690
691	/* N IS THE ORDER OF THE MATRIX. */
692
693	/* MB IS THE (HALF) BAND WIDTH OF THE MATRIX, DEFINED AS THE */
694	/* NUMBER OF ADJACENT DIAGONALS, INCLUDING THE PRINCIPAL */
695	/* DIAGONAL, REQUIRED TO SPECIFY THE NON-ZERO PORTION OF THE */
696	/* LOWER TRIANGLE OF THE MATRIX. */
697
698	/* A CONTAINS THE LOWER TRIANGLE OF THE SYMMETRIC BAND INPUT */
699	/* MATRIX STORED AS AN N BY MB ARRAY. ITS LOWEST SUBDIAGONAL */
700	/* IS STORED IN THE LAST N+1-MB POSITIONS OF THE FIRST COLUMN, */
701	/* ITS NEXT SUBDIAGONAL IN THE LAST N+2-MB POSITIONS OF THE */
702	/* SECOND COLUMN, FURTHER SUBDIAGONALS SIMILARLY, AND FINALLY */
703	/* ITS PRINCIPAL DIAGONAL IN THE N POSITIONS OF THE LAST COLUMN.
704	*/
705	/* CONTENTS OF STORAGES NOT PART OF THE MATRIX ARE ARBITRARY. */
706
707	/* MATZ SHOULD BE SET TO .TRUE. IF THE TRANSFORMATION MATRIX IS */
708	/* TO BE ACCUMULATED, AND TO .FALSE. OTHERWISE. */
709
710	/* ON OUTPUT */
711
712	/* A HAS BEEN DESTROYED, EXCEPT FOR ITS LAST TWO COLUMNS WHICH */
713	/* CONTAIN A COPY OF THE TRIDIAGONAL MATRIX. */
714
715	/* D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX. */
716
717	/* E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL */
718	/* MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS SET TO ZERO. */
719
720	/* E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. */
721	/* E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED. */
722
723	/* Z CONTAINS THE ORTHOGONAL TRANSFORMATION MATRIX PRODUCED IN */
724	/* THE REDUCTION IF MATZ HAS BEEN SET TO .TRUE. OTHERWISE, Z */
725	/* IS NOT REFERENCED. */
726
727	/* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
728	/* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
729	*/
730
731	/* THIS VERSION DATED AUGUST 1983. */
732
733	/* ------------------------------------------------------------------
734	*/
735
736	/* Parameter adjustments */
737	z_dim1 = *nm;
738	z_offset = z_dim1 + 1;
739	z -= z_offset;
740	--e2;
741	--e;
742	--d;
743	a_dim1 = *nm;
744	a_offset = a_dim1 + 1;
745	a -= a_offset;
746
747	/* Function Body */
748	dmin_ = 5.4210108624275222e-20;
749	dminrt = 2.3283064365386963e-10;
750	/* .......... INITIALIZE DIAGONAL SCALING MATRIX .......... */
751	i_1 = *n;
752	for (j = 1; j <= i_1; ++j) {
753	/* L30: */
754	d[j] = 1.;
755	}
756
757	if (! (*matz)) {
758	goto L60;
759	}
760
761	i_1 = *n;
762	for (j = 1; j <= i_1; ++j) {
763
764	i_2 = *n;
765	for (k = 1; k <= i_2; ++k) {
766	/* L40: */
767	z[j + k * z_dim1] = 0.;
768	}
769
770	z[j + j * z_dim1] = 1.;
771	/* L50: */
772	}
773
774	L60:
775	m1 = *mb - 1;
776	if ((i_1 = m1 - 1) < 0) {
777	goto L900;
778	} else if (i_1 == 0) {
779	goto L800;
780	} else {
781	goto L70;
782	}
783	L70:
784	n2 = *n - 2;
785
786	i_1 = n2;
787	for (k = 1; k <= i_1; ++k) {
788	/* Computing MIN */
789	i_2 = m1, i_3 = *n - k;
790	maxr = min(i_2,i_3);
791	/* .......... FOR R=MAXR STEP -1 UNTIL 2 DO -- .......... */
792	i_2 = maxr;
793	for (r1 = 2; r1 <= i_2; ++r1) {
794	r = maxr + 2 - r1;
795	kr = k + r;
796	mr = *mb - r;
797	g = a[kr + mr * a_dim1];
798	a[kr - 1 + a_dim1] = a[kr - 1 + (mr + 1) * a_dim1];
799	ugl = k;
800
801	i_3 = *n;
802	i_4 = m1;
803	for (j = kr; i_4 < 0 ? j >= i_3 : j <= i_3; j += i_4) {
804	j1 = j - 1;
805	j2 = j1 - 1;
806	if (g == 0.) {
807	goto L600;
808	}
809	b1 = a[j1 + a_dim1] / g;
810	b2 = b1 * d[j1] / d[j];
811	s2 = 1. / (b1 * b2 + 1.);
812	if (s2 >= .5) {
813	goto L450;
814	}
815	b1 = g / a[j1 + a_dim1];
816	b2 = b1 * d[j] / d[j1];
817	c2 = 1. - s2;
818	d[j1] = c2 * d[j1];
819	d[j] = c2 * d[j];
820	f1 = a[j + m1 * a_dim1] * 2.;
821	f2 = b1 * a[j1 + mb a_dim1];
822	a[j + m1 * a_dim1] = -b2 * (b1 * a[j + m1 * a_dim1] - a[j + *
823	mb * a_dim1]) - f2 + a[j + m1 * a_dim1];
824	a[j1 + mb a_dim1] = b2 * (b2 * a[j + mb a_dim1] + f1) +
825	a[j1 + mb a_dim1];
826	a[j + mb a_dim1] = b1 * (f2 - f1) + a[j + mb a_dim1];
827
828	i_5 = j2;
829	for (l = ugl; l <= i_5; ++l) {
830	i2 = *mb - j + l;
831	u = a[j1 + (i2 + 1) * a_dim1] + b2 * a[j + i2 * a_dim1];
832	a[j + i2 * a_dim1] = -b1 * a[j1 + (i2 + 1) * a_dim1] + a[
833	j + i2 * a_dim1];
834	a[j1 + (i2 + 1) * a_dim1] = u;
835	/* L200: */
836	}
837
838	ugl = j;
839	a[j1 + a_dim1] += b2 * g;
840	if (j == *n) {
841	goto L350;
842	}
843	/* Computing MIN */
844	i_5 = m1, i_6 = *n - j1;
845	maxl = min(i_5,i_6);
846
847	i_5 = maxl;
848	for (l = 2; l <= i_5; ++l) {
849	i1 = j1 + l;
850	i2 = *mb - l;
851	u = a[i1 + i2 * a_dim1] + b2 * a[i1 + (i2 + 1) * a_dim1];
852	a[i1 + (i2 + 1) * a_dim1] = -b1 * a[i1 + i2 * a_dim1] + a[
853	i1 + (i2 + 1) * a_dim1];
854	a[i1 + i2 * a_dim1] = u;
855	/* L300: */
856	}
857
858	i1 = j + m1;
859	if (i1 > *n) {
860	goto L350;
861	}
862	g = b2 * a[i1 + a_dim1];
863	L350:
864	if (! (*matz)) {
865	goto L500;
866	}
867
868	i_5 = *n;
869	for (l = 1; l <= i_5; ++l) {
870	u = z[l + j1 * z_dim1] + b2 * z[l + j * z_dim1];
871	z[l + j * z_dim1] = -b1 * z[l + j1 * z_dim1] + z[l + j *
872	z_dim1];
873	z[l + j1 * z_dim1] = u;
874	/* L400: */
875	}
876
877	goto L500;
878
879	L450:
880	u = d[j1];
881	d[j1] = s2 * d[j];
882	d[j] = s2 * u;
883	f1 = a[j + m1 * a_dim1] * 2.;
884	f2 = b1 * a[j + mb a_dim1];
885	u = b1 * (f2 - f1) + a[j1 + mb a_dim1];
886	a[j + m1 * a_dim1] = b2 * (b1 * a[j + m1 * a_dim1] - a[j1 + *
887	mb * a_dim1]) + f2 - a[j + m1 * a_dim1];
888	a[j1 + mb a_dim1] = b2 * (b2 * a[j1 + mb a_dim1] + f1)
889	+ a[j + mb a_dim1];
890	a[j + mb a_dim1] = u;
891
892	i_5 = j2;
893	for (l = ugl; l <= i_5; ++l) {
894	i2 = *mb - j + l;
895	u = b2 * a[j1 + (i2 + 1) * a_dim1] + a[j + i2 * a_dim1];
896	a[j + i2 * a_dim1] = -a[j1 + (i2 + 1) * a_dim1] + b1 * a[
897	j + i2 * a_dim1];
898	a[j1 + (i2 + 1) * a_dim1] = u;
899	/* L460: */
900	}
901
902	ugl = j;
903	a[j1 + a_dim1] = b2 * a[j1 + a_dim1] + g;
904	if (j == *n) {
905	goto L480;
906	}
907	/* Computing MIN */
908	i_5 = m1, i_6 = *n - j1;
909	maxl = min(i_5,i_6);
910
911	i_5 = maxl;
912	for (l = 2; l <= i_5; ++l) {
913	i1 = j1 + l;
914	i2 = *mb - l;
915	u = b2 * a[i1 + i2 * a_dim1] + a[i1 + (i2 + 1) * a_dim1];
916	a[i1 + (i2 + 1) * a_dim1] = -a[i1 + i2 * a_dim1] + b1 * a[
917	i1 + (i2 + 1) * a_dim1];
918	a[i1 + i2 * a_dim1] = u;
919	/* L470: */
920	}
921
922	i1 = j + m1;
923	if (i1 > *n) {
924	goto L480;
925	}
926	g = a[i1 + a_dim1];
927	a[i1 + a_dim1] = b1 * a[i1 + a_dim1];
928	L480:
929	if (! (*matz)) {
930	goto L500;
931	}
932
933	i_5 = *n;
934	for (l = 1; l <= i_5; ++l) {
935	u = b2 * z[l + j1 * z_dim1] + z[l + j * z_dim1];
936	z[l + j * z_dim1] = -z[l + j1 * z_dim1] + b1 * z[l + j *
937	z_dim1];
938	z[l + j1 * z_dim1] = u;
939	/* L490: */
940	}
941
942	L500:
943	;
944	}
945
946	L600:
947	;
948	}
949
950	if (k % 64 != 0) {
951	goto L700;
952	}
953	/* .......... RESCALE TO AVOID UNDERFLOW OR OVERFLOW .......... */
954	i_2 = *n;
955	for (j = k; j <= i_2; ++j) {
956	if (d[j] >= dmin_) {
957	goto L650;
958	}
959	/* Computing MAX */
960	i_4 = 1, i_3 = *mb + 1 - j;
961	maxl = max(i_4,i_3);
962
963	i_4 = m1;
964	for (l = maxl; l <= i_4; ++l) {
965	/* L610: */
966	a[j + l * a_dim1] = dminrt * a[j + l * a_dim1];
967	}
968
969	if (j == *n) {
970	goto L630;
971	}
972	/* Computing MIN */
973	i_4 = m1, i_3 = *n - j;
974	maxl = min(i_4,i_3);
975
976	i_4 = maxl;
977	for (l = 1; l <= i_4; ++l) {
978	i1 = j + l;
979	i2 = *mb - l;
980	a[i1 + i2 * a_dim1] = dminrt * a[i1 + i2 * a_dim1];
981	/* L620: */
982	}
983
984	L630:
985	if (! (*matz)) {
986	goto L645;
987	}
988
989	i_4 = *n;
990	for (l = 1; l <= i_4; ++l) {
991	/* L640: */
992	z[l + j * z_dim1] = dminrt * z[l + j * z_dim1];
993	}
994
995	L645:
996	a[j + mb a_dim1] = dmin_ * a[j + mb a_dim1];
997	d[j] /= dmin_;
998	L650:
999	;
1000	}
1001
1002	L700:
1003	;
1004	}
1005	/* .......... FORM SQUARE ROOT OF SCALING MATRIX .......... */
1006	L800:
1007	i_1 = *n;
1008	for (j = 2; j <= i_1; ++j) {
1009	/* L810: */
1010	e[j] = sqrt(d[j]);
1011	}
1012
1013	if (! (*matz)) {
1014	goto L840;
1015	}
1016
1017	i_1 = *n;
1018	for (j = 1; j <= i_1; ++j) {
1019
1020	i_2 = *n;
1021	for (k = 2; k <= i_2; ++k) {
1022	/* L820: */
1023	z[j + k * z_dim1] = e[k] * z[j + k * z_dim1];
1024	}
1025
1026	/* L830: */
1027	}
1028
1029	L840:
1030	u = 1.;
1031
1032	i_1 = *n;
1033	for (j = 2; j <= i_1; ++j) {
1034	a[j + m1 * a_dim1] = u * e[j] * a[j + m1 * a_dim1];
1035	u = e[j];
1036	/* Computing 2nd power */
1037	d_1 = a[j + m1 * a_dim1];
1038	e2[j] = d_1 * d_1;
1039	a[j + mb a_dim1] = d[j] * a[j + mb a_dim1];
1040	d[j] = a[j + mb a_dim1];
1041	e[j] = a[j + m1 * a_dim1];
1042	/* L850: */
1043	}
1044
1045	d[1] = a[mb a_dim1 + 1];
1046	e[1] = 0.;
1047	e2[1] = 0.;
1048	goto L1001;
1049
1050	L900:
1051	i_1 = *n;
1052	for (j = 1; j <= i_1; ++j) {
1053	d[j] = a[j + mb a_dim1];
1054	e[j] = 0.;
1055	e2[j] = 0.;
1056	/* L950: */
1057	}
1058
1059	L1001:
1060	return 0;
1061	} /* bandr_ */
1062
1063	/* Subroutine / int bandv_(integer nm, integer n, integer mbw, doublereal
1064	a, doublereal e21, integer m, doublereal w, doublereal *z,
1065	integer ierr, integer /nv/, doublereal rv, doublereal rv6)
1066	{
1067	/* System generated locals */
1068	integer a_dim1, a_offset, z_dim1, z_offset, i_1, i_2, i_3, i_4, i_5;
1069	doublereal d_1;
1070
1071	/* Builtin functions */
1072	double sqrt(doublereal), d_sign(doublereal , doublereal );
1073
1074	/* Local variables */
1075	static integer maxj, maxk;
1076	static doublereal norm;
1077	static integer i, j, k, r;
1078	static doublereal u, v, order;
1079	static integer group, m1;
1080	static doublereal x0, x1;
1081	static integer mb, m21, ii, ij, jj, kj;
1082	static doublereal uk, xu;
1083	extern doublereal pythag_(doublereal , doublereal ), epslon_(doublereal
1084	*);
1085	static integer ij1, kj1, its;
1086	static doublereal eps2, eps3, eps4;
1087
1088
1089
1090	/* THIS SUBROUTINE FINDS THOSE EIGENVECTORS OF A REAL SYMMETRIC */
1091	/* BAND MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES, USING INVERSE
1092	*/
1093	/* ITERATION. THE SUBROUTINE MAY ALSO BE USED TO SOLVE SYSTEMS */
1094	/* OF LINEAR EQUATIONS WITH A SYMMETRIC OR NON-SYMMETRIC BAND */
1095	/* COEFFICIENT MATRIX. */
1096
1097	/* ON INPUT */
1098
1099	/* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
1100	/* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
1101	/* DIMENSION STATEMENT. */
1102
1103	/* N IS THE ORDER OF THE MATRIX. */
1104
1105	/* MBW IS THE NUMBER OF COLUMNS OF THE ARRAY A USED TO STORE THE */
1106	/* BAND MATRIX. IF THE MATRIX IS SYMMETRIC, MBW IS ITS (HALF) */
1107	/* BAND WIDTH, DENOTED MB AND DEFINED AS THE NUMBER OF ADJACENT
1108	*/
1109	/* DIAGONALS, INCLUDING THE PRINCIPAL DIAGONAL, REQUIRED TO */
1110	/* SPECIFY THE NON-ZERO PORTION OF THE LOWER TRIANGLE OF THE */
1111	/* MATRIX. IF THE SUBROUTINE IS BEING USED TO SOLVE SYSTEMS */
1112	/* OF LINEAR EQUATIONS AND THE COEFFICIENT MATRIX IS NOT */
1113	/* SYMMETRIC, IT MUST HOWEVER HAVE THE SAME NUMBER OF ADJACENT */
1114	/* DIAGONALS ABOVE THE MAIN DIAGONAL AS BELOW, AND IN THIS */
1115	/* CASE, MBW=2MB-1. /
1116
1117	/* A CONTAINS THE LOWER TRIANGLE OF THE SYMMETRIC BAND INPUT */
1118	/* MATRIX STORED AS AN N BY MB ARRAY. ITS LOWEST SUBDIAGONAL */
1119	/* IS STORED IN THE LAST N+1-MB POSITIONS OF THE FIRST COLUMN, */
1120	/* ITS NEXT SUBDIAGONAL IN THE LAST N+2-MB POSITIONS OF THE */
1121	/* SECOND COLUMN, FURTHER SUBDIAGONALS SIMILARLY, AND FINALLY */
1122	/* ITS PRINCIPAL DIAGONAL IN THE N POSITIONS OF COLUMN MB. */
1123	/* IF THE SUBROUTINE IS BEING USED TO SOLVE SYSTEMS OF LINEAR */
1124	/* EQUATIONS AND THE COEFFICIENT MATRIX IS NOT SYMMETRIC, A IS */
1125	/* N BY 2MB-1 INSTEAD WITH LOWER TRIANGLE AS ABOVE AND WITH /
1126	/* ITS FIRST SUPERDIAGONAL STORED IN THE FIRST N-1 POSITIONS OF
1127	*/
1128	/* COLUMN MB+1, ITS SECOND SUPERDIAGONAL IN THE FIRST N-2 */
1129	/* POSITIONS OF COLUMN MB+2, FURTHER SUPERDIAGONALS SIMILARLY, */
1130	/* AND FINALLY ITS HIGHEST SUPERDIAGONAL IN THE FIRST N+1-MB */
1131	/* POSITIONS OF THE LAST COLUMN. */
1132	/* CONTENTS OF STORAGES NOT PART OF THE MATRIX ARE ARBITRARY. */
1133
1134	/* E21 SPECIFIES THE ORDERING OF THE EIGENVALUES AND CONTAINS */
1135	/* 0.0D0 IF THE EIGENVALUES ARE IN ASCENDING ORDER, OR */
1136	/* 2.0D0 IF THE EIGENVALUES ARE IN DESCENDING ORDER. */
1137	/* IF THE SUBROUTINE IS BEING USED TO SOLVE SYSTEMS OF LINEAR */
1138	/* EQUATIONS, E21 SHOULD BE SET TO 1.0D0 IF THE COEFFICIENT */
1139	/* MATRIX IS SYMMETRIC AND TO -1.0D0 IF NOT. */
1140
1141	/* M IS THE NUMBER OF SPECIFIED EIGENVALUES OR THE NUMBER OF */
1142	/* SYSTEMS OF LINEAR EQUATIONS. */
1143
1144	/* W CONTAINS THE M EIGENVALUES IN ASCENDING OR DESCENDING ORDER.
1145	*/
1146	/* IF THE SUBROUTINE IS BEING USED TO SOLVE SYSTEMS OF LINEAR */
1147	/* EQUATIONS (A-W(R)I)X(R)=B(R), WHERE I IS THE IDENTITY */
1148	/* MATRIX, W(R) SHOULD BE SET ACCORDINGLY, FOR R=1,2,...,M. */
1149
1150	/* Z CONTAINS THE CONSTANT MATRIX COLUMNS (B(R),R=1,2,...,M), IF */
1151	/* THE SUBROUTINE IS USED TO SOLVE SYSTEMS OF LINEAR EQUATIONS.
1152	*/
1153
1154	/* NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER RV */
1155	/* AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT. */
1156
1157	/* ON OUTPUT */
1158
1159	/* A AND W ARE UNALTERED. */
1160
1161	/* Z CONTAINS THE ASSOCIATED SET OF ORTHOGONAL EIGENVECTORS. */
1162	/* ANY VECTOR WHICH FAILS TO CONVERGE IS SET TO ZERO. IF THE */
1163	/* SUBROUTINE IS USED TO SOLVE SYSTEMS OF LINEAR EQUATIONS, */
1164	/* Z CONTAINS THE SOLUTION MATRIX COLUMNS (X(R),R=1,2,...,M). */
1165
1166	/* IERR IS SET TO */
1167	/* ZERO FOR NORMAL RETURN, */
1168	/* -R IF THE EIGENVECTOR CORRESPONDING TO THE R-TH */
1169	/* EIGENVALUE FAILS TO CONVERGE, OR IF THE R-TH */
1170	/* SYSTEM OF LINEAR EQUATIONS IS NEARLY SINGULAR. */
1171
1172	/* RV AND RV6 ARE TEMPORARY STORAGE ARRAYS. NOTE THAT RV IS */
1173	/* OF DIMENSION AT LEAST N(2MB-1). IF THE SUBROUTINE */
1174	/* IS BEING USED TO SOLVE SYSTEMS OF LINEAR EQUATIONS, THE */
1175	/* DETERMINANT (UP TO SIGN) OF A-W(M)I IS AVAILABLE, UPON /
1176	/* RETURN, AS THE PRODUCT OF THE FIRST N ELEMENTS OF RV. */
1177
1178	/* CALLS PYTHAG FOR DSQRT(AA + BB) . */
1179
1180	/* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
1181	/* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
1182	*/
1183
1184	/* THIS VERSION DATED AUGUST 1983. */
1185
1186	/* ------------------------------------------------------------------
1187	*/
1188
1189	/* Parameter adjustments */
1190	--rv6;
1191	a_dim1 = *nm;
1192	a_offset = a_dim1 + 1;
1193	a -= a_offset;
1194	z_dim1 = *nm;
1195	z_offset = z_dim1 + 1;
1196	z -= z_offset;
1197	--w;
1198	--rv;
1199
1200	/* Function Body */
1201	*ierr = 0;
1202	if (*m == 0) {
1203	goto L1001;
1204	}
1205	mb = *mbw;
1206	if (*e21 < 0.) {
1207	mb = (*mbw + 1) / 2;
1208	}
1209	m1 = mb - 1;
1210	m21 = m1 + mb;
1211	order = 1. - abs(*e21);
1212	/* .......... FIND VECTORS BY INVERSE ITERATION .......... */
1213	i_1 = *m;
1214	for (r = 1; r <= i_1; ++r) {
1215	its = 1;
1216	x1 = w[r];
1217	if (r != 1) {
1218	goto L100;
1219	}
1220	/* .......... COMPUTE NORM OF MATRIX .......... */
1221	norm = 0.;
1222
1223	i_2 = mb;
1224	for (j = 1; j <= i_2; ++j) {
1225	jj = mb + 1 - j;
1226	kj = jj + m1;
1227	ij = 1;
1228	v = 0.;
1229
1230	i_3 = *n;
1231	for (i = jj; i <= i_3; ++i) {
1232	v += (d_1 = a[i + j * a_dim1], abs(d_1));
1233	if (*e21 >= 0.) {
1234	goto L40;
1235	}
1236	v += (d_1 = a[ij + kj * a_dim1], abs(d_1));
1237	++ij;
1238	L40:
1239	;
1240	}
1241
1242	norm = max(norm,v);
1243	/* L60: */
1244	}
1245
1246	if (*e21 < 0.) {
1247	norm *= .5;
1248	}
1249	/* .......... EPS2 IS THE CRITERION FOR GROUPING, */
1250	/* EPS3 REPLACES ZERO PIVOTS AND EQUAL */
1251	/* ROOTS ARE MODIFIED BY EPS3, */
1252	/* EPS4 IS TAKEN VERY SMALL TO AVOID OVERFLOW .........
1253	. */
1254	if (norm == 0.) {
1255	norm = 1.;
1256	}
1257	eps2 = norm * .001 * abs(order);
1258	eps3 = epslon_(&norm);
1259	uk = (doublereal) (*n);
1260	uk = sqrt(uk);
1261	eps4 = uk * eps3;
1262	L80:
1263	group = 0;
1264	goto L120;
1265	/* .......... LOOK FOR CLOSE OR COINCIDENT ROOTS .......... */
1266	L100:
1267	if ((d_1 = x1 - x0, abs(d_1)) >= eps2) {
1268	goto L80;
1269	}
1270	++group;
1271	if (order * (x1 - x0) <= 0.) {
1272	x1 = x0 + order * eps3;
1273	}
1274	/* .......... EXPAND MATRIX, SUBTRACT EIGENVALUE, */
1275	/* AND INITIALIZE VECTOR .......... */
1276	L120:
1277	i_2 = *n;
1278	for (i = 1; i <= i_2; ++i) {
1279	/* Computing MIN */
1280	i_3 = 0, i_4 = i - m1;
1281	ij = i + min(i_3,i_4) * *n;
1282	kj = ij + mb * *n;
1283	ij1 = kj + m1 * *n;
1284	if (m1 == 0) {
1285	goto L180;
1286	}
1287
1288	i_3 = m1;
1289	for (j = 1; j <= i_3; ++j) {
1290	if (ij > m1) {
1291	goto L125;
1292	}
1293	if (ij > 0) {
1294	goto L130;
1295	}
1296	rv[ij1] = 0.;
1297	ij1 += *n;
1298	goto L130;
1299	L125:
1300	rv[ij] = a[i + j * a_dim1];
1301	L130:
1302	ij += *n;
1303	ii = i + j;
1304	if (ii > *n) {
1305	goto L150;
1306	}
1307	jj = mb - j;
1308	if (*e21 >= 0.) {
1309	goto L140;
1310	}
1311	ii = i;
1312	jj = mb + j;
1313	L140:
1314	rv[kj] = a[ii + jj * a_dim1];
1315	kj += *n;
1316	L150:
1317	;
1318	}
1319
1320	L180:
1321	rv[ij] = a[i + mb * a_dim1] - x1;
1322	rv6[i] = eps4;
1323	if (order == 0.) {
1324	rv6[i] = z[i + r * z_dim1];
1325	}
1326	/* L200: */
1327	}
1328
1329	if (m1 == 0) {
1330	goto L600;
1331	}
1332	/* .......... ELIMINATION WITH INTERCHANGES .......... */
1333	i_2 = *n;
1334	for (i = 1; i <= i_2; ++i) {
1335	ii = i + 1;
1336	/* Computing MIN */
1337	i_3 = i + m1 - 1;
1338	maxk = min(i_3,*n);
1339	/* Computing MIN */
1340	i_3 = *n - i, i_4 = m21 - 2;
1341	maxj = min(i_3,i_4) * *n;
1342
1343	i_3 = maxk;
1344	for (k = i; k <= i_3; ++k) {
1345	kj1 = k;
1346	j = kj1 + *n;
1347	jj = j + maxj;
1348
1349	i_4 = jj;
1350	i_5 = *n;
1351	for (kj = j; i_5 < 0 ? kj >= i_4 : kj <= i_4; kj += i_5) {
1352	rv[kj1] = rv[kj];
1353	kj1 = kj;
1354	/* L340: */
1355	}
1356
1357	rv[kj1] = 0.;
1358	/* L360: */
1359	}
1360
1361	if (i == *n) {
1362	goto L580;
1363	}
1364	u = 0.;
1365	/* Computing MIN */
1366	i_3 = i + m1;
1367	maxk = min(i_3,*n);
1368	/* Computing MIN */
1369	i_3 = *n - ii, i_5 = m21 - 2;
1370	maxj = min(i_3,i_5) * *n;
1371
1372	i_3 = maxk;
1373	for (j = i; j <= i_3; ++j) {
1374	if ((d_1 = rv[j], abs(d_1)) < abs(u)) {
1375	goto L450;
1376	}
1377	u = rv[j];
1378	k = j;
1379	L450:
1380	;
1381	}
1382
1383	j = i + *n;
1384	jj = j + maxj;
1385	if (k == i) {
1386	goto L520;
1387	}
1388	kj = k;
1389
1390	i_3 = jj;
1391	i_5 = *n;
1392	for (ij = i; i_5 < 0 ? ij >= i_3 : ij <= i_3; ij += i_5) {
1393	v = rv[ij];
1394	rv[ij] = rv[kj];
1395	rv[kj] = v;
1396	kj += *n;
1397	/* L500: */
1398	}
1399
1400	if (order != 0.) {
1401	goto L520;
1402	}
1403	v = rv6[i];
1404	rv6[i] = rv6[k];
1405	rv6[k] = v;
1406	L520:
1407	if (u == 0.) {
1408	goto L580;
1409	}
1410
1411	i_5 = maxk;
1412	for (k = ii; k <= i_5; ++k) {
1413	v = rv[k] / u;
1414	kj = k;
1415
1416	i_3 = jj;
1417	i_4 = *n;
1418	for (ij = j; i_4 < 0 ? ij >= i_3 : ij <= i_3; ij += i_4) {
1419	kj += *n;
1420	rv[kj] -= v * rv[ij];
1421	/* L540: */
1422	}
1423
1424	if (order == 0.) {
1425	rv6[k] -= v * rv6[i];
1426	}
1427	/* L560: */
1428	}
1429
1430	L580:
1431	;
1432	}
1433	/* .......... BACK SUBSTITUTION */
1434	/* FOR I=N STEP -1 UNTIL 1 DO -- .......... */
1435	L600:
1436	i_2 = *n;
1437	for (ii = 1; ii <= i_2; ++ii) {
1438	i = *n + 1 - ii;
1439	maxj = min(ii,m21);
1440	if (maxj == 1) {
1441	goto L620;
1442	}
1443	ij1 = i;
1444	j = ij1 + *n;
1445	jj = j + (maxj - 2) * *n;
1446
1447	i_5 = jj;
1448	i_4 = *n;
1449	for (ij = j; i_4 < 0 ? ij >= i_5 : ij <= i_5; ij += i_4) {
1450	++ij1;
1451	rv6[i] -= rv[ij] * rv6[ij1];
1452	/* L610: */
1453	}
1454
1455	L620:
1456	v = rv[i];
1457	if (abs(v) >= eps3) {
1458	goto L625;
1459	}
1460	/* .......... SET ERROR -- NEARLY SINGULAR LINEAR SYSTEM .....
1461	..... */
1462	if (order == 0.) {
1463	*ierr = -r;
1464	}
1465	v = d_sign(&eps3, &v);
1466	L625:
1467	rv6[i] /= v;
1468	/* L630: */
1469	}
1470
1471	xu = 1.;
1472	if (order == 0.) {
1473	goto L870;
1474	}
1475	/* .......... ORTHOGONALIZE WITH RESPECT TO PREVIOUS */
1476	/* MEMBERS OF GROUP .......... */
1477	if (group == 0) {
1478	goto L700;
1479	}
1480
1481	i_2 = group;
1482	for (jj = 1; jj <= i_2; ++jj) {
1483	j = r - group - 1 + jj;
1484	xu = 0.;
1485
1486	i_4 = *n;
1487	for (i = 1; i <= i_4; ++i) {
1488	/* L640: */
1489	xu += rv6[i] * z[i + j * z_dim1];
1490	}
1491
1492	i_4 = *n;
1493	for (i = 1; i <= i_4; ++i) {
1494	/* L660: */
1495	rv6[i] -= xu * z[i + j * z_dim1];
1496	}
1497
1498	/* L680: */
1499	}
1500
1501	L700:
1502	norm = 0.;
1503
1504	i_2 = *n;
1505	for (i = 1; i <= i_2; ++i) {
1506	/* L720: */
1507	norm += (d_1 = rv6[i], abs(d_1));
1508	}
1509
1510	if (norm >= .1) {
1511	goto L840;
1512	}
1513	/* .......... IN-LINE PROCEDURE FOR CHOOSING */
1514	/* A NEW STARTING VECTOR .......... */
1515	if (its >= *n) {
1516	goto L830;
1517	}
1518	++its;
1519	xu = eps4 / (uk + 1.);
1520	rv6[1] = eps4;
1521
1522	i_2 = *n;
1523	for (i = 2; i <= i_2; ++i) {
1524	/* L760: */
1525	rv6[i] = xu;
1526	}
1527
1528	rv6[its] -= eps4 * uk;
1529	goto L600;
1530	/* .......... SET ERROR -- NON-CONVERGED EIGENVECTOR .......... */
1531	L830:
1532	*ierr = -r;
1533	xu = 0.;
1534	goto L870;
1535	/* .......... NORMALIZE SO THAT SUM OF SQUARES IS */
1536	/* 1 AND EXPAND TO FULL ORDER .......... */
1537	L840:
1538	u = 0.;
1539
1540	i_2 = *n;
1541	for (i = 1; i <= i_2; ++i) {
1542	/* L860: */
1543	u = pythag_(&u, &rv6[i]);
1544	}
1545
1546	xu = 1. / u;
1547
1548	L870:
1549	i_2 = *n;
1550	for (i = 1; i <= i_2; ++i) {
1551	/* L900: */
1552	z[i + r * z_dim1] = rv6[i] * xu;
1553	}
1554
1555	x0 = x1;
1556	/* L920: */
1557	}
1558
1559	L1001:
1560	return 0;
1561	} /* bandv_ */
1562
1563	/* Subroutine / int bisect_(integer n, doublereal eps1, doublereal d,
1564	doublereal e, doublereal e2, doublereal lb, doublereal ub,
1565	integer mm, integer m, doublereal w, integer ind, integer *ierr,
1566	doublereal rv4, doublereal rv5)
1567	{
1568	/* System generated locals */
1569	integer i_1, i_2;
1570	doublereal d_1, d_2, d_3;
1571
1572	/* Local variables */
1573	static integer i, j, k, l, p, q, r, s;
1574	static doublereal u, v;
1575	static integer m1, m2;
1576	static doublereal t1, t2, x0, x1;
1577	static integer ii;
1578	static doublereal xu;
1579	extern doublereal epslon_(doublereal *);
1580	static integer isturm, tag;
1581	static doublereal tst1, tst2;
1582
1583
1584
1585	/* THIS SUBROUTINE IS A TRANSLATION OF THE BISECTION TECHNIQUE */
1586	/* IN THE ALGOL PROCEDURE TRISTURM BY PETERS AND WILKINSON. */
1587	/* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971). */
1588
1589	/* THIS SUBROUTINE FINDS THOSE EIGENVALUES OF A TRIDIAGONAL */
1590	/* SYMMETRIC MATRIX WHICH LIE IN A SPECIFIED INTERVAL, */
1591	/* USING BISECTION. */
1592
1593	/* ON INPUT */
1594
1595	/* N IS THE ORDER OF THE MATRIX. */
1596
1597	/* EPS1 IS AN ABSOLUTE ERROR TOLERANCE FOR THE COMPUTED */
1598	/* EIGENVALUES. IF THE INPUT EPS1 IS NON-POSITIVE, */
1599	/* IT IS RESET FOR EACH SUBMATRIX TO A DEFAULT VALUE, */
1600	/* NAMELY, MINUS THE PRODUCT OF THE RELATIVE MACHINE */
1601	/* PRECISION AND THE 1-NORM OF THE SUBMATRIX. */
1602
1603	/* D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. */
1604
1605	/* E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX */
1606	/* IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY. */
1607
1608	/* E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. */
1609	/* E2(1) IS ARBITRARY. */
1610
1611	/* LB AND UB DEFINE THE INTERVAL TO BE SEARCHED FOR EIGENVALUES. */
1612	/* IF LB IS NOT LESS THAN UB, NO EIGENVALUES WILL BE FOUND. */
1613
1614	/* MM SHOULD BE SET TO AN UPPER BOUND FOR THE NUMBER OF */
1615	/* EIGENVALUES IN THE INTERVAL. WARNING. IF MORE THAN */
1616	/* MM EIGENVALUES ARE DETERMINED TO LIE IN THE INTERVAL, */
1617	/* AN ERROR RETURN IS MADE WITH NO EIGENVALUES FOUND. */
1618
1619	/* ON OUTPUT */
1620
1621	/* EPS1 IS UNALTERED UNLESS IT HAS BEEN RESET TO ITS */
1622	/* (LAST) DEFAULT VALUE. */
1623
1624	/* D AND E ARE UNALTERED. */
1625
1626	/* ELEMENTS OF E2, CORRESPONDING TO ELEMENTS OF E REGARDED */
1627	/* AS NEGLIGIBLE, HAVE BEEN REPLACED BY ZERO CAUSING THE */
1628	/* MATRIX TO SPLIT INTO A DIRECT SUM OF SUBMATRICES. */
1629	/* E2(1) IS ALSO SET TO ZERO. */
1630
1631	/* M IS THE NUMBER OF EIGENVALUES DETERMINED TO LIE IN (LB,UB). */
1632
1633	/* W CONTAINS THE M EIGENVALUES IN ASCENDING ORDER. */
1634
1635	/* IND CONTAINS IN ITS FIRST M POSITIONS THE SUBMATRIX INDICES */
1636	/* ASSOCIATED WITH THE CORRESPONDING EIGENVALUES IN W -- */
1637	/* 1 FOR EIGENVALUES BELONGING TO THE FIRST SUBMATRIX FROM */
1638	/* THE TOP, 2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC..
1639	*/
1640
1641	/* IERR IS SET TO */
1642	/* ZERO FOR NORMAL RETURN, */
1643	/* 3N+1 IF M EXCEEDS MM. /
1644
1645	/* RV4 AND RV5 ARE TEMPORARY STORAGE ARRAYS. */
1646
1647	/* THE ALGOL PROCEDURE STURMCNT CONTAINED IN TRISTURM */
1648	/* APPEARS IN BISECT IN-LINE. */
1649
1650	/* NOTE THAT SUBROUTINE TQL1 OR IMTQL1 IS GENERALLY FASTER THAN */
1651	/* BISECT, IF MORE THAN N/4 EIGENVALUES ARE TO BE FOUND. */
1652
1653	/* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
1654	/* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
1655	*/
1656
1657	/* THIS VERSION DATED AUGUST 1983. */
1658
1659	/* ------------------------------------------------------------------
1660	*/
1661
1662	/* Parameter adjustments */
1663	--rv5;
1664	--rv4;
1665	--e2;
1666	--e;
1667	--d;
1668	--ind;
1669	--w;
1670
1671	/* Function Body */
1672	*ierr = 0;
1673	tag = 0;
1674	t1 = *lb;
1675	t2 = *ub;
1676	/* .......... LOOK FOR SMALL SUB-DIAGONAL ENTRIES .......... */
1677	i_1 = *n;
1678	for (i = 1; i <= i_1; ++i) {
1679	if (i == 1) {
1680	goto L20;
1681	}
1682	tst1 = (d_1 = d[i], abs(d_1)) + (d_2 = d[i - 1], abs(d_2));
1683	tst2 = tst1 + (d_1 = e[i], abs(d_1));
1684	if (tst2 > tst1) {
1685	goto L40;
1686	}
1687	L20:
1688	e2[i] = 0.;
1689	L40:
1690	;
1691	}
1692	/* .......... DETERMINE THE NUMBER OF EIGENVALUES */
1693	/* IN THE INTERVAL .......... */
1694	p = 1;
1695	q = *n;
1696	x1 = *ub;
1697	isturm = 1;
1698	goto L320;
1699	L60:
1700	*m = s;
1701	x1 = *lb;
1702	isturm = 2;
1703	goto L320;
1704	L80:
1705	*m -= s;
1706	if (m > mm) {
1707	goto L980;
1708	}
1709	q = 0;
1710	r = 0;
1711	/* .......... ESTABLISH AND PROCESS NEXT SUBMATRIX, REFINING */
1712	/* INTERVAL BY THE GERSCHGORIN BOUNDS .......... */
1713	L100:
1714	if (r == *m) {
1715	goto L1001;
1716	}
1717	++tag;
1718	p = q + 1;
1719	xu = d[p];
1720	x0 = d[p];
1721	u = 0.;
1722
1723	i_1 = *n;
1724	for (q = p; q <= i_1; ++q) {
1725	x1 = u;
1726	u = 0.;
1727	v = 0.;
1728	if (q == *n) {
1729	goto L110;
1730	}
1731	u = (d_1 = e[q + 1], abs(d_1));
1732	v = e2[q + 1];
1733	L110:
1734	/* Computing MIN */
1735	d_1 = d[q] - (x1 + u);
1736	xu = min(d_1,xu);
1737	/* Computing MAX */
1738	d_1 = d[q] + (x1 + u);
1739	x0 = max(d_1,x0);
1740	if (v == 0.) {
1741	goto L140;
1742	}
1743	/* L120: */
1744	}
1745
1746	L140:
1747	/* Computing MAX */
1748	d_2 = abs(xu), d_3 = abs(x0);
1749	d_1 = max(d_2,d_3);
1750	x1 = epslon_(&d_1);
1751	if (*eps1 <= 0.) {
1752	*eps1 = -x1;
1753	}
1754	if (p != q) {
1755	goto L180;
1756	}
1757	/* .......... CHECK FOR ISOLATED ROOT WITHIN INTERVAL .......... */
1758	if (t1 > d[p] \|\| d[p] >= t2) {
1759	goto L940;
1760	}
1761	m1 = p;
1762	m2 = p;
1763	rv5[p] = d[p];
1764	goto L900;
1765	L180:
1766	x1 *= q - p + 1;
1767	/* Computing MAX */
1768	d_1 = t1, d_2 = xu - x1;
1769	*lb = max(d_1,d_2);
1770	/* Computing MIN */
1771	d_1 = t2, d_2 = x0 + x1;
1772	*ub = min(d_1,d_2);
1773	x1 = *lb;
1774	isturm = 3;
1775	goto L320;
1776	L200:
1777	m1 = s + 1;
1778	x1 = *ub;
1779	isturm = 4;
1780	goto L320;
1781	L220:
1782	m2 = s;
1783	if (m1 > m2) {
1784	goto L940;
1785	}
1786	/* .......... FIND ROOTS BY BISECTION .......... */
1787	x0 = *ub;
1788	isturm = 5;
1789
1790	i_1 = m2;
1791	for (i = m1; i <= i_1; ++i) {
1792	rv5[i] = *ub;
1793	rv4[i] = *lb;
1794	/* L240: */
1795	}
1796	/* .......... LOOP FOR K-TH EIGENVALUE */
1797	/* FOR K=M2 STEP -1 UNTIL M1 DO -- */
1798	/* (-DO- NOT USED TO LEGALIZE -COMPUTED GO TO-) ..........
1799	*/
1800	k = m2;
1801	L250:
1802	xu = *lb;
1803	/* .......... FOR I=K STEP -1 UNTIL M1 DO -- .......... */
1804	i_1 = k;
1805	for (ii = m1; ii <= i_1; ++ii) {
1806	i = m1 + k - ii;
1807	if (xu >= rv4[i]) {
1808	goto L260;
1809	}
1810	xu = rv4[i];
1811	goto L280;
1812	L260:
1813	;
1814	}
1815
1816	L280:
1817	if (x0 > rv5[k]) {
1818	x0 = rv5[k];
1819	}
1820	/* .......... NEXT BISECTION STEP .......... */
1821	L300:
1822	x1 = (xu + x0) * .5;
1823	if (x0 - xu <= abs(*eps1)) {
1824	goto L420;
1825	}
1826	tst1 = (abs(xu) + abs(x0)) * 2.;
1827	tst2 = tst1 + (x0 - xu);
1828	if (tst2 == tst1) {
1829	goto L420;
1830	}
1831	/* .......... IN-LINE PROCEDURE FOR STURM SEQUENCE .......... */
1832	L320:
1833	s = p - 1;
1834	u = 1.;
1835
1836	i_1 = q;
1837	for (i = p; i <= i_1; ++i) {
1838	if (u != 0.) {
1839	goto L325;
1840	}
1841	v = (d_1 = e[i], abs(d_1)) / epslon_(&c_b141);
1842	if (e2[i] == 0.) {
1843	v = 0.;
1844	}
1845	goto L330;
1846	L325:
1847	v = e2[i] / u;
1848	L330:
1849	u = d[i] - x1 - v;
1850	if (u < 0.) {
1851	++s;
1852	}
1853	/* L340: */
1854	}
1855
1856	switch (isturm) {
1857	case 1: goto L60;
1858	case 2: goto L80;
1859	case 3: goto L200;
1860	case 4: goto L220;
1861	case 5: goto L360;
1862	}
1863	/* .......... REFINE INTERVALS .......... */
1864	L360:
1865	if (s >= k) {
1866	goto L400;
1867	}
1868	xu = x1;
1869	if (s >= m1) {
1870	goto L380;
1871	}
1872	rv4[m1] = x1;
1873	goto L300;
1874	L380:
1875	rv4[s + 1] = x1;
1876	if (rv5[s] > x1) {
1877	rv5[s] = x1;
1878	}
1879	goto L300;
1880	L400:
1881	x0 = x1;
1882	goto L300;
1883	/* .......... K-TH EIGENVALUE FOUND .......... */
1884	L420:
1885	rv5[k] = x1;
1886	--k;
1887	if (k >= m1) {
1888	goto L250;
1889	}
1890	/* .......... ORDER EIGENVALUES TAGGED WITH THEIR */
1891	/* SUBMATRIX ASSOCIATIONS .......... */
1892	L900:
1893	s = r;
1894	r = r + m2 - m1 + 1;
1895	j = 1;
1896	k = m1;
1897
1898	i_1 = r;
1899	for (l = 1; l <= i_1; ++l) {
1900	if (j > s) {
1901	goto L910;
1902	}
1903	if (k > m2) {
1904	goto L940;
1905	}
1906	if (rv5[k] >= w[l]) {
1907	goto L915;
1908	}
1909
1910	i_2 = s;
1911	for (ii = j; ii <= i_2; ++ii) {
1912	i = l + s - ii;
1913	w[i + 1] = w[i];
1914	ind[i + 1] = ind[i];
1915	/* L905: */
1916	}
1917
1918	L910:
1919	w[l] = rv5[k];
1920	ind[l] = tag;
1921	++k;
1922	goto L920;
1923	L915:
1924	++j;
1925	L920:
1926	;
1927	}
1928
1929	L940:
1930	if (q < *n) {
1931	goto L100;
1932	}
1933	goto L1001;
1934	/* .......... SET ERROR -- UNDERESTIMATE OF NUMBER OF */
1935	/* EIGENVALUES IN INTERVAL .......... */
1936	L980:
1937	ierr = n * 3 + 1;
1938	L1001:
1939	*lb = t1;
1940	*ub = t2;
1941	return 0;
1942	} /* bisect_ */
1943
1944	/* Subroutine / int bqr_(integer nm, integer n, integer mb, doublereal *a,
1945	doublereal t, doublereal r, integer ierr, integer /nv/, doublereal
1946	*rv)
1947	{
1948	/* System generated locals */
1949	integer a_dim1, a_offset, i_1, i_2, i_3;
1950	doublereal d_1;
1951
1952	/* Builtin functions */
1953	double d_sign(doublereal , doublereal ), sqrt(doublereal);
1954
1955	/* Local variables */
1956	static doublereal f, g;
1957	static integer i, j, k, l, m;
1958	static doublereal q, s, scale;
1959	static integer imult, m1, m2, m3, m4, m21, m31, ii, ik, jk, kj, jm, kk,
1960	km, ll, mk, mn, ni, mz;
1961	extern doublereal pythag_(doublereal , doublereal );
1962	static integer kj1, its;
1963	static doublereal tst1, tst2;
1964
1965
1966
1967	/* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BQR, */
1968	/* NUM. MATH. 16, 85-92(1970) BY MARTIN, REINSCH, AND WILKINSON. */
1969	/* HANDBOOK FOR AUTO. COMP., VOL II-LINEAR ALGEBRA, 266-272(1971). */
1970
1971	/* THIS SUBROUTINE FINDS THE EIGENVALUE OF SMALLEST (USUALLY) */
1972	/* MAGNITUDE OF A REAL SYMMETRIC BAND MATRIX USING THE */
1973	/* QR ALGORITHM WITH SHIFTS OF ORIGIN. CONSECUTIVE CALLS */
1974	/* CAN BE MADE TO FIND FURTHER EIGENVALUES. */
1975
1976	/* ON INPUT */
1977
1978	/* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
1979	/* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
1980	/* DIMENSION STATEMENT. */
1981
1982	/* N IS THE ORDER OF THE MATRIX. */
1983
1984	/* MB IS THE (HALF) BAND WIDTH OF THE MATRIX, DEFINED AS THE */
1985	/* NUMBER OF ADJACENT DIAGONALS, INCLUDING THE PRINCIPAL */
1986	/* DIAGONAL, REQUIRED TO SPECIFY THE NON-ZERO PORTION OF THE */
1987	/* LOWER TRIANGLE OF THE MATRIX. */
1988
1989	/* A CONTAINS THE LOWER TRIANGLE OF THE SYMMETRIC BAND INPUT */
1990	/* MATRIX STORED AS AN N BY MB ARRAY. ITS LOWEST SUBDIAGONAL */
1991	/* IS STORED IN THE LAST N+1-MB POSITIONS OF THE FIRST COLUMN, */
1992	/* ITS NEXT SUBDIAGONAL IN THE LAST N+2-MB POSITIONS OF THE */
1993	/* SECOND COLUMN, FURTHER SUBDIAGONALS SIMILARLY, AND FINALLY */
1994	/* ITS PRINCIPAL DIAGONAL IN THE N POSITIONS OF THE LAST COLUMN.
1995	*/
1996	/* CONTENTS OF STORAGES NOT PART OF THE MATRIX ARE ARBITRARY. */
1997	/* ON A SUBSEQUENT CALL, ITS OUTPUT CONTENTS FROM THE PREVIOUS */
1998	/* CALL SHOULD BE PASSED. */
1999
2000	/* T SPECIFIES THE SHIFT (OF EIGENVALUES) APPLIED TO THE DIAGONAL
2001	*/
2002	/* OF A IN FORMING THE INPUT MATRIX. WHAT IS ACTUALLY DETERMINED
2003	*/
2004	/* IS THE EIGENVALUE OF A+TI (I IS THE IDENTITY MATRIX) NEAREST
2005	*/
2006	/* TO T. ON A SUBSEQUENT CALL, THE OUTPUT VALUE OF T FROM THE */
2007	/* PREVIOUS CALL SHOULD BE PASSED IF THE NEXT NEAREST EIGENVALUE
2008	*/
2009	/* IS SOUGHT. */
2010
2011	/* R SHOULD BE SPECIFIED AS ZERO ON THE FIRST CALL, AND AS ITS */
2012	/* OUTPUT VALUE FROM THE PREVIOUS CALL ON A SUBSEQUENT CALL. */
2013	/* IT IS USED TO DETERMINE WHEN THE LAST ROW AND COLUMN OF */
2014	/* THE TRANSFORMED BAND MATRIX CAN BE REGARDED AS NEGLIGIBLE. */
2015
2016	/* NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER RV */
2017	/* AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT. */
2018
2019	/* ON OUTPUT */
2020
2021	/* A CONTAINS THE TRANSFORMED BAND MATRIX. THE MATRIX A+TI */
2022	/* DERIVED FROM THE OUTPUT PARAMETERS IS SIMILAR TO THE */
2023	/* INPUT A+TI TO WITHIN ROUNDING ERRORS. ITS LAST ROW AND */
2024	/* COLUMN ARE NULL (IF IERR IS ZERO). */
2025
2026	/* T CONTAINS THE COMPUTED EIGENVALUE OF A+TI (IF IERR IS ZERO). */
2027
2028	/* R CONTAINS THE MAXIMUM OF ITS INPUT VALUE AND THE NORM OF THE */
2029	/* LAST COLUMN OF THE INPUT MATRIX A. */
2030
2031	/* IERR IS SET TO */
2032	/* ZERO FOR NORMAL RETURN, */
2033	/* N IF THE EIGENVALUE HAS NOT BEEN */
2034	/* DETERMINED AFTER 30 ITERATIONS. */
2035
2036	/* RV IS A TEMPORARY STORAGE ARRAY OF DIMENSION AT LEAST */
2037	/* (2MB2+4MB-3). THE FIRST (3MB-2) LOCATIONS CORRESPOND /
2038	/* TO THE ALGOL ARRAY B, THE NEXT (2*MB-1) LOCATIONS CORRESPOND
2039	*/
2040	/* TO THE ALGOL ARRAY H, AND THE FINAL (2MB2-MB) LOCATIONS /
2041	/* CORRESPOND TO THE MB BY (2MB-1) ALGOL ARRAY U. /
2042
2043	/* NOTE. FOR A SUBSEQUENT CALL, N SHOULD BE REPLACED BY N-1, BUT */
2044	/* MB SHOULD NOT BE ALTERED EVEN WHEN IT EXCEEDS THE CURRENT N. */
2045
2046	/* CALLS PYTHAG FOR DSQRT(AA + BB) . */
2047
2048	/* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
2049	/* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
2050	*/
2051
2052	/* THIS VERSION DATED AUGUST 1983. */
2053
2054	/* ------------------------------------------------------------------
2055	*/
2056
2057	/* Parameter adjustments */
2058	a_dim1 = *nm;
2059	a_offset = a_dim1 + 1;
2060	a -= a_offset;
2061	--rv;
2062
2063	/* Function Body */
2064	*ierr = 0;
2065	m1 = min(mb,n);
2066	m = m1 - 1;
2067	m2 = m + m;
2068	m21 = m2 + 1;
2069	m3 = m21 + m;
2070	m31 = m3 + 1;
2071	m4 = m31 + m2;
2072	mn = m + *n;
2073	mz = *mb - m1;
2074	its = 0;
2075	/* .......... TEST FOR CONVERGENCE .......... */
2076	L40:
2077	g = a[n + mb * a_dim1];
2078	if (m == 0) {
2079	goto L360;
2080	}
2081	f = 0.;
2082
2083	i_1 = m;
2084	for (k = 1; k <= i_1; ++k) {
2085	mk = k + mz;
2086	f += (d_1 = a[n + mk a_dim1], abs(d_1));
2087	/* L50: */
2088	}
2089
2090	if (its == 0 && f > *r) {
2091	*r = f;
2092	}
2093	tst1 = *r;
2094	tst2 = tst1 + f;
2095	if (tst2 <= tst1) {
2096	goto L360;
2097	}
2098	if (its == 30) {
2099	goto L1000;
2100	}
2101	++its;
2102	/* .......... FORM SHIFT FROM BOTTOM 2 BY 2 MINOR .......... */
2103	if (f > r .25 && its < 5) {
2104	goto L90;
2105	}
2106	f = a[n + (mb - 1) * a_dim1];
2107	if (f == 0.) {
2108	goto L70;
2109	}
2110	q = (a[n - 1 + mb * a_dim1] - g) / (f * 2.);
2111	s = pythag_(&q, &c_b141);
2112	g -= f / (q + d_sign(&s, &q));
2113	L70:
2114	*t += g;
2115
2116	i_1 = *n;
2117	for (i = 1; i <= i_1; ++i) {
2118	/* L80: */
2119	a[i + mb a_dim1] -= g;
2120	}
2121
2122	L90:
2123	i_1 = m4;
2124	for (k = m31; k <= i_1; ++k) {
2125	/* L100: */
2126	rv[k] = 0.;
2127	}
2128
2129	i_1 = mn;
2130	for (ii = 1; ii <= i_1; ++ii) {
2131	i = ii - m;
2132	ni = *n - ii;
2133	if (ni < 0) {
2134	goto L230;
2135	}
2136	/* .......... FORM COLUMN OF SHIFTED MATRIX A-GI .......... /
2137	/* Computing MAX */
2138	i_2 = 1, i_3 = 2 - i;
2139	l = max(i_2,i_3);
2140
2141	i_2 = m3;
2142	for (k = 1; k <= i_2; ++k) {
2143	/* L110: */
2144	rv[k] = 0.;
2145	}
2146
2147	i_2 = m1;
2148	for (k = l; k <= i_2; ++k) {
2149	km = k + m;
2150	mk = k + mz;
2151	rv[km] = a[ii + mk * a_dim1];
2152	/* L120: */
2153	}
2154
2155	ll = min(m,ni);
2156	if (ll == 0) {
2157	goto L135;
2158	}
2159
2160	i_2 = ll;
2161	for (k = 1; k <= i_2; ++k) {
2162	km = k + m21;
2163	ik = ii + k;
2164	mk = *mb - k;
2165	rv[km] = a[ik + mk * a_dim1];
2166	/* L130: */
2167	}
2168	/* .......... PRE-MULTIPLY WITH HOUSEHOLDER REFLECTIONS ..........
2169	*/
2170	L135:
2171	ll = m2;
2172	imult = 0;
2173	/* .......... MULTIPLICATION PROCEDURE .......... */
2174	L140:
2175	kj = m4 - m1;
2176
2177	i_2 = ll;
2178	for (j = 1; j <= i_2; ++j) {
2179	kj += m1;
2180	jm = j + m3;
2181	if (rv[jm] == 0.) {
2182	goto L170;
2183	}
2184	f = 0.;
2185
2186	i_3 = m1;
2187	for (k = 1; k <= i_3; ++k) {
2188	++kj;
2189	jk = j + k - 1;
2190	f += rv[kj] * rv[jk];
2191	/* L150: */
2192	}
2193
2194	f /= rv[jm];
2195	kj -= m1;
2196
2197	i_3 = m1;
2198	for (k = 1; k <= i_3; ++k) {
2199	++kj;
2200	jk = j + k - 1;
2201	rv[jk] -= rv[kj] * f;
2202	/* L160: */
2203	}
2204
2205	kj -= m1;
2206	L170:
2207	;
2208	}
2209
2210	if (imult != 0) {
2211	goto L280;
2212	}
2213	/* .......... HOUSEHOLDER REFLECTION .......... */
2214	f = rv[m21];
2215	s = 0.;
2216	rv[m4] = 0.;
2217	scale = 0.;
2218
2219	i_2 = m3;
2220	for (k = m21; k <= i_2; ++k) {
2221	/* L180: */
2222	scale += (d_1 = rv[k], abs(d_1));
2223	}
2224
2225	if (scale == 0.) {
2226	goto L210;
2227	}
2228
2229	i_2 = m3;
2230	for (k = m21; k <= i_2; ++k) {
2231	/* L190: */
2232	/* Computing 2nd power */
2233	d_1 = rv[k] / scale;
2234	s += d_1 * d_1;
2235	}
2236
2237	s = scale * scale * s;
2238	d_1 = sqrt(s);
2239	g = -d_sign(&d_1, &f);
2240	rv[m21] = g;
2241	rv[m4] = s - f * g;
2242	kj = m4 + m2 * m1 + 1;
2243	rv[kj] = f - g;
2244
2245	i_2 = m1;
2246	for (k = 2; k <= i_2; ++k) {
2247	++kj;
2248	km = k + m2;
2249	rv[kj] = rv[km];
2250	/* L200: */
2251	}
2252	/* .......... SAVE COLUMN OF TRIANGULAR FACTOR R .......... */
2253	L210:
2254	i_2 = m1;
2255	for (k = l; k <= i_2; ++k) {
2256	km = k + m;
2257	mk = k + mz;
2258	a[ii + mk * a_dim1] = rv[km];
2259	/* L220: */
2260	}
2261
2262	L230:
2263	/* Computing MAX */
2264	i_2 = 1, i_3 = m1 + 1 - i;
2265	l = max(i_2,i_3);
2266	if (i <= 0) {
2267	goto L300;
2268	}
2269	/* .......... PERFORM ADDITIONAL STEPS .......... */
2270	i_2 = m21;
2271	for (k = 1; k <= i_2; ++k) {
2272	/* L240: */
2273	rv[k] = 0.;
2274	}
2275
2276	/* Computing MIN */
2277	i_2 = m1, i_3 = ni + m1;
2278	ll = min(i_2,i_3);
2279	/* .......... GET ROW OF TRIANGULAR FACTOR R .......... */
2280	i_2 = ll;
2281	for (kk = 1; kk <= i_2; ++kk) {
2282	k = kk - 1;
2283	km = k + m1;
2284	ik = i + k;
2285	mk = *mb - k;
2286	rv[km] = a[ik + mk * a_dim1];
2287	/* L250: */
2288	}
2289	/* .......... POST-MULTIPLY WITH HOUSEHOLDER REFLECTIONS .........
2290	. */
2291	ll = m1;
2292	imult = 1;
2293	goto L140;
2294	/* .......... STORE COLUMN OF NEW A MATRIX .......... */
2295	L280:
2296	i_2 = m1;
2297	for (k = l; k <= i_2; ++k) {
2298	mk = k + mz;
2299	a[i + mk * a_dim1] = rv[k];
2300	/* L290: */
2301	}
2302	/* .......... UPDATE HOUSEHOLDER REFLECTIONS .......... */
2303	L300:
2304	if (l > 1) {
2305	--l;
2306	}
2307	kj1 = m4 + l * m1;
2308
2309	i_2 = m2;
2310	for (j = l; j <= i_2; ++j) {
2311	jm = j + m3;
2312	rv[jm] = rv[jm + 1];
2313
2314	i_3 = m1;
2315	for (k = 1; k <= i_3; ++k) {
2316	++kj1;
2317	kj = kj1 - m1;
2318	rv[kj] = rv[kj1];
2319	/* L320: */
2320	}
2321	}
2322
2323	/* L350: */
2324	}
2325
2326	goto L40;
2327	/* .......... CONVERGENCE .......... */
2328	L360:
2329	*t += g;
2330
2331	i_1 = *n;
2332	for (i = 1; i <= i_1; ++i) {
2333	/* L380: */
2334	a[i + mb a_dim1] -= g;
2335	}
2336
2337	i_1 = m1;
2338	for (k = 1; k <= i_1; ++k) {
2339	mk = k + mz;
2340	a[n + mk a_dim1] = 0.;
2341	/* L400: */
2342	}
2343
2344	goto L1001;
2345	/* .......... SET ERROR -- NO CONVERGENCE TO */
2346	/* EIGENVALUE AFTER 30 ITERATIONS .......... */
2347	L1000:
2348	ierr = n;
2349	L1001:
2350	return 0;
2351	} /* bqr_ */
2352
2353	/* Subroutine / int cbabk2_(integer nm, integer n, integer low, integer *
2354	igh, doublereal scale, integer m, doublereal zr, doublereal zi)
2355	{
2356	/* System generated locals */
2357	integer zr_dim1, zr_offset, zi_dim1, zi_offset, i_1, i_2;
2358
2359	/* Local variables */
2360	static integer i, j, k;
2361	static doublereal s;
2362	static integer ii;
2363
2364
2365
2366	/* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE */
2367	/* CBABK2, WHICH IS A COMPLEX VERSION OF BALBAK, */
2368	/* NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH. */
2369	/* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971). */
2370
2371	/* THIS SUBROUTINE FORMS THE EIGENVECTORS OF A COMPLEX GENERAL */
2372	/* MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING */
2373	/* BALANCED MATRIX DETERMINED BY CBAL. */
2374
2375	/* ON INPUT */
2376
2377	/* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
2378	/* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
2379	/* DIMENSION STATEMENT. */
2380
2381	/* N IS THE ORDER OF THE MATRIX. */
2382
2383	/* LOW AND IGH ARE INTEGERS DETERMINED BY CBAL. */
2384
2385	/* SCALE CONTAINS INFORMATION DETERMINING THE PERMUTATIONS */
2386	/* AND SCALING FACTORS USED BY CBAL. */
2387
2388	/* M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED. */
2389
2390	/* ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, */
2391	/* RESPECTIVELY, OF THE EIGENVECTORS TO BE */
2392	/* BACK TRANSFORMED IN THEIR FIRST M COLUMNS. */
2393
2394	/* ON OUTPUT */
2395
2396	/* ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, */
2397	/* RESPECTIVELY, OF THE TRANSFORMED EIGENVECTORS */
2398	/* IN THEIR FIRST M COLUMNS. */
2399
2400	/* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
2401	/* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
2402	*/
2403
2404	/* THIS VERSION DATED AUGUST 1983. */
2405
2406	/* ------------------------------------------------------------------
2407	*/
2408
2409	/* Parameter adjustments */
2410	--scale;
2411	zi_dim1 = *nm;
2412	zi_offset = zi_dim1 + 1;
2413	zi -= zi_offset;
2414	zr_dim1 = *nm;
2415	zr_offset = zr_dim1 + 1;
2416	zr -= zr_offset;
2417
2418	/* Function Body */
2419	if (*m == 0) {
2420	goto L200;
2421	}
2422	if (igh == low) {
2423	goto L120;
2424	}
2425
2426	i_1 = *igh;
2427	for (i = *low; i <= i_1; ++i) {
2428	s = scale[i];
2429	/* .......... LEFT HAND EIGENVECTORS ARE BACK TRANSFORMED */
2430	/* IF THE FOREGOING STATEMENT IS REPLACED BY */
2431	/* S=1.0D0/SCALE(I). .......... */
2432	i_2 = *m;
2433	for (j = 1; j <= i_2; ++j) {
2434	zr[i + j * zr_dim1] *= s;
2435	zi[i + j * zi_dim1] *= s;
2436	/* L100: */
2437	}
2438
2439	/* L110: */
2440	}
2441	/* .......... FOR I=LOW-1 STEP -1 UNTIL 1, */
2442	/* IGH+1 STEP 1 UNTIL N DO -- .......... */
2443	L120:
2444	i_1 = *n;
2445	for (ii = 1; ii <= i_1; ++ii) {
2446	i = ii;
2447	if (i >= low && i <= igh) {
2448	goto L140;
2449	}
2450	if (i < *low) {
2451	i = *low - ii;
2452	}
2453	k = (integer) scale[i];
2454	if (k == i) {
2455	goto L140;
2456	}
2457
2458	i_2 = *m;
2459	for (j = 1; j <= i_2; ++j) {
2460	s = zr[i + j * zr_dim1];
2461	zr[i + j * zr_dim1] = zr[k + j * zr_dim1];
2462	zr[k + j * zr_dim1] = s;
2463	s = zi[i + j * zi_dim1];
2464	zi[i + j * zi_dim1] = zi[k + j * zi_dim1];
2465	zi[k + j * zi_dim1] = s;
2466	/* L130: */
2467	}
2468
2469	L140:
2470	;
2471	}
2472
2473	L200:
2474	return 0;
2475	} /* cbabk2_ */
2476
2477	/* Subroutine / int cbal_(integer nm, integer n, doublereal ar,
2478	doublereal ai, integer low, integer igh, doublereal scale)
2479	{
2480	/* System generated locals */
2481	integer ar_dim1, ar_offset, ai_dim1, ai_offset, i_1, i_2;
2482	doublereal d_1, d_2;
2483
2484	/* Local variables */
2485	static integer iexc;
2486	static doublereal c, f, g;
2487	static integer i, j, k, l, m;
2488	static doublereal r, s, radix, b2;
2489	static integer jj;
2490	static logical noconv;
2491
2492
2493
2494	/* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE */
2495	/* CBALANCE, WHICH IS A COMPLEX VERSION OF BALANCE, */
2496	/* NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH. */
2497	/* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971). */
2498
2499	/* THIS SUBROUTINE BALANCES A COMPLEX MATRIX AND ISOLATES */
2500	/* EIGENVALUES WHENEVER POSSIBLE. */
2501
2502	/* ON INPUT */
2503
2504	/* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
2505	/* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
2506	/* DIMENSION STATEMENT. */
2507
2508	/* N IS THE ORDER OF THE MATRIX. */
2509
2510	/* AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, */
2511	/* RESPECTIVELY, OF THE COMPLEX MATRIX TO BE BALANCED. */
2512
2513	/* ON OUTPUT */
2514
2515	/* AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, */
2516	/* RESPECTIVELY, OF THE BALANCED MATRIX. */
2517
2518	/* LOW AND IGH ARE TWO INTEGERS SUCH THAT AR(I,J) AND AI(I,J) */
2519	/* ARE EQUAL TO ZERO IF */
2520	/* (1) I IS GREATER THAN J AND */
2521	/* (2) J=1,...,LOW-1 OR I=IGH+1,...,N. */
2522
2523	/* SCALE CONTAINS INFORMATION DETERMINING THE */
2524	/* PERMUTATIONS AND SCALING FACTORS USED. */
2525
2526	/* SUPPOSE THAT THE PRINCIPAL SUBMATRIX IN ROWS LOW THROUGH IGH */
2527	/* HAS BEEN BALANCED, THAT P(J) DENOTES THE INDEX INTERCHANGED */
2528	/* WITH J DURING THE PERMUTATION STEP, AND THAT THE ELEMENTS */
2529	/* OF THE DIAGONAL MATRIX USED ARE DENOTED BY D(I,J). THEN */
2530	/* SCALE(J) = P(J), FOR J = 1,...,LOW-1 */
2531	/* = D(J,J) J = LOW,...,IGH */
2532	/* = P(J) J = IGH+1,...,N. */
2533	/* THE ORDER IN WHICH THE INTERCHANGES ARE MADE IS N TO IGH+1, */
2534	/* THEN 1 TO LOW-1. */
2535
2536	/* NOTE THAT 1 IS RETURNED FOR IGH IF IGH IS ZERO FORMALLY. */
2537
2538	/* THE ALGOL PROCEDURE EXC CONTAINED IN CBALANCE APPEARS IN */
2539	/* CBAL IN LINE. (NOTE THAT THE ALGOL ROLES OF IDENTIFIERS */
2540	/* K,L HAVE BEEN REVERSED.) */
2541
2542	/* ARITHMETIC IS REAL THROUGHOUT. */
2543
2544	/* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
2545	/* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
2546	*/
2547
2548	/* THIS VERSION DATED AUGUST 1983. */
2549
2550	/* ------------------------------------------------------------------
2551	*/
2552
2553	/* Parameter adjustments */
2554	--scale;
2555	ai_dim1 = *nm;
2556	ai_offset = ai_dim1 + 1;
2557	ai -= ai_offset;
2558	ar_dim1 = *nm;
2559	ar_offset = ar_dim1 + 1;
2560	ar -= ar_offset;
2561
2562	/* Function Body */
2563	radix = 16.;
2564
2565	b2 = radix * radix;
2566	k = 1;
2567	l = *n;
2568	goto L100;
2569	/* .......... IN-LINE PROCEDURE FOR ROW AND */
2570	/* COLUMN EXCHANGE .......... */
2571	L20:
2572	scale[m] = (doublereal) j;
2573	if (j == m) {
2574	goto L50;
2575	}
2576
2577	i_1 = l;
2578	for (i = 1; i <= i_1; ++i) {
2579	f = ar[i + j * ar_dim1];
2580	ar[i + j * ar_dim1] = ar[i + m * ar_dim1];
2581	ar[i + m * ar_dim1] = f;
2582	f = ai[i + j * ai_dim1];
2583	ai[i + j * ai_dim1] = ai[i + m * ai_dim1];
2584	ai[i + m * ai_dim1] = f;
2585	/* L30: */
2586	}
2587
2588	i_1 = *n;
2589	for (i = k; i <= i_1; ++i) {
2590	f = ar[j + i * ar_dim1];
2591	ar[j + i * ar_dim1] = ar[m + i * ar_dim1];
2592	ar[m + i * ar_dim1] = f;
2593	f = ai[j + i * ai_dim1];
2594	ai[j + i * ai_dim1] = ai[m + i * ai_dim1];
2595	ai[m + i * ai_dim1] = f;
2596	/* L40: */
2597	}
2598
2599	L50:
2600	switch (iexc) {
2601	case 1: goto L80;
2602	case 2: goto L130;
2603	}
2604	/* .......... SEARCH FOR ROWS ISOLATING AN EIGENVALUE */
2605	/* AND PUSH THEM DOWN .......... */
2606	L80:
2607	if (l == 1) {
2608	goto L280;
2609	}
2610	--l;
2611	/* .......... FOR J=L STEP -1 UNTIL 1 DO -- .......... */
2612	L100:
2613	i_1 = l;
2614	for (jj = 1; jj <= i_1; ++jj) {
2615	j = l + 1 - jj;
2616
2617	i_2 = l;
2618	for (i = 1; i <= i_2; ++i) {
2619	if (i == j) {
2620	goto L110;
2621	}
2622	if (ar[j + i * ar_dim1] != 0. \|\| ai[j + i * ai_dim1] != 0.) {
2623	goto L120;
2624	}
2625	L110:
2626	;
2627	}
2628
2629	m = l;
2630	iexc = 1;
2631	goto L20;
2632	L120:
2633	;
2634	}
2635
2636	goto L140;
2637	/* .......... SEARCH FOR COLUMNS ISOLATING AN EIGENVALUE */
2638	/* AND PUSH THEM LEFT .......... */
2639	L130:
2640	++k;
2641
2642	L140:
2643	i_1 = l;
2644	for (j = k; j <= i_1; ++j) {
2645
2646	i_2 = l;
2647	for (i = k; i <= i_2; ++i) {
2648	if (i == j) {
2649	goto L150;
2650	}
2651	if (ar[i + j * ar_dim1] != 0. \|\| ai[i + j * ai_dim1] != 0.) {
2652	goto L170;
2653	}
2654	L150:
2655	;
2656	}
2657
2658	m = k;
2659	iexc = 2;
2660	goto L20;
2661	L170:
2662	;
2663	}
2664	/* .......... NOW BALANCE THE SUBMATRIX IN ROWS K TO L .......... */
2665	i_1 = l;
2666	for (i = k; i <= i_1; ++i) {
2667	/* L180: */
2668	scale[i] = 1.;
2669	}
2670	/* .......... ITERATIVE LOOP FOR NORM REDUCTION .......... */
2671	L190:
2672	noconv = FALSE_;
2673
2674	i_1 = l;
2675	for (i = k; i <= i_1; ++i) {
2676	c = 0.;
2677	r = 0.;
2678
2679	i_2 = l;
2680	for (j = k; j <= i_2; ++j) {
2681	if (j == i) {
2682	goto L200;
2683	}
2684	c = c + (d_1 = ar[j + i * ar_dim1], abs(d_1)) + (d_2 = ai[j +
2685	i * ai_dim1], abs(d_2));
2686	r = r + (d_1 = ar[i + j * ar_dim1], abs(d_1)) + (d_2 = ai[i +
2687	j * ai_dim1], abs(d_2));
2688	L200:
2689	;
2690	}
2691	/* .......... GUARD AGAINST ZERO C OR R DUE TO UNDERFLOW .........
2692	. */
2693	if (c == 0. \|\| r == 0.) {
2694	goto L270;
2695	}
2696	g = r / radix;
2697	f = 1.;
2698	s = c + r;
2699	L210:
2700	if (c >= g) {
2701	goto L220;
2702	}
2703	f *= radix;
2704	c *= b2;
2705	goto L210;
2706	L220:
2707	g = r * radix;
2708	L230:
2709	if (c < g) {
2710	goto L240;
2711	}
2712	f /= radix;
2713	c /= b2;
2714	goto L230;
2715	/* .......... NOW BALANCE .......... */
2716	L240:
2717	if ((c + r) / f >= s * .95) {
2718	goto L270;
2719	}
2720	g = 1. / f;
2721	scale[i] *= f;
2722	noconv = TRUE_;
2723
2724	i_2 = *n;
2725	for (j = k; j <= i_2; ++j) {
2726	ar[i + j * ar_dim1] *= g;
2727	ai[i + j * ai_dim1] *= g;
2728	/* L250: */
2729	}
2730
2731	i_2 = l;
2732	for (j = 1; j <= i_2; ++j) {
2733	ar[j + i * ar_dim1] *= f;
2734	ai[j + i * ai_dim1] *= f;
2735	/* L260: */
2736	}
2737
2738	L270:
2739	;
2740	}
2741
2742	if (noconv) {
2743	goto L190;
2744	}
2745
2746	L280:
2747	*low = k;
2748	*igh = l;
2749	return 0;
2750	} /* cbal_ */
2751
2752	/* Subroutine / int cg_(integer nm, integer n, doublereal ar, doublereal *
2753	ai, doublereal wr, doublereal wi, integer matz, doublereal zr,
2754	doublereal zi, doublereal fv1, doublereal fv2, doublereal fv3,
2755	integer *ierr)
2756	{
2757	/* System generated locals */
2758	integer ar_dim1, ar_offset, ai_dim1, ai_offset, zr_dim1, zr_offset,
2759	zi_dim1, zi_offset;
2760
2761	/* Local variables */
2762	extern /* Subroutine / int cbal_(integer , integer , doublereal ,
2763	doublereal , integer , integer , doublereal ), corth_(integer
2764	, integer , integer , integer , doublereal , doublereal ,
2765	doublereal , doublereal ), comqr_(integer , integer , integer
2766	, integer , doublereal , doublereal , doublereal *,
2767	doublereal , integer ), cbabk2_(integer , integer , integer *,
2768	integer , doublereal , integer , doublereal , doublereal *),
2769	comqr2_(integer , integer , integer , integer , doublereal *,
2770	doublereal , doublereal , doublereal , doublereal ,
2771	doublereal , doublereal , doublereal , integer );
2772	static integer is1, is2;
2773
2774
2775
2776	/* THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF */
2777	/* SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK) */
2778	/* TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED) */
2779	/* OF A COMPLEX GENERAL MATRIX. */
2780
2781	/* ON INPUT */
2782
2783	/* NM MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL */
2784	/* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
2785	/* DIMENSION STATEMENT. */
2786
2787	/* N IS THE ORDER OF THE MATRIX A=(AR,AI). */
2788
2789	/* AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, */
2790	/* RESPECTIVELY, OF THE COMPLEX GENERAL MATRIX. */
2791
2792	/* MATZ IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF */
2793	/* ONLY EIGENVALUES ARE DESIRED. OTHERWISE IT IS SET TO */
2794	/* ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS. */
2795
2796	/* ON OUTPUT */
2797
2798	/* WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, */
2799	/* RESPECTIVELY, OF THE EIGENVALUES. */
2800
2801	/* ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, */
2802	/* RESPECTIVELY, OF THE EIGENVECTORS IF MATZ IS NOT ZERO. */
2803
2804	/* IERR IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR */
2805	/* COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR COMQR */
2806	/* AND COMQR2. THE NORMAL COMPLETION CODE IS ZERO. */
2807
2808	/* FV1, FV2, AND FV3 ARE TEMPORARY STORAGE ARRAYS. */
2809
2810	/* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
2811	/* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
2812	*/
2813
2814	/* THIS VERSION DATED AUGUST 1983. */
2815
2816	/* ------------------------------------------------------------------
2817	*/
2818
2819	/* Parameter adjustments */
2820	--fv3;
2821	--fv2;
2822	--fv1;
2823	zi_dim1 = *nm;
2824	zi_offset = zi_dim1 + 1;
2825	zi -= zi_offset;
2826	zr_dim1 = *nm;
2827	zr_offset = zr_dim1 + 1;
2828	zr -= zr_offset;
2829	--wi;
2830	--wr;
2831	ai_dim1 = *nm;
2832	ai_offset = ai_dim1 + 1;
2833	ai -= ai_offset;
2834	ar_dim1 = *nm;
2835	ar_offset = ar_dim1 + 1;
2836	ar -= ar_offset;
2837
2838	/* Function Body */
2839	if (n <= nm) {
2840	goto L10;
2841	}
2842	ierr = n * 10;
2843	goto L50;
2844
2845	L10:
2846	cbal_(nm, n, &ar[ar_offset], &ai[ai_offset], &is1, &is2, &fv1[1]);
2847	corth_(nm, n, &is1, &is2, &ar[ar_offset], &ai[ai_offset], &fv2[1], &fv3[1]
2848	);
2849	if (*matz != 0) {
2850	goto L20;
2851	}
2852	/* .......... FIND EIGENVALUES ONLY .......... */
2853	comqr_(nm, n, &is1, &is2, &ar[ar_offset], &ai[ai_offset], &wr[1], &wi[1],
2854	ierr);
2855	goto L50;
2856	/* .......... FIND BOTH EIGENVALUES AND EIGENVECTORS .......... */
2857	L20:
2858	comqr2_(nm, n, &is1, &is2, &fv2[1], &fv3[1], &ar[ar_offset], &ai[
2859	ai_offset], &wr[1], &wi[1], &zr[zr_offset], &zi[zi_offset], ierr);
2860	if (*ierr != 0) {
2861	goto L50;
2862	}
2863	cbabk2_(nm, n, &is1, &is2, &fv1[1], n, &zr[zr_offset], &zi[zi_offset]);
2864	L50:
2865	return 0;
2866	} /* cg_ */
2867
2868	/* Subroutine / int ch_(integer nm, integer n, doublereal ar, doublereal *
2869	ai, doublereal w, integer matz, doublereal zr, doublereal zi,
2870	doublereal fv1, doublereal fv2, doublereal fm1, integer ierr)
2871	{
2872	/* System generated locals */
2873	integer ar_dim1, ar_offset, ai_dim1, ai_offset, zr_dim1, zr_offset,
2874	zi_dim1, zi_offset, i_1, i_2;
2875
2876	/* Local variables */
2877	static integer i, j;
2878	extern /* Subroutine / int htridi_(integer , integer , doublereal ,
2879	doublereal , doublereal , doublereal , doublereal ,
2880	doublereal ), htribk_(integer , integer , doublereal ,
2881	doublereal , doublereal , integer , doublereal , doublereal *)
2882	, tqlrat_(integer , doublereal , doublereal , integer ),
2883	tql2_(integer , integer , doublereal , doublereal ,
2884	doublereal , integer );
2885
2886
2887
2888	/* THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF */
2889	/* SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK) */
2890	/* TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED) */
2891	/* OF A COMPLEX HERMITIAN MATRIX. */
2892
2893	/* ON INPUT */
2894
2895	/* NM MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL */
2896	/* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
2897	/* DIMENSION STATEMENT. */
2898
2899	/* N IS THE ORDER OF THE MATRIX A=(AR,AI). */
2900
2901	/* AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, */
2902	/* RESPECTIVELY, OF THE COMPLEX HERMITIAN MATRIX. */
2903
2904	/* MATZ IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF */
2905	/* ONLY EIGENVALUES ARE DESIRED. OTHERWISE IT IS SET TO */
2906	/* ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS. */
2907
2908	/* ON OUTPUT */
2909
2910	/* W CONTAINS THE EIGENVALUES IN ASCENDING ORDER. */
2911
2912	/* ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, */
2913	/* RESPECTIVELY, OF THE EIGENVECTORS IF MATZ IS NOT ZERO. */
2914
2915	/* IERR IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR */
2916	/* COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR TQLRAT */
2917	/* AND TQL2. THE NORMAL COMPLETION CODE IS ZERO. */
2918
2919	/* FV1, FV2, AND FM1 ARE TEMPORARY STORAGE ARRAYS. */
2920
2921	/* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
2922	/* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
2923	*/
2924
2925	/* THIS VERSION DATED AUGUST 1983. */
2926
2927	/* ------------------------------------------------------------------
2928	*/
2929
2930	/* Parameter adjustments */
2931	fm1 -= 3;
2932	--fv2;
2933	--fv1;
2934	zi_dim1 = *nm;
2935	zi_offset = zi_dim1 + 1;
2936	zi -= zi_offset;
2937	zr_dim1 = *nm;
2938	zr_offset = zr_dim1 + 1;
2939	zr -= zr_offset;
2940	--w;
2941	ai_dim1 = *nm;
2942	ai_offset = ai_dim1 + 1;
2943	ai -= ai_offset;
2944	ar_dim1 = *nm;
2945	ar_offset = ar_dim1 + 1;
2946	ar -= ar_offset;
2947
2948	/* Function Body */
2949	if (n <= nm) {
2950	goto L10;
2951	}
2952	ierr = n * 10;
2953	goto L50;
2954
2955	L10:
2956	htridi_(nm, n, &ar[ar_offset], &ai[ai_offset], &w[1], &fv1[1], &fv2[1], &
2957	fm1[3]);
2958	if (*matz != 0) {
2959	goto L20;
2960	}
2961	/* .......... FIND EIGENVALUES ONLY .......... */
2962	tqlrat_(n, &w[1], &fv2[1], ierr);
2963	goto L50;
2964	/* .......... FIND BOTH EIGENVALUES AND EIGENVECTORS .......... */
2965	L20:
2966	i_1 = *n;
2967	for (i = 1; i <= i_1; ++i) {
2968
2969	i_2 = *n;
2970	for (j = 1; j <= i_2; ++j) {
2971	zr[j + i * zr_dim1] = 0.;
2972	/* L30: */
2973	}
2974
2975	zr[i + i * zr_dim1] = 1.;
2976	/* L40: */
2977	}
2978
2979	tql2_(nm, n, &w[1], &fv1[1], &zr[zr_offset], ierr);
2980	if (*ierr != 0) {
2981	goto L50;
2982	}
2983	htribk_(nm, n, &ar[ar_offset], &ai[ai_offset], &fm1[3], n, &zr[zr_offset],
2984	&zi[zi_offset]);
2985	L50:
2986	return 0;
2987	} /* ch_ */
2988
2989	/* Subroutine / int cinvit_(integer nm, integer n, doublereal ar,
2990	doublereal ai, doublereal wr, doublereal wi, logical select,
2991	integer mm, integer m, doublereal zr, doublereal zi, integer *
2992	ierr, doublereal rm1, doublereal rm2, doublereal rv1, doublereal
2993	rv2)
2994	{
2995	/* System generated locals */
2996	integer ar_dim1, ar_offset, ai_dim1, ai_offset, zr_dim1, zr_offset,
2997	zi_dim1, zi_offset, rm1_dim1, rm1_offset, rm2_dim1, rm2_offset,
2998	i_1, i_2, i_3;
2999	doublereal d_1, d_2;
3000
3001	/* Builtin functions */
3002	double sqrt(doublereal);
3003
3004	/* Local variables */
3005	extern /* Subroutine / int cdiv_(doublereal , doublereal , doublereal
3006	, doublereal , doublereal , doublereal *);
3007	static doublereal norm;
3008	static integer i, j, k, s;
3009	static doublereal x, y, normv;
3010	static integer ii;
3011	static doublereal ilambd;
3012	static integer mp, uk;
3013	static doublereal rlambd;
3014	extern doublereal pythag_(doublereal , doublereal ), epslon_(doublereal
3015	*);
3016	static integer km1, ip1;
3017	static doublereal growto, ukroot;
3018	static integer its;
3019	static doublereal eps3;
3020
3021
3022
3023	/* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE CX INVIT */
3024	/* BY PETERS AND WILKINSON. */
3025	/* HANDBOOK FOR AUTO. COMP. VOL.II-LINEAR ALGEBRA, 418-439(1971). */
3026
3027	/* THIS SUBROUTINE FINDS THOSE EIGENVECTORS OF A COMPLEX UPPER */
3028	/* HESSENBERG MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES, */
3029	/* USING INVERSE ITERATION. */
3030
3031	/* ON INPUT */
3032
3033	/* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
3034	/* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
3035	/* DIMENSION STATEMENT. */
3036
3037	/* N IS THE ORDER OF THE MATRIX. */
3038
3039	/* AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, */
3040	/* RESPECTIVELY, OF THE HESSENBERG MATRIX. */
3041
3042	/* WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, RESPECTIVELY, */
3043	/* OF THE EIGENVALUES OF THE MATRIX. THE EIGENVALUES MUST BE */
3044	/* STORED IN A MANNER IDENTICAL TO THAT OF SUBROUTINE COMLR, */
3045	/* WHICH RECOGNIZES POSSIBLE SPLITTING OF THE MATRIX. */
3046
3047	/* SELECT SPECIFIES THE EIGENVECTORS TO BE FOUND. THE */
3048	/* EIGENVECTOR CORRESPONDING TO THE J-TH EIGENVALUE IS */
3049	/* SPECIFIED BY SETTING SELECT(J) TO .TRUE.. */
3050
3051	/* MM SHOULD BE SET TO AN UPPER BOUND FOR THE NUMBER OF */
3052	/* EIGENVECTORS TO BE FOUND. */
3053
3054	/* ON OUTPUT */
3055
3056	/* AR, AI, WI, AND SELECT ARE UNALTERED. */
3057
3058	/* WR MAY HAVE BEEN ALTERED SINCE CLOSE EIGENVALUES ARE PERTURBED
3059	*/
3060	/* SLIGHTLY IN SEARCHING FOR INDEPENDENT EIGENVECTORS. */
3061
3062	/* M IS THE NUMBER OF EIGENVECTORS ACTUALLY FOUND. */
3063
3064	/* ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, RESPECTIVELY, */
3065	/* OF THE EIGENVECTORS. THE EIGENVECTORS ARE NORMALIZED */
3066	/* SO THAT THE COMPONENT OF LARGEST MAGNITUDE IS 1. */
3067	/* ANY VECTOR WHICH FAILS THE ACCEPTANCE TEST IS SET TO ZERO. */
3068
3069	/* IERR IS SET TO */
3070	/* ZERO FOR NORMAL RETURN, */
3071	/* -(2*N+1) IF MORE THAN MM EIGENVECTORS HAVE BEEN SPECIFIED,
3072	*/
3073	/* -K IF THE ITERATION CORRESPONDING TO THE K-TH */
3074	/* VALUE FAILS, */
3075	/* -(N+K) IF BOTH ERROR SITUATIONS OCCUR. */
3076
3077	/* RM1, RM2, RV1, AND RV2 ARE TEMPORARY STORAGE ARRAYS. */
3078
3079	/* THE ALGOL PROCEDURE GUESSVEC APPEARS IN CINVIT IN LINE. */
3080
3081	/* CALLS CDIV FOR COMPLEX DIVISION. */
3082	/* CALLS PYTHAG FOR DSQRT(AA + BB) . */
3083
3084	/* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
3085	/* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
3086	*/
3087
3088	/* THIS VERSION DATED AUGUST 1983. */
3089
3090	/* ------------------------------------------------------------------
3091	*/
3092
3093	/* Parameter adjustments */
3094	--rv2;
3095	--rv1;
3096	rm2_dim1 = *n;
3097	rm2_offset = rm2_dim1 + 1;
3098	rm2 -= rm2_offset;
3099	rm1_dim1 = *n;
3100	rm1_offset = rm1_dim1 + 1;
3101	rm1 -= rm1_offset;
3102	--select;
3103	--wi;
3104	--wr;
3105	ai_dim1 = *nm;
3106	ai_offset = ai_dim1 + 1;
3107	ai -= ai_offset;
3108	ar_dim1 = *nm;
3109	ar_offset = ar_dim1 + 1;
3110	ar -= ar_offset;
3111	zi_dim1 = *nm;
3112	zi_offset = zi_dim1 + 1;
3113	zi -= zi_offset;
3114	zr_dim1 = *nm;
3115	zr_offset = zr_dim1 + 1;
3116	zr -= zr_offset;
3117
3118	/* Function Body */
3119	*ierr = 0;
3120	uk = 0;
3121	s = 1;
3122
3123	i_1 = *n;
3124	for (k = 1; k <= i_1; ++k) {
3125	if (! select[k]) {
3126	goto L980;
3127	}
3128	if (s > *mm) {
3129	goto L1000;
3130	}
3131	if (uk >= k) {
3132	goto L200;
3133	}
3134	/* .......... CHECK FOR POSSIBLE SPLITTING .......... */
3135	i_2 = *n;
3136	for (uk = k; uk <= i_2; ++uk) {
3137	if (uk == *n) {
3138	goto L140;
3139	}
3140	if (ar[uk + 1 + uk * ar_dim1] == 0. && ai[uk + 1 + uk * ai_dim1]
3141	== 0.) {
3142	goto L140;
3143	}
3144	/* L120: */
3145	}
3146	/* .......... COMPUTE INFINITY NORM OF LEADING UK BY UK */
3147	/* (HESSENBERG) MATRIX .......... */
3148	L140:
3149	norm = 0.;
3150	mp = 1;
3151
3152	i_2 = uk;
3153	for (i = 1; i <= i_2; ++i) {
3154	x = 0.;
3155
3156	i_3 = uk;
3157	for (j = mp; j <= i_3; ++j) {
3158	/* L160: */
3159	x += pythag_(&ar[i + j * ar_dim1], &ai[i + j * ai_dim1]);
3160	}
3161
3162	if (x > norm) {
3163	norm = x;
3164	}
3165	mp = i;
3166	/* L180: */
3167	}
3168	/* .......... EPS3 REPLACES ZERO PIVOT IN DECOMPOSITION */
3169	/* AND CLOSE ROOTS ARE MODIFIED BY EPS3 .......... */
3170	if (norm == 0.) {
3171	norm = 1.;
3172	}
3173	eps3 = epslon_(&norm);
3174	/* .......... GROWTO IS THE CRITERION FOR GROWTH .......... */
3175	ukroot = (doublereal) uk;
3176	ukroot = sqrt(ukroot);
3177	growto = .1 / ukroot;
3178	L200:
3179	rlambd = wr[k];
3180	ilambd = wi[k];
3181	if (k == 1) {
3182	goto L280;
3183	}
3184	km1 = k - 1;
3185	goto L240;
3186	/* .......... PERTURB EIGENVALUE IF IT IS CLOSE */
3187	/* TO ANY PREVIOUS EIGENVALUE .......... */
3188	L220:
3189	rlambd += eps3;
3190	/* .......... FOR I=K-1 STEP -1 UNTIL 1 DO -- .......... */
3191	L240:
3192	i_2 = km1;
3193	for (ii = 1; ii <= i_2; ++ii) {
3194	i = k - ii;
3195	if (select[i] && (d_1 = wr[i] - rlambd, abs(d_1)) < eps3 && (
3196	d_2 = wi[i] - ilambd, abs(d_2)) < eps3) {
3197	goto L220;
3198	}
3199	/* L260: */
3200	}
3201
3202	wr[k] = rlambd;
3203	/* .......... FORM UPPER HESSENBERG (AR,AI)-(RLAMBD,ILAMBD)I /
3204	/* AND INITIAL COMPLEX VECTOR .......... */
3205	L280:
3206	mp = 1;
3207
3208	i_2 = uk;
3209	for (i = 1; i <= i_2; ++i) {
3210
3211	i_3 = uk;
3212	for (j = mp; j <= i_3; ++j) {
3213	rm1[i + j * rm1_dim1] = ar[i + j * ar_dim1];
3214	rm2[i + j * rm2_dim1] = ai[i + j * ai_dim1];
3215	/* L300: */
3216	}
3217
3218	rm1[i + i * rm1_dim1] -= rlambd;
3219	rm2[i + i * rm2_dim1] -= ilambd;
3220	mp = i;
3221	rv1[i] = eps3;
3222	/* L320: */
3223	}
3224	/* .......... TRIANGULAR DECOMPOSITION WITH INTERCHANGES, */
3225	/* REPLACING ZERO PIVOTS BY EPS3 .......... */
3226	if (uk == 1) {
3227	goto L420;
3228	}
3229
3230	i_2 = uk;
3231	for (i = 2; i <= i_2; ++i) {
3232	mp = i - 1;
3233	if (pythag_(&rm1[i + mp * rm1_dim1], &rm2[i + mp * rm2_dim1]) <=
3234	pythag_(&rm1[mp + mp * rm1_dim1], &rm2[mp + mp * rm2_dim1]
3235	)) {
3236	goto L360;
3237	}
3238
3239	i_3 = uk;
3240	for (j = mp; j <= i_3; ++j) {
3241	y = rm1[i + j * rm1_dim1];
3242	rm1[i + j * rm1_dim1] = rm1[mp + j * rm1_dim1];
3243	rm1[mp + j * rm1_dim1] = y;
3244	y = rm2[i + j * rm2_dim1];
3245	rm2[i + j * rm2_dim1] = rm2[mp + j * rm2_dim1];
3246	rm2[mp + j * rm2_dim1] = y;
3247	/* L340: */
3248	}
3249
3250	L360:
3251	if (rm1[mp + mp * rm1_dim1] == 0. && rm2[mp + mp * rm2_dim1] ==
3252	0.) {
3253	rm1[mp + mp * rm1_dim1] = eps3;
3254	}
3255	cdiv_(&rm1[i + mp * rm1_dim1], &rm2[i + mp * rm2_dim1], &rm1[mp +
3256	mp * rm1_dim1], &rm2[mp + mp * rm2_dim1], &x, &y);
3257	if (x == 0. && y == 0.) {
3258	goto L400;
3259	}
3260
3261	i_3 = uk;
3262	for (j = i; j <= i_3; ++j) {
3263	rm1[i + j * rm1_dim1] = rm1[i + j * rm1_dim1] - x * rm1[mp +
3264	j * rm1_dim1] + y * rm2[mp + j * rm2_dim1];
3265	rm2[i + j * rm2_dim1] = rm2[i + j * rm2_dim1] - x * rm2[mp +
3266	j * rm2_dim1] - y * rm1[mp + j * rm1_dim1];
3267	/* L380: */
3268	}
3269
3270	L400:
3271	;
3272	}
3273
3274	L420:
3275	if (rm1[uk + uk * rm1_dim1] == 0. && rm2[uk + uk * rm2_dim1] == 0.) {
3276	rm1[uk + uk * rm1_dim1] = eps3;
3277	}
3278	its = 0;
3279	/* .......... BACK SUBSTITUTION */
3280	/* FOR I=UK STEP -1 UNTIL 1 DO -- .......... */
3281	L660:
3282	i_2 = uk;
3283	for (ii = 1; ii <= i_2; ++ii) {
3284	i = uk + 1 - ii;
3285	x = rv1[i];
3286	y = 0.;
3287	if (i == uk) {
3288	goto L700;
3289	}
3290	ip1 = i + 1;
3291
3292	i_3 = uk;
3293	for (j = ip1; j <= i_3; ++j) {
3294	x = x - rm1[i + j * rm1_dim1] * rv1[j] + rm2[i + j * rm2_dim1]
3295	* rv2[j];
3296	y = y - rm1[i + j * rm1_dim1] * rv2[j] - rm2[i + j * rm2_dim1]
3297	* rv1[j];
3298	/* L680: */
3299	}
3300
3301	L700:
3302	cdiv_(&x, &y, &rm1[i + i * rm1_dim1], &rm2[i + i * rm2_dim1], &
3303	rv1[i], &rv2[i]);
3304	/* L720: */
3305	}
3306	/* .......... ACCEPTANCE TEST FOR EIGENVECTOR */
3307	/* AND NORMALIZATION .......... */
3308	++its;
3309	norm = 0.;
3310	normv = 0.;
3311
3312	i_2 = uk;
3313	for (i = 1; i <= i_2; ++i) {
3314	x = pythag_(&rv1[i], &rv2[i]);
3315	if (normv >= x) {
3316	goto L760;
3317	}
3318	normv = x;
3319	j = i;
3320	L760:
3321	norm += x;
3322	/* L780: */
3323	}
3324
3325	if (norm < growto) {
3326	goto L840;
3327	}
3328	/* .......... ACCEPT VECTOR .......... */
3329	x = rv1[j];
3330	y = rv2[j];
3331
3332	i_2 = uk;
3333	for (i = 1; i <= i_2; ++i) {
3334	cdiv_(&rv1[i], &rv2[i], &x, &y, &zr[i + s * zr_dim1], &zi[i + s *
3335	zi_dim1]);
3336	/* L820: */
3337	}
3338
3339	if (uk == *n) {
3340	goto L940;
3341	}
3342	j = uk + 1;
3343	goto L900;
3344	/* .......... IN-LINE PROCEDURE FOR CHOOSING */
3345	/* A NEW STARTING VECTOR .......... */
3346	L840:
3347	if (its >= uk) {
3348	goto L880;
3349	}
3350	x = ukroot;
3351	y = eps3 / (x + 1.);
3352	rv1[1] = eps3;
3353
3354	i_2 = uk;
3355	for (i = 2; i <= i_2; ++i) {
3356	/* L860: */
3357	rv1[i] = y;
3358	}
3359
3360	j = uk - its + 1;
3361	rv1[j] -= eps3 * x;
3362	goto L660;
3363	/* .......... SET ERROR -- UNACCEPTED EIGENVECTOR .......... */
3364	L880:
3365	j = 1;
3366	*ierr = -k;
3367	/* .......... SET REMAINING VECTOR COMPONENTS TO ZERO ..........
3368	*/
3369	L900:
3370	i_2 = *n;
3371	for (i = j; i <= i_2; ++i) {
3372	zr[i + s * zr_dim1] = 0.;
3373	zi[i + s * zi_dim1] = 0.;
3374	/* L920: */
3375	}
3376
3377	L940:
3378	++s;
3379	L980:
3380	;
3381	}
3382
3383	goto L1001;
3384	/* .......... SET ERROR -- UNDERESTIMATE OF EIGENVECTOR */
3385	/* SPACE REQUIRED .......... */
3386	L1000:
3387	if (*ierr != 0) {
3388	ierr -= n;
3389	}
3390	if (*ierr == 0) {
3391	ierr = -((n << 1) + 1);
3392	}
3393	L1001:
3394	*m = s - 1;
3395	return 0;
3396	} /* cinvit_ */
3397
3398	/* Subroutine / int combak_(integer nm, integer low, integer igh,
3399	doublereal ar, doublereal ai, integer int_, integer m,
3400	doublereal zr, doublereal zi)
3401	{
3402	/* System generated locals */
3403	integer ar_dim1, ar_offset, ai_dim1, ai_offset, zr_dim1, zr_offset,
3404	zi_dim1, zi_offset, i_1, i_2, i_3;
3405
3406	/* Local variables */
3407	static integer i, j, la, mm, mp;
3408	static doublereal xi, xr;
3409	static integer kp1, mp1;
3410
3411
3412
3413	/* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE COMBAK, */
3414	/* NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON. */
3415	/* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971). */
3416
3417	/* THIS SUBROUTINE FORMS THE EIGENVECTORS OF A COMPLEX GENERAL */
3418	/* MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING */
3419	/* UPPER HESSENBERG MATRIX DETERMINED BY COMHES. */
3420
3421	/* ON INPUT */
3422
3423	/* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
3424	/* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
3425	/* DIMENSION STATEMENT. */
3426
3427	/* LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
3428	/* SUBROUTINE CBAL. IF CBAL HAS NOT BEEN USED, */
3429	/* SET LOW=1 AND IGH EQUAL TO THE ORDER OF THE MATRIX. */
3430
3431	/* AR AND AI CONTAIN THE MULTIPLIERS WHICH WERE USED IN THE */
3432	/* REDUCTION BY COMHES IN THEIR LOWER TRIANGLES */
3433	/* BELOW THE SUBDIAGONAL. */
3434
3435	/* INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS */
3436	/* INTERCHANGED IN THE REDUCTION BY COMHES. */
3437	/* ONLY ELEMENTS LOW THROUGH IGH ARE USED. */
3438
3439	/* M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED. */
3440
3441	/* ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, */
3442	/* RESPECTIVELY, OF THE EIGENVECTORS TO BE */
3443	/* BACK TRANSFORMED IN THEIR FIRST M COLUMNS. */
3444
3445	/* ON OUTPUT */
3446
3447	/* ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, */
3448	/* RESPECTIVELY, OF THE TRANSFORMED EIGENVECTORS */
3449	/* IN THEIR FIRST M COLUMNS. */
3450
3451	/* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
3452	/* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
3453	*/
3454
3455	/* THIS VERSION DATED AUGUST 1983. */
3456
3457	/* ------------------------------------------------------------------
3458	*/
3459
3460	/* Parameter adjustments */
3461	--int_;
3462	ai_dim1 = *nm;
3463	ai_offset = ai_dim1 + 1;
3464	ai -= ai_offset;
3465	ar_dim1 = *nm;
3466	ar_offset = ar_dim1 + 1;
3467	ar -= ar_offset;
3468	zi_dim1 = *nm;
3469	zi_offset = zi_dim1 + 1;
3470	zi -= zi_offset;
3471	zr_dim1 = *nm;
3472	zr_offset = zr_dim1 + 1;
3473	zr -= zr_offset;
3474
3475	/* Function Body */
3476	if (*m == 0) {
3477	goto L200;
3478	}
3479	la = *igh - 1;
3480	kp1 = *low + 1;
3481	if (la < kp1) {
3482	goto L200;
3483	}
3484	/* .......... FOR MP=IGH-1 STEP -1 UNTIL LOW+1 DO -- .......... */
3485	i_1 = la;
3486	for (mm = kp1; mm <= i_1; ++mm) {
3487	mp = low + igh - mm;
3488	mp1 = mp + 1;
3489
3490	i_2 = *igh;
3491	for (i = mp1; i <= i_2; ++i) {
3492	xr = ar[i + (mp - 1) * ar_dim1];
3493	xi = ai[i + (mp - 1) * ai_dim1];
3494	if (xr == 0. && xi == 0.) {
3495	goto L110;
3496	}
3497
3498	i_3 = *m;
3499	for (j = 1; j <= i_3; ++j) {
3500	zr[i + j * zr_dim1] = zr[i + j * zr_dim1] + xr * zr[mp + j *
3501	zr_dim1] - xi * zi[mp + j * zi_dim1];
3502	zi[i + j * zi_dim1] = zi[i + j * zi_dim1] + xr * zi[mp + j *
3503	zi_dim1] + xi * zr[mp + j * zr_dim1];
3504	/* L100: */
3505	}
3506
3507	L110:
3508	;
3509	}
3510
3511	i = int_[mp];
3512	if (i == mp) {
3513	goto L140;
3514	}
3515
3516	i_2 = *m;
3517	for (j = 1; j <= i_2; ++j) {
3518	xr = zr[i + j * zr_dim1];
3519	zr[i + j * zr_dim1] = zr[mp + j * zr_dim1];
3520	zr[mp + j * zr_dim1] = xr;
3521	xi = zi[i + j * zi_dim1];
3522	zi[i + j * zi_dim1] = zi[mp + j * zi_dim1];
3523	zi[mp + j * zi_dim1] = xi;
3524	/* L130: */
3525	}
3526
3527	L140:
3528	;
3529	}
3530
3531	L200:
3532	return 0;
3533	} /* combak_ */
3534
3535	/* Subroutine / int comhes_(integer nm, integer n, integer low, integer *
3536	igh, doublereal ar, doublereal ai, integer *int_)
3537	{
3538	/* System generated locals */
3539	integer ar_dim1, ar_offset, ai_dim1, ai_offset, i_1, i_2, i_3;
3540	doublereal d_1, d_2;
3541
3542	/* Local variables */
3543	extern /* Subroutine / int cdiv_(doublereal , doublereal , doublereal
3544	, doublereal , doublereal , doublereal *);
3545	static integer i, j, m, la;
3546	static doublereal xi, yi, xr, yr;
3547	static integer mm1, kp1, mp1;
3548
3549
3550
3551	/* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE COMHES, */
3552	/* NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON. */
3553	/* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971). */
3554
3555	/* GIVEN A COMPLEX GENERAL MATRIX, THIS SUBROUTINE */
3556	/* REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS */
3557	/* LOW THROUGH IGH TO UPPER HESSENBERG FORM BY */
3558	/* STABILIZED ELEMENTARY SIMILARITY TRANSFORMATIONS. */
3559
3560	/* ON INPUT */
3561
3562	/* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
3563	/* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
3564	/* DIMENSION STATEMENT. */
3565
3566	/* N IS THE ORDER OF THE MATRIX. */
3567
3568	/* LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
3569	/* SUBROUTINE CBAL. IF CBAL HAS NOT BEEN USED, */
3570	/* SET LOW=1, IGH=N. */
3571
3572	/* AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, */
3573	/* RESPECTIVELY, OF THE COMPLEX INPUT MATRIX. */
3574
3575	/* ON OUTPUT */
3576
3577	/* AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, */
3578	/* RESPECTIVELY, OF THE HESSENBERG MATRIX. THE */
3579	/* MULTIPLIERS WHICH WERE USED IN THE REDUCTION */
3580	/* ARE STORED IN THE REMAINING TRIANGLES UNDER THE */
3581	/* HESSENBERG MATRIX. */
3582
3583	/* INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS */
3584	/* INTERCHANGED IN THE REDUCTION. */
3585	/* ONLY ELEMENTS LOW THROUGH IGH ARE USED. */
3586
3587	/* CALLS CDIV FOR COMPLEX DIVISION. */
3588
3589	/* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
3590	/* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
3591	*/
3592
3593	/* THIS VERSION DATED AUGUST 1983. */
3594
3595	/* ------------------------------------------------------------------
3596	*/
3597
3598	/* Parameter adjustments */
3599	ai_dim1 = *nm;
3600	ai_offset = ai_dim1 + 1;
3601	ai -= ai_offset;
3602	ar_dim1 = *nm;
3603	ar_offset = ar_dim1 + 1;
3604	ar -= ar_offset;
3605	--int_;
3606
3607	/* Function Body */
3608	la = *igh - 1;
3609	kp1 = *low + 1;
3610	if (la < kp1) {
3611	goto L200;
3612	}
3613
3614	i_1 = la;
3615	for (m = kp1; m <= i_1; ++m) {
3616	mm1 = m - 1;
3617	xr = 0.;
3618	xi = 0.;
3619	i = m;
3620
3621	i_2 = *igh;
3622	for (j = m; j <= i_2; ++j) {
3623	if ((d_1 = ar[j + mm1 * ar_dim1], abs(d_1)) + (d_2 = ai[j +
3624	mm1 * ai_dim1], abs(d_2)) <= abs(xr) + abs(xi)) {
3625	goto L100;
3626	}
3627	xr = ar[j + mm1 * ar_dim1];
3628	xi = ai[j + mm1 * ai_dim1];
3629	i = j;
3630	L100:
3631	;
3632	}
3633
3634	int_[m] = i;
3635	if (i == m) {
3636	goto L130;
3637	}
3638	/* .......... INTERCHANGE ROWS AND COLUMNS OF AR AND AI ..........
3639	*/
3640	i_2 = *n;
3641	for (j = mm1; j <= i_2; ++j) {
3642	yr = ar[i + j * ar_dim1];
3643	ar[i + j * ar_dim1] = ar[m + j * ar_dim1];
3644	ar[m + j * ar_dim1] = yr;
3645	yi = ai[i + j * ai_dim1];
3646	ai[i + j * ai_dim1] = ai[m + j * ai_dim1];
3647	ai[m + j * ai_dim1] = yi;
3648	/* L110: */
3649	}
3650
3651	i_2 = *igh;
3652	for (j = 1; j <= i_2; ++j) {
3653	yr = ar[j + i * ar_dim1];
3654	ar[j + i * ar_dim1] = ar[j + m * ar_dim1];
3655	ar[j + m * ar_dim1] = yr;
3656	yi = ai[j + i * ai_dim1];
3657	ai[j + i * ai_dim1] = ai[j + m * ai_dim1];
3658	ai[j + m * ai_dim1] = yi;
3659	/* L120: */
3660	}
3661	/* .......... END INTERCHANGE .......... */
3662	L130:
3663	if (xr == 0. && xi == 0.) {
3664	goto L180;
3665	}
3666	mp1 = m + 1;
3667
3668	i_2 = *igh;
3669	for (i = mp1; i <= i_2; ++i) {
3670	yr = ar[i + mm1 * ar_dim1];
3671	yi = ai[i + mm1 * ai_dim1];
3672	if (yr == 0. && yi == 0.) {
3673	goto L160;
3674	}
3675	cdiv_(&yr, &yi, &xr, &xi, &yr, &yi);
3676	ar[i + mm1 * ar_dim1] = yr;
3677	ai[i + mm1 * ai_dim1] = yi;
3678
3679	i_3 = *n;
3680	for (j = m; j <= i_3; ++j) {
3681	ar[i + j * ar_dim1] = ar[i + j * ar_dim1] - yr * ar[m + j *
3682	ar_dim1] + yi * ai[m + j * ai_dim1];
3683	ai[i + j * ai_dim1] = ai[i + j * ai_dim1] - yr * ai[m + j *
3684	ai_dim1] - yi * ar[m + j * ar_dim1];
3685	/* L140: */
3686	}
3687
3688	i_3 = *igh;
3689	for (j = 1; j <= i_3; ++j) {
3690	ar[j + m * ar_dim1] = ar[j + m * ar_dim1] + yr * ar[j + i *
3691	ar_dim1] - yi * ai[j + i * ai_dim1];
3692	ai[j + m * ai_dim1] = ai[j + m * ai_dim1] + yr * ai[j + i *
3693	ai_dim1] + yi * ar[j + i * ar_dim1];
3694	/* L150: */
3695	}
3696
3697	L160:
3698	;
3699	}
3700
3701	L180:
3702	;
3703	}
3704
3705	L200:
3706	return 0;
3707	} /* comhes_ */
3708
3709	/* Subroutine / int comlr_(integer nm, integer n, integer low, integer *
3710	igh, doublereal hr, doublereal hi, doublereal wr, doublereal wi,
3711	integer *ierr)
3712	{
3713	/* System generated locals */
3714	integer hr_dim1, hr_offset, hi_dim1, hi_offset, i_1, i_2;
3715	doublereal d_1, d_2, d_3, d_4;
3716
3717	/* Local variables */
3718	extern /* Subroutine / int cdiv_(doublereal , doublereal , doublereal
3719	, doublereal , doublereal , doublereal *);
3720	static integer i, j, l, m, en, ll, mm;
3721	static doublereal si, ti, xi, yi, sr, tr, xr, yr;
3722	static integer im1;
3723	extern /* Subroutine / int csroot_(doublereal , doublereal *,
3724	doublereal , doublereal );
3725	static integer mp1, itn, its;
3726	static doublereal zzi, zzr;
3727	static integer enm1;
3728	static doublereal tst1, tst2;
3729
3730
3731
3732	/* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE COMLR, */
3733	/* NUM. MATH. 12, 369-376(1968) BY MARTIN AND WILKINSON. */
3734	/* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 396-403(1971). */
3735
3736	/* THIS SUBROUTINE FINDS THE EIGENVALUES OF A COMPLEX */
3737	/* UPPER HESSENBERG MATRIX BY THE MODIFIED LR METHOD. */
3738
3739	/* ON INPUT */
3740
3741	/* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
3742	/* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
3743	/* DIMENSION STATEMENT. */
3744
3745	/* N IS THE ORDER OF THE MATRIX. */
3746
3747	/* LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
3748	/* SUBROUTINE CBAL. IF CBAL HAS NOT BEEN USED, */
3749	/* SET LOW=1, IGH=N. */
3750
3751	/* HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS, */
3752	/* RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX. */
3753	/* THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN THE */
3754	/* MULTIPLIERS WHICH WERE USED IN THE REDUCTION BY COMHES, */
3755	/* IF PERFORMED. */
3756
3757	/* ON OUTPUT */
3758
3759	/* THE UPPER HESSENBERG PORTIONS OF HR AND HI HAVE BEEN */
3760	/* DESTROYED. THEREFORE, THEY MUST BE SAVED BEFORE */
3761	/* CALLING COMLR IF SUBSEQUENT CALCULATION OF */
3762	/* EIGENVECTORS IS TO BE PERFORMED. */
3763
3764	/* WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, */
3765	/* RESPECTIVELY, OF THE EIGENVALUES. IF AN ERROR */
3766	/* EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT */
3767	/* FOR INDICES IERR+1,...,N. */
3768
3769	/* IERR IS SET TO */
3770	/* ZERO FOR NORMAL RETURN, */
3771	/* J IF THE LIMIT OF 30N ITERATIONS IS EXHAUSTED /
3772	/* WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. */
3773
3774	/* CALLS CDIV FOR COMPLEX DIVISION. */
3775	/* CALLS CSROOT FOR COMPLEX SQUARE ROOT. */
3776
3777	/* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
3778	/* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
3779	*/
3780
3781	/* THIS VERSION DATED AUGUST 1983. */
3782
3783	/* ------------------------------------------------------------------
3784	*/
3785
3786	/* Parameter adjustments */
3787	--wi;
3788	--wr;
3789	hi_dim1 = *nm;
3790	hi_offset = hi_dim1 + 1;
3791	hi -= hi_offset;
3792	hr_dim1 = *nm;
3793	hr_offset = hr_dim1 + 1;
3794	hr -= hr_offset;
3795
3796	/* Function Body */
3797	*ierr = 0;
3798	/* .......... STORE ROOTS ISOLATED BY CBAL .......... */
3799	i_1 = *n;
3800	for (i = 1; i <= i_1; ++i) {
3801	if (i >= low && i <= igh) {
3802	goto L200;
3803	}
3804	wr[i] = hr[i + i * hr_dim1];
3805	wi[i] = hi[i + i * hi_dim1];
3806	L200:
3807	;
3808	}
3809
3810	en = *igh;
3811	tr = 0.;
3812	ti = 0.;
3813	itn = n 30;
3814	/* .......... SEARCH FOR NEXT EIGENVALUE .......... */
3815	L220:
3816	if (en < *low) {
3817	goto L1001;
3818	}
3819	its = 0;
3820	enm1 = en - 1;
3821	/* .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT */
3822	/* FOR L=EN STEP -1 UNTIL LOW D0 -- .......... */
3823	L240:
3824	i_1 = en;
3825	for (ll = *low; ll <= i_1; ++ll) {
3826	l = en + *low - ll;
3827	if (l == *low) {
3828	goto L300;
3829	}
3830	tst1 = (d_1 = hr[l - 1 + (l - 1) * hr_dim1], abs(d_1)) + (d_2 = hi[
3831	l - 1 + (l - 1) * hi_dim1], abs(d_2)) + (d_3 = hr[l + l *
3832	hr_dim1], abs(d_3)) + (d_4 = hi[l + l * hi_dim1], abs(d_4))
3833	;
3834	tst2 = tst1 + (d_1 = hr[l + (l - 1) * hr_dim1], abs(d_1)) + (d_2 =
3835	hi[l + (l - 1) * hi_dim1], abs(d_2));
3836	if (tst2 == tst1) {
3837	goto L300;
3838	}
3839	/* L260: */
3840	}
3841	/* .......... FORM SHIFT .......... */
3842	L300:
3843	if (l == en) {
3844	goto L660;
3845	}
3846	if (itn == 0) {
3847	goto L1000;
3848	}
3849	if (its == 10 \|\| its == 20) {
3850	goto L320;
3851	}
3852	sr = hr[en + en * hr_dim1];
3853	si = hi[en + en * hi_dim1];
3854	xr = hr[enm1 + en * hr_dim1] * hr[en + enm1 * hr_dim1] - hi[enm1 + en *
3855	hi_dim1] * hi[en + enm1 * hi_dim1];
3856	xi = hr[enm1 + en * hr_dim1] * hi[en + enm1 * hi_dim1] + hi[enm1 + en *
3857	hi_dim1] * hr[en + enm1 * hr_dim1];
3858	if (xr == 0. && xi == 0.) {
3859	goto L340;
3860	}
3861	yr = (hr[enm1 + enm1 * hr_dim1] - sr) / 2.;
3862	yi = (hi[enm1 + enm1 * hi_dim1] - si) / 2.;
3863	/* Computing 2nd power */
3864	d_2 = yr;
3865	/* Computing 2nd power */
3866	d_3 = yi;
3867	d_1 = d_2 * d_2 - d_3 * d_3 + xr;
3868	d_4 = yr * 2. * yi + xi;
3869	csroot_(&d_1, &d_4, &zzr, &zzi);
3870	if (yr * zzr + yi * zzi >= 0.) {
3871	goto L310;
3872	}
3873	zzr = -zzr;
3874	zzi = -zzi;
3875	L310:
3876	d_1 = yr + zzr;
3877	d_2 = yi + zzi;
3878	cdiv_(&xr, &xi, &d_1, &d_2, &xr, &xi);
3879	sr -= xr;
3880	si -= xi;
3881	goto L340;
3882	/* .......... FORM EXCEPTIONAL SHIFT .......... */
3883	L320:
3884	sr = (d_1 = hr[en + enm1 * hr_dim1], abs(d_1)) + (d_2 = hr[enm1 + (en
3885	- 2) * hr_dim1], abs(d_2));
3886	si = (d_1 = hi[en + enm1 * hi_dim1], abs(d_1)) + (d_2 = hi[enm1 + (en
3887	- 2) * hi_dim1], abs(d_2));
3888
3889	L340:
3890	i_1 = en;
3891	for (i = *low; i <= i_1; ++i) {
3892	hr[i + i * hr_dim1] -= sr;
3893	hi[i + i * hi_dim1] -= si;
3894	/* L360: */
3895	}
3896
3897	tr += sr;
3898	ti += si;
3899	++its;
3900	--itn;
3901	/* .......... LOOK FOR TWO CONSECUTIVE SMALL */
3902	/* SUB-DIAGONAL ELEMENTS .......... */
3903	xr = (d_1 = hr[enm1 + enm1 * hr_dim1], abs(d_1)) + (d_2 = hi[enm1 +
3904	enm1 * hi_dim1], abs(d_2));
3905	yr = (d_1 = hr[en + enm1 * hr_dim1], abs(d_1)) + (d_2 = hi[en + enm1 *
3906	hi_dim1], abs(d_2));
3907	zzr = (d_1 = hr[en + en * hr_dim1], abs(d_1)) + (d_2 = hi[en + en *
3908	hi_dim1], abs(d_2));
3909	/* .......... FOR M=EN-1 STEP -1 UNTIL L DO -- .......... */
3910	i_1 = enm1;
3911	for (mm = l; mm <= i_1; ++mm) {
3912	m = enm1 + l - mm;
3913	if (m == l) {
3914	goto L420;
3915	}
3916	yi = yr;
3917	yr = (d_1 = hr[m + (m - 1) * hr_dim1], abs(d_1)) + (d_2 = hi[m + (
3918	m - 1) * hi_dim1], abs(d_2));
3919	xi = zzr;
3920	zzr = xr;
3921	xr = (d_1 = hr[m - 1 + (m - 1) * hr_dim1], abs(d_1)) + (d_2 = hi[m
3922	- 1 + (m - 1) * hi_dim1], abs(d_2));
3923	tst1 = zzr / yi * (zzr + xr + xi);
3924	tst2 = tst1 + yr;
3925	if (tst2 == tst1) {
3926	goto L420;
3927	}
3928	/* L380: */
3929	}
3930	/* .......... TRIANGULAR DECOMPOSITION H=LR .......... /
3931	L420:
3932	mp1 = m + 1;
3933
3934	i_1 = en;
3935	for (i = mp1; i <= i_1; ++i) {
3936	im1 = i - 1;
3937	xr = hr[im1 + im1 * hr_dim1];
3938	xi = hi[im1 + im1 * hi_dim1];
3939	yr = hr[i + im1 * hr_dim1];
3940	yi = hi[i + im1 * hi_dim1];
3941	if (abs(xr) + abs(xi) >= abs(yr) + abs(yi)) {
3942	goto L460;
3943	}
3944	/* .......... INTERCHANGE ROWS OF HR AND HI .......... */
3945	i_2 = en;
3946	for (j = im1; j <= i_2; ++j) {
3947	zzr = hr[im1 + j * hr_dim1];
3948	hr[im1 + j * hr_dim1] = hr[i + j * hr_dim1];
3949	hr[i + j * hr_dim1] = zzr;
3950	zzi = hi[im1 + j * hi_dim1];
3951	hi[im1 + j * hi_dim1] = hi[i + j * hi_dim1];
3952	hi[i + j * hi_dim1] = zzi;
3953	/* L440: */
3954	}
3955
3956	cdiv_(&xr, &xi, &yr, &yi, &zzr, &zzi);
3957	wr[i] = 1.;
3958	goto L480;
3959	L460:
3960	cdiv_(&yr, &yi, &xr, &xi, &zzr, &zzi);
3961	wr[i] = -1.;
3962	L480:
3963	hr[i + im1 * hr_dim1] = zzr;
3964	hi[i + im1 * hi_dim1] = zzi;
3965
3966	i_2 = en;
3967	for (j = i; j <= i_2; ++j) {
3968	hr[i + j * hr_dim1] = hr[i + j * hr_dim1] - zzr * hr[im1 + j *
3969	hr_dim1] + zzi * hi[im1 + j * hi_dim1];
3970	hi[i + j * hi_dim1] = hi[i + j * hi_dim1] - zzr * hi[im1 + j *
3971	hi_dim1] - zzi * hr[im1 + j * hr_dim1];
3972	/* L500: */
3973	}
3974
3975	/* L520: */
3976	}
3977	/* .......... COMPOSITION RL=H .......... /
3978	i_1 = en;
3979	for (j = mp1; j <= i_1; ++j) {
3980	xr = hr[j + (j - 1) * hr_dim1];
3981	xi = hi[j + (j - 1) * hi_dim1];
3982	hr[j + (j - 1) * hr_dim1] = 0.;
3983	hi[j + (j - 1) * hi_dim1] = 0.;
3984	/* .......... INTERCHANGE COLUMNS OF HR AND HI, */
3985	/* IF NECESSARY .......... */
3986	if (wr[j] <= 0.) {
3987	goto L580;
3988	}
3989
3990	i_2 = j;
3991	for (i = l; i <= i_2; ++i) {
3992	zzr = hr[i + (j - 1) * hr_dim1];
3993	hr[i + (j - 1) * hr_dim1] = hr[i + j * hr_dim1];
3994	hr[i + j * hr_dim1] = zzr;
3995	zzi = hi[i + (j - 1) * hi_dim1];
3996	hi[i + (j - 1) * hi_dim1] = hi[i + j * hi_dim1];
3997	hi[i + j * hi_dim1] = zzi;
3998	/* L540: */
3999	}
4000
4001	L580:
4002	i_2 = j;
4003	for (i = l; i <= i_2; ++i) {
4004	hr[i + (j - 1) * hr_dim1] = hr[i + (j - 1) * hr_dim1] + xr * hr[i
4005	+ j * hr_dim1] - xi * hi[i + j * hi_dim1];
4006	hi[i + (j - 1) * hi_dim1] = hi[i + (j - 1) * hi_dim1] + xr * hi[i
4007	+ j * hi_dim1] + xi * hr[i + j * hr_dim1];
4008	/* L600: */
4009	}
4010
4011	/* L640: */
4012	}
4013
4014	goto L240;
4015	/* .......... A ROOT FOUND .......... */
4016	L660:
4017	wr[en] = hr[en + en * hr_dim1] + tr;
4018	wi[en] = hi[en + en * hi_dim1] + ti;
4019	en = enm1;
4020	goto L220;
4021	/* .......... SET ERROR -- ALL EIGENVALUES HAVE NOT */
4022	/* CONVERGED AFTER 30N ITERATIONS .......... /
4023	L1000:
4024	*ierr = en;
4025	L1001:
4026	return 0;
4027	} /* comlr_ */
4028
4029	/* Subroutine / int comlr2_(integer nm, integer n, integer low, integer *
4030	igh, integer int_, doublereal hr, doublereal hi, doublereal wr,
4031	doublereal wi, doublereal zr, doublereal zi, integer ierr)
4032	{
4033	/* System generated locals */
4034	integer hr_dim1, hr_offset, hi_dim1, hi_offset, zr_dim1, zr_offset,
4035	zi_dim1, zi_offset, i_1, i_2, i_3;
4036	doublereal d_1, d_2, d_3, d_4;
4037
4038	/* Local variables */
4039	static integer iend;
4040	extern /* Subroutine / int cdiv_(doublereal , doublereal , doublereal
4041	, doublereal , doublereal , doublereal *);
4042	static doublereal norm;
4043	static integer i, j, k, l, m, ii, en, jj, ll, mm, nn;
4044	static doublereal si, ti, xi, yi, sr, tr, xr, yr;
4045	static integer im1;
4046	extern /* Subroutine / int csroot_(doublereal , doublereal *,
4047	doublereal , doublereal );
4048	static integer ip1, mp1, itn, its;
4049	static doublereal zzi, zzr;
4050	static integer enm1;
4051	static doublereal tst1, tst2;
4052
4053
4054
4055	/* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE COMLR2, */
4056	/* NUM. MATH. 16, 181-204(1970) BY PETERS AND WILKINSON. */
4057	/* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971). */
4058
4059	/* THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS */
4060	/* OF A COMPLEX UPPER HESSENBERG MATRIX BY THE MODIFIED LR */
4061	/* METHOD. THE EIGENVECTORS OF A COMPLEX GENERAL MATRIX */
4062	/* CAN ALSO BE FOUND IF COMHES HAS BEEN USED TO REDUCE */
4063	/* THIS GENERAL MATRIX TO HESSENBERG FORM. */
4064
4065	/* ON INPUT */
4066
4067	/* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
4068	/* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
4069	/* DIMENSION STATEMENT. */
4070
4071	/* N IS THE ORDER OF THE MATRIX. */
4072
4073	/* LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
4074	/* SUBROUTINE CBAL. IF CBAL HAS NOT BEEN USED, */
4075	/* SET LOW=1, IGH=N. */
4076
4077	/* INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS INTERCHANGED */
4078	/* IN THE REDUCTION BY COMHES, IF PERFORMED. ONLY ELEMENTS */
4079	/* LOW THROUGH IGH ARE USED. IF THE EIGENVECTORS OF THE HESSEN-
4080	*/
4081	/* BERG MATRIX ARE DESIRED, SET INT(J)=J FOR THESE ELEMENTS. */
4082
4083	/* HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS, */
4084	/* RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX. */
4085	/* THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN THE */
4086	/* MULTIPLIERS WHICH WERE USED IN THE REDUCTION BY COMHES, */
4087	/* IF PERFORMED. IF THE EIGENVECTORS OF THE HESSENBERG */
4088	/* MATRIX ARE DESIRED, THESE ELEMENTS MUST BE SET TO ZERO. */
4089
4090	/* ON OUTPUT */
4091
4092	/* THE UPPER HESSENBERG PORTIONS OF HR AND HI HAVE BEEN */
4093	/* DESTROYED, BUT THE LOCATION HR(1,1) CONTAINS THE NORM */
4094	/* OF THE TRIANGULARIZED MATRIX. */
4095
4096	/* WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, */
4097	/* RESPECTIVELY, OF THE EIGENVALUES. IF AN ERROR */
4098	/* EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT */
4099	/* FOR INDICES IERR+1,...,N. */
4100
4101	/* ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, */
4102	/* RESPECTIVELY, OF THE EIGENVECTORS. THE EIGENVECTORS */
4103	/* ARE UNNORMALIZED. IF AN ERROR EXIT IS MADE, NONE OF */
4104	/* THE EIGENVECTORS HAS BEEN FOUND. */
4105
4106	/* IERR IS SET TO */
4107	/* ZERO FOR NORMAL RETURN, */
4108	/* J IF THE LIMIT OF 30N ITERATIONS IS EXHAUSTED /
4109	/* WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. */
4110
4111
4112	/* CALLS CDIV FOR COMPLEX DIVISION. */
4113	/* CALLS CSROOT FOR COMPLEX SQUARE ROOT. */
4114
4115	/* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
4116	/* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
4117	*/
4118
4119	/* THIS VERSION DATED AUGUST 1983. */
4120
4121	/* ------------------------------------------------------------------
4122	*/
4123
4124	/* Parameter adjustments */
4125	zi_dim1 = *nm;
4126	zi_offset = zi_dim1 + 1;
4127	zi -= zi_offset;
4128	zr_dim1 = *nm;
4129	zr_offset = zr_dim1 + 1;
4130	zr -= zr_offset;
4131	--wi;
4132	--wr;
4133	hi_dim1 = *nm;
4134	hi_offset = hi_dim1 + 1;
4135	hi -= hi_offset;
4136	hr_dim1 = *nm;
4137	hr_offset = hr_dim1 + 1;
4138	hr -= hr_offset;
4139	--int_;
4140
4141	/* Function Body */
4142	*ierr = 0;
4143	/* .......... INITIALIZE EIGENVECTOR MATRIX .......... */
4144	i_1 = *n;
4145	for (i = 1; i <= i_1; ++i) {
4146
4147	i_2 = *n;
4148	for (j = 1; j <= i_2; ++j) {
4149	zr[i + j * zr_dim1] = 0.;
4150	zi[i + j * zi_dim1] = 0.;
4151	if (i == j) {
4152	zr[i + j * zr_dim1] = 1.;
4153	}
4154	/* L100: */
4155	}
4156	}
4157	/* .......... FORM THE MATRIX OF ACCUMULATED TRANSFORMATIONS */
4158	/* FROM THE INFORMATION LEFT BY COMHES .......... */
4159	iend = igh - low - 1;
4160	if (iend <= 0) {
4161	goto L180;
4162	}
4163	/* .......... FOR I=IGH-1 STEP -1 UNTIL LOW+1 DO -- .......... */
4164	i_2 = iend;
4165	for (ii = 1; ii <= i_2; ++ii) {
4166	i = *igh - ii;
4167	ip1 = i + 1;
4168
4169	i_1 = *igh;
4170	for (k = ip1; k <= i_1; ++k) {
4171	zr[k + i * zr_dim1] = hr[k + (i - 1) * hr_dim1];
4172	zi[k + i * zi_dim1] = hi[k + (i - 1) * hi_dim1];
4173	/* L120: */
4174	}
4175
4176	j = int_[i];
4177	if (i == j) {
4178	goto L160;
4179	}
4180
4181	i_1 = *igh;
4182	for (k = i; k <= i_1; ++k) {
4183	zr[i + k * zr_dim1] = zr[j + k * zr_dim1];
4184	zi[i + k * zi_dim1] = zi[j + k * zi_dim1];
4185	zr[j + k * zr_dim1] = 0.;
4186	zi[j + k * zi_dim1] = 0.;
4187	/* L140: */
4188	}
4189
4190	zr[j + i * zr_dim1] = 1.;
4191	L160:
4192	;
4193	}
4194	/* .......... STORE ROOTS ISOLATED BY CBAL .......... */
4195	L180:
4196	i_2 = *n;
4197	for (i = 1; i <= i_2; ++i) {
4198	if (i >= low && i <= igh) {
4199	goto L200;
4200	}
4201	wr[i] = hr[i + i * hr_dim1];
4202	wi[i] = hi[i + i * hi_dim1];
4203	L200:
4204	;
4205	}
4206
4207	en = *igh;
4208	tr = 0.;
4209	ti = 0.;
4210	itn = n 30;
4211	/* .......... SEARCH FOR NEXT EIGENVALUE .......... */
4212	L220:
4213	if (en < *low) {
4214	goto L680;
4215	}
4216	its = 0;
4217	enm1 = en - 1;
4218	/* .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT */
4219	/* FOR L=EN STEP -1 UNTIL LOW DO -- .......... */
4220	L240:
4221	i_2 = en;
4222	for (ll = *low; ll <= i_2; ++ll) {
4223	l = en + *low - ll;
4224	if (l == *low) {
4225	goto L300;
4226	}
4227	tst1 = (d_1 = hr[l - 1 + (l - 1) * hr_dim1], abs(d_1)) + (d_2 = hi[
4228	l - 1 + (l - 1) * hi_dim1], abs(d_2)) + (d_3 = hr[l + l *
4229	hr_dim1], abs(d_3)) + (d_4 = hi[l + l * hi_dim1], abs(d_4))
4230	;
4231	tst2 = tst1 + (d_1 = hr[l + (l - 1) * hr_dim1], abs(d_1)) + (d_2 =
4232	hi[l + (l - 1) * hi_dim1], abs(d_2));
4233	if (tst2 == tst1) {
4234	goto L300;
4235	}
4236	/* L260: */
4237	}
4238	/* .......... FORM SHIFT .......... */
4239	L300:
4240	if (l == en) {
4241	goto L660;
4242	}
4243	if (itn == 0) {
4244	goto L1000;
4245	}
4246	if (its == 10 \|\| its == 20) {
4247	goto L320;
4248	}
4249	sr = hr[en + en * hr_dim1];
4250	si = hi[en + en * hi_dim1];
4251	xr = hr[enm1 + en * hr_dim1] * hr[en + enm1 * hr_dim1] - hi[enm1 + en *
4252	hi_dim1] * hi[en + enm1 * hi_dim1];
4253	xi = hr[enm1 + en * hr_dim1] * hi[en + enm1 * hi_dim1] + hi[enm1 + en *
4254	hi_dim1] * hr[en + enm1 * hr_dim1];
4255	if (xr == 0. && xi == 0.) {
4256	goto L340;
4257	}
4258	yr = (hr[enm1 + enm1 * hr_dim1] - sr) / 2.;
4259	yi = (hi[enm1 + enm1 * hi_dim1] - si) / 2.;
4260	/* Computing 2nd power */
4261	d_2 = yr;
4262	/* Computing 2nd power */
4263	d_3 = yi;
4264	d_1 = d_2 * d_2 - d_3 * d_3 + xr;
4265	d_4 = yr * 2. * yi + xi;
4266	csroot_(&d_1, &d_4, &zzr, &zzi);
4267	if (yr * zzr + yi * zzi >= 0.) {
4268	goto L310;
4269	}
4270	zzr = -zzr;
4271	zzi = -zzi;
4272	L310:
4273	d_1 = yr + zzr;
4274	d_2 = yi + zzi;
4275	cdiv_(&xr, &xi, &d_1, &d_2, &xr, &xi);
4276	sr -= xr;
4277	si -= xi;
4278	goto L340;
4279	/* .......... FORM EXCEPTIONAL SHIFT .......... */
4280	L320:
4281	sr = (d_1 = hr[en + enm1 * hr_dim1], abs(d_1)) + (d_2 = hr[enm1 + (en
4282	- 2) * hr_dim1], abs(d_2));
4283	si = (d_1 = hi[en + enm1 * hi_dim1], abs(d_1)) + (d_2 = hi[enm1 + (en
4284	- 2) * hi_dim1], abs(d_2));
4285
4286	L340:
4287	i_2 = en;
4288	for (i = *low; i <= i_2; ++i) {
4289	hr[i + i * hr_dim1] -= sr;
4290	hi[i + i * hi_dim1] -= si;
4291	/* L360: */
4292	}
4293
4294	tr += sr;
4295	ti += si;
4296	++its;
4297	--itn;
4298	/* .......... LOOK FOR TWO CONSECUTIVE SMALL */
4299	/* SUB-DIAGONAL ELEMENTS .......... */
4300	xr = (d_1 = hr[enm1 + enm1 * hr_dim1], abs(d_1)) + (d_2 = hi[enm1 +
4301	enm1 * hi_dim1], abs(d_2));
4302	yr = (d_1 = hr[en + enm1 * hr_dim1], abs(d_1)) + (d_2 = hi[en + enm1 *
4303	hi_dim1], abs(d_2));
4304	zzr = (d_1 = hr[en + en * hr_dim1], abs(d_1)) + (d_2 = hi[en + en *
4305	hi_dim1], abs(d_2));
4306	/* .......... FOR M=EN-1 STEP -1 UNTIL L DO -- .......... */
4307	i_2 = enm1;
4308	for (mm = l; mm <= i_2; ++mm) {
4309	m = enm1 + l - mm;
4310	if (m == l) {
4311	goto L420;
4312	}
4313	yi = yr;
4314	yr = (d_1 = hr[m + (m - 1) * hr_dim1], abs(d_1)) + (d_2 = hi[m + (
4315	m - 1) * hi_dim1], abs(d_2));
4316	xi = zzr;
4317	zzr = xr;
4318	xr = (d_1 = hr[m - 1 + (m - 1) * hr_dim1], abs(d_1)) + (d_2 = hi[m
4319	- 1 + (m - 1) * hi_dim1], abs(d_2));
4320	tst1 = zzr / yi * (zzr + xr + xi);
4321	tst2 = tst1 + yr;
4322	if (tst2 == tst1) {
4323	goto L420;
4324	}
4325	/* L380: */
4326	}
4327	/* .......... TRIANGULAR DECOMPOSITION H=LR .......... /
4328	L420:
4329	mp1 = m + 1;
4330
4331	i_2 = en;
4332	for (i = mp1; i <= i_2; ++i) {
4333	im1 = i - 1;
4334	xr = hr[im1 + im1 * hr_dim1];
4335	xi = hi[im1 + im1 * hi_dim1];
4336	yr = hr[i + im1 * hr_dim1];
4337	yi = hi[i + im1 * hi_dim1];
4338	if (abs(xr) + abs(xi) >= abs(yr) + abs(yi)) {
4339	goto L460;
4340	}
4341	/* .......... INTERCHANGE ROWS OF HR AND HI .......... */
4342	i_1 = *n;
4343	for (j = im1; j <= i_1; ++j) {
4344	zzr = hr[im1 + j * hr_dim1];
4345	hr[im1 + j * hr_dim1] = hr[i + j * hr_dim1];
4346	hr[i + j * hr_dim1] = zzr;
4347	zzi = hi[im1 + j * hi_dim1];
4348	hi[im1 + j * hi_dim1] = hi[i + j * hi_dim1];
4349	hi[i + j * hi_dim1] = zzi;
4350	/* L440: */
4351	}
4352
4353	cdiv_(&xr, &xi, &yr, &yi, &zzr, &zzi);
4354	wr[i] = 1.;
4355	goto L480;
4356	L460:
4357	cdiv_(&yr, &yi, &xr, &xi, &zzr, &zzi);
4358	wr[i] = -1.;
4359	L480:
4360	hr[i + im1 * hr_dim1] = zzr;
4361	hi[i + im1 * hi_dim1] = zzi;
4362
4363	i_1 = *n;
4364	for (j = i; j <= i_1; ++j) {
4365	hr[i + j * hr_dim1] = hr[i + j * hr_dim1] - zzr * hr[im1 + j *
4366	hr_dim1] + zzi * hi[im1 + j * hi_dim1];
4367	hi[i + j * hi_dim1] = hi[i + j * hi_dim1] - zzr * hi[im1 + j *
4368	hi_dim1] - zzi * hr[im1 + j * hr_dim1];
4369	/* L500: */
4370	}
4371
4372	/* L520: */
4373	}
4374	/* .......... COMPOSITION RL=H .......... /
4375	i_2 = en;
4376	for (j = mp1; j <= i_2; ++j) {
4377	xr = hr[j + (j - 1) * hr_dim1];
4378	xi = hi[j + (j - 1) * hi_dim1];
4379	hr[j + (j - 1) * hr_dim1] = 0.;
4380	hi[j + (j - 1) * hi_dim1] = 0.;
4381	/* .......... INTERCHANGE COLUMNS OF HR, HI, ZR, AND ZI, */
4382	/* IF NECESSARY .......... */
4383	if (wr[j] <= 0.) {
4384	goto L580;
4385	}
4386
4387	i_1 = j;
4388	for (i = 1; i <= i_1; ++i) {
4389	zzr = hr[i + (j - 1) * hr_dim1];
4390	hr[i + (j - 1) * hr_dim1] = hr[i + j * hr_dim1];
4391	hr[i + j * hr_dim1] = zzr;
4392	zzi = hi[i + (j - 1) * hi_dim1];
4393	hi[i + (j - 1) * hi_dim1] = hi[i + j * hi_dim1];
4394	hi[i + j * hi_dim1] = zzi;
4395	/* L540: */
4396	}
4397
4398	i_1 = *igh;
4399	for (i = *low; i <= i_1; ++i) {
4400	zzr = zr[i + (j - 1) * zr_dim1];
4401	zr[i + (j - 1) * zr_dim1] = zr[i + j * zr_dim1];
4402	zr[i + j * zr_dim1] = zzr;
4403	zzi = zi[i + (j - 1) * zi_dim1];
4404	zi[i + (j - 1) * zi_dim1] = zi[i + j * zi_dim1];
4405	zi[i + j * zi_dim1] = zzi;
4406	/* L560: */
4407	}
4408
4409	L580:
4410	i_1 = j;
4411	for (i = 1; i <= i_1; ++i) {
4412	hr[i + (j - 1) * hr_dim1] = hr[i + (j - 1) * hr_dim1] + xr * hr[i
4413	+ j * hr_dim1] - xi * hi[i + j * hi_dim1];
4414	hi[i + (j - 1) * hi_dim1] = hi[i + (j - 1) * hi_dim1] + xr * hi[i
4415	+ j * hi_dim1] + xi * hr[i + j * hr_dim1];
4416	/* L600: */
4417	}
4418	/* .......... ACCUMULATE TRANSFORMATIONS .......... */
4419	i_1 = *igh;
4420	for (i = *low; i <= i_1; ++i) {
4421	zr[i + (j - 1) * zr_dim1] = zr[i + (j - 1) * zr_dim1] + xr * zr[i
4422	+ j * zr_dim1] - xi * zi[i + j * zi_dim1];
4423	zi[i + (j - 1) * zi_dim1] = zi[i + (j - 1) * zi_dim1] + xr * zi[i
4424	+ j * zi_dim1] + xi * zr[i + j * zr_dim1];
4425	/* L620: */
4426	}
4427
4428	/* L640: */
4429	}
4430
4431	goto L240;
4432	/* .......... A ROOT FOUND .......... */
4433	L660:
4434	hr[en + en * hr_dim1] += tr;
4435	wr[en] = hr[en + en * hr_dim1];
4436	hi[en + en * hi_dim1] += ti;
4437	wi[en] = hi[en + en * hi_dim1];
4438	en = enm1;
4439	goto L220;
4440	/* .......... ALL ROOTS FOUND. BACKSUBSTITUTE TO FIND */
4441	/* VECTORS OF UPPER TRIANGULAR FORM .......... */
4442	L680:
4443	norm = 0.;
4444
4445	i_2 = *n;
4446	for (i = 1; i <= i_2; ++i) {
4447
4448	i_1 = *n;
4449	for (j = i; j <= i_1; ++j) {
4450	tr = (d_1 = hr[i + j * hr_dim1], abs(d_1)) + (d_2 = hi[i + j *
4451	hi_dim1], abs(d_2));
4452	if (tr > norm) {
4453	norm = tr;
4454	}
4455	/* L720: */
4456	}
4457	}
4458
4459	hr[hr_dim1 + 1] = norm;
4460	if (*n == 1 \|\| norm == 0.) {
4461	goto L1001;
4462	}
4463	/* .......... FOR EN=N STEP -1 UNTIL 2 DO -- .......... */
4464	i_1 = *n;
4465	for (nn = 2; nn <= i_1; ++nn) {
4466	en = *n + 2 - nn;
4467	xr = wr[en];
4468	xi = wi[en];
4469	hr[en + en * hr_dim1] = 1.;
4470	hi[en + en * hi_dim1] = 0.;
4471	enm1 = en - 1;
4472	/* .......... FOR I=EN-1 STEP -1 UNTIL 1 DO -- .......... */
4473	i_2 = enm1;
4474	for (ii = 1; ii <= i_2; ++ii) {
4475	i = en - ii;
4476	zzr = 0.;
4477	zzi = 0.;
4478	ip1 = i + 1;
4479
4480	i_3 = en;
4481	for (j = ip1; j <= i_3; ++j) {
4482	zzr = zzr + hr[i + j * hr_dim1] * hr[j + en * hr_dim1] - hi[i
4483	+ j * hi_dim1] * hi[j + en * hi_dim1];
4484	zzi = zzi + hr[i + j * hr_dim1] * hi[j + en * hi_dim1] + hi[i
4485	+ j * hi_dim1] * hr[j + en * hr_dim1];
4486	/* L740: */
4487	}
4488
4489	yr = xr - wr[i];
4490	yi = xi - wi[i];
4491	if (yr != 0. \|\| yi != 0.) {
4492	goto L765;
4493	}
4494	tst1 = norm;
4495	yr = tst1;
4496	L760:
4497	yr *= .01;
4498	tst2 = norm + yr;
4499	if (tst2 > tst1) {
4500	goto L760;
4501	}
4502	L765:
4503	cdiv_(&zzr, &zzi, &yr, &yi, &hr[i + en * hr_dim1], &hi[i + en *
4504	hi_dim1]);
4505	/* .......... OVERFLOW CONTROL .......... */
4506	tr = (d_1 = hr[i + en * hr_dim1], abs(d_1)) + (d_2 = hi[i + en
4507	* hi_dim1], abs(d_2));
4508	if (tr == 0.) {
4509	goto L780;
4510	}
4511	tst1 = tr;
4512	tst2 = tst1 + 1. / tst1;
4513	if (tst2 > tst1) {
4514	goto L780;
4515	}
4516	i_3 = en;
4517	for (j = i; j <= i_3; ++j) {
4518	hr[j + en * hr_dim1] /= tr;
4519	hi[j + en * hi_dim1] /= tr;
4520	/* L770: */
4521	}
4522
4523	L780:
4524	;
4525	}
4526
4527	/* L800: */
4528	}
4529	/* .......... END BACKSUBSTITUTION .......... */
4530	enm1 = *n - 1;
4531	/* .......... VECTORS OF ISOLATED ROOTS .......... */
4532	i_1 = enm1;
4533	for (i = 1; i <= i_1; ++i) {
4534	if (i >= low && i <= igh) {
4535	goto L840;
4536	}
4537	ip1 = i + 1;
4538
4539	i_2 = *n;
4540	for (j = ip1; j <= i_2; ++j) {
4541	zr[i + j * zr_dim1] = hr[i + j * hr_dim1];
4542	zi[i + j * zi_dim1] = hi[i + j * hi_dim1];
4543	/* L820: */
4544	}
4545
4546	L840:
4547	;
4548	}
4549	/* .......... MULTIPLY BY TRANSFORMATION MATRIX TO GIVE */
4550	/* VECTORS OF ORIGINAL FULL MATRIX. */
4551	/* FOR J=N STEP -1 UNTIL LOW+1 DO -- .......... */
4552	i_1 = enm1;
4553	for (jj = *low; jj <= i_1; ++jj) {
4554	j = n + low - jj;
4555	m = min(j,*igh);
4556
4557	i_2 = *igh;
4558	for (i = *low; i <= i_2; ++i) {
4559	zzr = 0.;
4560	zzi = 0.;
4561
4562	i_3 = m;
4563	for (k = *low; k <= i_3; ++k) {
4564	zzr = zzr + zr[i + k * zr_dim1] * hr[k + j * hr_dim1] - zi[i
4565	+ k * zi_dim1] * hi[k + j * hi_dim1];
4566	zzi = zzi + zr[i + k * zr_dim1] * hi[k + j * hi_dim1] + zi[i
4567	+ k * zi_dim1] * hr[k + j * hr_dim1];
4568	/* L860: */
4569	}
4570
4571	zr[i + j * zr_dim1] = zzr;
4572	zi[i + j * zi_dim1] = zzi;
4573	/* L880: */
4574	}
4575	}
4576
4577	goto L1001;
4578	/* .......... SET ERROR -- ALL EIGENVALUES HAVE NOT */
4579	/* CONVERGED AFTER 30N ITERATIONS .......... /
4580	L1000:
4581	*ierr = en;
4582	L1001:
4583	return 0;
4584	} /* comlr2_ */
4585
4586	/* Subroutine / int comqr_(integer nm, integer n, integer low, integer *
4587	igh, doublereal hr, doublereal hi, doublereal wr, doublereal wi,
4588	integer *ierr)
4589	{
4590	/* System generated locals */
4591	integer hr_dim1, hr_offset, hi_dim1, hi_offset, i_1, i_2;
4592	doublereal d_1, d_2, d_3, d_4;
4593
4594	/* Local variables */
4595	extern /* Subroutine / int cdiv_(doublereal , doublereal , doublereal
4596	, doublereal , doublereal , doublereal *);
4597	static doublereal norm;
4598	static integer i, j, l, en, ll;
4599	static doublereal si, ti, xi, yi, sr, tr, xr, yr;
4600	extern doublereal pythag_(doublereal , doublereal );
4601	extern /* Subroutine / int csroot_(doublereal , doublereal *,
4602	doublereal , doublereal );
4603	static integer lp1, itn, its;
4604	static doublereal zzi, zzr;
4605	static integer enm1;
4606	static doublereal tst1, tst2;
4607
4608
4609
4610	/* THIS SUBROUTINE IS A TRANSLATION OF A UNITARY ANALOGUE OF THE */
4611	/* ALGOL PROCEDURE COMLR, NUM. MATH. 12, 369-376(1968) BY MARTIN */
4612	/* AND WILKINSON. */
4613	/* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 396-403(1971). */
4614	/* THE UNITARY ANALOGUE SUBSTITUTES THE QR ALGORITHM OF FRANCIS */
4615	/* (COMP. JOUR. 4, 332-345(1962)) FOR THE LR ALGORITHM. */
4616
4617	/* THIS SUBROUTINE FINDS THE EIGENVALUES OF A COMPLEX */
4618	/* UPPER HESSENBERG MATRIX BY THE QR METHOD. */
4619
4620	/* ON INPUT */
4621
4622	/* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
4623	/* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
4624	/* DIMENSION STATEMENT. */
4625
4626	/* N IS THE ORDER OF THE MATRIX. */
4627
4628	/* LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
4629	/* SUBROUTINE CBAL. IF CBAL HAS NOT BEEN USED, */
4630	/* SET LOW=1, IGH=N. */
4631
4632	/* HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS, */
4633	/* RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX. */
4634	/* THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN */
4635	/* INFORMATION ABOUT THE UNITARY TRANSFORMATIONS USED IN */
4636	/* THE REDUCTION BY CORTH, IF PERFORMED. */
4637
4638	/* ON OUTPUT */
4639
4640	/* THE UPPER HESSENBERG PORTIONS OF HR AND HI HAVE BEEN */
4641	/* DESTROYED. THEREFORE, THEY MUST BE SAVED BEFORE */
4642	/* CALLING COMQR IF SUBSEQUENT CALCULATION OF */
4643	/* EIGENVECTORS IS TO BE PERFORMED. */
4644
4645	/* WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, */
4646	/* RESPECTIVELY, OF THE EIGENVALUES. IF AN ERROR */
4647	/* EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT */
4648	/* FOR INDICES IERR+1,...,N. */
4649
4650	/* IERR IS SET TO */
4651	/* ZERO FOR NORMAL RETURN, */
4652	/* J IF THE LIMIT OF 30N ITERATIONS IS EXHAUSTED /
4653	/* WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. */
4654
4655	/* CALLS CDIV FOR COMPLEX DIVISION. */
4656	/* CALLS CSROOT FOR COMPLEX SQUARE ROOT. */
4657	/* CALLS PYTHAG FOR DSQRT(AA + BB) . */
4658
4659	/* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
4660	/* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
4661	*/
4662
4663	/* THIS VERSION DATED AUGUST 1983. */
4664
4665	/* ------------------------------------------------------------------
4666	*/
4667
4668	/* Parameter adjustments */
4669	--wi;
4670	--wr;
4671	hi_dim1 = *nm;
4672	hi_offset = hi_dim1 + 1;
4673	hi -= hi_offset;
4674	hr_dim1 = *nm;
4675	hr_offset = hr_dim1 + 1;
4676	hr -= hr_offset;
4677
4678	/* Function Body */
4679	*ierr = 0;
4680	if (low == igh) {
4681	goto L180;
4682	}
4683	/* .......... CREATE REAL SUBDIAGONAL ELEMENTS .......... */
4684	l = *low + 1;
4685
4686	i_1 = *igh;
4687	for (i = l; i <= i_1; ++i) {
4688	/* Computing MIN */
4689	i_2 = i + 1;
4690	ll = min(i_2,*igh);
4691	if (hi[i + (i - 1) * hi_dim1] == 0.) {
4692	goto L170;
4693	}
4694	norm = pythag_(&hr[i + (i - 1) * hr_dim1], &hi[i + (i - 1) * hi_dim1])
4695	;
4696	yr = hr[i + (i - 1) * hr_dim1] / norm;
4697	yi = hi[i + (i - 1) * hi_dim1] / norm;
4698	hr[i + (i - 1) * hr_dim1] = norm;
4699	hi[i + (i - 1) * hi_dim1] = 0.;
4700
4701	i_2 = *igh;
4702	for (j = i; j <= i_2; ++j) {
4703	si = yr * hi[i + j * hi_dim1] - yi * hr[i + j * hr_dim1];
4704	hr[i + j * hr_dim1] = yr * hr[i + j * hr_dim1] + yi * hi[i + j *
4705	hi_dim1];
4706	hi[i + j * hi_dim1] = si;
4707	/* L155: */
4708	}
4709
4710	i_2 = ll;
4711	for (j = *low; j <= i_2; ++j) {
4712	si = yr * hi[j + i * hi_dim1] + yi * hr[j + i * hr_dim1];
4713	hr[j + i * hr_dim1] = yr * hr[j + i * hr_dim1] - yi * hi[j + i *
4714	hi_dim1];
4715	hi[j + i * hi_dim1] = si;
4716	/* L160: */
4717	}
4718
4719	L170:
4720	;
4721	}
4722	/* .......... STORE ROOTS ISOLATED BY CBAL .......... */
4723	L180:
4724	i_1 = *n;
4725	for (i = 1; i <= i_1; ++i) {
4726	if (i >= low && i <= igh) {
4727	goto L200;
4728	}
4729	wr[i] = hr[i + i * hr_dim1];
4730	wi[i] = hi[i + i * hi_dim1];
4731	L200:
4732	;
4733	}
4734
4735	en = *igh;
4736	tr = 0.;
4737	ti = 0.;
4738	itn = n 30;
4739	/* .......... SEARCH FOR NEXT EIGENVALUE .......... */
4740	L220:
4741	if (en < *low) {
4742	goto L1001;
4743	}
4744	its = 0;
4745	enm1 = en - 1;
4746	/* .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT */
4747	/* FOR L=EN STEP -1 UNTIL LOW D0 -- .......... */
4748	L240:
4749	i_1 = en;
4750	for (ll = *low; ll <= i_1; ++ll) {
4751	l = en + *low - ll;
4752	if (l == *low) {
4753	goto L300;
4754	}
4755	tst1 = (d_1 = hr[l - 1 + (l - 1) * hr_dim1], abs(d_1)) + (d_2 = hi[
4756	l - 1 + (l - 1) * hi_dim1], abs(d_2)) + (d_3 = hr[l + l *
4757	hr_dim1], abs(d_3)) + (d_4 = hi[l + l * hi_dim1], abs(d_4))
4758	;
4759	tst2 = tst1 + (d_1 = hr[l + (l - 1) * hr_dim1], abs(d_1));
4760	if (tst2 == tst1) {
4761	goto L300;
4762	}
4763	/* L260: */
4764	}
4765	/* .......... FORM SHIFT .......... */
4766	L300:
4767	if (l == en) {
4768	goto L660;
4769	}
4770	if (itn == 0) {
4771	goto L1000;
4772	}
4773	if (its == 10 \|\| its == 20) {
4774	goto L320;
4775	}
4776	sr = hr[en + en * hr_dim1];
4777	si = hi[en + en * hi_dim1];
4778	xr = hr[enm1 + en * hr_dim1] * hr[en + enm1 * hr_dim1];
4779	xi = hi[enm1 + en * hi_dim1] * hr[en + enm1 * hr_dim1];
4780	if (xr == 0. && xi == 0.) {
4781	goto L340;
4782	}
4783	yr = (hr[enm1 + enm1 * hr_dim1] - sr) / 2.;
4784	yi = (hi[enm1 + enm1 * hi_dim1] - si) / 2.;
4785	/* Computing 2nd power */
4786	d_2 = yr;
4787	/* Computing 2nd power */
4788	d_3 = yi;
4789	d_1 = d_2 * d_2 - d_3 * d_3 + xr;
4790	d_4 = yr * 2. * yi + xi;
4791	csroot_(&d_1, &d_4, &zzr, &zzi);
4792	if (yr * zzr + yi * zzi >= 0.) {
4793	goto L310;
4794	}
4795	zzr = -zzr;
4796	zzi = -zzi;
4797	L310:
4798	d_1 = yr + zzr;
4799	d_2 = yi + zzi;
4800	cdiv_(&xr, &xi, &d_1, &d_2, &xr, &xi);
4801	sr -= xr;
4802	si -= xi;
4803	goto L340;
4804	/* .......... FORM EXCEPTIONAL SHIFT .......... */
4805	L320:
4806	sr = (d_1 = hr[en + enm1 * hr_dim1], abs(d_1)) + (d_2 = hr[enm1 + (en
4807	- 2) * hr_dim1], abs(d_2));
4808	si = 0.;
4809
4810	L340:
4811	i_1 = en;
4812	for (i = *low; i <= i_1; ++i) {
4813	hr[i + i * hr_dim1] -= sr;
4814	hi[i + i * hi_dim1] -= si;
4815	/* L360: */
4816	}
4817
4818	tr += sr;
4819	ti += si;
4820	++its;
4821	--itn;
4822	/* .......... REDUCE TO TRIANGLE (ROWS) .......... */
4823	lp1 = l + 1;
4824
4825	i_1 = en;
4826	for (i = lp1; i <= i_1; ++i) {
4827	sr = hr[i + (i - 1) * hr_dim1];
4828	hr[i + (i - 1) * hr_dim1] = 0.;
4829	d_1 = pythag_(&hr[i - 1 + (i - 1) * hr_dim1], &hi[i - 1 + (i - 1) *
4830	hi_dim1]);
4831	norm = pythag_(&d_1, &sr);
4832	xr = hr[i - 1 + (i - 1) * hr_dim1] / norm;
4833	wr[i - 1] = xr;
4834	xi = hi[i - 1 + (i - 1) * hi_dim1] / norm;
4835	wi[i - 1] = xi;
4836	hr[i - 1 + (i - 1) * hr_dim1] = norm;
4837	hi[i - 1 + (i - 1) * hi_dim1] = 0.;
4838	hi[i + (i - 1) * hi_dim1] = sr / norm;
4839
4840	i_2 = en;
4841	for (j = i; j <= i_2; ++j) {
4842	yr = hr[i - 1 + j * hr_dim1];
4843	yi = hi[i - 1 + j * hi_dim1];
4844	zzr = hr[i + j * hr_dim1];
4845	zzi = hi[i + j * hi_dim1];
4846	hr[i - 1 + j * hr_dim1] = xr * yr + xi * yi + hi[i + (i - 1) *
4847	hi_dim1] * zzr;
4848	hi[i - 1 + j * hi_dim1] = xr * yi - xi * yr + hi[i + (i - 1) *
4849	hi_dim1] * zzi;
4850	hr[i + j * hr_dim1] = xr * zzr - xi * zzi - hi[i + (i - 1) *
4851	hi_dim1] * yr;
4852	hi[i + j * hi_dim1] = xr * zzi + xi * zzr - hi[i + (i - 1) *
4853	hi_dim1] * yi;
4854	/* L490: */
4855	}
4856
4857	/* L500: */
4858	}
4859
4860	si = hi[en + en * hi_dim1];
4861	if (si == 0.) {
4862	goto L540;
4863	}
4864	norm = pythag_(&hr[en + en * hr_dim1], &si);
4865	sr = hr[en + en * hr_dim1] / norm;
4866	si /= norm;
4867	hr[en + en * hr_dim1] = norm;
4868	hi[en + en * hi_dim1] = 0.;
4869	/* .......... INVERSE OPERATION (COLUMNS) .......... */
4870	L540:
4871	i_1 = en;
4872	for (j = lp1; j <= i_1; ++j) {
4873	xr = wr[j - 1];
4874	xi = wi[j - 1];
4875
4876	i_2 = j;
4877	for (i = l; i <= i_2; ++i) {
4878	yr = hr[i + (j - 1) * hr_dim1];
4879	yi = 0.;
4880	zzr = hr[i + j * hr_dim1];
4881	zzi = hi[i + j * hi_dim1];
4882	if (i == j) {
4883	goto L560;
4884	}
4885	yi = hi[i + (j - 1) * hi_dim1];
4886	hi[i + (j - 1) * hi_dim1] = xr * yi + xi * yr + hi[j + (j - 1) *
4887	hi_dim1] * zzi;
4888	L560:
4889	hr[i + (j - 1) * hr_dim1] = xr * yr - xi * yi + hi[j + (j - 1) *
4890	hi_dim1] * zzr;
4891	hr[i + j * hr_dim1] = xr * zzr + xi * zzi - hi[j + (j - 1) *
4892	hi_dim1] * yr;
4893	hi[i + j * hi_dim1] = xr * zzi - xi * zzr - hi[j + (j - 1) *
4894	hi_dim1] * yi;
4895	/* L580: */
4896	}
4897
4898	/* L600: */
4899	}
4900
4901	if (si == 0.) {
4902	goto L240;
4903	}
4904
4905	i_1 = en;
4906	for (i = l; i <= i_1; ++i) {
4907	yr = hr[i + en * hr_dim1];
4908	yi = hi[i + en * hi_dim1];
4909	hr[i + en * hr_dim1] = sr * yr - si * yi;
4910	hi[i + en * hi_dim1] = sr * yi + si * yr;
4911	/* L630: */
4912	}
4913
4914	goto L240;
4915	/* .......... A ROOT FOUND .......... */
4916	L660:
4917	wr[en] = hr[en + en * hr_dim1] + tr;
4918	wi[en] = hi[en + en * hi_dim1] + ti;
4919	en = enm1;
4920	goto L220;
4921	/* .......... SET ERROR -- ALL EIGENVALUES HAVE NOT */
4922	/* CONVERGED AFTER 30N ITERATIONS .......... /
4923	L1000:
4924	*ierr = en;
4925	L1001:
4926	return 0;
4927	} /* comqr_ */
4928
4929	/* Subroutine / int comqr2_(integer nm, integer n, integer low, integer *
4930	igh, doublereal ortr, doublereal orti, doublereal hr, doublereal
4931	hi, doublereal wr, doublereal wi, doublereal zr, doublereal zi,
4932	integer *ierr)
4933	{
4934	/* System generated locals */
4935	integer hr_dim1, hr_offset, hi_dim1, hi_offset, zr_dim1, zr_offset,
4936	zi_dim1, zi_offset, i_1, i_2, i_3;
4937	doublereal d_1, d_2, d_3, d_4;
4938
4939	/* Local variables */
4940	static integer iend;
4941	extern /* Subroutine / int cdiv_(doublereal , doublereal , doublereal
4942	, doublereal , doublereal , doublereal *);
4943	static doublereal norm;
4944	static integer i, j, k, l, m, ii, en, jj, ll, nn;
4945	static doublereal si, ti, xi, yi, sr, tr, xr, yr;
4946	extern doublereal pythag_(doublereal , doublereal );
4947	extern /* Subroutine / int csroot_(doublereal , doublereal *,
4948	doublereal , doublereal );
4949	static integer ip1, lp1, itn, its;
4950	static doublereal zzi, zzr;
4951	static integer enm1;
4952	static doublereal tst1, tst2;
4953
4954
4955
4956	/* THIS SUBROUTINE IS A TRANSLATION OF A UNITARY ANALOGUE OF THE */
4957	/* ALGOL PROCEDURE COMLR2, NUM. MATH. 16, 181-204(1970) BY PETERS */
4958	/* AND WILKINSON. */
4959	/* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971). */
4960	/* THE UNITARY ANALOGUE SUBSTITUTES THE QR ALGORITHM OF FRANCIS */
4961	/* (COMP. JOUR. 4, 332-345(1962)) FOR THE LR ALGORITHM. */
4962
4963	/* THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS */
4964	/* OF A COMPLEX UPPER HESSENBERG MATRIX BY THE QR */
4965	/* METHOD. THE EIGENVECTORS OF A COMPLEX GENERAL MATRIX */
4966	/* CAN ALSO BE FOUND IF CORTH HAS BEEN USED TO REDUCE */
4967	/* THIS GENERAL MATRIX TO HESSENBERG FORM. */
4968
4969	/* ON INPUT */
4970
4971	/* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
4972	/* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
4973	/* DIMENSION STATEMENT. */
4974
4975	/* N IS THE ORDER OF THE MATRIX. */
4976
4977	/* LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
4978	/* SUBROUTINE CBAL. IF CBAL HAS NOT BEEN USED, */
4979	/* SET LOW=1, IGH=N. */
4980
4981	/* ORTR AND ORTI CONTAIN INFORMATION ABOUT THE UNITARY TRANS- */
4982	/* FORMATIONS USED IN THE REDUCTION BY CORTH, IF PERFORMED. */
4983	/* ONLY ELEMENTS LOW THROUGH IGH ARE USED. IF THE EIGENVECTORS
4984	*/
4985	/* OF THE HESSENBERG MATRIX ARE DESIRED, SET ORTR(J) AND */
4986	/* ORTI(J) TO 0.0D0 FOR THESE ELEMENTS. */
4987
4988	/* HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS, */
4989	/* RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX. */
4990	/* THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN FURTHER */
4991	/* INFORMATION ABOUT THE TRANSFORMATIONS WHICH WERE USED IN THE
4992	*/
4993	/* REDUCTION BY CORTH, IF PERFORMED. IF THE EIGENVECTORS OF */
4994	/* THE HESSENBERG MATRIX ARE DESIRED, THESE ELEMENTS MAY BE */
4995	/* ARBITRARY. */
4996
4997	/* ON OUTPUT */
4998
4999	/* ORTR, ORTI, AND THE UPPER HESSENBERG PORTIONS OF HR AND HI */
5000	/* HAVE BEEN DESTROYED. */
5001
5002	/* WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, */
5003	/* RESPECTIVELY, OF THE EIGENVALUES. IF AN ERROR */
5004	/* EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT */
5005	/* FOR INDICES IERR+1,...,N. */
5006
5007	/* ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, */
5008	/* RESPECTIVELY, OF THE EIGENVECTORS. THE EIGENVECTORS */
5009	/* ARE UNNORMALIZED. IF AN ERROR EXIT IS MADE, NONE OF */
5010	/* THE EIGENVECTORS HAS BEEN FOUND. */
5011
5012	/* IERR IS SET TO */
5013	/* ZERO FOR NORMAL RETURN, */
5014	/* J IF THE LIMIT OF 30N ITERATIONS IS EXHAUSTED /
5015	/* WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. */
5016
5017	/* CALLS CDIV FOR COMPLEX DIVISION. */
5018	/* CALLS CSROOT FOR COMPLEX SQUARE ROOT. */
5019	/* CALLS PYTHAG FOR DSQRT(AA + BB) . */
5020
5021	/* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
5022	/* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
5023	*/
5024
5025	/* THIS VERSION DATED AUGUST 1983. */
5026
5027	/* ------------------------------------------------------------------
5028	*/
5029
5030	/* Parameter adjustments */
5031	zi_dim1 = *nm;
5032	zi_offset = zi_dim1 + 1;
5033	zi -= zi_offset;
5034	zr_dim1 = *nm;
5035	zr_offset = zr_dim1 + 1;
5036	zr -= zr_offset;
5037	--wi;
5038	--wr;
5039	hi_dim1 = *nm;
5040	hi_offset = hi_dim1 + 1;
5041	hi -= hi_offset;
5042	hr_dim1 = *nm;
5043	hr_offset = hr_dim1 + 1;
5044	hr -= hr_offset;
5045	--orti;
5046	--ortr;
5047
5048	/* Function Body */
5049	*ierr = 0;
5050	/* .......... INITIALIZE EIGENVECTOR MATRIX .......... */
5051	i_1 = *n;
5052	for (j = 1; j <= i_1; ++j) {
5053
5054	i_2 = *n;
5055	for (i = 1; i <= i_2; ++i) {
5056	zr[i + j * zr_dim1] = 0.;
5057	zi[i + j * zi_dim1] = 0.;
5058	/* L100: */
5059	}
5060	zr[j + j * zr_dim1] = 1.;
5061	/* L101: */
5062	}
5063	/* .......... FORM THE MATRIX OF ACCUMULATED TRANSFORMATIONS */
5064	/* FROM THE INFORMATION LEFT BY CORTH .......... */
5065	iend = igh - low - 1;
5066	if (iend < 0) {
5067	goto L180;
5068	} else if (iend == 0) {
5069	goto L150;
5070	} else {
5071	goto L105;
5072	}
5073	/* .......... FOR I=IGH-1 STEP -1 UNTIL LOW+1 DO -- .......... */
5074	L105:
5075	i_1 = iend;
5076	for (ii = 1; ii <= i_1; ++ii) {
5077	i = *igh - ii;
5078	if (ortr[i] == 0. && orti[i] == 0.) {
5079	goto L140;
5080	}
5081	if (hr[i + (i - 1) * hr_dim1] == 0. && hi[i + (i - 1) * hi_dim1] ==
5082	0.) {
5083	goto L140;
5084	}
5085	/* .......... NORM BELOW IS NEGATIVE OF H FORMED IN CORTH ........
5086	.. */
5087	norm = hr[i + (i - 1) * hr_dim1] * ortr[i] + hi[i + (i - 1) * hi_dim1]
5088	* orti[i];
5089	ip1 = i + 1;
5090
5091	i_2 = *igh;
5092	for (k = ip1; k <= i_2; ++k) {
5093	ortr[k] = hr[k + (i - 1) * hr_dim1];
5094	orti[k] = hi[k + (i - 1) * hi_dim1];
5095	/* L110: */
5096	}
5097
5098	i_2 = *igh;
5099	for (j = i; j <= i_2; ++j) {
5100	sr = 0.;
5101	si = 0.;
5102
5103	i_3 = *igh;
5104	for (k = i; k <= i_3; ++k) {
5105	sr = sr + ortr[k] * zr[k + j * zr_dim1] + orti[k] * zi[k + j *
5106	zi_dim1];
5107	si = si + ortr[k] * zi[k + j * zi_dim1] - orti[k] * zr[k + j *
5108	zr_dim1];
5109	/* L115: */
5110	}
5111
5112	sr /= norm;
5113	si /= norm;
5114
5115	i_3 = *igh;
5116	for (k = i; k <= i_3; ++k) {
5117	zr[k + j * zr_dim1] = zr[k + j * zr_dim1] + sr * ortr[k] - si
5118	* orti[k];
5119	zi[k + j * zi_dim1] = zi[k + j * zi_dim1] + sr * orti[k] + si
5120	* ortr[k];
5121	/* L120: */
5122	}
5123
5124	/* L130: */
5125	}
5126
5127	L140:
5128	;
5129	}
5130	/* .......... CREATE REAL SUBDIAGONAL ELEMENTS .......... */
5131	L150:
5132	l = *low + 1;
5133
5134	i_1 = *igh;
5135	for (i = l; i <= i_1; ++i) {
5136	/* Computing MIN */
5137	i_2 = i + 1;
5138	ll = min(i_2,*igh);
5139	if (hi[i + (i - 1) * hi_dim1] == 0.) {
5140	goto L170;
5141	}
5142	norm = pythag_(&hr[i + (i - 1) * hr_dim1], &hi[i + (i - 1) * hi_dim1])
5143	;
5144	yr = hr[i + (i - 1) * hr_dim1] / norm;
5145	yi = hi[i + (i - 1) * hi_dim1] / norm;
5146	hr[i + (i - 1) * hr_dim1] = norm;
5147	hi[i + (i - 1) * hi_dim1] = 0.;
5148
5149	i_2 = *n;
5150	for (j = i; j <= i_2; ++j) {
5151	si = yr * hi[i + j * hi_dim1] - yi * hr[i + j * hr_dim1];
5152	hr[i + j * hr_dim1] = yr * hr[i + j * hr_dim1] + yi * hi[i + j *
5153	hi_dim1];
5154	hi[i + j * hi_dim1] = si;
5155	/* L155: */
5156	}
5157
5158	i_2 = ll;
5159	for (j = 1; j <= i_2; ++j) {
5160	si = yr * hi[j + i * hi_dim1] + yi * hr[j + i * hr_dim1];
5161	hr[j + i * hr_dim1] = yr * hr[j + i * hr_dim1] - yi * hi[j + i *
5162	hi_dim1];
5163	hi[j + i * hi_dim1] = si;
5164	/* L160: */
5165	}
5166
5167	i_2 = *igh;
5168	for (j = *low; j <= i_2; ++j) {
5169	si = yr * zi[j + i * zi_dim1] + yi * zr[j + i * zr_dim1];
5170	zr[j + i * zr_dim1] = yr * zr[j + i * zr_dim1] - yi * zi[j + i *
5171	zi_dim1];
5172	zi[j + i * zi_dim1] = si;
5173	/* L165: */
5174	}
5175
5176	L170:
5177	;
5178	}
5179	/* .......... STORE ROOTS ISOLATED BY CBAL .......... */
5180	L180:
5181	i_1 = *n;
5182	for (i = 1; i <= i_1; ++i) {
5183	if (i >= low && i <= igh) {
5184	goto L200;
5185	}
5186	wr[i] = hr[i + i * hr_dim1];
5187	wi[i] = hi[i + i * hi_dim1];
5188	L200:
5189	;
5190	}
5191
5192	en = *igh;
5193	tr = 0.;
5194	ti = 0.;
5195	itn = n 30;
5196	/* .......... SEARCH FOR NEXT EIGENVALUE .......... */
5197	L220:
5198	if (en < *low) {
5199	goto L680;
5200	}
5201	its = 0;
5202	enm1 = en - 1;
5203	/* .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT */
5204	/* FOR L=EN STEP -1 UNTIL LOW DO -- .......... */
5205	L240:
5206	i_1 = en;
5207	for (ll = *low; ll <= i_1; ++ll) {
5208	l = en + *low - ll;
5209	if (l == *low) {
5210	goto L300;
5211	}
5212	tst1 = (d_1 = hr[l - 1 + (l - 1) * hr_dim1], abs(d_1)) + (d_2 = hi[
5213	l - 1 + (l - 1) * hi_dim1], abs(d_2)) + (d_3 = hr[l + l *
5214	hr_dim1], abs(d_3)) + (d_4 = hi[l + l * hi_dim1], abs(d_4))
5215	;
5216	tst2 = tst1 + (d_1 = hr[l + (l - 1) * hr_dim1], abs(d_1));
5217	if (tst2 == tst1) {
5218	goto L300;
5219	}
5220	/* L260: */
5221	}
5222	/* .......... FORM SHIFT .......... */
5223	L300:
5224	if (l == en) {
5225	goto L660;
5226	}
5227	if (itn == 0) {
5228	goto L1000;
5229	}
5230	if (its == 10 \|\| its == 20) {
5231	goto L320;
5232	}
5233	sr = hr[en + en * hr_dim1];
5234	si = hi[en + en * hi_dim1];
5235	xr = hr[enm1 + en * hr_dim1] * hr[en + enm1 * hr_dim1];
5236	xi = hi[enm1 + en * hi_dim1] * hr[en + enm1 * hr_dim1];
5237	if (xr == 0. && xi == 0.) {
5238	goto L340;
5239	}
5240	yr = (hr[enm1 + enm1 * hr_dim1] - sr) / 2.;
5241	yi = (hi[enm1 + enm1 * hi_dim1] - si) / 2.;
5242	/* Computing 2nd power */
5243	d_2 = yr;
5244	/* Computing 2nd power */
5245	d_3 = yi;
5246	d_1 = d_2 * d_2 - d_3 * d_3 + xr;
5247	d_4 = yr * 2. * yi + xi;
5248	csroot_(&d_1, &d_4, &zzr, &zzi);
5249	if (yr * zzr + yi * zzi >= 0.) {
5250	goto L310;
5251	}
5252	zzr = -zzr;
5253	zzi = -zzi;
5254	L310:
5255	d_1 = yr + zzr;
5256	d_2 = yi + zzi;
5257	cdiv_(&xr, &xi, &d_1, &d_2, &xr, &xi);
5258	sr -= xr;
5259	si -= xi;
5260	goto L340;
5261	/* .......... FORM EXCEPTIONAL SHIFT .......... */
5262	L320:
5263	sr = (d_1 = hr[en + enm1 * hr_dim1], abs(d_1)) + (d_2 = hr[enm1 + (en
5264	- 2) * hr_dim1], abs(d_2));
5265	si = 0.;
5266
5267	L340:
5268	i_1 = en;
5269	for (i = *low; i <= i_1; ++i) {
5270	hr[i + i * hr_dim1] -= sr;
5271	hi[i + i * hi_dim1] -= si;
5272	/* L360: */
5273	}
5274
5275	tr += sr;
5276	ti += si;
5277	++its;
5278	--itn;
5279	/* .......... REDUCE TO TRIANGLE (ROWS) .......... */
5280	lp1 = l + 1;
5281
5282	i_1 = en;
5283	for (i = lp1; i <= i_1; ++i) {
5284	sr = hr[i + (i - 1) * hr_dim1];
5285	hr[i + (i - 1) * hr_dim1] = 0.;
5286	d_1 = pythag_(&hr[i - 1 + (i - 1) * hr_dim1], &hi[i - 1 + (i - 1) *
5287	hi_dim1]);
5288	norm = pythag_(&d_1, &sr);
5289	xr = hr[i - 1 + (i - 1) * hr_dim1] / norm;
5290	wr[i - 1] = xr;
5291	xi = hi[i - 1 + (i - 1) * hi_dim1] / norm;
5292	wi[i - 1] = xi;
5293	hr[i - 1 + (i - 1) * hr_dim1] = norm;
5294	hi[i - 1 + (i - 1) * hi_dim1] = 0.;
5295	hi[i + (i - 1) * hi_dim1] = sr / norm;
5296
5297	i_2 = *n;
5298	for (j = i; j <= i_2; ++j) {
5299	yr = hr[i - 1 + j * hr_dim1];
5300	yi = hi[i - 1 + j * hi_dim1];
5301	zzr = hr[i + j * hr_dim1];
5302	zzi = hi[i + j * hi_dim1];
5303	hr[i - 1 + j * hr_dim1] = xr * yr + xi * yi + hi[i + (i - 1) *
5304	hi_dim1] * zzr;
5305	hi[i - 1 + j * hi_dim1] = xr * yi - xi * yr + hi[i + (i - 1) *
5306	hi_dim1] * zzi;
5307	hr[i + j * hr_dim1] = xr * zzr - xi * zzi - hi[i + (i - 1) *
5308	hi_dim1] * yr;
5309	hi[i + j * hi_dim1] = xr * zzi + xi * zzr - hi[i + (i - 1) *
5310	hi_dim1] * yi;
5311	/* L490: */
5312	}
5313
5314	/* L500: */
5315	}
5316
5317	si = hi[en + en * hi_dim1];
5318	if (si == 0.) {
5319	goto L540;
5320	}
5321	norm = pythag_(&hr[en + en * hr_dim1], &si);
5322	sr = hr[en + en * hr_dim1] / norm;
5323	si /= norm;
5324	hr[en + en * hr_dim1] = norm;
5325	hi[en + en * hi_dim1] = 0.;
5326	if (en == *n) {
5327	goto L540;
5328	}
5329	ip1 = en + 1;
5330
5331	i_1 = *n;
5332	for (j = ip1; j <= i_1; ++j) {
5333	yr = hr[en + j * hr_dim1];
5334	yi = hi[en + j * hi_dim1];
5335	hr[en + j * hr_dim1] = sr * yr + si * yi;
5336	hi[en + j * hi_dim1] = sr * yi - si * yr;
5337	/* L520: */
5338	}
5339	/* .......... INVERSE OPERATION (COLUMNS) .......... */
5340	L540:
5341	i_1 = en;
5342	for (j = lp1; j <= i_1; ++j) {
5343	xr = wr[j - 1];
5344	xi = wi[j - 1];
5345
5346	i_2 = j;
5347	for (i = 1; i <= i_2; ++i) {
5348	yr = hr[i + (j - 1) * hr_dim1];
5349	yi = 0.;
5350	zzr = hr[i + j * hr_dim1];
5351	zzi = hi[i + j * hi_dim1];
5352	if (i == j) {
5353	goto L560;
5354	}
5355	yi = hi[i + (j - 1) * hi_dim1];
5356	hi[i + (j - 1) * hi_dim1] = xr * yi + xi * yr + hi[j + (j - 1) *
5357	hi_dim1] * zzi;
5358	L560:
5359	hr[i + (j - 1) * hr_dim1] = xr * yr - xi * yi + hi[j + (j - 1) *
5360	hi_dim1] * zzr;
5361	hr[i + j * hr_dim1] = xr * zzr + xi * zzi - hi[j + (j - 1) *
5362	hi_dim1] * yr;
5363	hi[i + j * hi_dim1] = xr * zzi - xi * zzr - hi[j + (j - 1) *
5364	hi_dim1] * yi;
5365	/* L580: */
5366	}
5367
5368	i_2 = *igh;
5369	for (i = *low; i <= i_2; ++i) {
5370	yr = zr[i + (j - 1) * zr_dim1];
5371	yi = zi[i + (j - 1) * zi_dim1];
5372	zzr = zr[i + j * zr_dim1];
5373	zzi = zi[i + j * zi_dim1];
5374	zr[i + (j - 1) * zr_dim1] = xr * yr - xi * yi + hi[j + (j - 1) *
5375	hi_dim1] * zzr;
5376	zi[i + (j - 1) * zi_dim1] = xr * yi + xi * yr + hi[j + (j - 1) *
5377	hi_dim1] * zzi;
5378	zr[i + j * zr_dim1] = xr * zzr + xi * zzi - hi[j + (j - 1) *
5379	hi_dim1] * yr;
5380	zi[i + j * zi_dim1] = xr * zzi - xi * zzr - hi[j + (j - 1) *
5381	hi_dim1] * yi;
5382	/* L590: */
5383	}
5384
5385	/* L600: */
5386	}
5387
5388	if (si == 0.) {
5389	goto L240;
5390	}
5391
5392	i_1 = en;
5393	for (i = 1; i <= i_1; ++i) {
5394	yr = hr[i + en * hr_dim1];
5395	yi = hi[i + en * hi_dim1];
5396	hr[i + en * hr_dim1] = sr * yr - si * yi;
5397	hi[i + en * hi_dim1] = sr * yi + si * yr;
5398	/* L630: */
5399	}
5400
5401	i_1 = *igh;
5402	for (i = *low; i <= i_1; ++i) {
5403	yr = zr[i + en * zr_dim1];
5404	yi = zi[i + en * zi_dim1];
5405	zr[i + en * zr_dim1] = sr * yr - si * yi;
5406	zi[i + en * zi_dim1] = sr * yi + si * yr;
5407	/* L640: */
5408	}
5409
5410	goto L240;
5411	/* .......... A ROOT FOUND .......... */
5412	L660:
5413	hr[en + en * hr_dim1] += tr;
5414	wr[en] = hr[en + en * hr_dim1];
5415	hi[en + en * hi_dim1] += ti;
5416	wi[en] = hi[en + en * hi_dim1];
5417	en = enm1;
5418	goto L220;
5419	/* .......... ALL ROOTS FOUND. BACKSUBSTITUTE TO FIND */
5420	/* VECTORS OF UPPER TRIANGULAR FORM .......... */
5421	L680:
5422	norm = 0.;
5423
5424	i_1 = *n;
5425	for (i = 1; i <= i_1; ++i) {
5426
5427	i_2 = *n;
5428	for (j = i; j <= i_2; ++j) {
5429	tr = (d_1 = hr[i + j * hr_dim1], abs(d_1)) + (d_2 = hi[i + j *
5430	hi_dim1], abs(d_2));
5431	if (tr > norm) {
5432	norm = tr;
5433	}
5434	/* L720: */
5435	}
5436	}
5437
5438	if (*n == 1 \|\| norm == 0.) {
5439	goto L1001;
5440	}
5441	/* .......... FOR EN=N STEP -1 UNTIL 2 DO -- .......... */
5442	i_2 = *n;
5443	for (nn = 2; nn <= i_2; ++nn) {
5444	en = *n + 2 - nn;
5445	xr = wr[en];
5446	xi = wi[en];
5447	hr[en + en * hr_dim1] = 1.;
5448	hi[en + en * hi_dim1] = 0.;
5449	enm1 = en - 1;
5450	/* .......... FOR I=EN-1 STEP -1 UNTIL 1 DO -- .......... */
5451	i_1 = enm1;
5452	for (ii = 1; ii <= i_1; ++ii) {
5453	i = en - ii;
5454	zzr = 0.;
5455	zzi = 0.;
5456	ip1 = i + 1;
5457
5458	i_3 = en;
5459	for (j = ip1; j <= i_3; ++j) {
5460	zzr = zzr + hr[i + j * hr_dim1] * hr[j + en * hr_dim1] - hi[i
5461	+ j * hi_dim1] * hi[j + en * hi_dim1];
5462	zzi = zzi + hr[i + j * hr_dim1] * hi[j + en * hi_dim1] + hi[i
5463	+ j * hi_dim1] * hr[j + en * hr_dim1];
5464	/* L740: */
5465	}
5466
5467	yr = xr - wr[i];
5468	yi = xi - wi[i];
5469	if (yr != 0. \|\| yi != 0.) {
5470	goto L765;
5471	}
5472	tst1 = norm;
5473	yr = tst1;
5474	L760:
5475	yr *= .01;
5476	tst2 = norm + yr;
5477	if (tst2 > tst1) {
5478	goto L760;
5479	}
5480	L765:
5481	cdiv_(&zzr, &zzi, &yr, &yi, &hr[i + en * hr_dim1], &hi[i + en *
5482	hi_dim1]);
5483	/* .......... OVERFLOW CONTROL .......... */
5484	tr = (d_1 = hr[i + en * hr_dim1], abs(d_1)) + (d_2 = hi[i + en
5485	* hi_dim1], abs(d_2));
5486	if (tr == 0.) {
5487	goto L780;
5488	}
5489	tst1 = tr;
5490	tst2 = tst1 + 1. / tst1;
5491	if (tst2 > tst1) {
5492	goto L780;
5493	}
5494	i_3 = en;
5495	for (j = i; j <= i_3; ++j) {
5496	hr[j + en * hr_dim1] /= tr;
5497	hi[j + en * hi_dim1] /= tr;
5498	/* L770: */
5499	}
5500
5501	L780:
5502	;
5503	}
5504
5505	/* L800: */
5506	}
5507	/* .......... END BACKSUBSTITUTION .......... */
5508	enm1 = *n - 1;
5509	/* .......... VECTORS OF ISOLATED ROOTS .......... */
5510	i_2 = enm1;
5511	for (i = 1; i <= i_2; ++i) {
5512	if (i >= low && i <= igh) {
5513	goto L840;
5514	}
5515	ip1 = i + 1;
5516
5517	i_1 = *n;
5518	for (j = ip1; j <= i_1; ++j) {
5519	zr[i + j * zr_dim1] = hr[i + j * hr_dim1];
5520	zi[i + j * zi_dim1] = hi[i + j * hi_dim1];
5521	/* L820: */
5522	}
5523
5524	L840:
5525	;
5526	}
5527	/* .......... MULTIPLY BY TRANSFORMATION MATRIX TO GIVE */
5528	/* VECTORS OF ORIGINAL FULL MATRIX. */
5529	/* FOR J=N STEP -1 UNTIL LOW+1 DO -- .......... */
5530	i_2 = enm1;
5531	for (jj = *low; jj <= i_2; ++jj) {
5532	j = n + low - jj;
5533	m = min(j,*igh);
5534
5535	i_1 = *igh;
5536	for (i = *low; i <= i_1; ++i) {
5537	zzr = 0.;
5538	zzi = 0.;
5539
5540	i_3 = m;
5541	for (k = *low; k <= i_3; ++k) {
5542	zzr = zzr + zr[i + k * zr_dim1] * hr[k + j * hr_dim1] - zi[i
5543	+ k * zi_dim1] * hi[k + j * hi_dim1];
5544	zzi = zzi + zr[i + k * zr_dim1] * hi[k + j * hi_dim1] + zi[i
5545	+ k * zi_dim1] * hr[k + j * hr_dim1];
5546	/* L860: */
5547	}
5548
5549	zr[i + j * zr_dim1] = zzr;
5550	zi[i + j * zi_dim1] = zzi;
5551	/* L880: */
5552	}
5553	}
5554
5555	goto L1001;
5556	/* .......... SET ERROR -- ALL EIGENVALUES HAVE NOT */
5557	/* CONVERGED AFTER 30N ITERATIONS .......... /
5558	L1000:
5559	*ierr = en;
5560	L1001:
5561	return 0;
5562	} /* comqr2_ */
5563
5564	/* Subroutine / int cortb_(integer nm, integer low, integer igh,
5565	doublereal ar, doublereal ai, doublereal ortr, doublereal orti,
5566	integer m, doublereal zr, doublereal *zi)
5567	{
5568	/* System generated locals */
5569	integer ar_dim1, ar_offset, ai_dim1, ai_offset, zr_dim1, zr_offset,
5570	zi_dim1, zi_offset, i_1, i_2, i_3;
5571
5572	/* Local variables */
5573	static doublereal h;
5574	static integer i, j, la;
5575	static doublereal gi, gr;
5576	static integer mm, mp, kp1, mp1;
5577
5578
5579
5580	/* THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF */
5581	/* THE ALGOL PROCEDURE ORTBAK, NUM. MATH. 12, 349-368(1968) */
5582	/* BY MARTIN AND WILKINSON. */
5583	/* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971). */
5584
5585	/* THIS SUBROUTINE FORMS THE EIGENVECTORS OF A COMPLEX GENERAL */
5586	/* MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING */
5587	/* UPPER HESSENBERG MATRIX DETERMINED BY CORTH. */
5588
5589	/* ON INPUT */
5590
5591	/* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
5592	/* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
5593	/* DIMENSION STATEMENT. */
5594
5595	/* LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
5596	/* SUBROUTINE CBAL. IF CBAL HAS NOT BEEN USED, */
5597	/* SET LOW=1 AND IGH EQUAL TO THE ORDER OF THE MATRIX. */
5598
5599	/* AR AND AI CONTAIN INFORMATION ABOUT THE UNITARY */
5600	/* TRANSFORMATIONS USED IN THE REDUCTION BY CORTH */
5601	/* IN THEIR STRICT LOWER TRIANGLES. */
5602
5603	/* ORTR AND ORTI CONTAIN FURTHER INFORMATION ABOUT THE */
5604	/* TRANSFORMATIONS USED IN THE REDUCTION BY CORTH. */
5605	/* ONLY ELEMENTS LOW THROUGH IGH ARE USED. */
5606
5607	/* M IS THE NUMBER OF COLUMNS OF ZR AND ZI TO BE BACK TRANSFORMED.
5608	*/
5609
5610	/* ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, */
5611	/* RESPECTIVELY, OF THE EIGENVECTORS TO BE */
5612	/* BACK TRANSFORMED IN THEIR FIRST M COLUMNS. */
5613
5614	/* ON OUTPUT */
5615
5616	/* ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, */
5617	/* RESPECTIVELY, OF THE TRANSFORMED EIGENVECTORS */
5618	/* IN THEIR FIRST M COLUMNS. */
5619
5620	/* ORTR AND ORTI HAVE BEEN ALTERED. */
5621
5622	/* NOTE THAT CORTB PRESERVES VECTOR EUCLIDEAN NORMS. */
5623
5624	/* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
5625	/* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
5626	*/
5627
5628	/* THIS VERSION DATED AUGUST 1983. */
5629
5630	/* ------------------------------------------------------------------
5631	*/
5632
5633	/* Parameter adjustments */
5634	--orti;
5635	--ortr;
5636	ai_dim1 = *nm;
5637	ai_offset = ai_dim1 + 1;
5638	ai -= ai_offset;
5639	ar_dim1 = *nm;
5640	ar_offset = ar_dim1 + 1;
5641	ar -= ar_offset;
5642	zi_dim1 = *nm;
5643	zi_offset = zi_dim1 + 1;
5644	zi -= zi_offset;
5645	zr_dim1 = *nm;
5646	zr_offset = zr_dim1 + 1;
5647	zr -= zr_offset;
5648
5649	/* Function Body */
5650	if (*m == 0) {
5651	goto L200;
5652	}
5653	la = *igh - 1;
5654	kp1 = *low + 1;
5655	if (la < kp1) {
5656	goto L200;
5657	}
5658	/* .......... FOR MP=IGH-1 STEP -1 UNTIL LOW+1 DO -- .......... */
5659	i_1 = la;
5660	for (mm = kp1; mm <= i_1; ++mm) {
5661	mp = low + igh - mm;
5662	if (ar[mp + (mp - 1) * ar_dim1] == 0. && ai[mp + (mp - 1) * ai_dim1]
5663	== 0.) {
5664	goto L140;
5665	}
5666	/* .......... H BELOW IS NEGATIVE OF H FORMED IN CORTH ..........
5667	*/
5668	h = ar[mp + (mp - 1) * ar_dim1] * ortr[mp] + ai[mp + (mp - 1) *
5669	ai_dim1] * orti[mp];
5670	mp1 = mp + 1;
5671
5672	i_2 = *igh;
5673	for (i = mp1; i <= i_2; ++i) {
5674	ortr[i] = ar[i + (mp - 1) * ar_dim1];
5675	orti[i] = ai[i + (mp - 1) * ai_dim1];
5676	/* L100: */
5677	}
5678
5679	i_2 = *m;
5680	for (j = 1; j <= i_2; ++j) {
5681	gr = 0.;
5682	gi = 0.;
5683
5684	i_3 = *igh;
5685	for (i = mp; i <= i_3; ++i) {
5686	gr = gr + ortr[i] * zr[i + j * zr_dim1] + orti[i] * zi[i + j *
5687	zi_dim1];
5688	gi = gi + ortr[i] * zi[i + j * zi_dim1] - orti[i] * zr[i + j *
5689	zr_dim1];
5690	/* L110: */
5691	}
5692
5693	gr /= h;
5694	gi /= h;
5695
5696	i_3 = *igh;
5697	for (i = mp; i <= i_3; ++i) {
5698	zr[i + j * zr_dim1] = zr[i + j * zr_dim1] + gr * ortr[i] - gi
5699	* orti[i];
5700	zi[i + j * zi_dim1] = zi[i + j * zi_dim1] + gr * orti[i] + gi
5701	* ortr[i];
5702	/* L120: */
5703	}
5704
5705	/* L130: */
5706	}
5707
5708	L140:
5709	;
5710	}
5711
5712	L200:
5713	return 0;
5714	} /* cortb_ */
5715
5716	/* Subroutine / int corth_(integer nm, integer n, integer low, integer *
5717	igh, doublereal ar, doublereal ai, doublereal ortr, doublereal
5718	orti)
5719	{
5720	/* System generated locals */
5721	integer ar_dim1, ar_offset, ai_dim1, ai_offset, i_1, i_2, i_3;
5722	doublereal d_1, d_2;
5723
5724	/* Builtin functions */
5725	double sqrt(doublereal);
5726
5727	/* Local variables */
5728	static doublereal f, g, h;
5729	static integer i, j, m;
5730	static doublereal scale;
5731	static integer la;
5732	static doublereal fi;
5733	static integer ii, jj;
5734	static doublereal fr;
5735	static integer mp;
5736	extern doublereal pythag_(doublereal , doublereal );
5737	static integer kp1;
5738
5739
5740
5741	/* THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF */
5742	/* THE ALGOL PROCEDURE ORTHES, NUM. MATH. 12, 349-368(1968) */
5743	/* BY MARTIN AND WILKINSON. */
5744	/* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971). */
5745
5746	/* GIVEN A COMPLEX GENERAL MATRIX, THIS SUBROUTINE */
5747	/* REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS */
5748	/* LOW THROUGH IGH TO UPPER HESSENBERG FORM BY */
5749	/* UNITARY SIMILARITY TRANSFORMATIONS. */
5750
5751	/* ON INPUT */
5752
5753	/* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
5754	/* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
5755	/* DIMENSION STATEMENT. */
5756
5757	/* N IS THE ORDER OF THE MATRIX. */
5758
5759	/* LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
5760	/* SUBROUTINE CBAL. IF CBAL HAS NOT BEEN USED, */
5761	/* SET LOW=1, IGH=N. */
5762
5763	/* AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, */
5764	/* RESPECTIVELY, OF THE COMPLEX INPUT MATRIX. */
5765
5766	/* ON OUTPUT */
5767
5768	/* AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS, */
5769	/* RESPECTIVELY, OF THE HESSENBERG MATRIX. INFORMATION */
5770	/* ABOUT THE UNITARY TRANSFORMATIONS USED IN THE REDUCTION */
5771	/* IS STORED IN THE REMAINING TRIANGLES UNDER THE */
5772	/* HESSENBERG MATRIX. */
5773
5774	/* ORTR AND ORTI CONTAIN FURTHER INFORMATION ABOUT THE */
5775	/* TRANSFORMATIONS. ONLY ELEMENTS LOW THROUGH IGH ARE USED. */
5776
5777	/* CALLS PYTHAG FOR DSQRT(AA + BB) . */
5778
5779	/* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
5780	/* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
5781	*/
5782
5783	/* THIS VERSION DATED AUGUST 1983. */
5784
5785	/* ------------------------------------------------------------------
5786	*/
5787
5788	/* Parameter adjustments */
5789	ai_dim1 = *nm;
5790	ai_offset = ai_dim1 + 1;
5791	ai -= ai_offset;
5792	ar_dim1 = *nm;
5793	ar_offset = ar_dim1 + 1;
5794	ar -= ar_offset;
5795	--orti;
5796	--ortr;
5797
5798	/* Function Body */
5799	la = *igh - 1;
5800	kp1 = *low + 1;
5801	if (la < kp1) {
5802	goto L200;
5803	}
5804
5805	i_1 = la;
5806	for (m = kp1; m <= i_1; ++m) {
5807	h = 0.;
5808	ortr[m] = 0.;
5809	orti[m] = 0.;
5810	scale = 0.;
5811	/* .......... SCALE COLUMN (ALGOL TOL THEN NOT NEEDED) ..........
5812	*/
5813	i_2 = *igh;
5814	for (i = m; i <= i_2; ++i) {
5815	/* L90: */
5816	scale = scale + (d_1 = ar[i + (m - 1) * ar_dim1], abs(d_1)) + (
5817	d_2 = ai[i + (m - 1) * ai_dim1], abs(d_2));
5818	}
5819
5820	if (scale == 0.) {
5821	goto L180;
5822	}
5823	mp = m + *igh;
5824	/* .......... FOR I=IGH STEP -1 UNTIL M DO -- .......... */
5825	i_2 = *igh;
5826	for (ii = m; ii <= i_2; ++ii) {
5827	i = mp - ii;
5828	ortr[i] = ar[i + (m - 1) * ar_dim1] / scale;
5829	orti[i] = ai[i + (m - 1) * ai_dim1] / scale;
5830	h = h + ortr[i] * ortr[i] + orti[i] * orti[i];
5831	/* L100: */
5832	}
5833
5834	g = sqrt(h);
5835	f = pythag_(&ortr[m], &orti[m]);
5836	if (f == 0.) {
5837	goto L103;
5838	}
5839	h += f * g;
5840	g /= f;
5841	ortr[m] = (g + 1.) * ortr[m];
5842	orti[m] = (g + 1.) * orti[m];
5843	goto L105;
5844
5845	L103:
5846	ortr[m] = g;
5847	ar[m + (m - 1) * ar_dim1] = scale;
5848	/* .......... FORM (I-(UUT)/H) A .......... */
5849	L105:
5850	i_2 = *n;
5851	for (j = m; j <= i_2; ++j) {
5852	fr = 0.;
5853	fi = 0.;
5854	/* .......... FOR I=IGH STEP -1 UNTIL M DO -- .......... */
5855	i_3 = *igh;
5856	for (ii = m; ii <= i_3; ++ii) {
5857	i = mp - ii;
5858	fr = fr + ortr[i] * ar[i + j * ar_dim1] + orti[i] * ai[i + j *
5859	ai_dim1];
5860	fi = fi + ortr[i] * ai[i + j * ai_dim1] - orti[i] * ar[i + j *
5861	ar_dim1];
5862	/* L110: */
5863	}
5864
5865	fr /= h;
5866	fi /= h;
5867
5868	i_3 = *igh;
5869	for (i = m; i <= i_3; ++i) {
5870	ar[i + j * ar_dim1] = ar[i + j * ar_dim1] - fr * ortr[i] + fi
5871	* orti[i];
5872	ai[i + j * ai_dim1] = ai[i + j * ai_dim1] - fr * orti[i] - fi
5873	* ortr[i];
5874	/* L120: */
5875	}
5876
5877	/* L130: */
5878	}
5879	/* .......... FORM (I-(UUT)/H)A(I-(UUT)/H) .......... */
5880	i_2 = *igh;
5881	for (i = 1; i <= i_2; ++i) {
5882	fr = 0.;
5883	fi = 0.;
5884	/* .......... FOR J=IGH STEP -1 UNTIL M DO -- .......... */
5885	i_3 = *igh;
5886	for (jj = m; jj <= i_3; ++jj) {
5887	j = mp - jj;
5888	fr = fr + ortr[j] * ar[i + j * ar_dim1] - orti[j] * ai[i + j *
5889	ai_dim1];
5890	fi = fi + ortr[j] * ai[i + j * ai_dim1] + orti[j] * ar[i + j *
5891	ar_dim1];
5892	/* L140: */
5893	}
5894
5895	fr /= h;
5896	fi /= h;
5897
5898	i_3 = *igh;
5899	for (j = m; j <= i_3; ++j) {
5900	ar[i + j * ar_dim1] = ar[i + j * ar_dim1] - fr * ortr[j] - fi
5901	* orti[j];
5902	ai[i + j * ai_dim1] = ai[i + j * ai_dim1] + fr * orti[j] - fi
5903	* ortr[j];
5904	/* L150: */
5905	}
5906
5907	/* L160: */
5908	}
5909
5910	ortr[m] = scale * ortr[m];
5911	orti[m] = scale * orti[m];
5912	ar[m + (m - 1) * ar_dim1] = -g * ar[m + (m - 1) * ar_dim1];
5913	ai[m + (m - 1) * ai_dim1] = -g * ai[m + (m - 1) * ai_dim1];
5914	L180:
5915	;
5916	}
5917
5918	L200:
5919	return 0;
5920	} /* corth_ */
5921
5922	/* Subroutine / int elmbak_(integer nm, integer low, integer igh,
5923	doublereal a, integer int_, integer m, doublereal z)
5924	{
5925	/* System generated locals */
5926	integer a_dim1, a_offset, z_dim1, z_offset, i_1, i_2, i_3;
5927
5928	/* Local variables */
5929	static integer i, j;
5930	static doublereal x;
5931	static integer la, mm, mp, kp1, mp1;
5932
5933
5934
5935	/* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ELMBAK, */
5936	/* NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON. */
5937	/* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971). */
5938
5939	/* THIS SUBROUTINE FORMS THE EIGENVECTORS OF A REAL GENERAL */
5940	/* MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING */
5941	/* UPPER HESSENBERG MATRIX DETERMINED BY ELMHES. */
5942
5943	/* ON INPUT */
5944
5945	/* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
5946	/* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
5947	/* DIMENSION STATEMENT. */
5948
5949	/* LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
5950	/* SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED, */
5951	/* SET LOW=1 AND IGH EQUAL TO THE ORDER OF THE MATRIX. */
5952
5953	/* A CONTAINS THE MULTIPLIERS WHICH WERE USED IN THE */
5954	/* REDUCTION BY ELMHES IN ITS LOWER TRIANGLE */
5955	/* BELOW THE SUBDIAGONAL. */
5956
5957	/* INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS */
5958	/* INTERCHANGED IN THE REDUCTION BY ELMHES. */
5959	/* ONLY ELEMENTS LOW THROUGH IGH ARE USED. */
5960
5961	/* M IS THE NUMBER OF COLUMNS OF Z TO BE BACK TRANSFORMED. */
5962
5963	/* Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGEN- */
5964	/* VECTORS TO BE BACK TRANSFORMED IN ITS FIRST M COLUMNS. */
5965
5966	/* ON OUTPUT */
5967
5968	/* Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE */
5969	/* TRANSFORMED EIGENVECTORS IN ITS FIRST M COLUMNS. */
5970
5971	/* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
5972	/* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
5973	*/
5974
5975	/* THIS VERSION DATED AUGUST 1983. */
5976
5977	/* ------------------------------------------------------------------
5978	*/
5979
5980	/* Parameter adjustments */
5981	--int_;
5982	a_dim1 = *nm;
5983	a_offset = a_dim1 + 1;
5984	a -= a_offset;
5985	z_dim1 = *nm;
5986	z_offset = z_dim1 + 1;
5987	z -= z_offset;
5988
5989	/* Function Body */
5990	if (*m == 0) {
5991	goto L200;
5992	}
5993	la = *igh - 1;
5994	kp1 = *low + 1;
5995	if (la < kp1) {
5996	goto L200;
5997	}
5998	/* .......... FOR MP=IGH-1 STEP -1 UNTIL LOW+1 DO -- .......... */
5999	i_1 = la;
6000	for (mm = kp1; mm <= i_1; ++mm) {
6001	mp = low + igh - mm;
6002	mp1 = mp + 1;
6003
6004	i_2 = *igh;
6005	for (i = mp1; i <= i_2; ++i) {
6006	x = a[i + (mp - 1) * a_dim1];
6007	if (x == 0.) {
6008	goto L110;
6009	}
6010
6011	i_3 = *m;
6012	for (j = 1; j <= i_3; ++j) {
6013	/* L100: */
6014	z[i + j * z_dim1] += x * z[mp + j * z_dim1];
6015	}
6016
6017	L110:
6018	;
6019	}
6020
6021	i = int_[mp];
6022	if (i == mp) {
6023	goto L140;
6024	}
6025
6026	i_2 = *m;
6027	for (j = 1; j <= i_2; ++j) {
6028	x = z[i + j * z_dim1];
6029	z[i + j * z_dim1] = z[mp + j * z_dim1];
6030	z[mp + j * z_dim1] = x;
6031	/* L130: */
6032	}
6033
6034	L140:
6035	;
6036	}
6037
6038	L200:
6039	return 0;
6040	} /* elmbak_ */
6041
6042	/* Subroutine / int elmhes_(integer nm, integer n, integer low, integer *
6043	igh, doublereal a, integer int_)
6044	{
6045	/* System generated locals */
6046	integer a_dim1, a_offset, i_1, i_2, i_3;
6047	doublereal d_1;
6048
6049	/* Local variables */
6050	static integer i, j, m;
6051	static doublereal x, y;
6052	static integer la, mm1, kp1, mp1;
6053
6054
6055
6056	/* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ELMHES, */
6057	/* NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON. */
6058	/* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971). */
6059
6060	/* GIVEN A REAL GENERAL MATRIX, THIS SUBROUTINE */
6061	/* REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS */
6062	/* LOW THROUGH IGH TO UPPER HESSENBERG FORM BY */
6063	/* STABILIZED ELEMENTARY SIMILARITY TRANSFORMATIONS. */
6064
6065	/* ON INPUT */
6066
6067	/* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
6068	/* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
6069	/* DIMENSION STATEMENT. */
6070
6071	/* N IS THE ORDER OF THE MATRIX. */
6072
6073	/* LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING */
6074	/* SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED, */
6075	/* SET LOW=1, IGH=N. */
6076
6077	/* A CONTAINS THE INPUT MATRIX. */
6078
6079	/* ON OUTPUT */
6080
6081	/* A CONTAINS THE HESSENBERG MATRIX. THE MULTIPLIERS */
6082	/* WHICH WERE USED IN THE REDUCTION ARE STORED IN THE */
6083	/* REMAINING TRIANGLE UNDER THE HESSENBERG MATRIX. */
6084
6085	/* INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS */
6086	/* INTERCHANGED IN THE REDUCTION. */
6087	/* ONLY ELEMENTS LOW THROUGH IGH ARE USED. */
6088
6089	/* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */
6090	/* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
6091	*/
6092
6093	/* THIS VERSION DATED AUGUST 1983. */
6094
6095	/* ------------------------------------------------------------------
6096	*/
6097
6098	/* Parameter adjustments */
6099	a_dim1 = *nm;
6100	a_offset = a_dim1 + 1;
6101	a -= a_offset;
6102	--int_;
6103
6104	/* Function Body */
6105	la = *igh - 1;
6106	kp1 = *low + 1;
6107	if (la < kp1) {
6108	goto L200;
6109	}
6110
6111	i_1 = la;
6112	for (m = kp1; m <= i_1; ++m) {
6113	mm1 = m - 1;
6114	x = 0.;
6115	i = m;
6116
6117	i_2 = *igh;
6118	for (j = m; j <= i_2; ++j) {
6119	if ((d_1 = a[j + mm1 * a_dim1], abs(d_1)) <= abs(x)) {
6120	goto L100;
6121	}
6122	x = a[j + mm1 * a_dim1];
6123	i = j;
6124	L100:
6125	;
6126	}
6127
6128	int_[m] = i;
6129	if (i == m) {
6130	goto L130;
6131	}
6132	/* .......... INTERCHANGE ROWS AND COLUMNS OF A .......... */
6133	i_2 = *n;
6134	for (j = mm1; j <= i_2; ++j) {
6135	y = a[i + j * a_dim1];
6136	a[i + j * a_dim1] = a[m + j * a_dim1];
6137	a[m + j * a_dim1] = y;
6138	/* L110: */
6139	}
6140
6141	i_2 = *igh;
6142	for (j = 1; j <= i_2; ++j) {
6143	y = a[j + i * a_dim1];
6144	a[j + i * a_dim1] = a[j + m * a_dim1];
6145	a[j + m * a_dim1] = y;
6146	/* L120: */
6147	}
6148	/* .......... END INTERCHANGE .......... */
6149	L130:
6150	if (x == 0.) {
6151	goto L180;
6152	}
6153	mp1 = m + 1;
6154
6155	i_2 = *igh;
6156	for (i = mp1; i <= i_2; ++i) {
6157	y = a[i + mm1 * a_dim1];
6158	if (y == 0.) {
6159	goto L160;
6160	}
6161	y /= x;
6162	a[i + mm1 * a_dim1] = y;
6163
6164	i_3 = *n;
6165	for (j = m; j <= i_3; ++j) {
6166	/* L140: */
6167	a[i + j * a_dim1] -= y * a[m + j * a_dim1];
6168	}
6169
6170	i_3 = *igh;
6171	for (j = 1; j <= i_3; ++j) {
6172	/* L150: */
6173	a[j + m * a_dim1] += y * a[j + i * a_dim1];
6174	}
6175
6176	L160:
6177	;
6178	}
6179
6180	L180:
6181	;
6182	}
6183
6184	L200:
6185	return 0;
6186	} /* elmhes_ */
6187
6188	/* Subroutine / int eltran_(integer nm, integer n, integer low, integer *
6189	igh, doublereal a, integer int_, doublereal *z)
6190	{
6191	/* System generated locals */
6192	integer a_dim1, a_offset, z_dim1, z_offset, i_1, i_2;
6193
6194	/* Local variables */
6195	static integer i, j, kl, mm, mp, mp1;
6196
6197
6198
6199	/* THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ELMTRANS,
6200	*/
6201	/* NUM. MATH. 16, 181-204(1970) BY PETERS AND WILKINSON. */
6202	/* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971). */
6203
6204	/* THIS SUBROUTINE ACCUMULATES THE STABILIZED ELEMENTARY */
6205	/* SIMILARITY TRANSFORMATIONS USED IN THE REDUCTION OF A */
6206	/* REAL GENERAL MATRIX TO UPPER HESSENBERG FORM BY ELMHES. */
6207
6208	/* ON INPUT */
6209
6210	/* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */
6211	/* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */
6212	/* DIMENSION STATEMENT. */
6213
6214	/* N IS THE ORDER OF THE MATRIX. */
6215
6216	/* LOW AND IGH ARE INTEGERS D