Context Navigation

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.1.2.2591/ALGLIB-2.1.2.2591/hsschur.cs @ 2645

Visit:

Last change on this file since 2645 was 2645, checked in by mkommend, 14 years ago
extracted external libraries and adapted dependent plugins (ticket #837)
File size: 48.7 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 hsschur
32	{
33	/*************************************************************************
34	Subroutine performing the Schur decomposition of a matrix in upper
35	Hessenberg form using the QR algorithm with multiple shifts.
36
37	The source matrix H is represented as S'HS = T, where H - matrix in
38	upper Hessenberg form, S - orthogonal matrix (Schur vectors), T - upper
39	quasi-triangular matrix (with blocks of sizes 1x1 and 2x2 on the main
40	diagonal).
41
42	Input parameters:
43	H - matrix to be decomposed.
44	Array whose indexes range within [1..N, 1..N].
45	N - size of H, N>=0.
46
47
48	Output parameters:
49	H contains the matrix T.
50	Array whose indexes range within [1..N, 1..N].
51	All elements below the blocks on the main diagonal are equal
52	to 0.
53	S - contains Schur vectors.
54	Array whose indexes range within [1..N, 1..N].
55
56	Note 1:
57	The block structure of matrix T could be easily recognized: since all
58	the elements below the blocks are zeros, the elements a[i+1,i] which
59	are equal to 0 show the block border.
60
61	Note 2:
62	the algorithm performance depends on the value of the internal
63	parameter NS of InternalSchurDecomposition subroutine which defines
64	the number of shifts in the QR algorithm (analog of the block width
65	in block matrix algorithms in linear algebra). If you require maximum
66	performance on your machine, it is recommended to adjust this
67	parameter manually.
68
69	Result:
70	True, if the algorithm has converged and the parameters H and S contain
71	the result.
72	False, if the algorithm has not converged.
73
74	Algorithm implemented on the basis of subroutine DHSEQR (LAPACK 3.0 library).
75	*************************************************************************/
76	public static bool upperhessenbergschurdecomposition(ref double[,] h,
77	int n,
78	ref double[,] s)
79	{
80	bool result = new bool();
81	double[] wi = new double[0];
82	double[] wr = new double[0];
83	int info = 0;
84
85	internalschurdecomposition(ref h, n, 1, 2, ref wr, ref wi, ref s, ref info);
86	result = info==0;
87	return result;
88	}
89
90
91	public static void internalschurdecomposition(ref double[,] h,
92	int n,
93	int tneeded,
94	int zneeded,
95	ref double[] wr,
96	ref double[] wi,
97	ref double[,] z,
98	ref int info)
99	{
100	double[] work = new double[0];
101	int i = 0;
102	int i1 = 0;
103	int i2 = 0;
104	int ierr = 0;
105	int ii = 0;
106	int itemp = 0;
107	int itn = 0;
108	int its = 0;
109	int j = 0;
110	int k = 0;
111	int l = 0;
112	int maxb = 0;
113	int nr = 0;
114	int ns = 0;
115	int nv = 0;
116	double absw = 0;
117	double ovfl = 0;
118	double smlnum = 0;
119	double tau = 0;
120	double temp = 0;
121	double tst1 = 0;
122	double ulp = 0;
123	double unfl = 0;
124	double[,] s = new double[0,0];
125	double[] v = new double[0];
126	double[] vv = new double[0];
127	double[] workc1 = new double[0];
128	double[] works1 = new double[0];
129	double[] workv3 = new double[0];
130	double[] tmpwr = new double[0];
131	double[] tmpwi = new double[0];
132	bool initz = new bool();
133	bool wantt = new bool();
134	bool wantz = new bool();
135	double cnst = 0;
136	bool failflag = new bool();
137	int p1 = 0;
138	int p2 = 0;
139	double vt = 0;
140	int i_ = 0;
141	int i1_ = 0;
142
143
144	//
145	// Set the order of the multi-shift QR algorithm to be used.
146	// If you want to tune algorithm, change this values
147	//
148	ns = 12;
149	maxb = 50;
150
151	//
152	// Now 2 < NS <= MAXB < NH.
153	//
154	maxb = Math.Max(3, maxb);
155	ns = Math.Min(maxb, ns);
156
157	//
158	// Initialize
159	//
160	cnst = 1.5;
161	work = new double[Math.Max(n, 1)+1];
162	s = new double[ns+1, ns+1];
163	v = new double[ns+1+1];
164	vv = new double[ns+1+1];
165	wr = new double[Math.Max(n, 1)+1];
166	wi = new double[Math.Max(n, 1)+1];
167	workc1 = new double[1+1];
168	works1 = new double[1+1];
169	workv3 = new double[3+1];
170	tmpwr = new double[Math.Max(n, 1)+1];
171	tmpwi = new double[Math.Max(n, 1)+1];
172	System.Diagnostics.Debug.Assert(n>=0, "InternalSchurDecomposition: incorrect N!");
173	System.Diagnostics.Debug.Assert(tneeded==0 \| tneeded==1, "InternalSchurDecomposition: incorrect TNeeded!");
174	System.Diagnostics.Debug.Assert(zneeded==0 \| zneeded==1 \| zneeded==2, "InternalSchurDecomposition: incorrect ZNeeded!");
175	wantt = tneeded==1;
176	initz = zneeded==2;
177	wantz = zneeded!=0;
178	info = 0;
179
180	//
181	// Initialize Z, if necessary
182	//
183	if( initz )
184	{
185	z = new double[n+1, n+1];
186	for(i=1; i<=n; i++)
187	{
188	for(j=1; j<=n; j++)
189	{
190	if( i==j )
191	{
192	z[i,j] = 1;
193	}
194	else
195	{
196	z[i,j] = 0;
197	}
198	}
199	}
200	}
201
202	//
203	// Quick return if possible
204	//
205	if( n==0 )
206	{
207	return;
208	}
209	if( n==1 )
210	{
211	wr[1] = h[1,1];
212	wi[1] = 0;
213	return;
214	}
215
216	//
217	// Set rows and columns 1 to N to zero below the first
218	// subdiagonal.
219	//
220	for(j=1; j<=n-2; j++)
221	{
222	for(i=j+2; i<=n; i++)
223	{
224	h[i,j] = 0;
225	}
226	}
227
228	//
229	// Test if N is sufficiently small
230	//
231	if( ns<=2 \| ns>n \| maxb>=n )
232	{
233
234	//
235	// Use the standard double-shift algorithm
236	//
237	internalauxschur(wantt, wantz, n, 1, n, ref h, ref wr, ref wi, 1, n, ref z, ref work, ref workv3, ref workc1, ref works1, ref info);
238
239	//
240	// fill entries under diagonal blocks of T with zeros
241	//
242	if( wantt )
243	{
244	j = 1;
245	while( j<=n )
246	{
247	if( (double)(wi[j])==(double)(0) )
248	{
249	for(i=j+1; i<=n; i++)
250	{
251	h[i,j] = 0;
252	}
253	j = j+1;
254	}
255	else
256	{
257	for(i=j+2; i<=n; i++)
258	{
259	h[i,j] = 0;
260	h[i,j+1] = 0;
261	}
262	j = j+2;
263	}
264	}
265	}
266	return;
267	}
268	unfl = AP.Math.MinRealNumber;
269	ovfl = 1/unfl;
270	ulp = 2*AP.Math.MachineEpsilon;
271	smlnum = unfl*(n/ulp);
272
273	//
274	// I1 and I2 are the indices of the first row and last column of H
275	// to which transformations must be applied. If eigenvalues only are
276	// being computed, I1 and I2 are set inside the main loop.
277	//
278	i1 = 1;
279	i2 = n;
280
281	//
282	// ITN is the total number of multiple-shift QR iterations allowed.
283	//
284	itn = 30*n;
285
286	//
287	// The main loop begins here. I is the loop index and decreases from
288	// IHI to ILO in steps of at most MAXB. Each iteration of the loop
289	// works with the active submatrix in rows and columns L to I.
290	// Eigenvalues I+1 to IHI have already converged. Either L = ILO or
291	// H(L,L-1) is negligible so that the matrix splits.
292	//
293	i = n;
294	while( true )
295	{
296	l = 1;
297	if( i<1 )
298	{
299
300	//
301	// fill entries under diagonal blocks of T with zeros
302	//
303	if( wantt )
304	{
305	j = 1;
306	while( j<=n )
307	{
308	if( (double)(wi[j])==(double)(0) )
309	{
310	for(i=j+1; i<=n; i++)
311	{
312	h[i,j] = 0;
313	}
314	j = j+1;
315	}
316	else
317	{
318	for(i=j+2; i<=n; i++)
319	{
320	h[i,j] = 0;
321	h[i,j+1] = 0;
322	}
323	j = j+2;
324	}
325	}
326	}
327
328	//
329	// Exit
330	//
331	return;
332	}
333
334	//
335	// Perform multiple-shift QR iterations on rows and columns ILO to I
336	// until a submatrix of order at most MAXB splits off at the bottom
337	// because a subdiagonal element has become negligible.
338	//
339	failflag = true;
340	for(its=0; its<=itn; its++)
341	{
342
343	//
344	// Look for a single small subdiagonal element.
345	//
346	for(k=i; k>=l+1; k--)
347	{
348	tst1 = Math.Abs(h[k-1,k-1])+Math.Abs(h[k,k]);
349	if( (double)(tst1)==(double)(0) )
350	{
351	tst1 = blas.upperhessenberg1norm(ref h, l, i, l, i, ref work);
352	}
353	if( (double)(Math.Abs(h[k,k-1]))<=(double)(Math.Max(ulp*tst1, smlnum)) )
354	{
355	break;
356	}
357	}
358	l = k;
359	if( l>1 )
360	{
361
362	//
363	// H(L,L-1) is negligible.
364	//
365	h[l,l-1] = 0;
366	}
367
368	//
369	// Exit from loop if a submatrix of order <= MAXB has split off.
370	//
371	if( l>=i-maxb+1 )
372	{
373	failflag = false;
374	break;
375	}
376
377	//
378	// Now the active submatrix is in rows and columns L to I. If
379	// eigenvalues only are being computed, only the active submatrix
380	// need be transformed.
381	//
382	if( its==20 \| its==30 )
383	{
384
385	//
386	// Exceptional shifts.
387	//
388	for(ii=i-ns+1; ii<=i; ii++)
389	{
390	wr[ii] = cnst*(Math.Abs(h[ii,ii-1])+Math.Abs(h[ii,ii]));
391	wi[ii] = 0;
392	}
393	}
394	else
395	{
396
397	//
398	// Use eigenvalues of trailing submatrix of order NS as shifts.
399	//
400	blas.copymatrix(ref h, i-ns+1, i, i-ns+1, i, ref s, 1, ns, 1, ns);
401	internalauxschur(false, false, ns, 1, ns, ref s, ref tmpwr, ref tmpwi, 1, ns, ref z, ref work, ref workv3, ref workc1, ref works1, ref ierr);
402	for(p1=1; p1<=ns; p1++)
403	{
404	wr[i-ns+p1] = tmpwr[p1];
405	wi[i-ns+p1] = tmpwi[p1];
406	}
407	if( ierr>0 )
408	{
409
410	//
411	// If DLAHQR failed to compute all NS eigenvalues, use the
412	// unconverged diagonal elements as the remaining shifts.
413	//
414	for(ii=1; ii<=ierr; ii++)
415	{
416	wr[i-ns+ii] = s[ii,ii];
417	wi[i-ns+ii] = 0;
418	}
419	}
420	}
421
422	//
423	// Form the first column of (G-w(1)) (G-w(2)) . . . (G-w(ns))
424	// where G is the Hessenberg submatrix H(L:I,L:I) and w is
425	// the vector of shifts (stored in WR and WI). The result is
426	// stored in the local array V.
427	//
428	v[1] = 1;
429	for(ii=2; ii<=ns+1; ii++)
430	{
431	v[ii] = 0;
432	}
433	nv = 1;
434	for(j=i-ns+1; j<=i; j++)
435	{
436	if( (double)(wi[j])>=(double)(0) )
437	{
438	if( (double)(wi[j])==(double)(0) )
439	{
440
441	//
442	// real shift
443	//
444	p1 = nv+1;
445	for(i_=1; i_<=p1;i_++)
446	{
447	vv[i_] = v[i_];
448	}
449	blas.matrixvectormultiply(ref h, l, l+nv, l, l+nv-1, false, ref vv, 1, nv, 1.0, ref v, 1, nv+1, -wr[j]);
450	nv = nv+1;
451	}
452	else
453	{
454	if( (double)(wi[j])>(double)(0) )
455	{
456
457	//
458	// complex conjugate pair of shifts
459	//
460	p1 = nv+1;
461	for(i_=1; i_<=p1;i_++)
462	{
463	vv[i_] = v[i_];
464	}
465	blas.matrixvectormultiply(ref h, l, l+nv, l, l+nv-1, false, ref v, 1, nv, 1.0, ref vv, 1, nv+1, -(2*wr[j]));
466	itemp = blas.vectoridxabsmax(ref vv, 1, nv+1);
467	temp = 1/Math.Max(Math.Abs(vv[itemp]), smlnum);
468	p1 = nv+1;
469	for(i_=1; i_<=p1;i_++)
470	{
471	vv[i_] = temp*vv[i_];
472	}
473	absw = blas.pythag2(wr[j], wi[j]);
474	temp = tempabswabsw;
475	blas.matrixvectormultiply(ref h, l, l+nv+1, l, l+nv, false, ref vv, 1, nv+1, 1.0, ref v, 1, nv+2, temp);
476	nv = nv+2;
477	}
478	}
479
480	//
481	// Scale V(1:NV) so that max(abs(V(i))) = 1. If V is zero,
482	// reset it to the unit vector.
483	//
484	itemp = blas.vectoridxabsmax(ref v, 1, nv);
485	temp = Math.Abs(v[itemp]);
486	if( (double)(temp)==(double)(0) )
487	{
488	v[1] = 1;
489	for(ii=2; ii<=nv; ii++)
490	{
491	v[ii] = 0;
492	}
493	}
494	else
495	{
496	temp = Math.Max(temp, smlnum);
497	vt = 1/temp;
498	for(i_=1; i_<=nv;i_++)
499	{
500	v[i_] = vt*v[i_];
501	}
502	}
503	}
504	}
505
506	//
507	// Multiple-shift QR step
508	//
509	for(k=l; k<=i-1; k++)
510	{
511
512	//
513	// The first iteration of this loop determines a reflection G
514	// from the vector V and applies it from left and right to H,
515	// thus creating a nonzero bulge below the subdiagonal.
516	//
517	// Each subsequent iteration determines a reflection G to
518	// restore the Hessenberg form in the (K-1)th column, and thus
519	// chases the bulge one step toward the bottom of the active
520	// submatrix. NR is the order of G.
521	//
522	nr = Math.Min(ns+1, i-k+1);
523	if( k>l )
524	{
525	p1 = k-1;
526	p2 = k+nr-1;
527	i1_ = (k) - (1);
528	for(i_=1; i_<=nr;i_++)
529	{
530	v[i_] = h[i_+i1_,p1];
531	}
532	}
533	reflections.generatereflection(ref v, nr, ref tau);
534	if( k>l )
535	{
536	h[k,k-1] = v[1];
537	for(ii=k+1; ii<=i; ii++)
538	{
539	h[ii,k-1] = 0;
540	}
541	}
542	v[1] = 1;
543
544	//
545	// Apply G from the left to transform the rows of the matrix in
546	// columns K to I2.
547	//
548	reflections.applyreflectionfromtheleft(ref h, tau, ref v, k, k+nr-1, k, i2, ref work);
549
550	//
551	// Apply G from the right to transform the columns of the
552	// matrix in rows I1 to min(K+NR,I).
553	//
554	reflections.applyreflectionfromtheright(ref h, tau, ref v, i1, Math.Min(k+nr, i), k, k+nr-1, ref work);
555	if( wantz )
556	{
557
558	//
559	// Accumulate transformations in the matrix Z
560	//
561	reflections.applyreflectionfromtheright(ref z, tau, ref v, 1, n, k, k+nr-1, ref work);
562	}
563	}
564	}
565
566	//
567	// Failure to converge in remaining number of iterations
568	//
569	if( failflag )
570	{
571	info = i;
572	return;
573	}
574
575	//
576	// A submatrix of order <= MAXB in rows and columns L to I has split
577	// off. Use the double-shift QR algorithm to handle it.
578	//
579	internalauxschur(wantt, wantz, n, l, i, ref h, ref wr, ref wi, 1, n, ref z, ref work, ref workv3, ref workc1, ref works1, ref info);
580	if( info>0 )
581	{
582	return;
583	}
584
585	//
586	// Decrement number of remaining iterations, and return to start of
587	// the main loop with a new value of I.
588	//
589	itn = itn-its;
590	i = l-1;
591	}
592	}
593
594
595	private static void internalauxschur(bool wantt,
596	bool wantz,
597	int n,
598	int ilo,
599	int ihi,
600	ref double[,] h,
601	ref double[] wr,
602	ref double[] wi,
603	int iloz,
604	int ihiz,
605	ref double[,] z,
606	ref double[] work,
607	ref double[] workv3,
608	ref double[] workc1,
609	ref double[] works1,
610	ref int info)
611	{
612	int i = 0;
613	int i1 = 0;
614	int i2 = 0;
615	int itn = 0;
616	int its = 0;
617	int j = 0;
618	int k = 0;
619	int l = 0;
620	int m = 0;
621	int nh = 0;
622	int nr = 0;
623	int nz = 0;
624	double ave = 0;
625	double cs = 0;
626	double disc = 0;
627	double h00 = 0;
628	double h10 = 0;
629	double h11 = 0;
630	double h12 = 0;
631	double h21 = 0;
632	double h22 = 0;
633	double h33 = 0;
634	double h33s = 0;
635	double h43h34 = 0;
636	double h44 = 0;
637	double h44s = 0;
638	double ovfl = 0;
639	double s = 0;
640	double smlnum = 0;
641	double sn = 0;
642	double sum = 0;
643	double t1 = 0;
644	double t2 = 0;
645	double t3 = 0;
646	double tst1 = 0;
647	double unfl = 0;
648	double v1 = 0;
649	double v2 = 0;
650	double v3 = 0;
651	bool failflag = new bool();
652	double dat1 = 0;
653	double dat2 = 0;
654	int p1 = 0;
655	double him1im1 = 0;
656	double him1i = 0;
657	double hiim1 = 0;
658	double hii = 0;
659	double wrim1 = 0;
660	double wri = 0;
661	double wiim1 = 0;
662	double wii = 0;
663	double ulp = 0;
664
665	info = 0;
666	dat1 = 0.75;
667	dat2 = -0.4375;
668	ulp = AP.Math.MachineEpsilon;
669
670	//
671	// Quick return if possible
672	//
673	if( n==0 )
674	{
675	return;
676	}
677	if( ilo==ihi )
678	{
679	wr[ilo] = h[ilo,ilo];
680	wi[ilo] = 0;
681	return;
682	}
683	nh = ihi-ilo+1;
684	nz = ihiz-iloz+1;
685
686	//
687	// Set machine-dependent constants for the stopping criterion.
688	// If norm(H) <= sqrt(OVFL), overflow should not occur.
689	//
690	unfl = AP.Math.MinRealNumber;
691	ovfl = 1/unfl;
692	smlnum = unfl*(nh/ulp);
693
694	//
695	// I1 and I2 are the indices of the first row and last column of H
696	// to which transformations must be applied. If eigenvalues only are
697	// being computed, I1 and I2 are set inside the main loop.
698	//
699	i1 = 1;
700	i2 = n;
701
702	//
703	// ITN is the total number of QR iterations allowed.
704	//
705	itn = 30*nh;
706
707	//
708	// The main loop begins here. I is the loop index and decreases from
709	// IHI to ILO in steps of 1 or 2. Each iteration of the loop works
710	// with the active submatrix in rows and columns L to I.
711	// Eigenvalues I+1 to IHI have already converged. Either L = ILO or
712	// H(L,L-1) is negligible so that the matrix splits.
713	//
714	i = ihi;
715	while( true )
716	{
717	l = ilo;
718	if( i<ilo )
719	{
720	return;
721	}
722
723	//
724	// Perform QR iterations on rows and columns ILO to I until a
725	// submatrix of order 1 or 2 splits off at the bottom because a
726	// subdiagonal element has become negligible.
727	//
728	failflag = true;
729	for(its=0; its<=itn; its++)
730	{
731
732	//
733	// Look for a single small subdiagonal element.
734	//
735	for(k=i; k>=l+1; k--)
736	{
737	tst1 = Math.Abs(h[k-1,k-1])+Math.Abs(h[k,k]);
738	if( (double)(tst1)==(double)(0) )
739	{
740	tst1 = blas.upperhessenberg1norm(ref h, l, i, l, i, ref work);
741	}
742	if( (double)(Math.Abs(h[k,k-1]))<=(double)(Math.Max(ulp*tst1, smlnum)) )
743	{
744	break;
745	}
746	}
747	l = k;
748	if( l>ilo )
749	{
750
751	//
752	// H(L,L-1) is negligible
753	//
754	h[l,l-1] = 0;
755	}
756
757	//
758	// Exit from loop if a submatrix of order 1 or 2 has split off.
759	//
760	if( l>=i-1 )
761	{
762	failflag = false;
763	break;
764	}
765
766	//
767	// Now the active submatrix is in rows and columns L to I. If
768	// eigenvalues only are being computed, only the active submatrix
769	// need be transformed.
770	//
771	if( its==10 \| its==20 )
772	{
773
774	//
775	// Exceptional shift.
776	//
777	s = Math.Abs(h[i,i-1])+Math.Abs(h[i-1,i-2]);
778	h44 = dat1*s+h[i,i];
779	h33 = h44;
780	h43h34 = dat2ss;
781	}
782	else
783	{
784
785	//
786	// Prepare to use Francis' double shift
787	// (i.e. 2nd degree generalized Rayleigh quotient)
788	//
789	h44 = h[i,i];
790	h33 = h[i-1,i-1];
791	h43h34 = h[i,i-1]*h[i-1,i];
792	s = h[i-1,i-2]*h[i-1,i-2];
793	disc = (h33-h44)*0.5;
794	disc = disc*disc+h43h34;
795	if( (double)(disc)>(double)(0) )
796	{
797
798	//
799	// Real roots: use Wilkinson's shift twice
800	//
801	disc = Math.Sqrt(disc);
802	ave = 0.5*(h33+h44);
803	if( (double)(Math.Abs(h33)-Math.Abs(h44))>(double)(0) )
804	{
805	h33 = h33*h44-h43h34;
806	h44 = h33/(extschursign(disc, ave)+ave);
807	}
808	else
809	{
810	h44 = extschursign(disc, ave)+ave;
811	}
812	h33 = h44;
813	h43h34 = 0;
814	}
815	}
816
817	//
818	// Look for two consecutive small subdiagonal elements.
819	//
820	for(m=i-2; m>=l; m--)
821	{
822
823	//
824	// Determine the effect of starting the double-shift QR
825	// iteration at row M, and see if this would make H(M,M-1)
826	// negligible.
827	//
828	h11 = h[m,m];
829	h22 = h[m+1,m+1];
830	h21 = h[m+1,m];
831	h12 = h[m,m+1];
832	h44s = h44-h11;
833	h33s = h33-h11;
834	v1 = (h33s*h44s-h43h34)/h21+h12;
835	v2 = h22-h11-h33s-h44s;
836	v3 = h[m+2,m+1];
837	s = Math.Abs(v1)+Math.Abs(v2)+Math.Abs(v3);
838	v1 = v1/s;
839	v2 = v2/s;
840	v3 = v3/s;
841	workv3[1] = v1;
842	workv3[2] = v2;
843	workv3[3] = v3;
844	if( m==l )
845	{
846	break;
847	}
848	h00 = h[m-1,m-1];
849	h10 = h[m,m-1];
850	tst1 = Math.Abs(v1)*(Math.Abs(h00)+Math.Abs(h11)+Math.Abs(h22));
851	if( (double)(Math.Abs(h10)(Math.Abs(v2)+Math.Abs(v3)))<=(double)(ulptst1) )
852	{
853	break;
854	}
855	}
856
857	//
858	// Double-shift QR step
859	//
860	for(k=m; k<=i-1; k++)
861	{
862
863	//
864	// The first iteration of this loop determines a reflection G
865	// from the vector V and applies it from left and right to H,
866	// thus creating a nonzero bulge below the subdiagonal.
867	//
868	// Each subsequent iteration determines a reflection G to
869	// restore the Hessenberg form in the (K-1)th column, and thus
870	// chases the bulge one step toward the bottom of the active
871	// submatrix. NR is the order of G.
872	//
873	nr = Math.Min(3, i-k+1);
874	if( k>m )
875	{
876	for(p1=1; p1<=nr; p1++)
877	{
878	workv3[p1] = h[k+p1-1,k-1];
879	}
880	}
881	reflections.generatereflection(ref workv3, nr, ref t1);
882	if( k>m )
883	{
884	h[k,k-1] = workv3[1];
885	h[k+1,k-1] = 0;
886	if( k<i-1 )
887	{
888	h[k+2,k-1] = 0;
889	}
890	}
891	else
892	{
893	if( m>l )
894	{
895	h[k,k-1] = -h[k,k-1];
896	}
897	}
898	v2 = workv3[2];
899	t2 = t1*v2;
900	if( nr==3 )
901	{
902	v3 = workv3[3];
903	t3 = t1*v3;
904
905	//
906	// Apply G from the left to transform the rows of the matrix
907	// in columns K to I2.
908	//
909	for(j=k; j<=i2; j++)
910	{
911	sum = h[k,j]+v2h[k+1,j]+v3h[k+2,j];
912	h[k,j] = h[k,j]-sum*t1;
913	h[k+1,j] = h[k+1,j]-sum*t2;
914	h[k+2,j] = h[k+2,j]-sum*t3;
915	}
916
917	//
918	// Apply G from the right to transform the columns of the
919	// matrix in rows I1 to min(K+3,I).
920	//
921	for(j=i1; j<=Math.Min(k+3, i); j++)
922	{
923	sum = h[j,k]+v2h[j,k+1]+v3h[j,k+2];
924	h[j,k] = h[j,k]-sum*t1;
925	h[j,k+1] = h[j,k+1]-sum*t2;
926	h[j,k+2] = h[j,k+2]-sum*t3;
927	}
928	if( wantz )
929	{
930
931	//
932	// Accumulate transformations in the matrix Z
933	//
934	for(j=iloz; j<=ihiz; j++)
935	{
936	sum = z[j,k]+v2z[j,k+1]+v3z[j,k+2];
937	z[j,k] = z[j,k]-sum*t1;
938	z[j,k+1] = z[j,k+1]-sum*t2;
939	z[j,k+2] = z[j,k+2]-sum*t3;
940	}
941	}
942	}
943	else
944	{
945	if( nr==2 )
946	{
947
948	//
949	// Apply G from the left to transform the rows of the matrix
950	// in columns K to I2.
951	//
952	for(j=k; j<=i2; j++)
953	{
954	sum = h[k,j]+v2*h[k+1,j];
955	h[k,j] = h[k,j]-sum*t1;
956	h[k+1,j] = h[k+1,j]-sum*t2;
957	}
958
959	//
960	// Apply G from the right to transform the columns of the
961	// matrix in rows I1 to min(K+3,I).
962	//
963	for(j=i1; j<=i; j++)
964	{
965	sum = h[j,k]+v2*h[j,k+1];
966	h[j,k] = h[j,k]-sum*t1;
967	h[j,k+1] = h[j,k+1]-sum*t2;
968	}
969	if( wantz )
970	{
971
972	//
973	// Accumulate transformations in the matrix Z
974	//
975	for(j=iloz; j<=ihiz; j++)
976	{
977	sum = z[j,k]+v2*z[j,k+1];
978	z[j,k] = z[j,k]-sum*t1;
979	z[j,k+1] = z[j,k+1]-sum*t2;
980	}
981	}
982	}
983	}
984	}
985	}
986	if( failflag )
987	{
988
989	//
990	// Failure to converge in remaining number of iterations
991	//
992	info = i;
993	return;
994	}
995	if( l==i )
996	{
997
998	//
999	// H(I,I-1) is negligible: one eigenvalue has converged.
1000	//
1001	wr[i] = h[i,i];
1002	wi[i] = 0;
1003	}
1004	else
1005	{
1006	if( l==i-1 )
1007	{
1008
1009	//
1010	// H(I-1,I-2) is negligible: a pair of eigenvalues have converged.
1011	//
1012	// Transform the 2-by-2 submatrix to standard Schur form,
1013	// and compute and store the eigenvalues.
1014	//
1015	him1im1 = h[i-1,i-1];
1016	him1i = h[i-1,i];
1017	hiim1 = h[i,i-1];
1018	hii = h[i,i];
1019	aux2x2schur(ref him1im1, ref him1i, ref hiim1, ref hii, ref wrim1, ref wiim1, ref wri, ref wii, ref cs, ref sn);
1020	wr[i-1] = wrim1;
1021	wi[i-1] = wiim1;
1022	wr[i] = wri;
1023	wi[i] = wii;
1024	h[i-1,i-1] = him1im1;
1025	h[i-1,i] = him1i;
1026	h[i,i-1] = hiim1;
1027	h[i,i] = hii;
1028	if( wantt )
1029	{
1030
1031	//
1032	// Apply the transformation to the rest of H.
1033	//
1034	if( i2>i )
1035	{
1036	workc1[1] = cs;
1037	works1[1] = sn;
1038	rotations.applyrotationsfromtheleft(true, i-1, i, i+1, i2, ref workc1, ref works1, ref h, ref work);
1039	}
1040	workc1[1] = cs;
1041	works1[1] = sn;
1042	rotations.applyrotationsfromtheright(true, i1, i-2, i-1, i, ref workc1, ref works1, ref h, ref work);
1043	}
1044	if( wantz )
1045	{
1046
1047	//
1048	// Apply the transformation to Z.
1049	//
1050	workc1[1] = cs;
1051	works1[1] = sn;
1052	rotations.applyrotationsfromtheright(true, iloz, iloz+nz-1, i-1, i, ref workc1, ref works1, ref z, ref work);
1053	}
1054	}
1055	}
1056
1057	//
1058	// Decrement number of remaining iterations, and return to start of
1059	// the main loop with new value of I.
1060	//
1061	itn = itn-its;
1062	i = l-1;
1063	}
1064	}
1065
1066
1067	private static void aux2x2schur(ref double a,
1068	ref double b,
1069	ref double c,
1070	ref double d,
1071	ref double rt1r,
1072	ref double rt1i,
1073	ref double rt2r,
1074	ref double rt2i,
1075	ref double cs,
1076	ref double sn)
1077	{
1078	double multpl = 0;
1079	double aa = 0;
1080	double bb = 0;
1081	double bcmax = 0;
1082	double bcmis = 0;
1083	double cc = 0;
1084	double cs1 = 0;
1085	double dd = 0;
1086	double eps = 0;
1087	double p = 0;
1088	double sab = 0;
1089	double sac = 0;
1090	double scl = 0;
1091	double sigma = 0;
1092	double sn1 = 0;
1093	double tau = 0;
1094	double temp = 0;
1095	double z = 0;
1096
1097	multpl = 4.0;
1098	eps = AP.Math.MachineEpsilon;
1099	if( (double)(c)==(double)(0) )
1100	{
1101	cs = 1;
1102	sn = 0;
1103	}
1104	else
1105	{
1106	if( (double)(b)==(double)(0) )
1107	{
1108
1109	//
1110	// Swap rows and columns
1111	//
1112	cs = 0;
1113	sn = 1;
1114	temp = d;
1115	d = a;
1116	a = temp;
1117	b = -c;
1118	c = 0;
1119	}
1120	else
1121	{
1122	if( (double)(a-d)==(double)(0) & extschursigntoone(b)!=extschursigntoone(c) )
1123	{
1124	cs = 1;
1125	sn = 0;
1126	}
1127	else
1128	{
1129	temp = a-d;
1130	p = 0.5*temp;
1131	bcmax = Math.Max(Math.Abs(b), Math.Abs(c));
1132	bcmis = Math.Min(Math.Abs(b), Math.Abs(c))extschursigntoone(b)extschursigntoone(c);
1133	scl = Math.Max(Math.Abs(p), bcmax);
1134	z = p/sclp+bcmax/sclbcmis;
1135
1136	//
1137	// If Z is of the order of the machine accuracy, postpone the
1138	// decision on the nature of eigenvalues
1139	//
1140	if( (double)(z)>=(double)(multpl*eps) )
1141	{
1142
1143	//
1144	// Real eigenvalues. Compute A and D.
1145	//
1146	z = p+extschursign(Math.Sqrt(scl)*Math.Sqrt(z), p);
1147	a = d+z;
1148	d = d-bcmax/z*bcmis;
1149
1150	//
1151	// Compute B and the rotation matrix
1152	//
1153	tau = blas.pythag2(c, z);
1154	cs = z/tau;
1155	sn = c/tau;
1156	b = b-c;
1157	c = 0;
1158	}
1159	else
1160	{
1161
1162	//
1163	// Complex eigenvalues, or real (almost) equal eigenvalues.
1164	// Make diagonal elements equal.
1165	//
1166	sigma = b+c;
1167	tau = blas.pythag2(sigma, temp);
1168	cs = Math.Sqrt(0.5*(1+Math.Abs(sigma)/tau));
1169	sn = -(p/(taucs)extschursign(1, sigma));
1170
1171	//
1172	// Compute [ AA BB ] = [ A B ] [ CS -SN ]
1173	// [ CC DD ] [ C D ] [ SN CS ]
1174	//
1175	aa = acs+bsn;
1176	bb = -(asn)+bcs;
1177	cc = ccs+dsn;
1178	dd = -(csn)+dcs;
1179
1180	//
1181	// Compute [ A B ] = [ CS SN ] [ AA BB ]
1182	// [ C D ] [-SN CS ] [ CC DD ]
1183	//
1184	a = aacs+ccsn;
1185	b = bbcs+ddsn;
1186	c = -(aasn)+cccs;
1187	d = -(bbsn)+ddcs;
1188	temp = 0.5*(a+d);
1189	a = temp;
1190	d = temp;
1191	if( (double)(c)!=(double)(0) )
1192	{
1193	if( (double)(b)!=(double)(0) )
1194	{
1195	if( extschursigntoone(b)==extschursigntoone(c) )
1196	{
1197
1198	//
1199	// Real eigenvalues: reduce to upper triangular form
1200	//
1201	sab = Math.Sqrt(Math.Abs(b));
1202	sac = Math.Sqrt(Math.Abs(c));
1203	p = extschursign(sab*sac, c);
1204	tau = 1/Math.Sqrt(Math.Abs(b+c));
1205	a = temp+p;
1206	d = temp-p;
1207	b = b-c;
1208	c = 0;
1209	cs1 = sab*tau;
1210	sn1 = sac*tau;
1211	temp = cscs1-snsn1;
1212	sn = cssn1+sncs1;
1213	cs = temp;
1214	}
1215	}
1216	else
1217	{
1218	b = -c;
1219	c = 0;
1220	temp = cs;
1221	cs = -sn;
1222	sn = temp;
1223	}
1224	}
1225	}
1226	}
1227	}
1228	}
1229
1230	//
1231	// Store eigenvalues in (RT1R,RT1I) and (RT2R,RT2I).
1232	//
1233	rt1r = a;
1234	rt2r = d;
1235	if( (double)(c)==(double)(0) )
1236	{
1237	rt1i = 0;
1238	rt2i = 0;
1239	}
1240	else
1241	{
1242	rt1i = Math.Sqrt(Math.Abs(b))*Math.Sqrt(Math.Abs(c));
1243	rt2i = -rt1i;
1244	}
1245	}
1246
1247
1248	private static double extschursign(double a,
1249	double b)
1250	{
1251	double result = 0;
1252
1253	if( (double)(b)>=(double)(0) )
1254	{
1255	result = Math.Abs(a);
1256	}
1257	else
1258	{
1259	result = -Math.Abs(a);
1260	}
1261	return result;
1262	}
1263
1264
1265	private static int extschursigntoone(double b)
1266	{
1267	int result = 0;
1268
1269	if( (double)(b)>=(double)(0) )
1270	{
1271	result = 1;
1272	}
1273	else
1274	{
1275	result = -1;
1276	}
1277	return result;
1278	}
1279	}
1280	}

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

Download in other formats:

Update cookies preferences