Context Navigation

source: trunk/sources/ALGLIB/tdevd.cs @ 2615

Visit:

Last change on this file since 2615 was 2563, checked in by gkronber, 15 years ago
Updated ALGLIB to latest version. #751 (Plugin for for data-modeling with ANN (integrated into CEDMA))
File size: 43.9 KB

Rev	Line
[2563]	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 tdevd
	32	{
	33	/*************************************************************************
	34	Finding the eigenvalues and eigenvectors of a tridiagonal symmetric matrix
	35
	36	The algorithm finds the eigen pairs of a tridiagonal symmetric matrix by
	37	using an QL/QR algorithm with implicit shifts.
	38
	39	Input parameters:
	40	D - the main diagonal of a tridiagonal matrix.
	41	Array whose index ranges within [0..N-1].
	42	E - the secondary diagonal of a tridiagonal matrix.
	43	Array whose index ranges within [0..N-2].
	44	N - size of matrix A.
	45	ZNeeded - flag controlling whether the eigenvectors are needed or not.
	46	If ZNeeded is equal to:
	47	* 0, the eigenvectors are not needed;
	48	* 1, the eigenvectors of a tridiagonal matrix
	49	are multiplied by the square matrix Z. It is used if the
	50	tridiagonal matrix is obtained by the similarity
	51	transformation of a symmetric matrix;
	52	* 2, the eigenvectors of a tridiagonal matrix replace the
	53	square matrix Z;
	54	* 3, matrix Z contains the first row of the eigenvectors
	55	matrix.
	56	Z - if ZNeeded=1, Z contains the square matrix by which the
	57	eigenvectors are multiplied.
	58	Array whose indexes range within [0..N-1, 0..N-1].
	59
	60	Output parameters:
	61	D - eigenvalues in ascending order.
	62	Array whose index ranges within [0..N-1].
	63	Z - if ZNeeded is equal to:
	64	* 0, Z hasnt changed;
	65	* 1, Z contains the product of a given matrix (from the left)
	66	and the eigenvectors matrix (from the right);
	67	* 2, Z contains the eigenvectors.
	68	* 3, Z contains the first row of the eigenvectors matrix.
	69	If ZNeeded<3, Z is the array whose indexes range within [0..N-1, 0..N-1].
	70	In that case, the eigenvectors are stored in the matrix columns.
	71	If ZNeeded=3, Z is the array whose indexes range within [0..0, 0..N-1].
	72
	73	Result:
	74	True, if the algorithm has converged.
	75	False, if the algorithm hasn't converged.
	76
	77	-- LAPACK routine (version 3.0) --
	78	Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
	79	Courant Institute, Argonne National Lab, and Rice University
	80	September 30, 1994
	81	*************************************************************************/
	82	public static bool smatrixtdevd(ref double[] d,
	83	double[] e,
	84	int n,
	85	int zneeded,
	86	ref double[,] z)
	87	{
	88	bool result = new bool();
	89	double[] d1 = new double[0];
	90	double[] e1 = new double[0];
	91	double[,] z1 = new double[0,0];
	92	int i = 0;
	93	int i_ = 0;
	94	int i1_ = 0;
	95
	96	e = (double[])e.Clone();
	97
	98
	99	//
	100	// Prepare 1-based task
	101	//
	102	d1 = new double[n+1];
	103	e1 = new double[n+1];
	104	i1_ = (0) - (1);
	105	for(i_=1; i_<=n;i_++)
	106	{
	107	d1[i_] = d[i_+i1_];
	108	}
	109	if( n>1 )
	110	{
	111	i1_ = (0) - (1);
	112	for(i_=1; i_<=n-1;i_++)
	113	{
	114	e1[i_] = e[i_+i1_];
	115	}
	116	}
	117	if( zneeded==1 )
	118	{
	119	z1 = new double[n+1, n+1];
	120	for(i=1; i<=n; i++)
	121	{
	122	i1_ = (0) - (1);
	123	for(i_=1; i_<=n;i_++)
	124	{
	125	z1[i,i_] = z[i-1,i_+i1_];
	126	}
	127	}
	128	}
	129
	130	//
	131	// Solve 1-based task
	132	//
	133	result = tridiagonalevd(ref d1, e1, n, zneeded, ref z1);
	134	if( !result )
	135	{
	136	return result;
	137	}
	138
	139	//
	140	// Convert back to 0-based result
	141	//
	142	i1_ = (1) - (0);
	143	for(i_=0; i_<=n-1;i_++)
	144	{
	145	d[i_] = d1[i_+i1_];
	146	}
	147	if( zneeded!=0 )
	148	{
	149	if( zneeded==1 )
	150	{
	151	for(i=1; i<=n; i++)
	152	{
	153	i1_ = (1) - (0);
	154	for(i_=0; i_<=n-1;i_++)
	155	{
	156	z[i-1,i_] = z1[i,i_+i1_];
	157	}
	158	}
	159	return result;
	160	}
	161	if( zneeded==2 )
	162	{
	163	z = new double[n-1+1, n-1+1];
	164	for(i=1; i<=n; i++)
	165	{
	166	i1_ = (1) - (0);
	167	for(i_=0; i_<=n-1;i_++)
	168	{
	169	z[i-1,i_] = z1[i,i_+i1_];
	170	}
	171	}
	172	return result;
	173	}
	174	if( zneeded==3 )
	175	{
	176	z = new double[0+1, n-1+1];
	177	i1_ = (1) - (0);
	178	for(i_=0; i_<=n-1;i_++)
	179	{
	180	z[0,i_] = z1[1,i_+i1_];
	181	}
	182	return result;
	183	}
	184	System.Diagnostics.Debug.Assert(false, "SMatrixTDEVD: Incorrect ZNeeded!");
	185	}
	186	return result;
	187	}
	188
	189
	190	public static bool tridiagonalevd(ref double[] d,
	191	double[] e,
	192	int n,
	193	int zneeded,
	194	ref double[,] z)
	195	{
	196	bool result = new bool();
	197	int maxit = 0;
	198	int i = 0;
	199	int icompz = 0;
	200	int ii = 0;
	201	int iscale = 0;
	202	int j = 0;
	203	int jtot = 0;
	204	int k = 0;
	205	int t = 0;
	206	int l = 0;
	207	int l1 = 0;
	208	int lend = 0;
	209	int lendm1 = 0;
	210	int lendp1 = 0;
	211	int lendsv = 0;
	212	int lm1 = 0;
	213	int lsv = 0;
	214	int m = 0;
	215	int mm = 0;
	216	int mm1 = 0;
	217	int nm1 = 0;
	218	int nmaxit = 0;
	219	int tmpint = 0;
	220	double anorm = 0;
	221	double b = 0;
	222	double c = 0;
	223	double eps = 0;
	224	double eps2 = 0;
	225	double f = 0;
	226	double g = 0;
	227	double p = 0;
	228	double r = 0;
	229	double rt1 = 0;
	230	double rt2 = 0;
	231	double s = 0;
	232	double safmax = 0;
	233	double safmin = 0;
	234	double ssfmax = 0;
	235	double ssfmin = 0;
	236	double tst = 0;
	237	double tmp = 0;
	238	double[] work1 = new double[0];
	239	double[] work2 = new double[0];
	240	double[] workc = new double[0];
	241	double[] works = new double[0];
	242	double[] wtemp = new double[0];
	243	bool gotoflag = new bool();
	244	int zrows = 0;
	245	bool wastranspose = new bool();
	246	int i_ = 0;
	247
	248	e = (double[])e.Clone();
	249
	250	System.Diagnostics.Debug.Assert(zneeded>=0 & zneeded<=3, "TridiagonalEVD: Incorrent ZNeeded");
	251
	252	//
	253	// Quick return if possible
	254	//
	255	if( zneeded<0 \| zneeded>3 )
	256	{
	257	result = false;
	258	return result;
	259	}
	260	result = true;
	261	if( n==0 )
	262	{
	263	return result;
	264	}
	265	if( n==1 )
	266	{
	267	if( zneeded==2 \| zneeded==3 )
	268	{
	269	z = new double[1+1, 1+1];
	270	z[1,1] = 1;
	271	}
	272	return result;
	273	}
	274	maxit = 30;
	275
	276	//
	277	// Initialize arrays
	278	//
	279	wtemp = new double[n+1];
	280	work1 = new double[n-1+1];
	281	work2 = new double[n-1+1];
	282	workc = new double[n+1];
	283	works = new double[n+1];
	284
	285	//
	286	// Determine the unit roundoff and over/underflow thresholds.
	287	//
	288	eps = AP.Math.MachineEpsilon;
	289	eps2 = AP.Math.Sqr(eps);
	290	safmin = AP.Math.MinRealNumber;
	291	safmax = AP.Math.MaxRealNumber;
	292	ssfmax = Math.Sqrt(safmax)/3;
	293	ssfmin = Math.Sqrt(safmin)/eps2;
	294
	295	//
	296	// Prepare Z
	297	//
	298	// Here we are using transposition to get rid of column operations
	299	//
	300	//
	301	wastranspose = false;
	302	if( zneeded==0 )
	303	{
	304	zrows = 0;
	305	}
	306	if( zneeded==1 )
	307	{
	308	zrows = n;
	309	}
	310	if( zneeded==2 )
	311	{
	312	zrows = n;
	313	}
	314	if( zneeded==3 )
	315	{
	316	zrows = 1;
	317	}
	318	if( zneeded==1 )
	319	{
	320	wastranspose = true;
	321	blas.inplacetranspose(ref z, 1, n, 1, n, ref wtemp);
	322	}
	323	if( zneeded==2 )
	324	{
	325	wastranspose = true;
	326	z = new double[n+1, n+1];
	327	for(i=1; i<=n; i++)
	328	{
	329	for(j=1; j<=n; j++)
	330	{
	331	if( i==j )
	332	{
	333	z[i,j] = 1;
	334	}
	335	else
	336	{
	337	z[i,j] = 0;
	338	}
	339	}
	340	}
	341	}
	342	if( zneeded==3 )
	343	{
	344	wastranspose = false;
	345	z = new double[1+1, n+1];
	346	for(j=1; j<=n; j++)
	347	{
	348	if( j==1 )
	349	{
	350	z[1,j] = 1;
	351	}
	352	else
	353	{
	354	z[1,j] = 0;
	355	}
	356	}
	357	}
	358	nmaxit = n*maxit;
	359	jtot = 0;
	360
	361	//
	362	// Determine where the matrix splits and choose QL or QR iteration
	363	// for each block, according to whether top or bottom diagonal
	364	// element is smaller.
	365	//
	366	l1 = 1;
	367	nm1 = n-1;
	368	while( true )
	369	{
	370	if( l1>n )
	371	{
	372	break;
	373	}
	374	if( l1>1 )
	375	{
	376	e[l1-1] = 0;
	377	}
	378	gotoflag = false;
	379	if( l1<=nm1 )
	380	{
	381	for(m=l1; m<=nm1; m++)
	382	{
	383	tst = Math.Abs(e[m]);
	384	if( (double)(tst)==(double)(0) )
	385	{
	386	gotoflag = true;
	387	break;
	388	}
	389	if( (double)(tst)<=(double)(Math.Sqrt(Math.Abs(d[m]))Math.Sqrt(Math.Abs(d[m+1]))eps) )
	390	{
	391	e[m] = 0;
	392	gotoflag = true;
	393	break;
	394	}
	395	}
	396	}
	397	if( !gotoflag )
	398	{
	399	m = n;
	400	}
	401
	402	//
	403	// label 30:
	404	//
	405	l = l1;
	406	lsv = l;
	407	lend = m;
	408	lendsv = lend;
	409	l1 = m+1;
	410	if( lend==l )
	411	{
	412	continue;
	413	}
	414
	415	//
	416	// Scale submatrix in rows and columns L to LEND
	417	//
	418	if( l==lend )
	419	{
	420	anorm = Math.Abs(d[l]);
	421	}
	422	else
	423	{
	424	anorm = Math.Max(Math.Abs(d[l])+Math.Abs(e[l]), Math.Abs(e[lend-1])+Math.Abs(d[lend]));
	425	for(i=l+1; i<=lend-1; i++)
	426	{
	427	anorm = Math.Max(anorm, Math.Abs(d[i])+Math.Abs(e[i])+Math.Abs(e[i-1]));
	428	}
	429	}
	430	iscale = 0;
	431	if( (double)(anorm)==(double)(0) )
	432	{
	433	continue;
	434	}
	435	if( (double)(anorm)>(double)(ssfmax) )
	436	{
	437	iscale = 1;
	438	tmp = ssfmax/anorm;
	439	tmpint = lend-1;
	440	for(i_=l; i_<=lend;i_++)
	441	{
	442	d[i_] = tmp*d[i_];
	443	}
	444	for(i_=l; i_<=tmpint;i_++)
	445	{
	446	e[i_] = tmp*e[i_];
	447	}
	448	}
	449	if( (double)(anorm)<(double)(ssfmin) )
	450	{
	451	iscale = 2;
	452	tmp = ssfmin/anorm;
	453	tmpint = lend-1;
	454	for(i_=l; i_<=lend;i_++)
	455	{
	456	d[i_] = tmp*d[i_];
	457	}
	458	for(i_=l; i_<=tmpint;i_++)
	459	{
	460	e[i_] = tmp*e[i_];
	461	}
	462	}
	463
	464	//
	465	// Choose between QL and QR iteration
	466	//
	467	if( (double)(Math.Abs(d[lend]))<(double)(Math.Abs(d[l])) )
	468	{
	469	lend = lsv;
	470	l = lendsv;
	471	}
	472	if( lend>l )
	473	{
	474
	475	//
	476	// QL Iteration
	477	//
	478	// Look for small subdiagonal element.
	479	//
	480	while( true )
	481	{
	482	gotoflag = false;
	483	if( l!=lend )
	484	{
	485	lendm1 = lend-1;
	486	for(m=l; m<=lendm1; m++)
	487	{
	488	tst = AP.Math.Sqr(Math.Abs(e[m]));
	489	if( (double)(tst)<=(double)(eps2Math.Abs(d[m])Math.Abs(d[m+1])+safmin) )
	490	{
	491	gotoflag = true;
	492	break;
	493	}
	494	}
	495	}
	496	if( !gotoflag )
	497	{
	498	m = lend;
	499	}
	500	if( m<lend )
	501	{
	502	e[m] = 0;
	503	}
	504	p = d[l];
	505	if( m!=l )
	506	{
	507
	508	//
	509	// If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
	510	// to compute its eigensystem.
	511	//
	512	if( m==l+1 )
	513	{
	514	if( zneeded>0 )
	515	{
	516	tdevdev2(d[l], e[l], d[l+1], ref rt1, ref rt2, ref c, ref s);
	517	work1[l] = c;
	518	work2[l] = s;
	519	workc[1] = work1[l];
	520	works[1] = work2[l];
	521	if( !wastranspose )
	522	{
	523	rotations.applyrotationsfromtheright(false, 1, zrows, l, l+1, ref workc, ref works, ref z, ref wtemp);
	524	}
	525	else
	526	{
	527	rotations.applyrotationsfromtheleft(false, l, l+1, 1, zrows, ref workc, ref works, ref z, ref wtemp);
	528	}
	529	}
	530	else
	531	{
	532	tdevde2(d[l], e[l], d[l+1], ref rt1, ref rt2);
	533	}
	534	d[l] = rt1;
	535	d[l+1] = rt2;
	536	e[l] = 0;
	537	l = l+2;
	538	if( l<=lend )
	539	{
	540	continue;
	541	}
	542
	543	//
	544	// GOTO 140
	545	//
	546	break;
	547	}
	548	if( jtot==nmaxit )
	549	{
	550
	551	//
	552	// GOTO 140
	553	//
	554	break;
	555	}
	556	jtot = jtot+1;
	557
	558	//
	559	// Form shift.
	560	//
	561	g = (d[l+1]-p)/(2*e[l]);
	562	r = tdevdpythag(g, 1);
	563	g = d[m]-p+e[l]/(g+tdevdextsign(r, g));
	564	s = 1;
	565	c = 1;
	566	p = 0;
	567
	568	//
	569	// Inner loop
	570	//
	571	mm1 = m-1;
	572	for(i=mm1; i>=l; i--)
	573	{
	574	f = s*e[i];
	575	b = c*e[i];
	576	rotations.generaterotation(g, f, ref c, ref s, ref r);
	577	if( i!=m-1 )
	578	{
	579	e[i+1] = r;
	580	}
	581	g = d[i+1]-p;
	582	r = (d[i]-g)s+2c*b;
	583	p = s*r;
	584	d[i+1] = g+p;
	585	g = c*r-b;
	586
	587	//
	588	// If eigenvectors are desired, then save rotations.
	589	//
	590	if( zneeded>0 )
	591	{
	592	work1[i] = c;
	593	work2[i] = -s;
	594	}
	595	}
	596
	597	//
	598	// If eigenvectors are desired, then apply saved rotations.
	599	//
	600	if( zneeded>0 )
	601	{
	602	for(i=l; i<=m-1; i++)
	603	{
	604	workc[i-l+1] = work1[i];
	605	works[i-l+1] = work2[i];
	606	}
	607	if( !wastranspose )
	608	{
	609	rotations.applyrotationsfromtheright(false, 1, zrows, l, m, ref workc, ref works, ref z, ref wtemp);
	610	}
	611	else
	612	{
	613	rotations.applyrotationsfromtheleft(false, l, m, 1, zrows, ref workc, ref works, ref z, ref wtemp);
	614	}
	615	}
	616	d[l] = d[l]-p;
	617	e[l] = g;
	618	continue;
	619	}
	620
	621	//
	622	// Eigenvalue found.
	623	//
	624	d[l] = p;
	625	l = l+1;
	626	if( l<=lend )
	627	{
	628	continue;
	629	}
	630	break;
	631	}
	632	}
	633	else
	634	{
	635
	636	//
	637	// QR Iteration
	638	//
	639	// Look for small superdiagonal element.
	640	//
	641	while( true )
	642	{
	643	gotoflag = false;
	644	if( l!=lend )
	645	{
	646	lendp1 = lend+1;
	647	for(m=l; m>=lendp1; m--)
	648	{
	649	tst = AP.Math.Sqr(Math.Abs(e[m-1]));
	650	if( (double)(tst)<=(double)(eps2Math.Abs(d[m])Math.Abs(d[m-1])+safmin) )
	651	{
	652	gotoflag = true;
	653	break;
	654	}
	655	}
	656	}
	657	if( !gotoflag )
	658	{
	659	m = lend;
	660	}
	661	if( m>lend )
	662	{
	663	e[m-1] = 0;
	664	}
	665	p = d[l];
	666	if( m!=l )
	667	{
	668
	669	//
	670	// If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
	671	// to compute its eigensystem.
	672	//
	673	if( m==l-1 )
	674	{
	675	if( zneeded>0 )
	676	{
	677	tdevdev2(d[l-1], e[l-1], d[l], ref rt1, ref rt2, ref c, ref s);
	678	work1[m] = c;
	679	work2[m] = s;
	680	workc[1] = c;
	681	works[1] = s;
	682	if( !wastranspose )
	683	{
	684	rotations.applyrotationsfromtheright(true, 1, zrows, l-1, l, ref workc, ref works, ref z, ref wtemp);
	685	}
	686	else
	687	{
	688	rotations.applyrotationsfromtheleft(true, l-1, l, 1, zrows, ref workc, ref works, ref z, ref wtemp);
	689	}
	690	}
	691	else
	692	{
	693	tdevde2(d[l-1], e[l-1], d[l], ref rt1, ref rt2);
	694	}
	695	d[l-1] = rt1;
	696	d[l] = rt2;
	697	e[l-1] = 0;
	698	l = l-2;
	699	if( l>=lend )
	700	{
	701	continue;
	702	}
	703	break;
	704	}
	705	if( jtot==nmaxit )
	706	{
	707	break;
	708	}
	709	jtot = jtot+1;
	710
	711	//
	712	// Form shift.
	713	//
	714	g = (d[l-1]-p)/(2*e[l-1]);
	715	r = tdevdpythag(g, 1);
	716	g = d[m]-p+e[l-1]/(g+tdevdextsign(r, g));
	717	s = 1;
	718	c = 1;
	719	p = 0;
	720
	721	//
	722	// Inner loop
	723	//
	724	lm1 = l-1;
	725	for(i=m; i<=lm1; i++)
	726	{
	727	f = s*e[i];
	728	b = c*e[i];
	729	rotations.generaterotation(g, f, ref c, ref s, ref r);
	730	if( i!=m )
	731	{
	732	e[i-1] = r;
	733	}
	734	g = d[i]-p;
	735	r = (d[i+1]-g)s+2c*b;
	736	p = s*r;
	737	d[i] = g+p;
	738	g = c*r-b;
	739
	740	//
	741	// If eigenvectors are desired, then save rotations.
	742	//
	743	if( zneeded>0 )
	744	{
	745	work1[i] = c;
	746	work2[i] = s;
	747	}
	748	}
	749
	750	//
	751	// If eigenvectors are desired, then apply saved rotations.
	752	//
	753	if( zneeded>0 )
	754	{
	755	mm = l-m+1;
	756	for(i=m; i<=l-1; i++)
	757	{
	758	workc[i-m+1] = work1[i];
	759	works[i-m+1] = work2[i];
	760	}
	761	if( !wastranspose )
	762	{
	763	rotations.applyrotationsfromtheright(true, 1, zrows, m, l, ref workc, ref works, ref z, ref wtemp);
	764	}
	765	else
	766	{
	767	rotations.applyrotationsfromtheleft(true, m, l, 1, zrows, ref workc, ref works, ref z, ref wtemp);
	768	}
	769	}
	770	d[l] = d[l]-p;
	771	e[lm1] = g;
	772	continue;
	773	}
	774
	775	//
	776	// Eigenvalue found.
	777	//
	778	d[l] = p;
	779	l = l-1;
	780	if( l>=lend )
	781	{
	782	continue;
	783	}
	784	break;
	785	}
	786	}
	787
	788	//
	789	// Undo scaling if necessary
	790	//
	791	if( iscale==1 )
	792	{
	793	tmp = anorm/ssfmax;
	794	tmpint = lendsv-1;
	795	for(i_=lsv; i_<=lendsv;i_++)
	796	{
	797	d[i_] = tmp*d[i_];
	798	}
	799	for(i_=lsv; i_<=tmpint;i_++)
	800	{
	801	e[i_] = tmp*e[i_];
	802	}
	803	}
	804	if( iscale==2 )
	805	{
	806	tmp = anorm/ssfmin;
	807	tmpint = lendsv-1;
	808	for(i_=lsv; i_<=lendsv;i_++)
	809	{
	810	d[i_] = tmp*d[i_];
	811	}
	812	for(i_=lsv; i_<=tmpint;i_++)
	813	{
	814	e[i_] = tmp*e[i_];
	815	}
	816	}
	817
	818	//
	819	// Check for no convergence to an eigenvalue after a total
	820	// of N*MAXIT iterations.
	821	//
	822	if( jtot>=nmaxit )
	823	{
	824	result = false;
	825	if( wastranspose )
	826	{
	827	blas.inplacetranspose(ref z, 1, n, 1, n, ref wtemp);
	828	}
	829	return result;
	830	}
	831	}
	832
	833	//
	834	// Order eigenvalues and eigenvectors.
	835	//
	836	if( zneeded==0 )
	837	{
	838
	839	//
	840	// Sort
	841	//
	842	if( n==1 )
	843	{
	844	return result;
	845	}
	846	if( n==2 )
	847	{
	848	if( (double)(d[1])>(double)(d[2]) )
	849	{
	850	tmp = d[1];
	851	d[1] = d[2];
	852	d[2] = tmp;
	853	}
	854	return result;
	855	}
	856	i = 2;
	857	do
	858	{
	859	t = i;
	860	while( t!=1 )
	861	{
	862	k = t/2;
	863	if( (double)(d[k])>=(double)(d[t]) )
	864	{
	865	t = 1;
	866	}
	867	else
	868	{
	869	tmp = d[k];
	870	d[k] = d[t];
	871	d[t] = tmp;
	872	t = k;
	873	}
	874	}
	875	i = i+1;
	876	}
	877	while( i<=n );
	878	i = n-1;
	879	do
	880	{
	881	tmp = d[i+1];
	882	d[i+1] = d[1];
	883	d[+1] = tmp;
	884	t = 1;
	885	while( t!=0 )
	886	{
	887	k = 2*t;
	888	if( k>i )
	889	{
	890	t = 0;
	891	}
	892	else
	893	{
	894	if( k<i )
	895	{
	896	if( (double)(d[k+1])>(double)(d[k]) )
	897	{
	898	k = k+1;
	899	}
	900	}
	901	if( (double)(d[t])>=(double)(d[k]) )
	902	{
	903	t = 0;
	904	}
	905	else
	906	{
	907	tmp = d[k];
	908	d[k] = d[t];
	909	d[t] = tmp;
	910	t = k;
	911	}
	912	}
	913	}
	914	i = i-1;
	915	}
	916	while( i>=1 );
	917	}
	918	else
	919	{
	920
	921	//
	922	// Use Selection Sort to minimize swaps of eigenvectors
	923	//
	924	for(ii=2; ii<=n; ii++)
	925	{
	926	i = ii-1;
	927	k = i;
	928	p = d[i];
	929	for(j=ii; j<=n; j++)
	930	{
	931	if( (double)(d[j])<(double)(p) )
	932	{
	933	k = j;
	934	p = d[j];
	935	}
	936	}
	937	if( k!=i )
	938	{
	939	d[k] = d[i];
	940	d[i] = p;
	941	if( wastranspose )
	942	{
	943	for(i_=1; i_<=n;i_++)
	944	{
	945	wtemp[i_] = z[i,i_];
	946	}
	947	for(i_=1; i_<=n;i_++)
	948	{
	949	z[i,i_] = z[k,i_];
	950	}
	951	for(i_=1; i_<=n;i_++)
	952	{
	953	z[k,i_] = wtemp[i_];
	954	}
	955	}
	956	else
	957	{
	958	for(i_=1; i_<=zrows;i_++)
	959	{
	960	wtemp[i_] = z[i_,i];
	961	}
	962	for(i_=1; i_<=zrows;i_++)
	963	{
	964	z[i_,i] = z[i_,k];
	965	}
	966	for(i_=1; i_<=zrows;i_++)
	967	{
	968	z[i_,k] = wtemp[i_];
	969	}
	970	}
	971	}
	972	}
	973	if( wastranspose )
	974	{
	975	blas.inplacetranspose(ref z, 1, n, 1, n, ref wtemp);
	976	}
	977	}
	978	return result;
	979	}
	980
	981
	982	/*************************************************************************
	983	DLAE2 computes the eigenvalues of a 2-by-2 symmetric matrix
	984	[ A B ]
	985	[ B C ].
	986	On return, RT1 is the eigenvalue of larger absolute value, and RT2
	987	is the eigenvalue of smaller absolute value.
	988
	989	-- LAPACK auxiliary routine (version 3.0) --
	990	Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
	991	Courant Institute, Argonne National Lab, and Rice University
	992	October 31, 1992
	993	*************************************************************************/
	994	private static void tdevde2(double a,
	995	double b,
	996	double c,
	997	ref double rt1,
	998	ref double rt2)
	999	{
	1000	double ab = 0;
	1001	double acmn = 0;
	1002	double acmx = 0;
	1003	double adf = 0;
	1004	double df = 0;
	1005	double rt = 0;
	1006	double sm = 0;
	1007	double tb = 0;
	1008
	1009	sm = a+c;
	1010	df = a-c;
	1011	adf = Math.Abs(df);
	1012	tb = b+b;
	1013	ab = Math.Abs(tb);
	1014	if( (double)(Math.Abs(a))>(double)(Math.Abs(c)) )
	1015	{
	1016	acmx = a;
	1017	acmn = c;
	1018	}
	1019	else
	1020	{
	1021	acmx = c;
	1022	acmn = a;
	1023	}
	1024	if( (double)(adf)>(double)(ab) )
	1025	{
	1026	rt = adf*Math.Sqrt(1+AP.Math.Sqr(ab/adf));
	1027	}
	1028	else
	1029	{
	1030	if( (double)(adf)<(double)(ab) )
	1031	{
	1032	rt = ab*Math.Sqrt(1+AP.Math.Sqr(adf/ab));
	1033	}
	1034	else
	1035	{
	1036
	1037	//
	1038	// Includes case AB=ADF=0
	1039	//
	1040	rt = ab*Math.Sqrt(2);
	1041	}
	1042	}
	1043	if( (double)(sm)<(double)(0) )
	1044	{
	1045	rt1 = 0.5*(sm-rt);
	1046
	1047	//
	1048	// Order of execution important.
	1049	// To get fully accurate smaller eigenvalue,
	1050	// next line needs to be executed in higher precision.
	1051	//
	1052	rt2 = acmx/rt1acmn-b/rt1b;
	1053	}
	1054	else
	1055	{
	1056	if( (double)(sm)>(double)(0) )
	1057	{
	1058	rt1 = 0.5*(sm+rt);
	1059
	1060	//
	1061	// Order of execution important.
	1062	// To get fully accurate smaller eigenvalue,
	1063	// next line needs to be executed in higher precision.
	1064	//
	1065	rt2 = acmx/rt1acmn-b/rt1b;
	1066	}
	1067	else
	1068	{
	1069
	1070	//
	1071	// Includes case RT1 = RT2 = 0
	1072	//
	1073	rt1 = 0.5*rt;
	1074	rt2 = -(0.5*rt);
	1075	}
	1076	}
	1077	}
	1078
	1079
	1080	/*************************************************************************
	1081	DLAEV2 computes the eigendecomposition of a 2-by-2 symmetric matrix
	1082
	1083	[ A B ]
	1084	[ B C ].
	1085
	1086	On return, RT1 is the eigenvalue of larger absolute value, RT2 is the
	1087	eigenvalue of smaller absolute value, and (CS1,SN1) is the unit right
	1088	eigenvector for RT1, giving the decomposition
	1089
	1090	[ CS1 SN1 ] [ A B ] [ CS1 -SN1 ] = [ RT1 0 ]
	1091	[-SN1 CS1 ] [ B C ] [ SN1 CS1 ] [ 0 RT2 ].
	1092
	1093
	1094	-- LAPACK auxiliary routine (version 3.0) --
	1095	Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
	1096	Courant Institute, Argonne National Lab, and Rice University
	1097	October 31, 1992
	1098	*************************************************************************/
	1099	private static void tdevdev2(double a,
	1100	double b,
	1101	double c,
	1102	ref double rt1,
	1103	ref double rt2,
	1104	ref double cs1,
	1105	ref double sn1)
	1106	{
	1107	int sgn1 = 0;
	1108	int sgn2 = 0;
	1109	double ab = 0;
	1110	double acmn = 0;
	1111	double acmx = 0;
	1112	double acs = 0;
	1113	double adf = 0;
	1114	double cs = 0;
	1115	double ct = 0;
	1116	double df = 0;
	1117	double rt = 0;
	1118	double sm = 0;
	1119	double tb = 0;
	1120	double tn = 0;
	1121
	1122
	1123	//
	1124	// Compute the eigenvalues
	1125	//
	1126	sm = a+c;
	1127	df = a-c;
	1128	adf = Math.Abs(df);
	1129	tb = b+b;
	1130	ab = Math.Abs(tb);
	1131	if( (double)(Math.Abs(a))>(double)(Math.Abs(c)) )
	1132	{
	1133	acmx = a;
	1134	acmn = c;
	1135	}
	1136	else
	1137	{
	1138	acmx = c;
	1139	acmn = a;
	1140	}
	1141	if( (double)(adf)>(double)(ab) )
	1142	{
	1143	rt = adf*Math.Sqrt(1+AP.Math.Sqr(ab/adf));
	1144	}
	1145	else
	1146	{
	1147	if( (double)(adf)<(double)(ab) )
	1148	{
	1149	rt = ab*Math.Sqrt(1+AP.Math.Sqr(adf/ab));
	1150	}
	1151	else
	1152	{
	1153
	1154	//
	1155	// Includes case AB=ADF=0
	1156	//
	1157	rt = ab*Math.Sqrt(2);
	1158	}
	1159	}
	1160	if( (double)(sm)<(double)(0) )
	1161	{
	1162	rt1 = 0.5*(sm-rt);
	1163	sgn1 = -1;
	1164
	1165	//
	1166	// Order of execution important.
	1167	// To get fully accurate smaller eigenvalue,
	1168	// next line needs to be executed in higher precision.
	1169	//
	1170	rt2 = acmx/rt1acmn-b/rt1b;
	1171	}
	1172	else
	1173	{
	1174	if( (double)(sm)>(double)(0) )
	1175	{
	1176	rt1 = 0.5*(sm+rt);
	1177	sgn1 = 1;
	1178
	1179	//
	1180	// Order of execution important.
	1181	// To get fully accurate smaller eigenvalue,
	1182	// next line needs to be executed in higher precision.
	1183	//
	1184	rt2 = acmx/rt1acmn-b/rt1b;
	1185	}
	1186	else
	1187	{
	1188
	1189	//
	1190	// Includes case RT1 = RT2 = 0
	1191	//
	1192	rt1 = 0.5*rt;
	1193	rt2 = -(0.5*rt);
	1194	sgn1 = 1;
	1195	}
	1196	}
	1197
	1198	//
	1199	// Compute the eigenvector
	1200	//
	1201	if( (double)(df)>=(double)(0) )
	1202	{
	1203	cs = df+rt;
	1204	sgn2 = 1;
	1205	}
	1206	else
	1207	{
	1208	cs = df-rt;
	1209	sgn2 = -1;
	1210	}
	1211	acs = Math.Abs(cs);
	1212	if( (double)(acs)>(double)(ab) )
	1213	{
	1214	ct = -(tb/cs);
	1215	sn1 = 1/Math.Sqrt(1+ct*ct);
	1216	cs1 = ct*sn1;
	1217	}
	1218	else
	1219	{
	1220	if( (double)(ab)==(double)(0) )
	1221	{
	1222	cs1 = 1;
	1223	sn1 = 0;
	1224	}
	1225	else
	1226	{
	1227	tn = -(cs/tb);
	1228	cs1 = 1/Math.Sqrt(1+tn*tn);
	1229	sn1 = tn*cs1;
	1230	}
	1231	}
	1232	if( sgn1==sgn2 )
	1233	{
	1234	tn = cs1;
	1235	cs1 = -sn1;
	1236	sn1 = tn;
	1237	}
	1238	}
	1239
	1240
	1241	/*************************************************************************
	1242	Internal routine
	1243	*************************************************************************/
	1244	private static double tdevdpythag(double a,
	1245	double b)
	1246	{
	1247	double result = 0;
	1248
	1249	if( (double)(Math.Abs(a))<(double)(Math.Abs(b)) )
	1250	{
	1251	result = Math.Abs(b)*Math.Sqrt(1+AP.Math.Sqr(a/b));
	1252	}
	1253	else
	1254	{
	1255	result = Math.Abs(a)*Math.Sqrt(1+AP.Math.Sqr(b/a));
	1256	}
	1257	return result;
	1258	}
	1259
	1260
	1261	/*************************************************************************
	1262	Internal routine
	1263	*************************************************************************/
	1264	private static double tdevdextsign(double a,
	1265	double b)
	1266	{
	1267	double result = 0;
	1268
	1269	if( (double)(b)>=(double)(0) )
	1270	{
	1271	result = Math.Abs(a);
	1272	}
	1273	else
	1274	{
	1275	result = -Math.Abs(a);
	1276	}
	1277	return result;
	1278	}
	1279	}
	1280	}

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

Download in other formats:

Update cookies preferences