Context Navigation

source: trunk/sources/ALGLIB/hsschur.cs @ 2635

Visit:

Last change on this file since 2635 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: 48.7 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 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