Context Navigation

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.3.0/ALGLIB-2.3.0/tdevd.cs @ 2807

Visit:

Last change on this file since 2807 was 2806, checked in by gkronber, 15 years ago
Added plugin for new version of ALGLIB. #875 (Update ALGLIB sources)
File size: 43.8 KB

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

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

Download in other formats:

Update cookies preferences