Context Navigation

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.1.2.2591/ALGLIB-2.1.2.2591/tdevd.cs @ 3031

Visit:

Last change on this file since 3031 was 2645, checked in by mkommend, 15 years ago
extracted external libraries and adapted dependent plugins (ticket #837)
File size: 43.9 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 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