Context Navigation

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

Visit:

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

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

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

Download in other formats:

Update cookies preferences