Context Navigation

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.1.2.2591/ALGLIB-2.1.2.2591/tdbisinv.cs @ 3580

Visit:

Last change on this file since 3580 was 2645, checked in by mkommend, 15 years ago
extracted external libraries and adapted dependent plugins (ticket #837)
File size: 90.7 KB

Rev	Line
[2645]	1	/*************************************************************************
	2	Copyright (c) 1992-2007 The University of Tennessee. All rights reserved.
	3
	4	Contributors:
	5	* Sergey Bochkanov (ALGLIB project). Translation from FORTRAN to
	6	pseudocode.
	7
	8	See subroutines comments for additional copyrights.
	9
	10	>>> SOURCE LICENSE >>>
	11	This program is free software; you can redistribute it and/or modify
	12	it under the terms of the GNU General Public License as published by
	13	the Free Software Foundation (www.fsf.org); either version 2 of the
	14	License, or (at your option) any later version.
	15
	16	This program is distributed in the hope that it will be useful,
	17	but WITHOUT ANY WARRANTY; without even the implied warranty of
	18	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	19	GNU General Public License for more details.
	20
	21	A copy of the GNU General Public License is available at
	22	http://www.fsf.org/licensing/licenses
	23
	24	>>> END OF LICENSE >>>
	25	*************************************************************************/
	26
	27	using System;
	28
	29	namespace alglib
	30	{
	31	public class tdbisinv
	32	{
	33	/*************************************************************************
	34	Subroutine for finding the tridiagonal matrix eigenvalues/vectors in a
	35	given half-interval (A, B] by using bisection and inverse iteration.
	36
	37	Input parameters:
	38	D - the main diagonal of a tridiagonal matrix.
	39	Array whose index ranges within [0..N-1].
	40	E - the secondary diagonal of a tridiagonal matrix.
	41	Array whose index ranges within [0..N-2].
	42	N - size of matrix, N>=0.
	43	ZNeeded - flag controlling whether the eigenvectors are needed or not.
	44	If ZNeeded is equal to:
	45	* 0, the eigenvectors are not needed;
	46	* 1, the eigenvectors of a tridiagonal matrix are multiplied
	47	by the square matrix Z. It is used if the tridiagonal
	48	matrix is obtained by the similarity transformation
	49	of a symmetric matrix.
	50	* 2, the eigenvectors of a tridiagonal matrix replace matrix Z.
	51	A, B - half-interval (A, B] to search eigenvalues in.
	52	Z - if ZNeeded is equal to:
	53	* 0, Z isn't used and remains unchanged;
	54	* 1, Z contains the square matrix (array whose indexes range
	55	within [0..N-1, 0..N-1]) which reduces the given symmetric
	56	matrix to tridiagonal form;
	57	* 2, Z isn't used (but changed on the exit).
	58
	59	Output parameters:
	60	D - array of the eigenvalues found.
	61	Array whose index ranges within [0..M-1].
	62	M - number of eigenvalues found in the given half-interval (M>=0).
	63	Z - if ZNeeded is equal to:
	64	* 0, doesn't contain any information;
	65	* 1, contains the product of a given NxN matrix Z (from the
	66	left) and NxM matrix of the eigenvectors found (from the
	67	right). Array whose indexes range within [0..N-1, 0..M-1].
	68	* 2, contains the matrix of the eigenvectors found.
	69	Array whose indexes range within [0..N-1, 0..M-1].
	70
	71	Result:
	72
	73	True, if successful. In that case, M contains the number of eigenvalues
	74	in the given half-interval (could be equal to 0), D contains the eigenvalues,
	75	Z contains the eigenvectors (if needed).
	76	It should be noted that the subroutine changes the size of arrays D and Z.
	77
	78	False, if the bisection method subroutine wasn't able to find the
	79	eigenvalues in the given interval or if the inverse iteration subroutine
	80	wasn't able to find all the corresponding eigenvectors. In that case,
	81	the eigenvalues and eigenvectors are not returned, M is equal to 0.
	82
	83	-- ALGLIB --
	84	Copyright 31.03.2008 by Bochkanov Sergey
	85	*************************************************************************/
	86	public static bool smatrixtdevdr(ref double[] d,
	87	ref double[] e,
	88	int n,
	89	int zneeded,
	90	double a,
	91	double b,
	92	ref int m,
	93	ref double[,] z)
	94	{
	95	bool result = new bool();
	96	int errorcode = 0;
	97	int nsplit = 0;
	98	int i = 0;
	99	int j = 0;
	100	int k = 0;
	101	int cr = 0;
	102	int[] iblock = new int[0];
	103	int[] isplit = new int[0];
	104	int[] ifail = new int[0];
	105	double[] d1 = new double[0];
	106	double[] e1 = new double[0];
	107	double[] w = new double[0];
	108	double[,] z2 = new double[0,0];
	109	double[,] z3 = new double[0,0];
	110	double v = 0;
	111	int i_ = 0;
	112	int i1_ = 0;
	113
	114	System.Diagnostics.Debug.Assert(zneeded>=0 & zneeded<=2, "SMatrixTDEVDR: incorrect ZNeeded!");
	115
	116	//
	117	// Special cases
	118	//
	119	if( (double)(b)<=(double)(a) )
	120	{
	121	m = 0;
	122	result = true;
	123	return result;
	124	}
	125	if( n<=0 )
	126	{
	127	m = 0;
	128	result = true;
	129	return result;
	130	}
	131
	132	//
	133	// Copy D,E to D1, E1
	134	//
	135	d1 = new double[n+1];
	136	i1_ = (0) - (1);
	137	for(i_=1; i_<=n;i_++)
	138	{
	139	d1[i_] = d[i_+i1_];
	140	}
	141	if( n>1 )
	142	{
	143	e1 = new double[n-1+1];
	144	i1_ = (0) - (1);
	145	for(i_=1; i_<=n-1;i_++)
	146	{
	147	e1[i_] = e[i_+i1_];
	148	}
	149	}
	150
	151	//
	152	// No eigen vectors
	153	//
	154	if( zneeded==0 )
	155	{
	156	result = internalbisectioneigenvalues(d1, e1, n, 2, 1, a, b, 0, 0, -1, ref w, ref m, ref nsplit, ref iblock, ref isplit, ref errorcode);
	157	if( !result \| m==0 )
	158	{
	159	m = 0;
	160	return result;
	161	}
	162	d = new double[m-1+1];
	163	i1_ = (1) - (0);
	164	for(i_=0; i_<=m-1;i_++)
	165	{
	166	d[i_] = w[i_+i1_];
	167	}
	168	return result;
	169	}
	170
	171	//
	172	// Eigen vectors are multiplied by Z
	173	//
	174	if( zneeded==1 )
	175	{
	176
	177	//
	178	// Find eigen pairs
	179	//
	180	result = internalbisectioneigenvalues(d1, e1, n, 2, 2, a, b, 0, 0, -1, ref w, ref m, ref nsplit, ref iblock, ref isplit, ref errorcode);
	181	if( !result \| m==0 )
	182	{
	183	m = 0;
	184	return result;
	185	}
	186	internaldstein(n, ref d1, e1, m, w, ref iblock, ref isplit, ref z2, ref ifail, ref cr);
	187	if( cr!=0 )
	188	{
	189	m = 0;
	190	result = false;
	191	return result;
	192	}
	193
	194	//
	195	// Sort eigen values and vectors
	196	//
	197	for(i=1; i<=m; i++)
	198	{
	199	k = i;
	200	for(j=i; j<=m; j++)
	201	{
	202	if( (double)(w[j])<(double)(w[k]) )
	203	{
	204	k = j;
	205	}
	206	}
	207	v = w[i];
	208	w[i] = w[k];
	209	w[k] = v;
	210	for(j=1; j<=n; j++)
	211	{
	212	v = z2[j,i];
	213	z2[j,i] = z2[j,k];
	214	z2[j,k] = v;
	215	}
	216	}
	217
	218	//
	219	// Transform Z2 and overwrite Z
	220	//
	221	z3 = new double[m+1, n+1];
	222	for(i=1; i<=m; i++)
	223	{
	224	for(i_=1; i_<=n;i_++)
	225	{
	226	z3[i,i_] = z2[i_,i];
	227	}
	228	}
	229	for(i=1; i<=n; i++)
	230	{
	231	for(j=1; j<=m; j++)
	232	{
	233	i1_ = (1)-(0);
	234	v = 0.0;
	235	for(i_=0; i_<=n-1;i_++)
	236	{
	237	v += z[i-1,i_]*z3[j,i_+i1_];
	238	}
	239	z2[i,j] = v;
	240	}
	241	}
	242	z = new double[n-1+1, m-1+1];
	243	for(i=1; i<=m; i++)
	244	{
	245	i1_ = (1) - (0);
	246	for(i_=0; i_<=n-1;i_++)
	247	{
	248	z[i_,i-1] = z2[i_+i1_,i];
	249	}
	250	}
	251
	252	//
	253	// Store W
	254	//
	255	d = new double[m-1+1];
	256	for(i=1; i<=m; i++)
	257	{
	258	d[i-1] = w[i];
	259	}
	260	return result;
	261	}
	262
	263	//
	264	// Eigen vectors are stored in Z
	265	//
	266	if( zneeded==2 )
	267	{
	268
	269	//
	270	// Find eigen pairs
	271	//
	272	result = internalbisectioneigenvalues(d1, e1, n, 2, 2, a, b, 0, 0, -1, ref w, ref m, ref nsplit, ref iblock, ref isplit, ref errorcode);
	273	if( !result \| m==0 )
	274	{
	275	m = 0;
	276	return result;
	277	}
	278	internaldstein(n, ref d1, e1, m, w, ref iblock, ref isplit, ref z2, ref ifail, ref cr);
	279	if( cr!=0 )
	280	{
	281	m = 0;
	282	result = false;
	283	return result;
	284	}
	285
	286	//
	287	// Sort eigen values and vectors
	288	//
	289	for(i=1; i<=m; i++)
	290	{
	291	k = i;
	292	for(j=i; j<=m; j++)
	293	{
	294	if( (double)(w[j])<(double)(w[k]) )
	295	{
	296	k = j;
	297	}
	298	}
	299	v = w[i];
	300	w[i] = w[k];
	301	w[k] = v;
	302	for(j=1; j<=n; j++)
	303	{
	304	v = z2[j,i];
	305	z2[j,i] = z2[j,k];
	306	z2[j,k] = v;
	307	}
	308	}
	309
	310	//
	311	// Store W
	312	//
	313	d = new double[m-1+1];
	314	for(i=1; i<=m; i++)
	315	{
	316	d[i-1] = w[i];
	317	}
	318	z = new double[n-1+1, m-1+1];
	319	for(i=1; i<=m; i++)
	320	{
	321	i1_ = (1) - (0);
	322	for(i_=0; i_<=n-1;i_++)
	323	{
	324	z[i_,i-1] = z2[i_+i1_,i];
	325	}
	326	}
	327	return result;
	328	}
	329	result = false;
	330	return result;
	331	}
	332
	333
	334	/*************************************************************************
	335	Subroutine for finding tridiagonal matrix eigenvalues/vectors with given
	336	indexes (in ascending order) by using the bisection and inverse iteraion.
	337
	338	Input parameters:
	339	D - the main diagonal of a tridiagonal matrix.
	340	Array whose index ranges within [0..N-1].
	341	E - the secondary diagonal of a tridiagonal matrix.
	342	Array whose index ranges within [0..N-2].
	343	N - size of matrix. N>=0.
	344	ZNeeded - flag controlling whether the eigenvectors are needed or not.
	345	If ZNeeded is equal to:
	346	* 0, the eigenvectors are not needed;
	347	* 1, the eigenvectors of a tridiagonal matrix are multiplied
	348	by the square matrix Z. It is used if the
	349	tridiagonal matrix is obtained by the similarity transformation
	350	of a symmetric matrix.
	351	* 2, the eigenvectors of a tridiagonal matrix replace
	352	matrix Z.
	353	I1, I2 - index interval for searching (from I1 to I2).
	354	0 <= I1 <= I2 <= N-1.
	355	Z - if ZNeeded is equal to:
	356	* 0, Z isn't used and remains unchanged;
	357	* 1, Z contains the square matrix (array whose indexes range within [0..N-1, 0..N-1])
	358	which reduces the given symmetric matrix to tridiagonal form;
	359	* 2, Z isn't used (but changed on the exit).
	360
	361	Output parameters:
	362	D - array of the eigenvalues found.
	363	Array whose index ranges within [0..I2-I1].
	364	Z - if ZNeeded is equal to:
	365	* 0, doesn't contain any information;
	366	* 1, contains the product of a given NxN matrix Z (from the left) and
	367	Nx(I2-I1) matrix of the eigenvectors found (from the right).
	368	Array whose indexes range within [0..N-1, 0..I2-I1].
	369	* 2, contains the matrix of the eigenvalues found.
	370	Array whose indexes range within [0..N-1, 0..I2-I1].
	371
	372
	373	Result:
	374
	375	True, if successful. In that case, D contains the eigenvalues,
	376	Z contains the eigenvectors (if needed).
	377	It should be noted that the subroutine changes the size of arrays D and Z.
	378
	379	False, if the bisection method subroutine wasn't able to find the eigenvalues
	380	in the given interval or if the inverse iteration subroutine wasn't able
	381	to find all the corresponding eigenvectors. In that case, the eigenvalues
	382	and eigenvectors are not returned.
	383
	384	-- ALGLIB --
	385	Copyright 25.12.2005 by Bochkanov Sergey
	386	*************************************************************************/
	387	public static bool smatrixtdevdi(ref double[] d,
	388	ref double[] e,
	389	int n,
	390	int zneeded,
	391	int i1,
	392	int i2,
	393	ref double[,] z)
	394	{
	395	bool result = new bool();
	396	int errorcode = 0;
	397	int nsplit = 0;
	398	int i = 0;
	399	int j = 0;
	400	int k = 0;
	401	int m = 0;
	402	int cr = 0;
	403	int[] iblock = new int[0];
	404	int[] isplit = new int[0];
	405	int[] ifail = new int[0];
	406	double[] w = new double[0];
	407	double[] d1 = new double[0];
	408	double[] e1 = new double[0];
	409	double[,] z2 = new double[0,0];
	410	double[,] z3 = new double[0,0];
	411	double v = 0;
	412	int i_ = 0;
	413	int i1_ = 0;
	414
	415	System.Diagnostics.Debug.Assert(0<=i1 & i1<=i2 & i2<n, "SMatrixTDEVDI: incorrect I1/I2!");
	416
	417	//
	418	// Copy D,E to D1, E1
	419	//
	420	d1 = new double[n+1];
	421	i1_ = (0) - (1);
	422	for(i_=1; i_<=n;i_++)
	423	{
	424	d1[i_] = d[i_+i1_];
	425	}
	426	if( n>1 )
	427	{
	428	e1 = new double[n-1+1];
	429	i1_ = (0) - (1);
	430	for(i_=1; i_<=n-1;i_++)
	431	{
	432	e1[i_] = e[i_+i1_];
	433	}
	434	}
	435
	436	//
	437	// No eigen vectors
	438	//
	439	if( zneeded==0 )
	440	{
	441	result = internalbisectioneigenvalues(d1, e1, n, 3, 1, 0, 0, i1+1, i2+1, -1, ref w, ref m, ref nsplit, ref iblock, ref isplit, ref errorcode);
	442	if( !result )
	443	{
	444	return result;
	445	}
	446	if( m!=i2-i1+1 )
	447	{
	448	result = false;
	449	return result;
	450	}
	451	d = new double[m-1+1];
	452	for(i=1; i<=m; i++)
	453	{
	454	d[i-1] = w[i];
	455	}
	456	return result;
	457	}
	458
	459	//
	460	// Eigen vectors are multiplied by Z
	461	//
	462	if( zneeded==1 )
	463	{
	464
	465	//
	466	// Find eigen pairs
	467	//
	468	result = internalbisectioneigenvalues(d1, e1, n, 3, 2, 0, 0, i1+1, i2+1, -1, ref w, ref m, ref nsplit, ref iblock, ref isplit, ref errorcode);
	469	if( !result )
	470	{
	471	return result;
	472	}
	473	if( m!=i2-i1+1 )
	474	{
	475	result = false;
	476	return result;
	477	}
	478	internaldstein(n, ref d1, e1, m, w, ref iblock, ref isplit, ref z2, ref ifail, ref cr);
	479	if( cr!=0 )
	480	{
	481	result = false;
	482	return result;
	483	}
	484
	485	//
	486	// Sort eigen values and vectors
	487	//
	488	for(i=1; i<=m; i++)
	489	{
	490	k = i;
	491	for(j=i; j<=m; j++)
	492	{
	493	if( (double)(w[j])<(double)(w[k]) )
	494	{
	495	k = j;
	496	}
	497	}
	498	v = w[i];
	499	w[i] = w[k];
	500	w[k] = v;
	501	for(j=1; j<=n; j++)
	502	{
	503	v = z2[j,i];
	504	z2[j,i] = z2[j,k];
	505	z2[j,k] = v;
	506	}
	507	}
	508
	509	//
	510	// Transform Z2 and overwrite Z
	511	//
	512	z3 = new double[m+1, n+1];
	513	for(i=1; i<=m; i++)
	514	{
	515	for(i_=1; i_<=n;i_++)
	516	{
	517	z3[i,i_] = z2[i_,i];
	518	}
	519	}
	520	for(i=1; i<=n; i++)
	521	{
	522	for(j=1; j<=m; j++)
	523	{
	524	i1_ = (1)-(0);
	525	v = 0.0;
	526	for(i_=0; i_<=n-1;i_++)
	527	{
	528	v += z[i-1,i_]*z3[j,i_+i1_];
	529	}
	530	z2[i,j] = v;
	531	}
	532	}
	533	z = new double[n-1+1, m-1+1];
	534	for(i=1; i<=m; i++)
	535	{
	536	i1_ = (1) - (0);
	537	for(i_=0; i_<=n-1;i_++)
	538	{
	539	z[i_,i-1] = z2[i_+i1_,i];
	540	}
	541	}
	542
	543	//
	544	// Store W
	545	//
	546	d = new double[m-1+1];
	547	for(i=1; i<=m; i++)
	548	{
	549	d[i-1] = w[i];
	550	}
	551	return result;
	552	}
	553
	554	//
	555	// Eigen vectors are stored in Z
	556	//
	557	if( zneeded==2 )
	558	{
	559
	560	//
	561	// Find eigen pairs
	562	//
	563	result = internalbisectioneigenvalues(d1, e1, n, 3, 2, 0, 0, i1+1, i2+1, -1, ref w, ref m, ref nsplit, ref iblock, ref isplit, ref errorcode);
	564	if( !result )
	565	{
	566	return result;
	567	}
	568	if( m!=i2-i1+1 )
	569	{
	570	result = false;
	571	return result;
	572	}
	573	internaldstein(n, ref d1, e1, m, w, ref iblock, ref isplit, ref z2, ref ifail, ref cr);
	574	if( cr!=0 )
	575	{
	576	result = false;
	577	return result;
	578	}
	579
	580	//
	581	// Sort eigen values and vectors
	582	//
	583	for(i=1; i<=m; i++)
	584	{
	585	k = i;
	586	for(j=i; j<=m; j++)
	587	{
	588	if( (double)(w[j])<(double)(w[k]) )
	589	{
	590	k = j;
	591	}
	592	}
	593	v = w[i];
	594	w[i] = w[k];
	595	w[k] = v;
	596	for(j=1; j<=n; j++)
	597	{
	598	v = z2[j,i];
	599	z2[j,i] = z2[j,k];
	600	z2[j,k] = v;
	601	}
	602	}
	603
	604	//
	605	// Store Z
	606	//
	607	z = new double[n-1+1, m-1+1];
	608	for(i=1; i<=m; i++)
	609	{
	610	i1_ = (1) - (0);
	611	for(i_=0; i_<=n-1;i_++)
	612	{
	613	z[i_,i-1] = z2[i_+i1_,i];
	614	}
	615	}
	616
	617	//
	618	// Store W
	619	//
	620	d = new double[m-1+1];
	621	for(i=1; i<=m; i++)
	622	{
	623	d[i-1] = w[i];
	624	}
	625	return result;
	626	}
	627	result = false;
	628	return result;
	629	}
	630
	631
	632	public static bool tridiagonaleigenvaluesandvectorsininterval(ref double[] d,
	633	ref double[] e,
	634	int n,
	635	int zneeded,
	636	double a,
	637	double b,
	638	ref int m,
	639	ref double[,] z)
	640	{
	641	bool result = new bool();
	642	int errorcode = 0;
	643	int nsplit = 0;
	644	int i = 0;
	645	int j = 0;
	646	int k = 0;
	647	int cr = 0;
	648	int[] iblock = new int[0];
	649	int[] isplit = new int[0];
	650	int[] ifail = new int[0];
	651	double[] w = new double[0];
	652	double[,] z2 = new double[0,0];
	653	double[,] z3 = new double[0,0];
	654	double v = 0;
	655	int i_ = 0;
	656
	657
	658	//
	659	// No eigen vectors
	660	//
	661	if( zneeded==0 )
	662	{
	663	result = internalbisectioneigenvalues(d, e, n, 2, 1, a, b, 0, 0, -1, ref w, ref m, ref nsplit, ref iblock, ref isplit, ref errorcode);
	664	if( !result \| m==0 )
	665	{
	666	m = 0;
	667	return result;
	668	}
	669	d = new double[m+1];
	670	for(i=1; i<=m; i++)
	671	{
	672	d[i] = w[i];
	673	}
	674	return result;
	675	}
	676
	677	//
	678	// Eigen vectors are multiplied by Z
	679	//
	680	if( zneeded==1 )
	681	{
	682
	683	//
	684	// Find eigen pairs
	685	//
	686	result = internalbisectioneigenvalues(d, e, n, 2, 2, a, b, 0, 0, -1, ref w, ref m, ref nsplit, ref iblock, ref isplit, ref errorcode);
	687	if( !result \| m==0 )
	688	{
	689	m = 0;
	690	return result;
	691	}
	692	internaldstein(n, ref d, e, m, w, ref iblock, ref isplit, ref z2, ref ifail, ref cr);
	693	if( cr!=0 )
	694	{
	695	m = 0;
	696	result = false;
	697	return result;
	698	}
	699
	700	//
	701	// Sort eigen values and vectors
	702	//
	703	for(i=1; i<=m; i++)
	704	{
	705	k = i;
	706	for(j=i; j<=m; j++)
	707	{
	708	if( (double)(w[j])<(double)(w[k]) )
	709	{
	710	k = j;
	711	}
	712	}
	713	v = w[i];
	714	w[i] = w[k];
	715	w[k] = v;
	716	for(j=1; j<=n; j++)
	717	{
	718	v = z2[j,i];
	719	z2[j,i] = z2[j,k];
	720	z2[j,k] = v;
	721	}
	722	}
	723
	724	//
	725	// Transform Z2 and overwrite Z
	726	//
	727	z3 = new double[m+1, n+1];
	728	for(i=1; i<=m; i++)
	729	{
	730	for(i_=1; i_<=n;i_++)
	731	{
	732	z3[i,i_] = z2[i_,i];
	733	}
	734	}
	735	for(i=1; i<=n; i++)
	736	{
	737	for(j=1; j<=m; j++)
	738	{
	739	v = 0.0;
	740	for(i_=1; i_<=n;i_++)
	741	{
	742	v += z[i,i_]*z3[j,i_];
	743	}
	744	z2[i,j] = v;
	745	}
	746	}
	747	z = new double[n+1, m+1];
	748	for(i=1; i<=m; i++)
	749	{
	750	for(i_=1; i_<=n;i_++)
	751	{
	752	z[i_,i] = z2[i_,i];
	753	}
	754	}
	755
	756	//
	757	// Store W
	758	//
	759	d = new double[m+1];
	760	for(i=1; i<=m; i++)
	761	{
	762	d[i] = w[i];
	763	}
	764	return result;
	765	}
	766
	767	//
	768	// Eigen vectors are stored in Z
	769	//
	770	if( zneeded==2 )
	771	{
	772
	773	//
	774	// Find eigen pairs
	775	//
	776	result = internalbisectioneigenvalues(d, e, n, 2, 2, a, b, 0, 0, -1, ref w, ref m, ref nsplit, ref iblock, ref isplit, ref errorcode);
	777	if( !result \| m==0 )
	778	{
	779	m = 0;
	780	return result;
	781	}
	782	internaldstein(n, ref d, e, m, w, ref iblock, ref isplit, ref z, ref ifail, ref cr);
	783	if( cr!=0 )
	784	{
	785	m = 0;
	786	result = false;
	787	return result;
	788	}
	789
	790	//
	791	// Sort eigen values and vectors
	792	//
	793	for(i=1; i<=m; i++)
	794	{
	795	k = i;
	796	for(j=i; j<=m; j++)
	797	{
	798	if( (double)(w[j])<(double)(w[k]) )
	799	{
	800	k = j;
	801	}
	802	}
	803	v = w[i];
	804	w[i] = w[k];
	805	w[k] = v;
	806	for(j=1; j<=n; j++)
	807	{
	808	v = z[j,i];
	809	z[j,i] = z[j,k];
	810	z[j,k] = v;
	811	}
	812	}
	813
	814	//
	815	// Store W
	816	//
	817	d = new double[m+1];
	818	for(i=1; i<=m; i++)
	819	{
	820	d[i] = w[i];
	821	}
	822	return result;
	823	}
	824	result = false;
	825	return result;
	826	}
	827
	828
	829	public static bool tridiagonaleigenvaluesandvectorsbyindexes(ref double[] d,
	830	ref double[] e,
	831	int n,
	832	int zneeded,
	833	int i1,
	834	int i2,
	835	ref double[,] z)
	836	{
	837	bool result = new bool();
	838	int errorcode = 0;
	839	int nsplit = 0;
	840	int i = 0;
	841	int j = 0;
	842	int k = 0;
	843	int m = 0;
	844	int cr = 0;
	845	int[] iblock = new int[0];
	846	int[] isplit = new int[0];
	847	int[] ifail = new int[0];
	848	double[] w = new double[0];
	849	double[,] z2 = new double[0,0];
	850	double[,] z3 = new double[0,0];
	851	double v = 0;
	852	int i_ = 0;
	853
	854
	855	//
	856	// No eigen vectors
	857	//
	858	if( zneeded==0 )
	859	{
	860	result = internalbisectioneigenvalues(d, e, n, 3, 1, 0, 0, i1, i2, -1, ref w, ref m, ref nsplit, ref iblock, ref isplit, ref errorcode);
	861	if( !result )
	862	{
	863	return result;
	864	}
	865	d = new double[m+1];
	866	for(i=1; i<=m; i++)
	867	{
	868	d[i] = w[i];
	869	}
	870	return result;
	871	}
	872
	873	//
	874	// Eigen vectors are multiplied by Z
	875	//
	876	if( zneeded==1 )
	877	{
	878
	879	//
	880	// Find eigen pairs
	881	//
	882	result = internalbisectioneigenvalues(d, e, n, 3, 2, 0, 0, i1, i2, -1, ref w, ref m, ref nsplit, ref iblock, ref isplit, ref errorcode);
	883	if( !result )
	884	{
	885	return result;
	886	}
	887	internaldstein(n, ref d, e, m, w, ref iblock, ref isplit, ref z2, ref ifail, ref cr);
	888	if( cr!=0 )
	889	{
	890	result = false;
	891	return result;
	892	}
	893
	894	//
	895	// Sort eigen values and vectors
	896	//
	897	for(i=1; i<=m; i++)
	898	{
	899	k = i;
	900	for(j=i; j<=m; j++)
	901	{
	902	if( (double)(w[j])<(double)(w[k]) )
	903	{
	904	k = j;
	905	}
	906	}
	907	v = w[i];
	908	w[i] = w[k];
	909	w[k] = v;
	910	for(j=1; j<=n; j++)
	911	{
	912	v = z2[j,i];
	913	z2[j,i] = z2[j,k];
	914	z2[j,k] = v;
	915	}
	916	}
	917
	918	//
	919	// Transform Z2 and overwrite Z
	920	//
	921	z3 = new double[m+1, n+1];
	922	for(i=1; i<=m; i++)
	923	{
	924	for(i_=1; i_<=n;i_++)
	925	{
	926	z3[i,i_] = z2[i_,i];
	927	}
	928	}
	929	for(i=1; i<=n; i++)
	930	{
	931	for(j=1; j<=m; j++)
	932	{
	933	v = 0.0;
	934	for(i_=1; i_<=n;i_++)
	935	{
	936	v += z[i,i_]*z3[j,i_];
	937	}
	938	z2[i,j] = v;
	939	}
	940	}
	941	z = new double[n+1, m+1];
	942	for(i=1; i<=m; i++)
	943	{
	944	for(i_=1; i_<=n;i_++)
	945	{
	946	z[i_,i] = z2[i_,i];
	947	}
	948	}
	949
	950	//
	951	// Store W
	952	//
	953	d = new double[m+1];
	954	for(i=1; i<=m; i++)
	955	{
	956	d[i] = w[i];
	957	}
	958	return result;
	959	}
	960
	961	//
	962	// Eigen vectors are stored in Z
	963	//
	964	if( zneeded==2 )
	965	{
	966
	967	//
	968	// Find eigen pairs
	969	//
	970	result = internalbisectioneigenvalues(d, e, n, 3, 2, 0, 0, i1, i2, -1, ref w, ref m, ref nsplit, ref iblock, ref isplit, ref errorcode);
	971	if( !result )
	972	{
	973	return result;
	974	}
	975	internaldstein(n, ref d, e, m, w, ref iblock, ref isplit, ref z, ref ifail, ref cr);
	976	if( cr!=0 )
	977	{
	978	result = false;
	979	return result;
	980	}
	981
	982	//
	983	// Sort eigen values and vectors
	984	//
	985	for(i=1; i<=m; i++)
	986	{
	987	k = i;
	988	for(j=i; j<=m; j++)
	989	{
	990	if( (double)(w[j])<(double)(w[k]) )
	991	{
	992	k = j;
	993	}
	994	}
	995	v = w[i];
	996	w[i] = w[k];
	997	w[k] = v;
	998	for(j=1; j<=n; j++)
	999	{
	1000	v = z[j,i];
	1001	z[j,i] = z[j,k];
	1002	z[j,k] = v;
	1003	}
	1004	}
	1005
	1006	//
	1007	// Store W
	1008	//
	1009	d = new double[m+1];
	1010	for(i=1; i<=m; i++)
	1011	{
	1012	d[i] = w[i];
	1013	}
	1014	return result;
	1015	}
	1016	result = false;
	1017	return result;
	1018	}
	1019
	1020
	1021	public static bool internalbisectioneigenvalues(double[] d,
	1022	double[] e,
	1023	int n,
	1024	int irange,
	1025	int iorder,
	1026	double vl,
	1027	double vu,
	1028	int il,
	1029	int iu,
	1030	double abstol,
	1031	ref double[] w,
	1032	ref int m,
	1033	ref int nsplit,
	1034	ref int[] iblock,
	1035	ref int[] isplit,
	1036	ref int errorcode)
	1037	{
	1038	bool result = new bool();
	1039	double fudge = 0;
	1040	double relfac = 0;
	1041	bool ncnvrg = new bool();
	1042	bool toofew = new bool();
	1043	int ib = 0;
	1044	int ibegin = 0;
	1045	int idiscl = 0;
	1046	int idiscu = 0;
	1047	int ie = 0;
	1048	int iend = 0;
	1049	int iinfo = 0;
	1050	int im = 0;
	1051	int iin = 0;
	1052	int ioff = 0;
	1053	int iout = 0;
	1054	int itmax = 0;
	1055	int iw = 0;
	1056	int iwoff = 0;
	1057	int j = 0;
	1058	int itmp1 = 0;
	1059	int jb = 0;
	1060	int jdisc = 0;
	1061	int je = 0;
	1062	int nwl = 0;
	1063	int nwu = 0;
	1064	double atoli = 0;
	1065	double bnorm = 0;
	1066	double gl = 0;
	1067	double gu = 0;
	1068	double pivmin = 0;
	1069	double rtoli = 0;
	1070	double safemn = 0;
	1071	double tmp1 = 0;
	1072	double tmp2 = 0;
	1073	double tnorm = 0;
	1074	double ulp = 0;
	1075	double wkill = 0;
	1076	double wl = 0;
	1077	double wlu = 0;
	1078	double wu = 0;
	1079	double wul = 0;
	1080	double scalefactor = 0;
	1081	double t = 0;
	1082	int[] idumma = new int[0];
	1083	double[] work = new double[0];
	1084	int[] iwork = new int[0];
	1085	int[] ia1s2 = new int[0];
	1086	double[] ra1s2 = new double[0];
	1087	double[,] ra1s2x2 = new double[0,0];
	1088	int[,] ia1s2x2 = new int[0,0];
	1089	double[] ra1siin = new double[0];
	1090	double[] ra2siin = new double[0];
	1091	double[] ra3siin = new double[0];
	1092	double[] ra4siin = new double[0];
	1093	double[,] ra1siinx2 = new double[0,0];
	1094	int[,] ia1siinx2 = new int[0,0];
	1095	int[] iworkspace = new int[0];
	1096	double[] rworkspace = new double[0];
	1097	int tmpi = 0;
	1098
	1099	d = (double[])d.Clone();
	1100	e = (double[])e.Clone();
	1101
	1102
	1103	//
	1104	// Quick return if possible
	1105	//
	1106	m = 0;
	1107	if( n==0 )
	1108	{
	1109	result = true;
	1110	return result;
	1111	}
	1112
	1113	//
	1114	// Get machine constants
	1115	// NB is the minimum vector length for vector bisection, or 0
	1116	// if only scalar is to be done.
	1117	//
	1118	fudge = 2;
	1119	relfac = 2;
	1120	safemn = AP.Math.MinRealNumber;
	1121	ulp = 2*AP.Math.MachineEpsilon;
	1122	rtoli = ulp*relfac;
	1123	idumma = new int[1+1];
	1124	work = new double[4*n+1];
	1125	iwork = new int[3*n+1];
	1126	w = new double[n+1];
	1127	iblock = new int[n+1];
	1128	isplit = new int[n+1];
	1129	ia1s2 = new int[2+1];
	1130	ra1s2 = new double[2+1];
	1131	ra1s2x2 = new double[2+1, 2+1];
	1132	ia1s2x2 = new int[2+1, 2+1];
	1133	ra1siin = new double[n+1];
	1134	ra2siin = new double[n+1];
	1135	ra3siin = new double[n+1];
	1136	ra4siin = new double[n+1];
	1137	ra1siinx2 = new double[n+1, 2+1];
	1138	ia1siinx2 = new int[n+1, 2+1];
	1139	iworkspace = new int[n+1];
	1140	rworkspace = new double[n+1];
	1141
	1142	//
	1143	// Check for Errors
	1144	//
	1145	result = false;
	1146	errorcode = 0;
	1147	if( irange<=0 \| irange>=4 )
	1148	{
	1149	errorcode = -4;
	1150	}
	1151	if( iorder<=0 \| iorder>=3 )
	1152	{
	1153	errorcode = -5;
	1154	}
	1155	if( n<0 )
	1156	{
	1157	errorcode = -3;
	1158	}
	1159	if( irange==2 & (double)(vl)>=(double)(vu) )
	1160	{
	1161	errorcode = -6;
	1162	}
	1163	if( irange==3 & (il<1 \| il>Math.Max(1, n)) )
	1164	{
	1165	errorcode = -8;
	1166	}
	1167	if( irange==3 & (iu<Math.Min(n, il) \| iu>n) )
	1168	{
	1169	errorcode = -9;
	1170	}
	1171	if( errorcode!=0 )
	1172	{
	1173	return result;
	1174	}
	1175
	1176	//
	1177	// Initialize error flags
	1178	//
	1179	ncnvrg = false;
	1180	toofew = false;
	1181
	1182	//
	1183	// Simplifications:
	1184	//
	1185	if( irange==3 & il==1 & iu==n )
	1186	{
	1187	irange = 1;
	1188	}
	1189
	1190	//
	1191	// Special Case when N=1
	1192	//
	1193	if( n==1 )
	1194	{
	1195	nsplit = 1;
	1196	isplit[1] = 1;
	1197	if( irange==2 & ((double)(vl)>=(double)(d[1]) \| (double)(vu)<(double)(d[1])) )
	1198	{
	1199	m = 0;
	1200	}
	1201	else
	1202	{
	1203	w[1] = d[1];
	1204	iblock[1] = 1;
	1205	m = 1;
	1206	}
	1207	result = true;
	1208	return result;
	1209	}
	1210
	1211	//
	1212	// Scaling
	1213	//
	1214	t = Math.Abs(d[n]);
	1215	for(j=1; j<=n-1; j++)
	1216	{
	1217	t = Math.Max(t, Math.Abs(d[j]));
	1218	t = Math.Max(t, Math.Abs(e[j]));
	1219	}
	1220	scalefactor = 1;
	1221	if( (double)(t)!=(double)(0) )
	1222	{
	1223	if( (double)(t)>(double)(Math.Sqrt(Math.Sqrt(AP.Math.MinRealNumber))*Math.Sqrt(AP.Math.MaxRealNumber)) )
	1224	{
	1225	scalefactor = t;
	1226	}
	1227	if( (double)(t)<(double)(Math.Sqrt(Math.Sqrt(AP.Math.MaxRealNumber))*Math.Sqrt(AP.Math.MinRealNumber)) )
	1228	{
	1229	scalefactor = t;
	1230	}
	1231	for(j=1; j<=n-1; j++)
	1232	{
	1233	d[j] = d[j]/scalefactor;
	1234	e[j] = e[j]/scalefactor;
	1235	}
	1236	d[n] = d[n]/scalefactor;
	1237	}
	1238
	1239	//
	1240	// Compute Splitting Points
	1241	//
	1242	nsplit = 1;
	1243	work[n] = 0;
	1244	pivmin = 1;
	1245	for(j=2; j<=n; j++)
	1246	{
	1247	tmp1 = AP.Math.Sqr(e[j-1]);
	1248	if( (double)(Math.Abs(d[j]d[j-1])AP.Math.Sqr(ulp)+safemn)>(double)(tmp1) )
	1249	{
	1250	isplit[nsplit] = j-1;
	1251	nsplit = nsplit+1;
	1252	work[j-1] = 0;
	1253	}
	1254	else
	1255	{
	1256	work[j-1] = tmp1;
	1257	pivmin = Math.Max(pivmin, tmp1);
	1258	}
	1259	}
	1260	isplit[nsplit] = n;
	1261	pivmin = pivmin*safemn;
	1262
	1263	//
	1264	// Compute Interval and ATOLI
	1265	//
	1266	if( irange==3 )
	1267	{
	1268
	1269	//
	1270	// RANGE='I': Compute the interval containing eigenvalues
	1271	// IL through IU.
	1272	//
	1273	// Compute Gershgorin interval for entire (split) matrix
	1274	// and use it as the initial interval
	1275	//
	1276	gu = d[1];
	1277	gl = d[1];
	1278	tmp1 = 0;
	1279	for(j=1; j<=n-1; j++)
	1280	{
	1281	tmp2 = Math.Sqrt(work[j]);
	1282	gu = Math.Max(gu, d[j]+tmp1+tmp2);
	1283	gl = Math.Min(gl, d[j]-tmp1-tmp2);
	1284	tmp1 = tmp2;
	1285	}
	1286	gu = Math.Max(gu, d[n]+tmp1);
	1287	gl = Math.Min(gl, d[n]-tmp1);
	1288	tnorm = Math.Max(Math.Abs(gl), Math.Abs(gu));
	1289	gl = gl-fudgetnormulpn-fudge2*pivmin;
	1290	gu = gu+fudgetnormulpn+fudgepivmin;
	1291
	1292	//
	1293	// Compute Iteration parameters
	1294	//
	1295	itmax = (int)Math.Ceiling((Math.Log(tnorm+pivmin)-Math.Log(pivmin))/Math.Log(2))+2;
	1296	if( (double)(abstol)<=(double)(0) )
	1297	{
	1298	atoli = ulp*tnorm;
	1299	}
	1300	else
	1301	{
	1302	atoli = abstol;
	1303	}
	1304	work[n+1] = gl;
	1305	work[n+2] = gl;
	1306	work[n+3] = gu;
	1307	work[n+4] = gu;
	1308	work[n+5] = gl;
	1309	work[n+6] = gu;
	1310	iwork[1] = -1;
	1311	iwork[2] = -1;
	1312	iwork[3] = n+1;
	1313	iwork[4] = n+1;
	1314	iwork[5] = il-1;
	1315	iwork[6] = iu;
	1316
	1317	//
	1318	// Calling DLAEBZ
	1319	//
	1320	// DLAEBZ( 3, ITMAX, N, 2, 2, NB, ATOLI, RTOLI, PIVMIN, D, E,
	1321	// WORK, IWORK( 5 ), WORK( N+1 ), WORK( N+5 ), IOUT,
	1322	// IWORK, W, IBLOCK, IINFO )
	1323	//
	1324	ia1s2[1] = iwork[5];
	1325	ia1s2[2] = iwork[6];
	1326	ra1s2[1] = work[n+5];
	1327	ra1s2[2] = work[n+6];
	1328	ra1s2x2[1,1] = work[n+1];
	1329	ra1s2x2[2,1] = work[n+2];
	1330	ra1s2x2[1,2] = work[n+3];
	1331	ra1s2x2[2,2] = work[n+4];
	1332	ia1s2x2[1,1] = iwork[1];
	1333	ia1s2x2[2,1] = iwork[2];
	1334	ia1s2x2[1,2] = iwork[3];
	1335	ia1s2x2[2,2] = iwork[4];
	1336	internaldlaebz(3, itmax, n, 2, 2, atoli, rtoli, pivmin, ref d, ref e, ref work, ref ia1s2, ref ra1s2x2, ref ra1s2, ref iout, ref ia1s2x2, ref w, ref iblock, ref iinfo);
	1337	iwork[5] = ia1s2[1];
	1338	iwork[6] = ia1s2[2];
	1339	work[n+5] = ra1s2[1];
	1340	work[n+6] = ra1s2[2];
	1341	work[n+1] = ra1s2x2[1,1];
	1342	work[n+2] = ra1s2x2[2,1];
	1343	work[n+3] = ra1s2x2[1,2];
	1344	work[n+4] = ra1s2x2[2,2];
	1345	iwork[1] = ia1s2x2[1,1];
	1346	iwork[2] = ia1s2x2[2,1];
	1347	iwork[3] = ia1s2x2[1,2];
	1348	iwork[4] = ia1s2x2[2,2];
	1349	if( iwork[6]==iu )
	1350	{
	1351	wl = work[n+1];
	1352	wlu = work[n+3];
	1353	nwl = iwork[1];
	1354	wu = work[n+4];
	1355	wul = work[n+2];
	1356	nwu = iwork[4];
	1357	}
	1358	else
	1359	{
	1360	wl = work[n+2];
	1361	wlu = work[n+4];
	1362	nwl = iwork[2];
	1363	wu = work[n+3];
	1364	wul = work[n+1];
	1365	nwu = iwork[3];
	1366	}
	1367	if( nwl<0 \| nwl>=n \| nwu<1 \| nwu>n )
	1368	{
	1369	errorcode = 4;
	1370	result = false;
	1371	return result;
	1372	}
	1373	}
	1374	else
	1375	{
	1376
	1377	//
	1378	// RANGE='A' or 'V' -- Set ATOLI
	1379	//
	1380	tnorm = Math.Max(Math.Abs(d[1])+Math.Abs(e[1]), Math.Abs(d[n])+Math.Abs(e[n-1]));
	1381	for(j=2; j<=n-1; j++)
	1382	{
	1383	tnorm = Math.Max(tnorm, Math.Abs(d[j])+Math.Abs(e[j-1])+Math.Abs(e[j]));
	1384	}
	1385	if( (double)(abstol)<=(double)(0) )
	1386	{
	1387	atoli = ulp*tnorm;
	1388	}
	1389	else
	1390	{
	1391	atoli = abstol;
	1392	}
	1393	if( irange==2 )
	1394	{
	1395	wl = vl;
	1396	wu = vu;
	1397	}
	1398	else
	1399	{
	1400	wl = 0;
	1401	wu = 0;
	1402	}
	1403	}
	1404
	1405	//
	1406	// Find Eigenvalues -- Loop Over Blocks and recompute NWL and NWU.
	1407	// NWL accumulates the number of eigenvalues .le. WL,
	1408	// NWU accumulates the number of eigenvalues .le. WU
	1409	//
	1410	m = 0;
	1411	iend = 0;
	1412	errorcode = 0;
	1413	nwl = 0;
	1414	nwu = 0;
	1415	for(jb=1; jb<=nsplit; jb++)
	1416	{
	1417	ioff = iend;
	1418	ibegin = ioff+1;
	1419	iend = isplit[jb];
	1420	iin = iend-ioff;
	1421	if( iin==1 )
	1422	{
	1423
	1424	//
	1425	// Special Case -- IIN=1
	1426	//
	1427	if( irange==1 \| (double)(wl)>=(double)(d[ibegin]-pivmin) )
	1428	{
	1429	nwl = nwl+1;
	1430	}
	1431	if( irange==1 \| (double)(wu)>=(double)(d[ibegin]-pivmin) )
	1432	{
	1433	nwu = nwu+1;
	1434	}
	1435	if( irange==1 \| (double)(wl)<(double)(d[ibegin]-pivmin) & (double)(wu)>=(double)(d[ibegin]-pivmin) )
	1436	{
	1437	m = m+1;
	1438	w[m] = d[ibegin];
	1439	iblock[m] = jb;
	1440	}
	1441	}
	1442	else
	1443	{
	1444
	1445	//
	1446	// General Case -- IIN > 1
	1447	//
	1448	// Compute Gershgorin Interval
	1449	// and use it as the initial interval
	1450	//
	1451	gu = d[ibegin];
	1452	gl = d[ibegin];
	1453	tmp1 = 0;
	1454	for(j=ibegin; j<=iend-1; j++)
	1455	{
	1456	tmp2 = Math.Abs(e[j]);
	1457	gu = Math.Max(gu, d[j]+tmp1+tmp2);
	1458	gl = Math.Min(gl, d[j]-tmp1-tmp2);
	1459	tmp1 = tmp2;
	1460	}
	1461	gu = Math.Max(gu, d[iend]+tmp1);
	1462	gl = Math.Min(gl, d[iend]-tmp1);
	1463	bnorm = Math.Max(Math.Abs(gl), Math.Abs(gu));
	1464	gl = gl-fudgebnormulpiin-fudgepivmin;
	1465	gu = gu+fudgebnormulpiin+fudgepivmin;
	1466
	1467	//
	1468	// Compute ATOLI for the current submatrix
	1469	//
	1470	if( (double)(abstol)<=(double)(0) )
	1471	{
	1472	atoli = ulp*Math.Max(Math.Abs(gl), Math.Abs(gu));
	1473	}
	1474	else
	1475	{
	1476	atoli = abstol;
	1477	}
	1478	if( irange>1 )
	1479	{
	1480	if( (double)(gu)<(double)(wl) )
	1481	{
	1482	nwl = nwl+iin;
	1483	nwu = nwu+iin;
	1484	continue;
	1485	}
	1486	gl = Math.Max(gl, wl);
	1487	gu = Math.Min(gu, wu);
	1488	if( (double)(gl)>=(double)(gu) )
	1489	{
	1490	continue;
	1491	}
	1492	}
	1493
	1494	//
	1495	// Set Up Initial Interval
	1496	//
	1497	work[n+1] = gl;
	1498	work[n+iin+1] = gu;
	1499
	1500	//
	1501	// Calling DLAEBZ
	1502	//
	1503	// CALL DLAEBZ( 1, 0, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN,
	1504	// D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ),
	1505	// IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IM,
	1506	// IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO )
	1507	//
	1508	for(tmpi=1; tmpi<=iin; tmpi++)
	1509	{
	1510	ra1siin[tmpi] = d[ibegin-1+tmpi];
	1511	if( ibegin-1+tmpi<n )
	1512	{
	1513	ra2siin[tmpi] = e[ibegin-1+tmpi];
	1514	}
	1515	ra3siin[tmpi] = work[ibegin-1+tmpi];
	1516	ra1siinx2[tmpi,1] = work[n+tmpi];
	1517	ra1siinx2[tmpi,2] = work[n+tmpi+iin];
	1518	ra4siin[tmpi] = work[n+2*iin+tmpi];
	1519	rworkspace[tmpi] = w[m+tmpi];
	1520	iworkspace[tmpi] = iblock[m+tmpi];
	1521	ia1siinx2[tmpi,1] = iwork[tmpi];
	1522	ia1siinx2[tmpi,2] = iwork[tmpi+iin];
	1523	}
	1524	internaldlaebz(1, 0, iin, iin, 1, atoli, rtoli, pivmin, ref ra1siin, ref ra2siin, ref ra3siin, ref idumma, ref ra1siinx2, ref ra4siin, ref im, ref ia1siinx2, ref rworkspace, ref iworkspace, ref iinfo);
	1525	for(tmpi=1; tmpi<=iin; tmpi++)
	1526	{
	1527	work[n+tmpi] = ra1siinx2[tmpi,1];
	1528	work[n+tmpi+iin] = ra1siinx2[tmpi,2];
	1529	work[n+2*iin+tmpi] = ra4siin[tmpi];
	1530	w[m+tmpi] = rworkspace[tmpi];
	1531	iblock[m+tmpi] = iworkspace[tmpi];
	1532	iwork[tmpi] = ia1siinx2[tmpi,1];
	1533	iwork[tmpi+iin] = ia1siinx2[tmpi,2];
	1534	}
	1535	nwl = nwl+iwork[1];
	1536	nwu = nwu+iwork[iin+1];
	1537	iwoff = m-iwork[1];
	1538
	1539	//
	1540	// Compute Eigenvalues
	1541	//
	1542	itmax = (int)Math.Ceiling((Math.Log(gu-gl+pivmin)-Math.Log(pivmin))/Math.Log(2))+2;
	1543
	1544	//
	1545	// Calling DLAEBZ
	1546	//
	1547	//CALL DLAEBZ( 2, ITMAX, IN, IN, 1, NB, ATOLI, RTOLI, PIVMIN,
	1548	// D( IBEGIN ), E( IBEGIN ), WORK( IBEGIN ),
	1549	// IDUMMA, WORK( N+1 ), WORK( N+2*IN+1 ), IOUT,
	1550	// IWORK, W( M+1 ), IBLOCK( M+1 ), IINFO )
	1551	//
	1552	for(tmpi=1; tmpi<=iin; tmpi++)
	1553	{
	1554	ra1siin[tmpi] = d[ibegin-1+tmpi];
	1555	if( ibegin-1+tmpi<n )
	1556	{
	1557	ra2siin[tmpi] = e[ibegin-1+tmpi];
	1558	}
	1559	ra3siin[tmpi] = work[ibegin-1+tmpi];
	1560	ra1siinx2[tmpi,1] = work[n+tmpi];
	1561	ra1siinx2[tmpi,2] = work[n+tmpi+iin];
	1562	ra4siin[tmpi] = work[n+2*iin+tmpi];
	1563	rworkspace[tmpi] = w[m+tmpi];
	1564	iworkspace[tmpi] = iblock[m+tmpi];
	1565	ia1siinx2[tmpi,1] = iwork[tmpi];
	1566	ia1siinx2[tmpi,2] = iwork[tmpi+iin];
	1567	}
	1568	internaldlaebz(2, itmax, iin, iin, 1, atoli, rtoli, pivmin, ref ra1siin, ref ra2siin, ref ra3siin, ref idumma, ref ra1siinx2, ref ra4siin, ref iout, ref ia1siinx2, ref rworkspace, ref iworkspace, ref iinfo);
	1569	for(tmpi=1; tmpi<=iin; tmpi++)
	1570	{
	1571	work[n+tmpi] = ra1siinx2[tmpi,1];
	1572	work[n+tmpi+iin] = ra1siinx2[tmpi,2];
	1573	work[n+2*iin+tmpi] = ra4siin[tmpi];
	1574	w[m+tmpi] = rworkspace[tmpi];
	1575	iblock[m+tmpi] = iworkspace[tmpi];
	1576	iwork[tmpi] = ia1siinx2[tmpi,1];
	1577	iwork[tmpi+iin] = ia1siinx2[tmpi,2];
	1578	}
	1579
	1580	//
	1581	// Copy Eigenvalues Into W and IBLOCK
	1582	// Use -JB for block number for unconverged eigenvalues.
	1583	//
	1584	for(j=1; j<=iout; j++)
	1585	{
	1586	tmp1 = 0.5*(work[j+n]+work[j+iin+n]);
	1587
	1588	//
	1589	// Flag non-convergence.
	1590	//
	1591	if( j>iout-iinfo )
	1592	{
	1593	ncnvrg = true;
	1594	ib = -jb;
	1595	}
	1596	else
	1597	{
	1598	ib = jb;
	1599	}
	1600	for(je=iwork[j]+1+iwoff; je<=iwork[j+iin]+iwoff; je++)
	1601	{
	1602	w[je] = tmp1;
	1603	iblock[je] = ib;
	1604	}
	1605	}
	1606	m = m+im;
	1607	}
	1608	}
	1609
	1610	//
	1611	// If RANGE='I', then (WL,WU) contains eigenvalues NWL+1,...,NWU
	1612	// If NWL+1 < IL or NWU > IU, discard extra eigenvalues.
	1613	//
	1614	if( irange==3 )
	1615	{
	1616	im = 0;
	1617	idiscl = il-1-nwl;
	1618	idiscu = nwu-iu;
	1619	if( idiscl>0 \| idiscu>0 )
	1620	{
	1621	for(je=1; je<=m; je++)
	1622	{
	1623	if( (double)(w[je])<=(double)(wlu) & idiscl>0 )
	1624	{
	1625	idiscl = idiscl-1;
	1626	}
	1627	else
	1628	{
	1629	if( (double)(w[je])>=(double)(wul) & idiscu>0 )
	1630	{
	1631	idiscu = idiscu-1;
	1632	}
	1633	else
	1634	{
	1635	im = im+1;
	1636	w[im] = w[je];
	1637	iblock[im] = iblock[je];
	1638	}
	1639	}
	1640	}
	1641	m = im;
	1642	}
	1643	if( idiscl>0 \| idiscu>0 )
	1644	{
	1645
	1646	//
	1647	// Code to deal with effects of bad arithmetic:
	1648	// Some low eigenvalues to be discarded are not in (WL,WLU],
	1649	// or high eigenvalues to be discarded are not in (WUL,WU]
	1650	// so just kill off the smallest IDISCL/largest IDISCU
	1651	// eigenvalues, by simply finding the smallest/largest
	1652	// eigenvalue(s).
	1653	//
	1654	// (If N(w) is monotone non-decreasing, this should never
	1655	// happen.)
	1656	//
	1657	if( idiscl>0 )
	1658	{
	1659	wkill = wu;
	1660	for(jdisc=1; jdisc<=idiscl; jdisc++)
	1661	{
	1662	iw = 0;
	1663	for(je=1; je<=m; je++)
	1664	{
	1665	if( iblock[je]!=0 & ((double)(w[je])<(double)(wkill) \| iw==0) )
	1666	{
	1667	iw = je;
	1668	wkill = w[je];
	1669	}
	1670	}
	1671	iblock[iw] = 0;
	1672	}
	1673	}
	1674	if( idiscu>0 )
	1675	{
	1676	wkill = wl;
	1677	for(jdisc=1; jdisc<=idiscu; jdisc++)
	1678	{
	1679	iw = 0;
	1680	for(je=1; je<=m; je++)
	1681	{
	1682	if( iblock[je]!=0 & ((double)(w[je])>(double)(wkill) \| iw==0) )
	1683	{
	1684	iw = je;
	1685	wkill = w[je];
	1686	}
	1687	}
	1688	iblock[iw] = 0;
	1689	}
	1690	}
	1691	im = 0;
	1692	for(je=1; je<=m; je++)
	1693	{
	1694	if( iblock[je]!=0 )
	1695	{
	1696	im = im+1;
	1697	w[im] = w[je];
	1698	iblock[im] = iblock[je];
	1699	}
	1700	}
	1701	m = im;
	1702	}
	1703	if( idiscl<0 \| idiscu<0 )
	1704	{
	1705	toofew = true;
	1706	}
	1707	}
	1708
	1709	//
	1710	// If ORDER='B', do nothing -- the eigenvalues are already sorted
	1711	// by block.
	1712	// If ORDER='E', sort the eigenvalues from smallest to largest
	1713	//
	1714	if( iorder==1 & nsplit>1 )
	1715	{
	1716	for(je=1; je<=m-1; je++)
	1717	{
	1718	ie = 0;
	1719	tmp1 = w[je];
	1720	for(j=je+1; j<=m; j++)
	1721	{
	1722	if( (double)(w[j])<(double)(tmp1) )
	1723	{
	1724	ie = j;
	1725	tmp1 = w[j];
	1726	}
	1727	}
	1728	if( ie!=0 )
	1729	{
	1730	itmp1 = iblock[ie];
	1731	w[ie] = w[je];
	1732	iblock[ie] = iblock[je];
	1733	w[je] = tmp1;
	1734	iblock[je] = itmp1;
	1735	}
	1736	}
	1737	}
	1738	for(j=1; j<=m; j++)
	1739	{
	1740	w[j] = w[j]*scalefactor;
	1741	}
	1742	errorcode = 0;
	1743	if( ncnvrg )
	1744	{
	1745	errorcode = errorcode+1;
	1746	}
	1747	if( toofew )
	1748	{
	1749	errorcode = errorcode+2;
	1750	}
	1751	result = errorcode==0;
	1752	return result;
	1753	}
	1754
	1755
	1756	public static void internaldstein(int n,
	1757	ref double[] d,
	1758	double[] e,
	1759	int m,
	1760	double[] w,
	1761	ref int[] iblock,
	1762	ref int[] isplit,
	1763	ref double[,] z,
	1764	ref int[] ifail,
	1765	ref int info)
	1766	{
	1767	int maxits = 0;
	1768	int extra = 0;
	1769	int b1 = 0;
	1770	int blksiz = 0;
	1771	int bn = 0;
	1772	int gpind = 0;
	1773	int i = 0;
	1774	int iinfo = 0;
	1775	int its = 0;
	1776	int j = 0;
	1777	int j1 = 0;
	1778	int jblk = 0;
	1779	int jmax = 0;
	1780	int nblk = 0;
	1781	int nrmchk = 0;
	1782	double dtpcrt = 0;
	1783	double eps = 0;
	1784	double eps1 = 0;
	1785	double nrm = 0;
	1786	double onenrm = 0;
	1787	double ortol = 0;
	1788	double pertol = 0;
	1789	double scl = 0;
	1790	double sep = 0;
	1791	double tol = 0;
	1792	double xj = 0;
	1793	double xjm = 0;
	1794	double ztr = 0;
	1795	double[] work1 = new double[0];
	1796	double[] work2 = new double[0];
	1797	double[] work3 = new double[0];
	1798	double[] work4 = new double[0];
	1799	double[] work5 = new double[0];
	1800	int[] iwork = new int[0];
	1801	bool tmpcriterion = new bool();
	1802	int ti = 0;
	1803	int i1 = 0;
	1804	int i2 = 0;
	1805	double v = 0;
	1806	int i_ = 0;
	1807	int i1_ = 0;
	1808
	1809	e = (double[])e.Clone();
	1810	w = (double[])w.Clone();
	1811
	1812	maxits = 5;
	1813	extra = 2;
	1814	work1 = new double[Math.Max(n, 1)+1];
	1815	work2 = new double[Math.Max(n-1, 1)+1];
	1816	work3 = new double[Math.Max(n, 1)+1];
	1817	work4 = new double[Math.Max(n, 1)+1];
	1818	work5 = new double[Math.Max(n, 1)+1];
	1819	iwork = new int[Math.Max(n, 1)+1];
	1820	ifail = new int[Math.Max(m, 1)+1];
	1821	z = new double[Math.Max(n, 1)+1, Math.Max(m, 1)+1];
	1822
	1823	//
	1824	// Test the input parameters.
	1825	//
	1826	info = 0;
	1827	for(i=1; i<=m; i++)
	1828	{
	1829	ifail[i] = 0;
	1830	}
	1831	if( n<0 )
	1832	{
	1833	info = -1;
	1834	return;
	1835	}
	1836	if( m<0 \| m>n )
	1837	{
	1838	info = -4;
	1839	return;
	1840	}
	1841	for(j=2; j<=m; j++)
	1842	{
	1843	if( iblock[j]<iblock[j-1] )
	1844	{
	1845	info = -6;
	1846	break;
	1847	}
	1848	if( iblock[j]==iblock[j-1] & (double)(w[j])<(double)(w[j-1]) )
	1849	{
	1850	info = -5;
	1851	break;
	1852	}
	1853	}
	1854	if( info!=0 )
	1855	{
	1856	return;
	1857	}
	1858
	1859	//
	1860	// Quick return if possible
	1861	//
	1862	if( n==0 \| m==0 )
	1863	{
	1864	return;
	1865	}
	1866	if( n==1 )
	1867	{
	1868	z[1,1] = 1;
	1869	return;
	1870	}
	1871
	1872	//
	1873	// Some preparations
	1874	//
	1875	ti = n-1;
	1876	for(i_=1; i_<=ti;i_++)
	1877	{
	1878	work1[i_] = e[i_];
	1879	}
	1880	e = new double[n+1];
	1881	for(i_=1; i_<=ti;i_++)
	1882	{
	1883	e[i_] = work1[i_];
	1884	}
	1885	for(i_=1; i_<=m;i_++)
	1886	{
	1887	work1[i_] = w[i_];
	1888	}
	1889	w = new double[n+1];
	1890	for(i_=1; i_<=m;i_++)
	1891	{
	1892	w[i_] = work1[i_];
	1893	}
	1894
	1895	//
	1896	// Get machine constants.
	1897	//
	1898	eps = AP.Math.MachineEpsilon;
	1899
	1900	//
	1901	// Compute eigenvectors of matrix blocks.
	1902	//
	1903	j1 = 1;
	1904	for(nblk=1; nblk<=iblock[m]; nblk++)
	1905	{
	1906
	1907	//
	1908	// Find starting and ending indices of block nblk.
	1909	//
	1910	if( nblk==1 )
	1911	{
	1912	b1 = 1;
	1913	}
	1914	else
	1915	{
	1916	b1 = isplit[nblk-1]+1;
	1917	}
	1918	bn = isplit[nblk];
	1919	blksiz = bn-b1+1;
	1920	if( blksiz!=1 )
	1921	{
	1922
	1923	//
	1924	// Compute reorthogonalization criterion and stopping criterion.
	1925	//
	1926	gpind = b1;
	1927	onenrm = Math.Abs(d[b1])+Math.Abs(e[b1]);
	1928	onenrm = Math.Max(onenrm, Math.Abs(d[bn])+Math.Abs(e[bn-1]));
	1929	for(i=b1+1; i<=bn-1; i++)
	1930	{
	1931	onenrm = Math.Max(onenrm, Math.Abs(d[i])+Math.Abs(e[i-1])+Math.Abs(e[i]));
	1932	}
	1933	ortol = 0.001*onenrm;
	1934	dtpcrt = Math.Sqrt(0.1/blksiz);
	1935	}
	1936
	1937	//
	1938	// Loop through eigenvalues of block nblk.
	1939	//
	1940	jblk = 0;
	1941	for(j=j1; j<=m; j++)
	1942	{
	1943	if( iblock[j]!=nblk )
	1944	{
	1945	j1 = j;
	1946	break;
	1947	}
	1948	jblk = jblk+1;
	1949	xj = w[j];
	1950	if( blksiz==1 )
	1951	{
	1952
	1953	//
	1954	// Skip all the work if the block size is one.
	1955	//
	1956	work1[1] = 1;
	1957	}
	1958	else
	1959	{
	1960
	1961	//
	1962	// If eigenvalues j and j-1 are too close, add a relatively
	1963	// small perturbation.
	1964	//
	1965	if( jblk>1 )
	1966	{
	1967	eps1 = Math.Abs(eps*xj);
	1968	pertol = 10*eps1;
	1969	sep = xj-xjm;
	1970	if( (double)(sep)<(double)(pertol) )
	1971	{
	1972	xj = xjm+pertol;
	1973	}
	1974	}
	1975	its = 0;
	1976	nrmchk = 0;
	1977
	1978	//
	1979	// Get random starting vector.
	1980	//
	1981	for(ti=1; ti<=blksiz; ti++)
	1982	{
	1983	work1[ti] = 2*AP.Math.RandomReal()-1;
	1984	}
	1985
	1986	//
	1987	// Copy the matrix T so it won't be destroyed in factorization.
	1988	//
	1989	for(ti=1; ti<=blksiz-1; ti++)
	1990	{
	1991	work2[ti] = e[b1+ti-1];
	1992	work3[ti] = e[b1+ti-1];
	1993	work4[ti] = d[b1+ti-1];
	1994	}
	1995	work4[blksiz] = d[b1+blksiz-1];
	1996
	1997	//
	1998	// Compute LU factors with partial pivoting ( PT = LU )
	1999	//
	2000	tol = 0;
	2001	tdininternaldlagtf(blksiz, ref work4, xj, ref work2, ref work3, tol, ref work5, ref iwork, ref iinfo);
	2002
	2003	//
	2004	// Update iteration count.
	2005	//
	2006	do
	2007	{
	2008	its = its+1;
	2009	if( its>maxits )
	2010	{
	2011
	2012	//
	2013	// If stopping criterion was not satisfied, update info and
	2014	// store eigenvector number in array ifail.
	2015	//
	2016	info = info+1;
	2017	ifail[info] = j;
	2018	break;
	2019	}
	2020
	2021	//
	2022	// Normalize and scale the righthand side vector Pb.
	2023	//
	2024	v = 0;
	2025	for(ti=1; ti<=blksiz; ti++)
	2026	{
	2027	v = v+Math.Abs(work1[ti]);
	2028	}
	2029	scl = blksizonenrmMath.Max(eps, Math.Abs(work4[blksiz]))/v;
	2030	for(i_=1; i_<=blksiz;i_++)
	2031	{
	2032	work1[i_] = scl*work1[i_];
	2033	}
	2034
	2035	//
	2036	// Solve the system LU = Pb.
	2037	//
	2038	tdininternaldlagts(blksiz, ref work4, ref work2, ref work3, ref work5, ref iwork, ref work1, ref tol, ref iinfo);
	2039
	2040	//
	2041	// Reorthogonalize by modified Gram-Schmidt if eigenvalues are
	2042	// close enough.
	2043	//
	2044	if( jblk!=1 )
	2045	{
	2046	if( (double)(Math.Abs(xj-xjm))>(double)(ortol) )
	2047	{
	2048	gpind = j;
	2049	}
	2050	if( gpind!=j )
	2051	{
	2052	for(i=gpind; i<=j-1; i++)
	2053	{
	2054	i1 = b1;
	2055	i2 = b1+blksiz-1;
	2056	i1_ = (i1)-(1);
	2057	ztr = 0.0;
	2058	for(i_=1; i_<=blksiz;i_++)
	2059	{
	2060	ztr += work1[i_]*z[i_+i1_,i];
	2061	}
	2062	i1_ = (i1) - (1);
	2063	for(i_=1; i_<=blksiz;i_++)
	2064	{
	2065	work1[i_] = work1[i_] - ztr*z[i_+i1_,i];
	2066	}
	2067	}
	2068	}
	2069	}
	2070
	2071	//
	2072	// Check the infinity norm of the iterate.
	2073	//
	2074	jmax = blas.vectoridxabsmax(ref work1, 1, blksiz);
	2075	nrm = Math.Abs(work1[jmax]);
	2076
	2077	//
	2078	// Continue for additional iterations after norm reaches
	2079	// stopping criterion.
	2080	//
	2081	tmpcriterion = false;
	2082	if( (double)(nrm)<(double)(dtpcrt) )
	2083	{
	2084	tmpcriterion = true;
	2085	}
	2086	else
	2087	{
	2088	nrmchk = nrmchk+1;
	2089	if( nrmchk<extra+1 )
	2090	{
	2091	tmpcriterion = true;
	2092	}
	2093	}
	2094	}
	2095	while( tmpcriterion );
	2096
	2097	//
	2098	// Accept iterate as jth eigenvector.
	2099	//
	2100	scl = 1/blas.vectornorm2(ref work1, 1, blksiz);
	2101	jmax = blas.vectoridxabsmax(ref work1, 1, blksiz);
	2102	if( (double)(work1[jmax])<(double)(0) )
	2103	{
	2104	scl = -scl;
	2105	}
	2106	for(i_=1; i_<=blksiz;i_++)
	2107	{
	2108	work1[i_] = scl*work1[i_];
	2109	}
	2110	}
	2111	for(i=1; i<=n; i++)
	2112	{
	2113	z[i,j] = 0;
	2114	}
	2115	for(i=1; i<=blksiz; i++)
	2116	{
	2117	z[b1+i-1,j] = work1[i];
	2118	}
	2119
	2120	//
	2121	// Save the shift to check eigenvalue spacing at next
	2122	// iteration.
	2123	//
	2124	xjm = xj;
	2125	}
	2126	}
	2127	}
	2128
	2129
	2130	private static void tdininternaldlagtf(int n,
	2131	ref double[] a,
	2132	double lambda,
	2133	ref double[] b,
	2134	ref double[] c,
	2135	double tol,
	2136	ref double[] d,
	2137	ref int[] iin,
	2138	ref int info)
	2139	{
	2140	int k = 0;
	2141	double eps = 0;
	2142	double mult = 0;
	2143	double piv1 = 0;
	2144	double piv2 = 0;
	2145	double scale1 = 0;
	2146	double scale2 = 0;
	2147	double temp = 0;
	2148	double tl = 0;
	2149
	2150	info = 0;
	2151	if( n<0 )
	2152	{
	2153	info = -1;
	2154	return;
	2155	}
	2156	if( n==0 )
	2157	{
	2158	return;
	2159	}
	2160	a[1] = a[1]-lambda;
	2161	iin[n] = 0;
	2162	if( n==1 )
	2163	{
	2164	if( (double)(a[1])==(double)(0) )
	2165	{
	2166	iin[1] = 1;
	2167	}
	2168	return;
	2169	}
	2170	eps = AP.Math.MachineEpsilon;
	2171	tl = Math.Max(tol, eps);
	2172	scale1 = Math.Abs(a[1])+Math.Abs(b[1]);
	2173	for(k=1; k<=n-1; k++)
	2174	{
	2175	a[k+1] = a[k+1]-lambda;
	2176	scale2 = Math.Abs(c[k])+Math.Abs(a[k+1]);
	2177	if( k<n-1 )
	2178	{
	2179	scale2 = scale2+Math.Abs(b[k+1]);
	2180	}
	2181	if( (double)(a[k])==(double)(0) )
	2182	{
	2183	piv1 = 0;
	2184	}
	2185	else
	2186	{
	2187	piv1 = Math.Abs(a[k])/scale1;
	2188	}
	2189	if( (double)(c[k])==(double)(0) )
	2190	{
	2191	iin[k] = 0;
	2192	piv2 = 0;
	2193	scale1 = scale2;
	2194	if( k<n-1 )
	2195	{
	2196	d[k] = 0;
	2197	}
	2198	}
	2199	else
	2200	{
	2201	piv2 = Math.Abs(c[k])/scale2;
	2202	if( (double)(piv2)<=(double)(piv1) )
	2203	{
	2204	iin[k] = 0;
	2205	scale1 = scale2;
	2206	c[k] = c[k]/a[k];
	2207	a[k+1] = a[k+1]-c[k]*b[k];
	2208	if( k<n-1 )
	2209	{
	2210	d[k] = 0;
	2211	}
	2212	}
	2213	else
	2214	{
	2215	iin[k] = 1;
	2216	mult = a[k]/c[k];
	2217	a[k] = c[k];
	2218	temp = a[k+1];
	2219	a[k+1] = b[k]-mult*temp;
	2220	if( k<n-1 )
	2221	{
	2222	d[k] = b[k+1];
	2223	b[k+1] = -(mult*d[k]);
	2224	}
	2225	b[k] = temp;
	2226	c[k] = mult;
	2227	}
	2228	}
	2229	if( (double)(Math.Max(piv1, piv2))<=(double)(tl) & iin[n]==0 )
	2230	{
	2231	iin[n] = k;
	2232	}
	2233	}
	2234	if( (double)(Math.Abs(a[n]))<=(double)(scale1*tl) & iin[n]==0 )
	2235	{
	2236	iin[n] = n;
	2237	}
	2238	}
	2239
	2240
	2241	private static void tdininternaldlagts(int n,
	2242	ref double[] a,
	2243	ref double[] b,
	2244	ref double[] c,
	2245	ref double[] d,
	2246	ref int[] iin,
	2247	ref double[] y,
	2248	ref double tol,
	2249	ref int info)
	2250	{
	2251	int k = 0;
	2252	double absak = 0;
	2253	double ak = 0;
	2254	double bignum = 0;
	2255	double eps = 0;
	2256	double pert = 0;
	2257	double sfmin = 0;
	2258	double temp = 0;
	2259
	2260	info = 0;
	2261	if( n<0 )
	2262	{
	2263	info = -1;
	2264	return;
	2265	}
	2266	if( n==0 )
	2267	{
	2268	return;
	2269	}
	2270	eps = AP.Math.MachineEpsilon;
	2271	sfmin = AP.Math.MinRealNumber;
	2272	bignum = 1/sfmin;
	2273	if( (double)(tol)<=(double)(0) )
	2274	{
	2275	tol = Math.Abs(a[1]);
	2276	if( n>1 )
	2277	{
	2278	tol = Math.Max(tol, Math.Max(Math.Abs(a[2]), Math.Abs(b[1])));
	2279	}
	2280	for(k=3; k<=n; k++)
	2281	{
	2282	tol = Math.Max(tol, Math.Max(Math.Abs(a[k]), Math.Max(Math.Abs(b[k-1]), Math.Abs(d[k-2]))));
	2283	}
	2284	tol = tol*eps;
	2285	if( (double)(tol)==(double)(0) )
	2286	{
	2287	tol = eps;
	2288	}
	2289	}
	2290	for(k=2; k<=n; k++)
	2291	{
	2292	if( iin[k-1]==0 )
	2293	{
	2294	y[k] = y[k]-c[k-1]*y[k-1];
	2295	}
	2296	else
	2297	{
	2298	temp = y[k-1];
	2299	y[k-1] = y[k];
	2300	y[k] = temp-c[k-1]*y[k];
	2301	}
	2302	}
	2303	for(k=n; k>=1; k--)
	2304	{
	2305	if( k<=n-2 )
	2306	{
	2307	temp = y[k]-b[k]y[k+1]-d[k]y[k+2];
	2308	}
	2309	else
	2310	{
	2311	if( k==n-1 )
	2312	{
	2313	temp = y[k]-b[k]*y[k+1];
	2314	}
	2315	else
	2316	{
	2317	temp = y[k];
	2318	}
	2319	}
	2320	ak = a[k];
	2321	pert = Math.Abs(tol);
	2322	if( (double)(ak)<(double)(0) )
	2323	{
	2324	pert = -pert;
	2325	}
	2326	while( true )
	2327	{
	2328	absak = Math.Abs(ak);
	2329	if( (double)(absak)<(double)(1) )
	2330	{
	2331	if( (double)(absak)<(double)(sfmin) )
	2332	{
	2333	if( (double)(absak)==(double)(0) \| (double)(Math.Abs(temp)*sfmin)>(double)(absak) )
	2334	{
	2335	ak = ak+pert;
	2336	pert = 2*pert;
	2337	continue;
	2338	}
	2339	else
	2340	{
	2341	temp = temp*bignum;
	2342	ak = ak*bignum;
	2343	}
	2344	}
	2345	else
	2346	{
	2347	if( (double)(Math.Abs(temp))>(double)(absak*bignum) )
	2348	{
	2349	ak = ak+pert;
	2350	pert = 2*pert;
	2351	continue;
	2352	}
	2353	}
	2354	}
	2355	break;
	2356	}
	2357	y[k] = temp/ak;
	2358	}
	2359	}
	2360
	2361
	2362	private static void internaldlaebz(int ijob,
	2363	int nitmax,
	2364	int n,
	2365	int mmax,
	2366	int minp,
	2367	double abstol,
	2368	double reltol,
	2369	double pivmin,
	2370	ref double[] d,
	2371	ref double[] e,
	2372	ref double[] e2,
	2373	ref int[] nval,
	2374	ref double[,] ab,
	2375	ref double[] c,
	2376	ref int mout,
	2377	ref int[,] nab,
	2378	ref double[] work,
	2379	ref int[] iwork,
	2380	ref int info)
	2381	{
	2382	int itmp1 = 0;
	2383	int itmp2 = 0;
	2384	int j = 0;
	2385	int ji = 0;
	2386	int jit = 0;
	2387	int jp = 0;
	2388	int kf = 0;
	2389	int kfnew = 0;
	2390	int kl = 0;
	2391	int klnew = 0;
	2392	double tmp1 = 0;
	2393	double tmp2 = 0;
	2394
	2395	info = 0;
	2396	if( ijob<1 \| ijob>3 )
	2397	{
	2398	info = -1;
	2399	return;
	2400	}
	2401
	2402	//
	2403	// Initialize NAB
	2404	//
	2405	if( ijob==1 )
	2406	{
	2407
	2408	//
	2409	// Compute the number of eigenvalues in the initial intervals.
	2410	//
	2411	mout = 0;
	2412
	2413	//
	2414	//DIR$ NOVECTOR
	2415	//
	2416	for(ji=1; ji<=minp; ji++)
	2417	{
	2418	for(jp=1; jp<=2; jp++)
	2419	{
	2420	tmp1 = d[1]-ab[ji,jp];
	2421	if( (double)(Math.Abs(tmp1))<(double)(pivmin) )
	2422	{
	2423	tmp1 = -pivmin;
	2424	}
	2425	nab[ji,jp] = 0;
	2426	if( (double)(tmp1)<=(double)(0) )
	2427	{
	2428	nab[ji,jp] = 1;
	2429	}
	2430	for(j=2; j<=n; j++)
	2431	{
	2432	tmp1 = d[j]-e2[j-1]/tmp1-ab[ji,jp];
	2433	if( (double)(Math.Abs(tmp1))<(double)(pivmin) )
	2434	{
	2435	tmp1 = -pivmin;
	2436	}
	2437	if( (double)(tmp1)<=(double)(0) )
	2438	{
	2439	nab[ji,jp] = nab[ji,jp]+1;
	2440	}
	2441	}
	2442	}
	2443	mout = mout+nab[ji,2]-nab[ji,1];
	2444	}
	2445	return;
	2446	}
	2447
	2448	//
	2449	// Initialize for loop
	2450	//
	2451	// KF and KL have the following meaning:
	2452	// Intervals 1,...,KF-1 have converged.
	2453	// Intervals KF,...,KL still need to be refined.
	2454	//
	2455	kf = 1;
	2456	kl = minp;
	2457
	2458	//
	2459	// If IJOB=2, initialize C.
	2460	// If IJOB=3, use the user-supplied starting point.
	2461	//
	2462	if( ijob==2 )
	2463	{
	2464	for(ji=1; ji<=minp; ji++)
	2465	{
	2466	c[ji] = 0.5*(ab[ji,1]+ab[ji,2]);
	2467	}
	2468	}
	2469
	2470	//
	2471	// Iteration loop
	2472	//
	2473	for(jit=1; jit<=nitmax; jit++)
	2474	{
	2475
	2476	//
	2477	// Loop over intervals
	2478	//
	2479	//
	2480	// Serial Version of the loop
	2481	//
	2482	klnew = kl;
	2483	for(ji=kf; ji<=kl; ji++)
	2484	{
	2485
	2486	//
	2487	// Compute N(w), the number of eigenvalues less than w
	2488	//
	2489	tmp1 = c[ji];
	2490	tmp2 = d[1]-tmp1;
	2491	itmp1 = 0;
	2492	if( (double)(tmp2)<=(double)(pivmin) )
	2493	{
	2494	itmp1 = 1;
	2495	tmp2 = Math.Min(tmp2, -pivmin);
	2496	}
	2497
	2498	//
	2499	// A series of compiler directives to defeat vectorization
	2500	// for the next loop
	2501	//
	2502	//*$PL$ CMCHAR=' '
	2503	//CDIR$ NEXTSCALAR
	2504	//C$DIR SCALAR
	2505	//CDIR$ NEXT SCALAR
	2506	//CVD$L NOVECTOR
	2507	//CDEC$ NOVECTOR
	2508	//CVD$ NOVECTOR
	2509	//*VDIR NOVECTOR
	2510	//*VOCL LOOP,SCALAR
	2511	//CIBM PREFER SCALAR
	2512	//$PL$ CMCHAR=''
	2513	//
	2514	for(j=2; j<=n; j++)
	2515	{
	2516	tmp2 = d[j]-e2[j-1]/tmp2-tmp1;
	2517	if( (double)(tmp2)<=(double)(pivmin) )
	2518	{
	2519	itmp1 = itmp1+1;
	2520	tmp2 = Math.Min(tmp2, -pivmin);
	2521	}
	2522	}
	2523	if( ijob<=2 )
	2524	{
	2525
	2526	//
	2527	// IJOB=2: Choose all intervals containing eigenvalues.
	2528	//
	2529	// Insure that N(w) is monotone
	2530	//
	2531	itmp1 = Math.Min(nab[ji,2], Math.Max(nab[ji,1], itmp1));
	2532
	2533	//
	2534	// Update the Queue -- add intervals if both halves
	2535	// contain eigenvalues.
	2536	//
	2537	if( itmp1==nab[ji,2] )
	2538	{
	2539
	2540	//
	2541	// No eigenvalue in the upper interval:
	2542	// just use the lower interval.
	2543	//
	2544	ab[ji,2] = tmp1;
	2545	}
	2546	else
	2547	{
	2548	if( itmp1==nab[ji,1] )
	2549	{
	2550
	2551	//
	2552	// No eigenvalue in the lower interval:
	2553	// just use the upper interval.
	2554	//
	2555	ab[ji,1] = tmp1;
	2556	}
	2557	else
	2558	{
	2559	if( klnew<mmax )
	2560	{
	2561
	2562	//
	2563	// Eigenvalue in both intervals -- add upper to queue.
	2564	//
	2565	klnew = klnew+1;
	2566	ab[klnew,2] = ab[ji,2];
	2567	nab[klnew,2] = nab[ji,2];
	2568	ab[klnew,1] = tmp1;
	2569	nab[klnew,1] = itmp1;
	2570	ab[ji,2] = tmp1;
	2571	nab[ji,2] = itmp1;
	2572	}
	2573	else
	2574	{
	2575	info = mmax+1;
	2576	return;
	2577	}
	2578	}
	2579	}
	2580	}
	2581	else
	2582	{
	2583
	2584	//
	2585	// IJOB=3: Binary search. Keep only the interval
	2586	// containing w s.t. N(w) = NVAL
	2587	//
	2588	if( itmp1<=nval[ji] )
	2589	{
	2590	ab[ji,1] = tmp1;
	2591	nab[ji,1] = itmp1;
	2592	}
	2593	if( itmp1>=nval[ji] )
	2594	{
	2595	ab[ji,2] = tmp1;
	2596	nab[ji,2] = itmp1;
	2597	}
	2598	}
	2599	}
	2600	kl = klnew;
	2601
	2602	//
	2603	// Check for convergence
	2604	//
	2605	kfnew = kf;
	2606	for(ji=kf; ji<=kl; ji++)
	2607	{
	2608	tmp1 = Math.Abs(ab[ji,2]-ab[ji,1]);
	2609	tmp2 = Math.Max(Math.Abs(ab[ji,2]), Math.Abs(ab[ji,1]));
	2610	if( (double)(tmp1)<(double)(Math.Max(abstol, Math.Max(pivmin, reltol*tmp2))) \| nab[ji,1]>=nab[ji,2] )
	2611	{
	2612
	2613	//
	2614	// Converged -- Swap with position KFNEW,
	2615	// then increment KFNEW
	2616	//
	2617	if( ji>kfnew )
	2618	{
	2619	tmp1 = ab[ji,1];
	2620	tmp2 = ab[ji,2];
	2621	itmp1 = nab[ji,1];
	2622	itmp2 = nab[ji,2];
	2623	ab[ji,1] = ab[kfnew,1];
	2624	ab[ji,2] = ab[kfnew,2];
	2625	nab[ji,1] = nab[kfnew,1];
	2626	nab[ji,2] = nab[kfnew,2];
	2627	ab[kfnew,1] = tmp1;
	2628	ab[kfnew,2] = tmp2;
	2629	nab[kfnew,1] = itmp1;
	2630	nab[kfnew,2] = itmp2;
	2631	if( ijob==3 )
	2632	{
	2633	itmp1 = nval[ji];
	2634	nval[ji] = nval[kfnew];
	2635	nval[kfnew] = itmp1;
	2636	}
	2637	}
	2638	kfnew = kfnew+1;
	2639	}
	2640	}
	2641	kf = kfnew;
	2642
	2643	//
	2644	// Choose Midpoints
	2645	//
	2646	for(ji=kf; ji<=kl; ji++)
	2647	{
	2648	c[ji] = 0.5*(ab[ji,1]+ab[ji,2]);
	2649	}
	2650
	2651	//
	2652	// If no more intervals to refine, quit.
	2653	//
	2654	if( kf>kl )
	2655	{
	2656	break;
	2657	}
	2658	}
	2659
	2660	//
	2661	// Converged
	2662	//
	2663	info = Math.Max(kl+1-kf, 0);
	2664	mout = kl;
	2665	}
	2666	}
	2667	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Update cookies preferences