Context Navigation

source: trunk/sources/ALGLIB/nsevd.cs @ 2575

Visit:

Last change on this file since 2575 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: 88.9 KB

Rev	Line
[2563]	1	/*************************************************************************
	2	Copyright (c) 1992-2007 The University of Tennessee. All rights reserved.
	3
	4	Contributors:
	5	* Sergey Bochkanov (ALGLIB project). Translation from FORTRAN to
	6	pseudocode.
	7
	8	See subroutines comments for additional copyrights.
	9
	10	>>> SOURCE LICENSE >>>
	11	This program is free software; you can redistribute it and/or modify
	12	it under the terms of the GNU General Public License as published by
	13	the Free Software Foundation (www.fsf.org); either version 2 of the
	14	License, or (at your option) any later version.
	15
	16	This program is distributed in the hope that it will be useful,
	17	but WITHOUT ANY WARRANTY; without even the implied warranty of
	18	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	19	GNU General Public License for more details.
	20
	21	A copy of the GNU General Public License is available at
	22	http://www.fsf.org/licensing/licenses
	23
	24	>>> END OF LICENSE >>>
	25	*************************************************************************/
	26
	27	using System;
	28
	29	namespace alglib
	30	{
	31	public class nsevd
	32	{
	33	/*************************************************************************
	34	Finding eigenvalues and eigenvectors of a general matrix
	35
	36	The algorithm finds eigenvalues and eigenvectors of a general matrix by
	37	using the QR algorithm with multiple shifts. The algorithm can find
	38	eigenvalues and both left and right eigenvectors.
	39
	40	The right eigenvector is a vector x such that Ax = wx, and the left
	41	eigenvector is a vector y such that y'A = wy' (here y' implies a complex
	42	conjugate transposition of vector y).
	43
	44	Input parameters:
	45	A - matrix. Array whose indexes range within [0..N-1, 0..N-1].
	46	N - size of matrix A.
	47	VNeeded - flag controlling whether eigenvectors are needed or not.
	48	If VNeeded is equal to:
	49	* 0, eigenvectors are not returned;
	50	* 1, right eigenvectors are returned;
	51	* 2, left eigenvectors are returned;
	52	* 3, both left and right eigenvectors are returned.
	53
	54	Output parameters:
	55	WR - real parts of eigenvalues.
	56	Array whose index ranges within [0..N-1].
	57	WR - imaginary parts of eigenvalues.
	58	Array whose index ranges within [0..N-1].
	59	VL, VR - arrays of left and right eigenvectors (if they are needed).
	60	If WI[i]=0, the respective eigenvalue is a real number,
	61	and it corresponds to the column number I of matrices VL/VR.
	62	If WI[i]>0, we have a pair of complex conjugate numbers with
	63	positive and negative imaginary parts:
	64	the first eigenvalue WR[i] + sqrt(-1)*WI[i];
	65	the second eigenvalue WR[i+1] + sqrt(-1)*WI[i+1];
	66	WI[i]>0
	67	WI[i+1] = -WI[i] < 0
	68	In that case, the eigenvector corresponding to the first
	69	eigenvalue is located in i and i+1 columns of matrices
	70	VL/VR (the column number i contains the real part, and the
	71	column number i+1 contains the imaginary part), and the vector
	72	corresponding to the second eigenvalue is a complex conjugate to
	73	the first vector.
	74	Arrays whose indexes range within [0..N-1, 0..N-1].
	75
	76	Result:
	77	True, if the algorithm has converged.
	78	False, if the algorithm has not converged.
	79
	80	Note 1:
	81	Some users may ask the following question: what if WI[N-1]>0?
	82	WI[N] must contain an eigenvalue which is complex conjugate to the
	83	N-th eigenvalue, but the array has only size N?
	84	The answer is as follows: such a situation cannot occur because the
	85	algorithm finds a pairs of eigenvalues, therefore, if WI[i]>0, I is
	86	strictly less than N-1.
	87
	88	Note 2:
	89	The algorithm performance depends on the value of the internal parameter
	90	NS of the InternalSchurDecomposition subroutine which defines the number
	91	of shifts in the QR algorithm (similarly to the block width in block-matrix
	92	algorithms of linear algebra). If you require maximum performance
	93	on your machine, it is recommended to adjust this parameter manually.
	94
	95
	96	See also the InternalTREVC subroutine.
	97
	98	The algorithm is based on the LAPACK 3.0 library.
	99	*************************************************************************/
	100	public static bool rmatrixevd(double[,] a,
	101	int n,
	102	int vneeded,
	103	ref double[] wr,
	104	ref double[] wi,
	105	ref double[,] vl,
	106	ref double[,] vr)
	107	{
	108	bool result = new bool();
	109	double[,] a1 = new double[0,0];
	110	double[,] vl1 = new double[0,0];
	111	double[,] vr1 = new double[0,0];
	112	double[] wr1 = new double[0];
	113	double[] wi1 = new double[0];
	114	int i = 0;
	115	double mx = 0;
	116	int i_ = 0;
	117	int i1_ = 0;
	118
	119	a = (double[,])a.Clone();
	120
	121	System.Diagnostics.Debug.Assert(vneeded>=0 & vneeded<=3, "RMatrixEVD: incorrect VNeeded!");
	122	a1 = new double[n+1, n+1];
	123	for(i=1; i<=n; i++)
	124	{
	125	i1_ = (0) - (1);
	126	for(i_=1; i_<=n;i_++)
	127	{
	128	a1[i,i_] = a[i-1,i_+i1_];
	129	}
	130	}
	131	result = nonsymmetricevd(a1, n, vneeded, ref wr1, ref wi1, ref vl1, ref vr1);
	132	if( result )
	133	{
	134	wr = new double[n-1+1];
	135	wi = new double[n-1+1];
	136	i1_ = (1) - (0);
	137	for(i_=0; i_<=n-1;i_++)
	138	{
	139	wr[i_] = wr1[i_+i1_];
	140	}
	141	i1_ = (1) - (0);
	142	for(i_=0; i_<=n-1;i_++)
	143	{
	144	wi[i_] = wi1[i_+i1_];
	145	}
	146	if( vneeded==2 \| vneeded==3 )
	147	{
	148	vl = new double[n-1+1, n-1+1];
	149	for(i=0; i<=n-1; i++)
	150	{
	151	i1_ = (1) - (0);
	152	for(i_=0; i_<=n-1;i_++)
	153	{
	154	vl[i,i_] = vl1[i+1,i_+i1_];
	155	}
	156	}
	157	}
	158	if( vneeded==1 \| vneeded==3 )
	159	{
	160	vr = new double[n-1+1, n-1+1];
	161	for(i=0; i<=n-1; i++)
	162	{
	163	i1_ = (1) - (0);
	164	for(i_=0; i_<=n-1;i_++)
	165	{
	166	vr[i,i_] = vr1[i+1,i_+i1_];
	167	}
	168	}
	169	}
	170	}
	171	return result;
	172	}
	173
	174
	175	public static bool nonsymmetricevd(double[,] a,
	176	int n,
	177	int vneeded,
	178	ref double[] wr,
	179	ref double[] wi,
	180	ref double[,] vl,
	181	ref double[,] vr)
	182	{
	183	bool result = new bool();
	184	double[,] s = new double[0,0];
	185	double[] tau = new double[0];
	186	bool[] sel = new bool[0];
	187	int i = 0;
	188	int info = 0;
	189	int m = 0;
	190	int i_ = 0;
	191
	192	a = (double[,])a.Clone();
	193
	194	System.Diagnostics.Debug.Assert(vneeded>=0 & vneeded<=3, "NonSymmetricEVD: incorrect VNeeded!");
	195	if( vneeded==0 )
	196	{
	197
	198	//
	199	// Eigen values only
	200	//
	201	hessenberg.toupperhessenberg(ref a, n, ref tau);
	202	hsschur.internalschurdecomposition(ref a, n, 0, 0, ref wr, ref wi, ref s, ref info);
	203	result = info==0;
	204	return result;
	205	}
	206
	207	//
	208	// Eigen values and vectors
	209	//
	210	hessenberg.toupperhessenberg(ref a, n, ref tau);
	211	hessenberg.unpackqfromupperhessenberg(ref a, n, ref tau, ref s);
	212	hsschur.internalschurdecomposition(ref a, n, 1, 1, ref wr, ref wi, ref s, ref info);
	213	result = info==0;
	214	if( !result )
	215	{
	216	return result;
	217	}
	218	if( vneeded==1 \| vneeded==3 )
	219	{
	220	vr = new double[n+1, n+1];
	221	for(i=1; i<=n; i++)
	222	{
	223	for(i_=1; i_<=n;i_++)
	224	{
	225	vr[i,i_] = s[i,i_];
	226	}
	227	}
	228	}
	229	if( vneeded==2 \| vneeded==3 )
	230	{
	231	vl = new double[n+1, n+1];
	232	for(i=1; i<=n; i++)
	233	{
	234	for(i_=1; i_<=n;i_++)
	235	{
	236	vl[i,i_] = s[i,i_];
	237	}
	238	}
	239	}
	240	internaltrevc(ref a, n, vneeded, 1, sel, ref vl, ref vr, ref m, ref info);
	241	result = info==0;
	242	return result;
	243	}
	244
	245
	246	private static void internaltrevc(ref double[,] t,
	247	int n,
	248	int side,
	249	int howmny,
	250	bool[] vselect,
	251	ref double[,] vl,
	252	ref double[,] vr,
	253	ref int m,
	254	ref int info)
	255	{
	256	bool allv = new bool();
	257	bool bothv = new bool();
	258	bool leftv = new bool();
	259	bool over = new bool();
	260	bool pair = new bool();
	261	bool rightv = new bool();
	262	bool somev = new bool();
	263	int i = 0;
	264	int ierr = 0;
	265	int ii = 0;
	266	int ip = 0;
	267	int iis = 0;
	268	int j = 0;
	269	int j1 = 0;
	270	int j2 = 0;
	271	int jnxt = 0;
	272	int k = 0;
	273	int ki = 0;
	274	int n2 = 0;
	275	double beta = 0;
	276	double bignum = 0;
	277	double emax = 0;
	278	double ovfl = 0;
	279	double rec = 0;
	280	double remax = 0;
	281	double scl = 0;
	282	double smin = 0;
	283	double smlnum = 0;
	284	double ulp = 0;
	285	double unfl = 0;
	286	double vcrit = 0;
	287	double vmax = 0;
	288	double wi = 0;
	289	double wr = 0;
	290	double xnorm = 0;
	291	double[,] x = new double[0,0];
	292	double[] work = new double[0];
	293	double[] temp = new double[0];
	294	double[,] temp11 = new double[0,0];
	295	double[,] temp22 = new double[0,0];
	296	double[,] temp11b = new double[0,0];
	297	double[,] temp21b = new double[0,0];
	298	double[,] temp12b = new double[0,0];
	299	double[,] temp22b = new double[0,0];
	300	bool skipflag = new bool();
	301	int k1 = 0;
	302	int k2 = 0;
	303	int k3 = 0;
	304	int k4 = 0;
	305	double vt = 0;
	306	bool[] rswap4 = new bool[0];
	307	bool[] zswap4 = new bool[0];
	308	int[,] ipivot44 = new int[0,0];
	309	double[] civ4 = new double[0];
	310	double[] crv4 = new double[0];
	311	int i_ = 0;
	312	int i1_ = 0;
	313
	314	vselect = (bool[])vselect.Clone();
	315
	316	x = new double[2+1, 2+1];
	317	temp11 = new double[1+1, 1+1];
	318	temp11b = new double[1+1, 1+1];
	319	temp21b = new double[2+1, 1+1];
	320	temp12b = new double[1+1, 2+1];
	321	temp22b = new double[2+1, 2+1];
	322	temp22 = new double[2+1, 2+1];
	323	work = new double[3*n+1];
	324	temp = new double[n+1];
	325	rswap4 = new bool[4+1];
	326	zswap4 = new bool[4+1];
	327	ipivot44 = new int[4+1, 4+1];
	328	civ4 = new double[4+1];
	329	crv4 = new double[4+1];
	330	if( howmny!=1 )
	331	{
	332	if( side==1 \| side==3 )
	333	{
	334	vr = new double[n+1, n+1];
	335	}
	336	if( side==2 \| side==3 )
	337	{
	338	vl = new double[n+1, n+1];
	339	}
	340	}
	341
	342	//
	343	// Decode and test the input parameters
	344	//
	345	bothv = side==3;
	346	rightv = side==1 \| bothv;
	347	leftv = side==2 \| bothv;
	348	allv = howmny==2;
	349	over = howmny==1;
	350	somev = howmny==3;
	351	info = 0;
	352	if( n<0 )
	353	{
	354	info = -2;
	355	return;
	356	}
	357	if( !rightv & !leftv )
	358	{
	359	info = -3;
	360	return;
	361	}
	362	if( !allv & !over & !somev )
	363	{
	364	info = -4;
	365	return;
	366	}
	367
	368	//
	369	// Set M to the number of columns required to store the selected
	370	// eigenvectors, standardize the array SELECT if necessary, and
	371	// test MM.
	372	//
	373	if( somev )
	374	{
	375	m = 0;
	376	pair = false;
	377	for(j=1; j<=n; j++)
	378	{
	379	if( pair )
	380	{
	381	pair = false;
	382	vselect[j] = false;
	383	}
	384	else
	385	{
	386	if( j<n )
	387	{
	388	if( (double)(t[j+1,j])==(double)(0) )
	389	{
	390	if( vselect[j] )
	391	{
	392	m = m+1;
	393	}
	394	}
	395	else
	396	{
	397	pair = true;
	398	if( vselect[j] \| vselect[j+1] )
	399	{
	400	vselect[j] = true;
	401	m = m+2;
	402	}
	403	}
	404	}
	405	else
	406	{
	407	if( vselect[n] )
	408	{
	409	m = m+1;
	410	}
	411	}
	412	}
	413	}
	414	}
	415	else
	416	{
	417	m = n;
	418	}
	419
	420	//
	421	// Quick return if possible.
	422	//
	423	if( n==0 )
	424	{
	425	return;
	426	}
	427
	428	//
	429	// Set the constants to control overflow.
	430	//
	431	unfl = AP.Math.MinRealNumber;
	432	ovfl = 1/unfl;
	433	ulp = AP.Math.MachineEpsilon;
	434	smlnum = unfl*(n/ulp);
	435	bignum = (1-ulp)/smlnum;
	436
	437	//
	438	// Compute 1-norm of each column of strictly upper triangular
	439	// part of T to control overflow in triangular solver.
	440	//
	441	work[1] = 0;
	442	for(j=2; j<=n; j++)
	443	{
	444	work[j] = 0;
	445	for(i=1; i<=j-1; i++)
	446	{
	447	work[j] = work[j]+Math.Abs(t[i,j]);
	448	}
	449	}
	450
	451	//
	452	// Index IP is used to specify the real or complex eigenvalue:
	453	// IP = 0, real eigenvalue,
	454	// 1, first of conjugate complex pair: (wr,wi)
	455	// -1, second of conjugate complex pair: (wr,wi)
	456	//
	457	n2 = 2*n;
	458	if( rightv )
	459	{
	460
	461	//
	462	// Compute right eigenvectors.
	463	//
	464	ip = 0;
	465	iis = m;
	466	for(ki=n; ki>=1; ki--)
	467	{
	468	skipflag = false;
	469	if( ip==1 )
	470	{
	471	skipflag = true;
	472	}
	473	else
	474	{
	475	if( ki!=1 )
	476	{
	477	if( (double)(t[ki,ki-1])!=(double)(0) )
	478	{
	479	ip = -1;
	480	}
	481	}
	482	if( somev )
	483	{
	484	if( ip==0 )
	485	{
	486	if( !vselect[ki] )
	487	{
	488	skipflag = true;
	489	}
	490	}
	491	else
	492	{
	493	if( !vselect[ki-1] )
	494	{
	495	skipflag = true;
	496	}
	497	}
	498	}
	499	}
	500	if( !skipflag )
	501	{
	502
	503	//
	504	// Compute the KI-th eigenvalue (WR,WI).
	505	//
	506	wr = t[ki,ki];
	507	wi = 0;
	508	if( ip!=0 )
	509	{
	510	wi = Math.Sqrt(Math.Abs(t[ki,ki-1]))*Math.Sqrt(Math.Abs(t[ki-1,ki]));
	511	}
	512	smin = Math.Max(ulp*(Math.Abs(wr)+Math.Abs(wi)), smlnum);
	513	if( ip==0 )
	514	{
	515
	516	//
	517	// Real right eigenvector
	518	//
	519	work[ki+n] = 1;
	520
	521	//
	522	// Form right-hand side
	523	//
	524	for(k=1; k<=ki-1; k++)
	525	{
	526	work[k+n] = -t[k,ki];
	527	}
	528
	529	//
	530	// Solve the upper quasi-triangular system:
	531	// (T(1:KI-1,1:KI-1) - WR)X = SCALEWORK.
	532	//
	533	jnxt = ki-1;
	534	for(j=ki-1; j>=1; j--)
	535	{
	536	if( j>jnxt )
	537	{
	538	continue;
	539	}
	540	j1 = j;
	541	j2 = j;
	542	jnxt = j-1;
	543	if( j>1 )
	544	{
	545	if( (double)(t[j,j-1])!=(double)(0) )
	546	{
	547	j1 = j-1;
	548	jnxt = j-2;
	549	}
	550	}
	551	if( j1==j2 )
	552	{
	553
	554	//
	555	// 1-by-1 diagonal block
	556	//
	557	temp11[1,1] = t[j,j];
	558	temp11b[1,1] = work[j+n];
	559	internalhsevdlaln2(false, 1, 1, smin, 1, ref temp11, 1.0, 1.0, ref temp11b, wr, 0.0, ref rswap4, ref zswap4, ref ipivot44, ref civ4, ref crv4, ref x, ref scl, ref xnorm, ref ierr);
	560
	561	//
	562	// Scale X(1,1) to avoid overflow when updating
	563	// the right-hand side.
	564	//
	565	if( (double)(xnorm)>(double)(1) )
	566	{
	567	if( (double)(work[j])>(double)(bignum/xnorm) )
	568	{
	569	x[1,1] = x[1,1]/xnorm;
	570	scl = scl/xnorm;
	571	}
	572	}
	573
	574	//
	575	// Scale if necessary
	576	//
	577	if( (double)(scl)!=(double)(1) )
	578	{
	579	k1 = n+1;
	580	k2 = n+ki;
	581	for(i_=k1; i_<=k2;i_++)
	582	{
	583	work[i_] = scl*work[i_];
	584	}
	585	}
	586	work[j+n] = x[1,1];
	587
	588	//
	589	// Update right-hand side
	590	//
	591	k1 = 1+n;
	592	k2 = j-1+n;
	593	k3 = j-1;
	594	vt = -x[1,1];
	595	i1_ = (1) - (k1);
	596	for(i_=k1; i_<=k2;i_++)
	597	{
	598	work[i_] = work[i_] + vt*t[i_+i1_,j];
	599	}
	600	}
	601	else
	602	{
	603
	604	//
	605	// 2-by-2 diagonal block
	606	//
	607	temp22[1,1] = t[j-1,j-1];
	608	temp22[1,2] = t[j-1,j];
	609	temp22[2,1] = t[j,j-1];
	610	temp22[2,2] = t[j,j];
	611	temp21b[1,1] = work[j-1+n];
	612	temp21b[2,1] = work[j+n];
	613	internalhsevdlaln2(false, 2, 1, smin, 1.0, ref temp22, 1.0, 1.0, ref temp21b, wr, 0, ref rswap4, ref zswap4, ref ipivot44, ref civ4, ref crv4, ref x, ref scl, ref xnorm, ref ierr);
	614
	615	//
	616	// Scale X(1,1) and X(2,1) to avoid overflow when
	617	// updating the right-hand side.
	618	//
	619	if( (double)(xnorm)>(double)(1) )
	620	{
	621	beta = Math.Max(work[j-1], work[j]);
	622	if( (double)(beta)>(double)(bignum/xnorm) )
	623	{
	624	x[1,1] = x[1,1]/xnorm;
	625	x[2,1] = x[2,1]/xnorm;
	626	scl = scl/xnorm;
	627	}
	628	}
	629
	630	//
	631	// Scale if necessary
	632	//
	633	if( (double)(scl)!=(double)(1) )
	634	{
	635	k1 = 1+n;
	636	k2 = ki+n;
	637	for(i_=k1; i_<=k2;i_++)
	638	{
	639	work[i_] = scl*work[i_];
	640	}
	641	}
	642	work[j-1+n] = x[1,1];
	643	work[j+n] = x[2,1];
	644
	645	//
	646	// Update right-hand side
	647	//
	648	k1 = 1+n;
	649	k2 = j-2+n;
	650	k3 = j-2;
	651	k4 = j-1;
	652	vt = -x[1,1];
	653	i1_ = (1) - (k1);
	654	for(i_=k1; i_<=k2;i_++)
	655	{
	656	work[i_] = work[i_] + vt*t[i_+i1_,k4];
	657	}
	658	vt = -x[2,1];
	659	i1_ = (1) - (k1);
	660	for(i_=k1; i_<=k2;i_++)
	661	{
	662	work[i_] = work[i_] + vt*t[i_+i1_,j];
	663	}
	664	}
	665	}
	666
	667	//
	668	// Copy the vector x or Q*x to VR and normalize.
	669	//
	670	if( !over )
	671	{
	672	k1 = 1+n;
	673	k2 = ki+n;
	674	i1_ = (k1) - (1);
	675	for(i_=1; i_<=ki;i_++)
	676	{
	677	vr[i_,iis] = work[i_+i1_];
	678	}
	679	ii = blas.columnidxabsmax(ref vr, 1, ki, iis);
	680	remax = 1/Math.Abs(vr[ii,iis]);
	681	for(i_=1; i_<=ki;i_++)
	682	{
	683	vr[i_,iis] = remax*vr[i_,iis];
	684	}
	685	for(k=ki+1; k<=n; k++)
	686	{
	687	vr[k,iis] = 0;
	688	}
	689	}
	690	else
	691	{
	692	if( ki>1 )
	693	{
	694	for(i_=1; i_<=n;i_++)
	695	{
	696	temp[i_] = vr[i_,ki];
	697	}
	698	blas.matrixvectormultiply(ref vr, 1, n, 1, ki-1, false, ref work, 1+n, ki-1+n, 1.0, ref temp, 1, n, work[ki+n]);
	699	for(i_=1; i_<=n;i_++)
	700	{
	701	vr[i_,ki] = temp[i_];
	702	}
	703	}
	704	ii = blas.columnidxabsmax(ref vr, 1, n, ki);
	705	remax = 1/Math.Abs(vr[ii,ki]);
	706	for(i_=1; i_<=n;i_++)
	707	{
	708	vr[i_,ki] = remax*vr[i_,ki];
	709	}
	710	}
	711	}
	712	else
	713	{
	714
	715	//
	716	// Complex right eigenvector.
	717	//
	718	// Initial solve
	719	// [ (T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I* WI)]*X = 0.
	720	// [ (T(KI,KI-1) T(KI,KI) ) ]
	721	//
	722	if( (double)(Math.Abs(t[ki-1,ki]))>=(double)(Math.Abs(t[ki,ki-1])) )
	723	{
	724	work[ki-1+n] = 1;
	725	work[ki+n2] = wi/t[ki-1,ki];
	726	}
	727	else
	728	{
	729	work[ki-1+n] = -(wi/t[ki,ki-1]);
	730	work[ki+n2] = 1;
	731	}
	732	work[ki+n] = 0;
	733	work[ki-1+n2] = 0;
	734
	735	//
	736	// Form right-hand side
	737	//
	738	for(k=1; k<=ki-2; k++)
	739	{
	740	work[k+n] = -(work[ki-1+n]*t[k,ki-1]);
	741	work[k+n2] = -(work[ki+n2]*t[k,ki]);
	742	}
	743
	744	//
	745	// Solve upper quasi-triangular system:
	746	// (T(1:KI-2,1:KI-2) - (WR+iWI))X = SCALE(WORK+iWORK2)
	747	//
	748	jnxt = ki-2;
	749	for(j=ki-2; j>=1; j--)
	750	{
	751	if( j>jnxt )
	752	{
	753	continue;
	754	}
	755	j1 = j;
	756	j2 = j;
	757	jnxt = j-1;
	758	if( j>1 )
	759	{
	760	if( (double)(t[j,j-1])!=(double)(0) )
	761	{
	762	j1 = j-1;
	763	jnxt = j-2;
	764	}
	765	}
	766	if( j1==j2 )
	767	{
	768
	769	//
	770	// 1-by-1 diagonal block
	771	//
	772	temp11[1,1] = t[j,j];
	773	temp12b[1,1] = work[j+n];
	774	temp12b[1,2] = work[j+n+n];
	775	internalhsevdlaln2(false, 1, 2, smin, 1.0, ref temp11, 1.0, 1.0, ref temp12b, wr, wi, ref rswap4, ref zswap4, ref ipivot44, ref civ4, ref crv4, ref x, ref scl, ref xnorm, ref ierr);
	776
	777	//
	778	// Scale X(1,1) and X(1,2) to avoid overflow when
	779	// updating the right-hand side.
	780	//
	781	if( (double)(xnorm)>(double)(1) )
	782	{
	783	if( (double)(work[j])>(double)(bignum/xnorm) )
	784	{
	785	x[1,1] = x[1,1]/xnorm;
	786	x[1,2] = x[1,2]/xnorm;
	787	scl = scl/xnorm;
	788	}
	789	}
	790
	791	//
	792	// Scale if necessary
	793	//
	794	if( (double)(scl)!=(double)(1) )
	795	{
	796	k1 = 1+n;
	797	k2 = ki+n;
	798	for(i_=k1; i_<=k2;i_++)
	799	{
	800	work[i_] = scl*work[i_];
	801	}
	802	k1 = 1+n2;
	803	k2 = ki+n2;
	804	for(i_=k1; i_<=k2;i_++)
	805	{
	806	work[i_] = scl*work[i_];
	807	}
	808	}
	809	work[j+n] = x[1,1];
	810	work[j+n2] = x[1,2];
	811
	812	//
	813	// Update the right-hand side
	814	//
	815	k1 = 1+n;
	816	k2 = j-1+n;
	817	k3 = 1;
	818	k4 = j-1;
	819	vt = -x[1,1];
	820	i1_ = (k3) - (k1);
	821	for(i_=k1; i_<=k2;i_++)
	822	{
	823	work[i_] = work[i_] + vt*t[i_+i1_,j];
	824	}
	825	k1 = 1+n2;
	826	k2 = j-1+n2;
	827	k3 = 1;
	828	k4 = j-1;
	829	vt = -x[1,2];
	830	i1_ = (k3) - (k1);
	831	for(i_=k1; i_<=k2;i_++)
	832	{
	833	work[i_] = work[i_] + vt*t[i_+i1_,j];
	834	}
	835	}
	836	else
	837	{
	838
	839	//
	840	// 2-by-2 diagonal block
	841	//
	842	temp22[1,1] = t[j-1,j-1];
	843	temp22[1,2] = t[j-1,j];
	844	temp22[2,1] = t[j,j-1];
	845	temp22[2,2] = t[j,j];
	846	temp22b[1,1] = work[j-1+n];
	847	temp22b[1,2] = work[j-1+n+n];
	848	temp22b[2,1] = work[j+n];
	849	temp22b[2,2] = work[j+n+n];
	850	internalhsevdlaln2(false, 2, 2, smin, 1.0, ref temp22, 1.0, 1.0, ref temp22b, wr, wi, ref rswap4, ref zswap4, ref ipivot44, ref civ4, ref crv4, ref x, ref scl, ref xnorm, ref ierr);
	851
	852	//
	853	// Scale X to avoid overflow when updating
	854	// the right-hand side.
	855	//
	856	if( (double)(xnorm)>(double)(1) )
	857	{
	858	beta = Math.Max(work[j-1], work[j]);
	859	if( (double)(beta)>(double)(bignum/xnorm) )
	860	{
	861	rec = 1/xnorm;
	862	x[1,1] = x[1,1]*rec;
	863	x[1,2] = x[1,2]*rec;
	864	x[2,1] = x[2,1]*rec;
	865	x[2,2] = x[2,2]*rec;
	866	scl = scl*rec;
	867	}
	868	}
	869
	870	//
	871	// Scale if necessary
	872	//
	873	if( (double)(scl)!=(double)(1) )
	874	{
	875	for(i_=1+n; i_<=ki+n;i_++)
	876	{
	877	work[i_] = scl*work[i_];
	878	}
	879	for(i_=1+n2; i_<=ki+n2;i_++)
	880	{
	881	work[i_] = scl*work[i_];
	882	}
	883	}
	884	work[j-1+n] = x[1,1];
	885	work[j+n] = x[2,1];
	886	work[j-1+n2] = x[1,2];
	887	work[j+n2] = x[2,2];
	888
	889	//
	890	// Update the right-hand side
	891	//
	892	vt = -x[1,1];
	893	i1_ = (1) - (n+1);
	894	for(i_=n+1; i_<=n+j-2;i_++)
	895	{
	896	work[i_] = work[i_] + vt*t[i_+i1_,j-1];
	897	}
	898	vt = -x[2,1];
	899	i1_ = (1) - (n+1);
	900	for(i_=n+1; i_<=n+j-2;i_++)
	901	{
	902	work[i_] = work[i_] + vt*t[i_+i1_,j];
	903	}
	904	vt = -x[1,2];
	905	i1_ = (1) - (n2+1);
	906	for(i_=n2+1; i_<=n2+j-2;i_++)
	907	{
	908	work[i_] = work[i_] + vt*t[i_+i1_,j-1];
	909	}
	910	vt = -x[2,2];
	911	i1_ = (1) - (n2+1);
	912	for(i_=n2+1; i_<=n2+j-2;i_++)
	913	{
	914	work[i_] = work[i_] + vt*t[i_+i1_,j];
	915	}
	916	}
	917	}
	918
	919	//
	920	// Copy the vector x or Q*x to VR and normalize.
	921	//
	922	if( !over )
	923	{
	924	i1_ = (n+1) - (1);
	925	for(i_=1; i_<=ki;i_++)
	926	{
	927	vr[i_,iis-1] = work[i_+i1_];
	928	}
	929	i1_ = (n2+1) - (1);
	930	for(i_=1; i_<=ki;i_++)
	931	{
	932	vr[i_,iis] = work[i_+i1_];
	933	}
	934	emax = 0;
	935	for(k=1; k<=ki; k++)
	936	{
	937	emax = Math.Max(emax, Math.Abs(vr[k,iis-1])+Math.Abs(vr[k,iis]));
	938	}
	939	remax = 1/emax;
	940	for(i_=1; i_<=ki;i_++)
	941	{
	942	vr[i_,iis-1] = remax*vr[i_,iis-1];
	943	}
	944	for(i_=1; i_<=ki;i_++)
	945	{
	946	vr[i_,iis] = remax*vr[i_,iis];
	947	}
	948	for(k=ki+1; k<=n; k++)
	949	{
	950	vr[k,iis-1] = 0;
	951	vr[k,iis] = 0;
	952	}
	953	}
	954	else
	955	{
	956	if( ki>2 )
	957	{
	958	for(i_=1; i_<=n;i_++)
	959	{
	960	temp[i_] = vr[i_,ki-1];
	961	}
	962	blas.matrixvectormultiply(ref vr, 1, n, 1, ki-2, false, ref work, 1+n, ki-2+n, 1.0, ref temp, 1, n, work[ki-1+n]);
	963	for(i_=1; i_<=n;i_++)
	964	{
	965	vr[i_,ki-1] = temp[i_];
	966	}
	967	for(i_=1; i_<=n;i_++)
	968	{
	969	temp[i_] = vr[i_,ki];
	970	}
	971	blas.matrixvectormultiply(ref vr, 1, n, 1, ki-2, false, ref work, 1+n2, ki-2+n2, 1.0, ref temp, 1, n, work[ki+n2]);
	972	for(i_=1; i_<=n;i_++)
	973	{
	974	vr[i_,ki] = temp[i_];
	975	}
	976	}
	977	else
	978	{
	979	vt = work[ki-1+n];
	980	for(i_=1; i_<=n;i_++)
	981	{
	982	vr[i_,ki-1] = vt*vr[i_,ki-1];
	983	}
	984	vt = work[ki+n2];
	985	for(i_=1; i_<=n;i_++)
	986	{
	987	vr[i_,ki] = vt*vr[i_,ki];
	988	}
	989	}
	990	emax = 0;
	991	for(k=1; k<=n; k++)
	992	{
	993	emax = Math.Max(emax, Math.Abs(vr[k,ki-1])+Math.Abs(vr[k,ki]));
	994	}
	995	remax = 1/emax;
	996	for(i_=1; i_<=n;i_++)
	997	{
	998	vr[i_,ki-1] = remax*vr[i_,ki-1];
	999	}
	1000	for(i_=1; i_<=n;i_++)
	1001	{
	1002	vr[i_,ki] = remax*vr[i_,ki];
	1003	}
	1004	}
	1005	}
	1006	iis = iis-1;
	1007	if( ip!=0 )
	1008	{
	1009	iis = iis-1;
	1010	}
	1011	}
	1012	if( ip==1 )
	1013	{
	1014	ip = 0;
	1015	}
	1016	if( ip==-1 )
	1017	{
	1018	ip = 1;
	1019	}
	1020	}
	1021	}
	1022	if( leftv )
	1023	{
	1024
	1025	//
	1026	// Compute left eigenvectors.
	1027	//
	1028	ip = 0;
	1029	iis = 1;
	1030	for(ki=1; ki<=n; ki++)
	1031	{
	1032	skipflag = false;
	1033	if( ip==-1 )
	1034	{
	1035	skipflag = true;
	1036	}
	1037	else
	1038	{
	1039	if( ki!=n )
	1040	{
	1041	if( (double)(t[ki+1,ki])!=(double)(0) )
	1042	{
	1043	ip = 1;
	1044	}
	1045	}
	1046	if( somev )
	1047	{
	1048	if( !vselect[ki] )
	1049	{
	1050	skipflag = true;
	1051	}
	1052	}
	1053	}
	1054	if( !skipflag )
	1055	{
	1056
	1057	//
	1058	// Compute the KI-th eigenvalue (WR,WI).
	1059	//
	1060	wr = t[ki,ki];
	1061	wi = 0;
	1062	if( ip!=0 )
	1063	{
	1064	wi = Math.Sqrt(Math.Abs(t[ki,ki+1]))*Math.Sqrt(Math.Abs(t[ki+1,ki]));
	1065	}
	1066	smin = Math.Max(ulp*(Math.Abs(wr)+Math.Abs(wi)), smlnum);
	1067	if( ip==0 )
	1068	{
	1069
	1070	//
	1071	// Real left eigenvector.
	1072	//
	1073	work[ki+n] = 1;
	1074
	1075	//
	1076	// Form right-hand side
	1077	//
	1078	for(k=ki+1; k<=n; k++)
	1079	{
	1080	work[k+n] = -t[ki,k];
	1081	}
	1082
	1083	//
	1084	// Solve the quasi-triangular system:
	1085	// (T(KI+1:N,KI+1:N) - WR)'X = SCALEWORK
	1086	//
	1087	vmax = 1;
	1088	vcrit = bignum;
	1089	jnxt = ki+1;
	1090	for(j=ki+1; j<=n; j++)
	1091	{
	1092	if( j<jnxt )
	1093	{
	1094	continue;
	1095	}
	1096	j1 = j;
	1097	j2 = j;
	1098	jnxt = j+1;
	1099	if( j<n )
	1100	{
	1101	if( (double)(t[j+1,j])!=(double)(0) )
	1102	{
	1103	j2 = j+1;
	1104	jnxt = j+2;
	1105	}
	1106	}
	1107	if( j1==j2 )
	1108	{
	1109
	1110	//
	1111	// 1-by-1 diagonal block
	1112	//
	1113	// Scale if necessary to avoid overflow when forming
	1114	// the right-hand side.
	1115	//
	1116	if( (double)(work[j])>(double)(vcrit) )
	1117	{
	1118	rec = 1/vmax;
	1119	for(i_=ki+n; i_<=n+n;i_++)
	1120	{
	1121	work[i_] = rec*work[i_];
	1122	}
	1123	vmax = 1;
	1124	vcrit = bignum;
	1125	}
	1126	i1_ = (ki+1+n)-(ki+1);
	1127	vt = 0.0;
	1128	for(i_=ki+1; i_<=j-1;i_++)
	1129	{
	1130	vt += t[i_,j]*work[i_+i1_];
	1131	}
	1132	work[j+n] = work[j+n]-vt;
	1133
	1134	//
	1135	// Solve (T(J,J)-WR)'*X = WORK
	1136	//
	1137	temp11[1,1] = t[j,j];
	1138	temp11b[1,1] = work[j+n];
	1139	internalhsevdlaln2(false, 1, 1, smin, 1.0, ref temp11, 1.0, 1.0, ref temp11b, wr, 0, ref rswap4, ref zswap4, ref ipivot44, ref civ4, ref crv4, ref x, ref scl, ref xnorm, ref ierr);
	1140
	1141	//
	1142	// Scale if necessary
	1143	//
	1144	if( (double)(scl)!=(double)(1) )
	1145	{
	1146	for(i_=ki+n; i_<=n+n;i_++)
	1147	{
	1148	work[i_] = scl*work[i_];
	1149	}
	1150	}
	1151	work[j+n] = x[1,1];
	1152	vmax = Math.Max(Math.Abs(work[j+n]), vmax);
	1153	vcrit = bignum/vmax;
	1154	}
	1155	else
	1156	{
	1157
	1158	//
	1159	// 2-by-2 diagonal block
	1160	//
	1161	// Scale if necessary to avoid overflow when forming
	1162	// the right-hand side.
	1163	//
	1164	beta = Math.Max(work[j], work[j+1]);
	1165	if( (double)(beta)>(double)(vcrit) )
	1166	{
	1167	rec = 1/vmax;
	1168	for(i_=ki+n; i_<=n+n;i_++)
	1169	{
	1170	work[i_] = rec*work[i_];
	1171	}
	1172	vmax = 1;
	1173	vcrit = bignum;
	1174	}
	1175	i1_ = (ki+1+n)-(ki+1);
	1176	vt = 0.0;
	1177	for(i_=ki+1; i_<=j-1;i_++)
	1178	{
	1179	vt += t[i_,j]*work[i_+i1_];
	1180	}
	1181	work[j+n] = work[j+n]-vt;
	1182	i1_ = (ki+1+n)-(ki+1);
	1183	vt = 0.0;
	1184	for(i_=ki+1; i_<=j-1;i_++)
	1185	{
	1186	vt += t[i_,j+1]*work[i_+i1_];
	1187	}
	1188	work[j+1+n] = work[j+1+n]-vt;
	1189
	1190	//
	1191	// Solve
	1192	// [T(J,J)-WR T(J,J+1) ]'* X = SCALE*( WORK1 )
	1193	// [T(J+1,J) T(J+1,J+1)-WR] ( WORK2 )
	1194	//
	1195	temp22[1,1] = t[j,j];
	1196	temp22[1,2] = t[j,j+1];
	1197	temp22[2,1] = t[j+1,j];
	1198	temp22[2,2] = t[j+1,j+1];
	1199	temp21b[1,1] = work[j+n];
	1200	temp21b[2,1] = work[j+1+n];
	1201	internalhsevdlaln2(true, 2, 1, smin, 1.0, ref temp22, 1.0, 1.0, ref temp21b, wr, 0, ref rswap4, ref zswap4, ref ipivot44, ref civ4, ref crv4, ref x, ref scl, ref xnorm, ref ierr);
	1202
	1203	//
	1204	// Scale if necessary
	1205	//
	1206	if( (double)(scl)!=(double)(1) )
	1207	{
	1208	for(i_=ki+n; i_<=n+n;i_++)
	1209	{
	1210	work[i_] = scl*work[i_];
	1211	}
	1212	}
	1213	work[j+n] = x[1,1];
	1214	work[j+1+n] = x[2,1];
	1215	vmax = Math.Max(Math.Abs(work[j+n]), Math.Max(Math.Abs(work[j+1+n]), vmax));
	1216	vcrit = bignum/vmax;
	1217	}
	1218	}
	1219
	1220	//
	1221	// Copy the vector x or Q*x to VL and normalize.
	1222	//
	1223	if( !over )
	1224	{
	1225	i1_ = (ki+n) - (ki);
	1226	for(i_=ki; i_<=n;i_++)
	1227	{
	1228	vl[i_,iis] = work[i_+i1_];
	1229	}
	1230	ii = blas.columnidxabsmax(ref vl, ki, n, iis);
	1231	remax = 1/Math.Abs(vl[ii,iis]);
	1232	for(i_=ki; i_<=n;i_++)
	1233	{
	1234	vl[i_,iis] = remax*vl[i_,iis];
	1235	}
	1236	for(k=1; k<=ki-1; k++)
	1237	{
	1238	vl[k,iis] = 0;
	1239	}
	1240	}
	1241	else
	1242	{
	1243	if( ki<n )
	1244	{
	1245	for(i_=1; i_<=n;i_++)
	1246	{
	1247	temp[i_] = vl[i_,ki];
	1248	}
	1249	blas.matrixvectormultiply(ref vl, 1, n, ki+1, n, false, ref work, ki+1+n, n+n, 1.0, ref temp, 1, n, work[ki+n]);
	1250	for(i_=1; i_<=n;i_++)
	1251	{
	1252	vl[i_,ki] = temp[i_];
	1253	}
	1254	}
	1255	ii = blas.columnidxabsmax(ref vl, 1, n, ki);
	1256	remax = 1/Math.Abs(vl[ii,ki]);
	1257	for(i_=1; i_<=n;i_++)
	1258	{
	1259	vl[i_,ki] = remax*vl[i_,ki];
	1260	}
	1261	}
	1262	}
	1263	else
	1264	{
	1265
	1266	//
	1267	// Complex left eigenvector.
	1268	//
	1269	// Initial solve:
	1270	// ((T(KI,KI) T(KI,KI+1) )' - (WR - I* WI))*X = 0.
	1271	// ((T(KI+1,KI) T(KI+1,KI+1)) )
	1272	//
	1273	if( (double)(Math.Abs(t[ki,ki+1]))>=(double)(Math.Abs(t[ki+1,ki])) )
	1274	{
	1275	work[ki+n] = wi/t[ki,ki+1];
	1276	work[ki+1+n2] = 1;
	1277	}
	1278	else
	1279	{
	1280	work[ki+n] = 1;
	1281	work[ki+1+n2] = -(wi/t[ki+1,ki]);
	1282	}
	1283	work[ki+1+n] = 0;
	1284	work[ki+n2] = 0;
	1285
	1286	//
	1287	// Form right-hand side
	1288	//
	1289	for(k=ki+2; k<=n; k++)
	1290	{
	1291	work[k+n] = -(work[ki+n]*t[ki,k]);
	1292	work[k+n2] = -(work[ki+1+n2]*t[ki+1,k]);
	1293	}
	1294
	1295	//
	1296	// Solve complex quasi-triangular system:
	1297	// ( T(KI+2,N:KI+2,N) - (WR-iWI) )X = WORK1+i*WORK2
	1298	//
	1299	vmax = 1;
	1300	vcrit = bignum;
	1301	jnxt = ki+2;
	1302	for(j=ki+2; j<=n; j++)
	1303	{
	1304	if( j<jnxt )
	1305	{
	1306	continue;
	1307	}
	1308	j1 = j;
	1309	j2 = j;
	1310	jnxt = j+1;
	1311	if( j<n )
	1312	{
	1313	if( (double)(t[j+1,j])!=(double)(0) )
	1314	{
	1315	j2 = j+1;
	1316	jnxt = j+2;
	1317	}
	1318	}
	1319	if( j1==j2 )
	1320	{
	1321
	1322	//
	1323	// 1-by-1 diagonal block
	1324	//
	1325	// Scale if necessary to avoid overflow when
	1326	// forming the right-hand side elements.
	1327	//
	1328	if( (double)(work[j])>(double)(vcrit) )
	1329	{
	1330	rec = 1/vmax;
	1331	for(i_=ki+n; i_<=n+n;i_++)
	1332	{
	1333	work[i_] = rec*work[i_];
	1334	}
	1335	for(i_=ki+n2; i_<=n+n2;i_++)
	1336	{
	1337	work[i_] = rec*work[i_];
	1338	}
	1339	vmax = 1;
	1340	vcrit = bignum;
	1341	}
	1342	i1_ = (ki+2+n)-(ki+2);
	1343	vt = 0.0;
	1344	for(i_=ki+2; i_<=j-1;i_++)
	1345	{
	1346	vt += t[i_,j]*work[i_+i1_];
	1347	}
	1348	work[j+n] = work[j+n]-vt;
	1349	i1_ = (ki+2+n2)-(ki+2);
	1350	vt = 0.0;
	1351	for(i_=ki+2; i_<=j-1;i_++)
	1352	{
	1353	vt += t[i_,j]*work[i_+i1_];
	1354	}
	1355	work[j+n2] = work[j+n2]-vt;
	1356
	1357	//
	1358	// Solve (T(J,J)-(WR-iWI))(X11+iX12)= WK+IWK2
	1359	//
	1360	temp11[1,1] = t[j,j];
	1361	temp12b[1,1] = work[j+n];
	1362	temp12b[1,2] = work[j+n+n];
	1363	internalhsevdlaln2(false, 1, 2, smin, 1.0, ref temp11, 1.0, 1.0, ref temp12b, wr, -wi, ref rswap4, ref zswap4, ref ipivot44, ref civ4, ref crv4, ref x, ref scl, ref xnorm, ref ierr);
	1364
	1365	//
	1366	// Scale if necessary
	1367	//
	1368	if( (double)(scl)!=(double)(1) )
	1369	{
	1370	for(i_=ki+n; i_<=n+n;i_++)
	1371	{
	1372	work[i_] = scl*work[i_];
	1373	}
	1374	for(i_=ki+n2; i_<=n+n2;i_++)
	1375	{
	1376	work[i_] = scl*work[i_];
	1377	}
	1378	}
	1379	work[j+n] = x[1,1];
	1380	work[j+n2] = x[1,2];
	1381	vmax = Math.Max(Math.Abs(work[j+n]), Math.Max(Math.Abs(work[j+n2]), vmax));
	1382	vcrit = bignum/vmax;
	1383	}
	1384	else
	1385	{
	1386
	1387	//
	1388	// 2-by-2 diagonal block
	1389	//
	1390	// Scale if necessary to avoid overflow when forming
	1391	// the right-hand side elements.
	1392	//
	1393	beta = Math.Max(work[j], work[j+1]);
	1394	if( (double)(beta)>(double)(vcrit) )
	1395	{
	1396	rec = 1/vmax;
	1397	for(i_=ki+n; i_<=n+n;i_++)
	1398	{
	1399	work[i_] = rec*work[i_];
	1400	}
	1401	for(i_=ki+n2; i_<=n+n2;i_++)
	1402	{
	1403	work[i_] = rec*work[i_];
	1404	}
	1405	vmax = 1;
	1406	vcrit = bignum;
	1407	}
	1408	i1_ = (ki+2+n)-(ki+2);
	1409	vt = 0.0;
	1410	for(i_=ki+2; i_<=j-1;i_++)
	1411	{
	1412	vt += t[i_,j]*work[i_+i1_];
	1413	}
	1414	work[j+n] = work[j+n]-vt;
	1415	i1_ = (ki+2+n2)-(ki+2);
	1416	vt = 0.0;
	1417	for(i_=ki+2; i_<=j-1;i_++)
	1418	{
	1419	vt += t[i_,j]*work[i_+i1_];
	1420	}
	1421	work[j+n2] = work[j+n2]-vt;
	1422	i1_ = (ki+2+n)-(ki+2);
	1423	vt = 0.0;
	1424	for(i_=ki+2; i_<=j-1;i_++)
	1425	{
	1426	vt += t[i_,j+1]*work[i_+i1_];
	1427	}
	1428	work[j+1+n] = work[j+1+n]-vt;
	1429	i1_ = (ki+2+n2)-(ki+2);
	1430	vt = 0.0;
	1431	for(i_=ki+2; i_<=j-1;i_++)
	1432	{
	1433	vt += t[i_,j+1]*work[i_+i1_];
	1434	}
	1435	work[j+1+n2] = work[j+1+n2]-vt;
	1436
	1437	//
	1438	// Solve 2-by-2 complex linear equation
	1439	// ([T(j,j) T(j,j+1) ]'-(wr-iwi)I)X = SCALEB
	1440	// ([T(j+1,j) T(j+1,j+1)] )
	1441	//
	1442	temp22[1,1] = t[j,j];
	1443	temp22[1,2] = t[j,j+1];
	1444	temp22[2,1] = t[j+1,j];
	1445	temp22[2,2] = t[j+1,j+1];
	1446	temp22b[1,1] = work[j+n];
	1447	temp22b[1,2] = work[j+n+n];
	1448	temp22b[2,1] = work[j+1+n];
	1449	temp22b[2,2] = work[j+1+n+n];
	1450	internalhsevdlaln2(true, 2, 2, smin, 1.0, ref temp22, 1.0, 1.0, ref temp22b, wr, -wi, ref rswap4, ref zswap4, ref ipivot44, ref civ4, ref crv4, ref x, ref scl, ref xnorm, ref ierr);
	1451
	1452	//
	1453	// Scale if necessary
	1454	//
	1455	if( (double)(scl)!=(double)(1) )
	1456	{
	1457	for(i_=ki+n; i_<=n+n;i_++)
	1458	{
	1459	work[i_] = scl*work[i_];
	1460	}
	1461	for(i_=ki+n2; i_<=n+n2;i_++)
	1462	{
	1463	work[i_] = scl*work[i_];
	1464	}
	1465	}
	1466	work[j+n] = x[1,1];
	1467	work[j+n2] = x[1,2];
	1468	work[j+1+n] = x[2,1];
	1469	work[j+1+n2] = x[2,2];
	1470	vmax = Math.Max(Math.Abs(x[1,1]), vmax);
	1471	vmax = Math.Max(Math.Abs(x[1,2]), vmax);
	1472	vmax = Math.Max(Math.Abs(x[2,1]), vmax);
	1473	vmax = Math.Max(Math.Abs(x[2,2]), vmax);
	1474	vcrit = bignum/vmax;
	1475	}
	1476	}
	1477
	1478	//
	1479	// Copy the vector x or Q*x to VL and normalize.
	1480	//
	1481	if( !over )
	1482	{
	1483	i1_ = (ki+n) - (ki);
	1484	for(i_=ki; i_<=n;i_++)
	1485	{
	1486	vl[i_,iis] = work[i_+i1_];
	1487	}
	1488	i1_ = (ki+n2) - (ki);
	1489	for(i_=ki; i_<=n;i_++)
	1490	{
	1491	vl[i_,iis+1] = work[i_+i1_];
	1492	}
	1493	emax = 0;
	1494	for(k=ki; k<=n; k++)
	1495	{
	1496	emax = Math.Max(emax, Math.Abs(vl[k,iis])+Math.Abs(vl[k,iis+1]));
	1497	}
	1498	remax = 1/emax;
	1499	for(i_=ki; i_<=n;i_++)
	1500	{
	1501	vl[i_,iis] = remax*vl[i_,iis];
	1502	}
	1503	for(i_=ki; i_<=n;i_++)
	1504	{
	1505	vl[i_,iis+1] = remax*vl[i_,iis+1];
	1506	}
	1507	for(k=1; k<=ki-1; k++)
	1508	{
	1509	vl[k,iis] = 0;
	1510	vl[k,iis+1] = 0;
	1511	}
	1512	}
	1513	else
	1514	{
	1515	if( ki<n-1 )
	1516	{
	1517	for(i_=1; i_<=n;i_++)
	1518	{
	1519	temp[i_] = vl[i_,ki];
	1520	}
	1521	blas.matrixvectormultiply(ref vl, 1, n, ki+2, n, false, ref work, ki+2+n, n+n, 1.0, ref temp, 1, n, work[ki+n]);
	1522	for(i_=1; i_<=n;i_++)
	1523	{
	1524	vl[i_,ki] = temp[i_];
	1525	}
	1526	for(i_=1; i_<=n;i_++)
	1527	{
	1528	temp[i_] = vl[i_,ki+1];
	1529	}
	1530	blas.matrixvectormultiply(ref vl, 1, n, ki+2, n, false, ref work, ki+2+n2, n+n2, 1.0, ref temp, 1, n, work[ki+1+n2]);
	1531	for(i_=1; i_<=n;i_++)
	1532	{
	1533	vl[i_,ki+1] = temp[i_];
	1534	}
	1535	}
	1536	else
	1537	{
	1538	vt = work[ki+n];
	1539	for(i_=1; i_<=n;i_++)
	1540	{
	1541	vl[i_,ki] = vt*vl[i_,ki];
	1542	}
	1543	vt = work[ki+1+n2];
	1544	for(i_=1; i_<=n;i_++)
	1545	{
	1546	vl[i_,ki+1] = vt*vl[i_,ki+1];
	1547	}
	1548	}
	1549	emax = 0;
	1550	for(k=1; k<=n; k++)
	1551	{
	1552	emax = Math.Max(emax, Math.Abs(vl[k,ki])+Math.Abs(vl[k,ki+1]));
	1553	}
	1554	remax = 1/emax;
	1555	for(i_=1; i_<=n;i_++)
	1556	{
	1557	vl[i_,ki] = remax*vl[i_,ki];
	1558	}
	1559	for(i_=1; i_<=n;i_++)
	1560	{
	1561	vl[i_,ki+1] = remax*vl[i_,ki+1];
	1562	}
	1563	}
	1564	}
	1565	iis = iis+1;
	1566	if( ip!=0 )
	1567	{
	1568	iis = iis+1;
	1569	}
	1570	}
	1571	if( ip==-1 )
	1572	{
	1573	ip = 0;
	1574	}
	1575	if( ip==1 )
	1576	{
	1577	ip = -1;
	1578	}
	1579	}
	1580	}
	1581	}
	1582
	1583
	1584	private static void internalhsevdlaln2(bool ltrans,
	1585	int na,
	1586	int nw,
	1587	double smin,
	1588	double ca,
	1589	ref double[,] a,
	1590	double d1,
	1591	double d2,
	1592	ref double[,] b,
	1593	double wr,
	1594	double wi,
	1595	ref bool[] rswap4,
	1596	ref bool[] zswap4,
	1597	ref int[,] ipivot44,
	1598	ref double[] civ4,
	1599	ref double[] crv4,
	1600	ref double[,] x,
	1601	ref double scl,
	1602	ref double xnorm,
	1603	ref int info)
	1604	{
	1605	int icmax = 0;
	1606	int j = 0;
	1607	double bbnd = 0;
	1608	double bi1 = 0;
	1609	double bi2 = 0;
	1610	double bignum = 0;
	1611	double bnorm = 0;
	1612	double br1 = 0;
	1613	double br2 = 0;
	1614	double ci21 = 0;
	1615	double ci22 = 0;
	1616	double cmax = 0;
	1617	double cnorm = 0;
	1618	double cr21 = 0;
	1619	double cr22 = 0;
	1620	double csi = 0;
	1621	double csr = 0;
	1622	double li21 = 0;
	1623	double lr21 = 0;
	1624	double smini = 0;
	1625	double smlnum = 0;
	1626	double temp = 0;
	1627	double u22abs = 0;
	1628	double ui11 = 0;
	1629	double ui11r = 0;
	1630	double ui12 = 0;
	1631	double ui12s = 0;
	1632	double ui22 = 0;
	1633	double ur11 = 0;
	1634	double ur11r = 0;
	1635	double ur12 = 0;
	1636	double ur12s = 0;
	1637	double ur22 = 0;
	1638	double xi1 = 0;
	1639	double xi2 = 0;
	1640	double xr1 = 0;
	1641	double xr2 = 0;
	1642	double tmp1 = 0;
	1643	double tmp2 = 0;
	1644
	1645	zswap4[1] = false;
	1646	zswap4[2] = false;
	1647	zswap4[3] = true;
	1648	zswap4[4] = true;
	1649	rswap4[1] = false;
	1650	rswap4[2] = true;
	1651	rswap4[3] = false;
	1652	rswap4[4] = true;
	1653	ipivot44[1,1] = 1;
	1654	ipivot44[2,1] = 2;
	1655	ipivot44[3,1] = 3;
	1656	ipivot44[4,1] = 4;
	1657	ipivot44[1,2] = 2;
	1658	ipivot44[2,2] = 1;
	1659	ipivot44[3,2] = 4;
	1660	ipivot44[4,2] = 3;
	1661	ipivot44[1,3] = 3;
	1662	ipivot44[2,3] = 4;
	1663	ipivot44[3,3] = 1;
	1664	ipivot44[4,3] = 2;
	1665	ipivot44[1,4] = 4;
	1666	ipivot44[2,4] = 3;
	1667	ipivot44[3,4] = 2;
	1668	ipivot44[4,4] = 1;
	1669	smlnum = 2*AP.Math.MinRealNumber;
	1670	bignum = 1/smlnum;
	1671	smini = Math.Max(smin, smlnum);
	1672
	1673	//
	1674	// Don't check for input errors
	1675	//
	1676	info = 0;
	1677
	1678	//
	1679	// Standard Initializations
	1680	//
	1681	scl = 1;
	1682	if( na==1 )
	1683	{
	1684
	1685	//
	1686	// 1 x 1 (i.e., scalar) system C X = B
	1687	//
	1688	if( nw==1 )
	1689	{
	1690
	1691	//
	1692	// Real 1x1 system.
	1693	//
	1694	// C = ca A - w D
	1695	//
	1696	csr = caa[1,1]-wrd1;
	1697	cnorm = Math.Abs(csr);
	1698
	1699	//
	1700	// If \| C \| < SMINI, use C = SMINI
	1701	//
	1702	if( (double)(cnorm)<(double)(smini) )
	1703	{
	1704	csr = smini;
	1705	cnorm = smini;
	1706	info = 1;
	1707	}
	1708
	1709	//
	1710	// Check scaling for X = B / C
	1711	//
	1712	bnorm = Math.Abs(b[1,1]);
	1713	if( (double)(cnorm)<(double)(1) & (double)(bnorm)>(double)(1) )
	1714	{
	1715	if( (double)(bnorm)>(double)(bignum*cnorm) )
	1716	{
	1717	scl = 1/bnorm;
	1718	}
	1719	}
	1720
	1721	//
	1722	// Compute X
	1723	//
	1724	x[1,1] = b[1,1]*scl/csr;
	1725	xnorm = Math.Abs(x[1,1]);
	1726	}
	1727	else
	1728	{
	1729
	1730	//
	1731	// Complex 1x1 system (w is complex)
	1732	//
	1733	// C = ca A - w D
	1734	//
	1735	csr = caa[1,1]-wrd1;
	1736	csi = -(wi*d1);
	1737	cnorm = Math.Abs(csr)+Math.Abs(csi);
	1738
	1739	//
	1740	// If \| C \| < SMINI, use C = SMINI
	1741	//
	1742	if( (double)(cnorm)<(double)(smini) )
	1743	{
	1744	csr = smini;
	1745	csi = 0;
	1746	cnorm = smini;
	1747	info = 1;
	1748	}
	1749
	1750	//
	1751	// Check scaling for X = B / C
	1752	//
	1753	bnorm = Math.Abs(b[1,1])+Math.Abs(b[1,2]);
	1754	if( (double)(cnorm)<(double)(1) & (double)(bnorm)>(double)(1) )
	1755	{
	1756	if( (double)(bnorm)>(double)(bignum*cnorm) )
	1757	{
	1758	scl = 1/bnorm;
	1759	}
	1760	}
	1761
	1762	//
	1763	// Compute X
	1764	//
	1765	internalhsevdladiv(sclb[1,1], sclb[1,2], csr, csi, ref tmp1, ref tmp2);
	1766	x[1,1] = tmp1;
	1767	x[1,2] = tmp2;
	1768	xnorm = Math.Abs(x[1,1])+Math.Abs(x[1,2]);
	1769	}
	1770	}
	1771	else
	1772	{
	1773
	1774	//
	1775	// 2x2 System
	1776	//
	1777	// Compute the real part of C = ca A - w D (or ca A' - w D )
	1778	//
	1779	crv4[1+0] = caa[1,1]-wrd1;
	1780	crv4[2+2] = caa[2,2]-wrd2;
	1781	if( ltrans )
	1782	{
	1783	crv4[1+2] = ca*a[2,1];
	1784	crv4[2+0] = ca*a[1,2];
	1785	}
	1786	else
	1787	{
	1788	crv4[2+0] = ca*a[2,1];
	1789	crv4[1+2] = ca*a[1,2];
	1790	}
	1791	if( nw==1 )
	1792	{
	1793
	1794	//
	1795	// Real 2x2 system (w is real)
	1796	//
	1797	// Find the largest element in C
	1798	//
	1799	cmax = 0;
	1800	icmax = 0;
	1801	for(j=1; j<=4; j++)
	1802	{
	1803	if( (double)(Math.Abs(crv4[j]))>(double)(cmax) )
	1804	{
	1805	cmax = Math.Abs(crv4[j]);
	1806	icmax = j;
	1807	}
	1808	}
	1809
	1810	//
	1811	// If norm(C) < SMINI, use SMINI*identity.
	1812	//
	1813	if( (double)(cmax)<(double)(smini) )
	1814	{
	1815	bnorm = Math.Max(Math.Abs(b[1,1]), Math.Abs(b[2,1]));
	1816	if( (double)(smini)<(double)(1) & (double)(bnorm)>(double)(1) )
	1817	{
	1818	if( (double)(bnorm)>(double)(bignum*smini) )
	1819	{
	1820	scl = 1/bnorm;
	1821	}
	1822	}
	1823	temp = scl/smini;
	1824	x[1,1] = temp*b[1,1];
	1825	x[2,1] = temp*b[2,1];
	1826	xnorm = temp*bnorm;
	1827	info = 1;
	1828	return;
	1829	}
	1830
	1831	//
	1832	// Gaussian elimination with complete pivoting.
	1833	//
	1834	ur11 = crv4[icmax];
	1835	cr21 = crv4[ipivot44[2,icmax]];
	1836	ur12 = crv4[ipivot44[3,icmax]];
	1837	cr22 = crv4[ipivot44[4,icmax]];
	1838	ur11r = 1/ur11;
	1839	lr21 = ur11r*cr21;
	1840	ur22 = cr22-ur12*lr21;
	1841
	1842	//
	1843	// If smaller pivot < SMINI, use SMINI
	1844	//
	1845	if( (double)(Math.Abs(ur22))<(double)(smini) )
	1846	{
	1847	ur22 = smini;
	1848	info = 1;
	1849	}
	1850	if( rswap4[icmax] )
	1851	{
	1852	br1 = b[2,1];
	1853	br2 = b[1,1];
	1854	}
	1855	else
	1856	{
	1857	br1 = b[1,1];
	1858	br2 = b[2,1];
	1859	}
	1860	br2 = br2-lr21*br1;
	1861	bbnd = Math.Max(Math.Abs(br1(ur22ur11r)), Math.Abs(br2));
	1862	if( (double)(bbnd)>(double)(1) & (double)(Math.Abs(ur22))<(double)(1) )
	1863	{
	1864	if( (double)(bbnd)>=(double)(bignum*Math.Abs(ur22)) )
	1865	{
	1866	scl = 1/bbnd;
	1867	}
	1868	}
	1869	xr2 = br2*scl/ur22;
	1870	xr1 = sclbr1ur11r-xr2(ur11rur12);
	1871	if( zswap4[icmax] )
	1872	{
	1873	x[1,1] = xr2;
	1874	x[2,1] = xr1;
	1875	}
	1876	else
	1877	{
	1878	x[1,1] = xr1;
	1879	x[2,1] = xr2;
	1880	}
	1881	xnorm = Math.Max(Math.Abs(xr1), Math.Abs(xr2));
	1882
	1883	//
	1884	// Further scaling if norm(A) norm(X) > overflow
	1885	//
	1886	if( (double)(xnorm)>(double)(1) & (double)(cmax)>(double)(1) )
	1887	{
	1888	if( (double)(xnorm)>(double)(bignum/cmax) )
	1889	{
	1890	temp = cmax/bignum;
	1891	x[1,1] = temp*x[1,1];
	1892	x[2,1] = temp*x[2,1];
	1893	xnorm = temp*xnorm;
	1894	scl = temp*scl;
	1895	}
	1896	}
	1897	}
	1898	else
	1899	{
	1900
	1901	//
	1902	// Complex 2x2 system (w is complex)
	1903	//
	1904	// Find the largest element in C
	1905	//
	1906	civ4[1+0] = -(wi*d1);
	1907	civ4[2+0] = 0;
	1908	civ4[1+2] = 0;
	1909	civ4[2+2] = -(wi*d2);
	1910	cmax = 0;
	1911	icmax = 0;
	1912	for(j=1; j<=4; j++)
	1913	{
	1914	if( (double)(Math.Abs(crv4[j])+Math.Abs(civ4[j]))>(double)(cmax) )
	1915	{
	1916	cmax = Math.Abs(crv4[j])+Math.Abs(civ4[j]);
	1917	icmax = j;
	1918	}
	1919	}
	1920
	1921	//
	1922	// If norm(C) < SMINI, use SMINI*identity.
	1923	//
	1924	if( (double)(cmax)<(double)(smini) )
	1925	{
	1926	bnorm = Math.Max(Math.Abs(b[1,1])+Math.Abs(b[1,2]), Math.Abs(b[2,1])+Math.Abs(b[2,2]));
	1927	if( (double)(smini)<(double)(1) & (double)(bnorm)>(double)(1) )
	1928	{
	1929	if( (double)(bnorm)>(double)(bignum*smini) )
	1930	{
	1931	scl = 1/bnorm;
	1932	}
	1933	}
	1934	temp = scl/smini;
	1935	x[1,1] = temp*b[1,1];
	1936	x[2,1] = temp*b[2,1];
	1937	x[1,2] = temp*b[1,2];
	1938	x[2,2] = temp*b[2,2];
	1939	xnorm = temp*bnorm;
	1940	info = 1;
	1941	return;
	1942	}
	1943
	1944	//
	1945	// Gaussian elimination with complete pivoting.
	1946	//
	1947	ur11 = crv4[icmax];
	1948	ui11 = civ4[icmax];
	1949	cr21 = crv4[ipivot44[2,icmax]];
	1950	ci21 = civ4[ipivot44[2,icmax]];
	1951	ur12 = crv4[ipivot44[3,icmax]];
	1952	ui12 = civ4[ipivot44[3,icmax]];
	1953	cr22 = crv4[ipivot44[4,icmax]];
	1954	ci22 = civ4[ipivot44[4,icmax]];
	1955	if( icmax==1 \| icmax==4 )
	1956	{
	1957
	1958	//
	1959	// Code when off-diagonals of pivoted C are real
	1960	//
	1961	if( (double)(Math.Abs(ur11))>(double)(Math.Abs(ui11)) )
	1962	{
	1963	temp = ui11/ur11;
	1964	ur11r = 1/(ur11*(1+AP.Math.Sqr(temp)));
	1965	ui11r = -(temp*ur11r);
	1966	}
	1967	else
	1968	{
	1969	temp = ur11/ui11;
	1970	ui11r = -(1/(ui11*(1+AP.Math.Sqr(temp))));
	1971	ur11r = -(temp*ui11r);
	1972	}
	1973	lr21 = cr21*ur11r;
	1974	li21 = cr21*ui11r;
	1975	ur12s = ur12*ur11r;
	1976	ui12s = ur12*ui11r;
	1977	ur22 = cr22-ur12*lr21;
	1978	ui22 = ci22-ur12*li21;
	1979	}
	1980	else
	1981	{
	1982
	1983	//
	1984	// Code when diagonals of pivoted C are real
	1985	//
	1986	ur11r = 1/ur11;
	1987	ui11r = 0;
	1988	lr21 = cr21*ur11r;
	1989	li21 = ci21*ur11r;
	1990	ur12s = ur12*ur11r;
	1991	ui12s = ui12*ur11r;
	1992	ur22 = cr22-ur12lr21+ui12li21;
	1993	ui22 = -(ur12li21)-ui12lr21;
	1994	}
	1995	u22abs = Math.Abs(ur22)+Math.Abs(ui22);
	1996
	1997	//
	1998	// If smaller pivot < SMINI, use SMINI
	1999	//
	2000	if( (double)(u22abs)<(double)(smini) )
	2001	{
	2002	ur22 = smini;
	2003	ui22 = 0;
	2004	info = 1;
	2005	}
	2006	if( rswap4[icmax] )
	2007	{
	2008	br2 = b[1,1];
	2009	br1 = b[2,1];
	2010	bi2 = b[1,2];
	2011	bi1 = b[2,2];
	2012	}
	2013	else
	2014	{
	2015	br1 = b[1,1];
	2016	br2 = b[2,1];
	2017	bi1 = b[1,2];
	2018	bi2 = b[2,2];
	2019	}
	2020	br2 = br2-lr21br1+li21bi1;
	2021	bi2 = bi2-li21br1-lr21bi1;
	2022	bbnd = Math.Max((Math.Abs(br1)+Math.Abs(bi1))(u22abs(Math.Abs(ur11r)+Math.Abs(ui11r))), Math.Abs(br2)+Math.Abs(bi2));
	2023	if( (double)(bbnd)>(double)(1) & (double)(u22abs)<(double)(1) )
	2024	{
	2025	if( (double)(bbnd)>=(double)(bignum*u22abs) )
	2026	{
	2027	scl = 1/bbnd;
	2028	br1 = scl*br1;
	2029	bi1 = scl*bi1;
	2030	br2 = scl*br2;
	2031	bi2 = scl*bi2;
	2032	}
	2033	}
	2034	internalhsevdladiv(br2, bi2, ur22, ui22, ref xr2, ref xi2);
	2035	xr1 = ur11rbr1-ui11rbi1-ur12sxr2+ui12sxi2;
	2036	xi1 = ui11rbr1+ur11rbi1-ui12sxr2-ur12sxi2;
	2037	if( zswap4[icmax] )
	2038	{
	2039	x[1,1] = xr2;
	2040	x[2,1] = xr1;
	2041	x[1,2] = xi2;
	2042	x[2,2] = xi1;
	2043	}
	2044	else
	2045	{
	2046	x[1,1] = xr1;
	2047	x[2,1] = xr2;
	2048	x[1,2] = xi1;
	2049	x[2,2] = xi2;
	2050	}
	2051	xnorm = Math.Max(Math.Abs(xr1)+Math.Abs(xi1), Math.Abs(xr2)+Math.Abs(xi2));
	2052
	2053	//
	2054	// Further scaling if norm(A) norm(X) > overflow
	2055	//
	2056	if( (double)(xnorm)>(double)(1) & (double)(cmax)>(double)(1) )
	2057	{
	2058	if( (double)(xnorm)>(double)(bignum/cmax) )
	2059	{
	2060	temp = cmax/bignum;
	2061	x[1,1] = temp*x[1,1];
	2062	x[2,1] = temp*x[2,1];
	2063	x[1,2] = temp*x[1,2];
	2064	x[2,2] = temp*x[2,2];
	2065	xnorm = temp*xnorm;
	2066	scl = temp*scl;
	2067	}
	2068	}
	2069	}
	2070	}
	2071	}
	2072
	2073
	2074	private static void internalhsevdladiv(double a,
	2075	double b,
	2076	double c,
	2077	double d,
	2078	ref double p,
	2079	ref double q)
	2080	{
	2081	double e = 0;
	2082	double f = 0;
	2083
	2084	if( (double)(Math.Abs(d))<(double)(Math.Abs(c)) )
	2085	{
	2086	e = d/c;
	2087	f = c+d*e;
	2088	p = (a+b*e)/f;
	2089	q = (b-a*e)/f;
	2090	}
	2091	else
	2092	{
	2093	e = c/d;
	2094	f = d+c*e;
	2095	p = (b+a*e)/f;
	2096	q = (-a+b*e)/f;
	2097	}
	2098	}
	2099	}
	2100	}

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

Download in other formats:

Update cookies preferences