Context Navigation

source: trunk/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.3.0/ALGLIB-2.3.0/matgen.cs @ 3021

Visit:

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

Rev	Line
[2806]	1	/*************************************************************************
	2	Copyright (c) 2007, Sergey Bochkanov (ALGLIB project).
	3
	4	>>> SOURCE LICENSE >>>
	5	This program is free software; you can redistribute it and/or modify
	6	it under the terms of the GNU General Public License as published by
	7	the Free Software Foundation (www.fsf.org); either version 2 of the
	8	License, or (at your option) any later version.
	9
	10	This program is distributed in the hope that it will be useful,
	11	but WITHOUT ANY WARRANTY; without even the implied warranty of
	12	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
	13	GNU General Public License for more details.
	14
	15	A copy of the GNU General Public License is available at
	16	http://www.fsf.org/licensing/licenses
	17
	18	>>> END OF LICENSE >>>
	19	*************************************************************************/
	20
	21	using System;
	22
	23	namespace alglib
	24	{
	25	public class matgen
	26	{
	27	/*************************************************************************
	28	Generation of a random uniformly distributed (Haar) orthogonal matrix
	29
	30	INPUT PARAMETERS:
	31	N - matrix size, N>=1
	32
	33	OUTPUT PARAMETERS:
	34	A - orthogonal NxN matrix, array[0..N-1,0..N-1]
	35
	36	-- ALGLIB routine --
	37	04.12.2009
	38	Bochkanov Sergey
	39	*************************************************************************/
	40	public static void rmatrixrndorthogonal(int n,
	41	ref double[,] a)
	42	{
	43	int i = 0;
	44	int j = 0;
	45
	46	System.Diagnostics.Debug.Assert(n>=1, "RMatrixRndOrthogonal: N<1!");
	47	a = new double[n-1+1, n-1+1];
	48	for(i=0; i<=n-1; i++)
	49	{
	50	for(j=0; j<=n-1; j++)
	51	{
	52	if( i==j )
	53	{
	54	a[i,j] = 1;
	55	}
	56	else
	57	{
	58	a[i,j] = 0;
	59	}
	60	}
	61	}
	62	rmatrixrndorthogonalfromtheright(ref a, n, n);
	63	}
	64
	65
	66	/*************************************************************************
	67	Generation of random NxN matrix with given condition number and norm2(A)=1
	68
	69	INPUT PARAMETERS:
	70	N - matrix size
	71	C - condition number (in 2-norm)
	72
	73	OUTPUT PARAMETERS:
	74	A - random matrix with norm2(A)=1 and cond(A)=C
	75
	76	-- ALGLIB routine --
	77	04.12.2009
	78	Bochkanov Sergey
	79	*************************************************************************/
	80	public static void rmatrixrndcond(int n,
	81	double c,
	82	ref double[,] a)
	83	{
	84	int i = 0;
	85	int j = 0;
	86	double l1 = 0;
	87	double l2 = 0;
	88
	89	System.Diagnostics.Debug.Assert(n>=1 & (double)(c)>=(double)(1), "RMatrixRndCond: N<1 or C<1!");
	90	a = new double[n-1+1, n-1+1];
	91	if( n==1 )
	92	{
	93
	94	//
	95	// special case
	96	//
	97	a[0,0] = 2*AP.Math.RandomInteger(2)-1;
	98	return;
	99	}
	100	l1 = 0;
	101	l2 = Math.Log(1/c);
	102	for(i=0; i<=n-1; i++)
	103	{
	104	for(j=0; j<=n-1; j++)
	105	{
	106	a[i,j] = 0;
	107	}
	108	}
	109	a[0,0] = Math.Exp(l1);
	110	for(i=1; i<=n-2; i++)
	111	{
	112	a[i,i] = Math.Exp(AP.Math.RandomReal()*(l2-l1)+l1);
	113	}
	114	a[n-1,n-1] = Math.Exp(l2);
	115	rmatrixrndorthogonalfromtheleft(ref a, n, n);
	116	rmatrixrndorthogonalfromtheright(ref a, n, n);
	117	}
	118
	119
	120	/*************************************************************************
	121	Generation of a random Haar distributed orthogonal complex matrix
	122
	123	INPUT PARAMETERS:
	124	N - matrix size, N>=1
	125
	126	OUTPUT PARAMETERS:
	127	A - orthogonal NxN matrix, array[0..N-1,0..N-1]
	128
	129	-- ALGLIB routine --
	130	04.12.2009
	131	Bochkanov Sergey
	132	*************************************************************************/
	133	public static void cmatrixrndorthogonal(int n,
	134	ref AP.Complex[,] a)
	135	{
	136	int i = 0;
	137	int j = 0;
	138
	139	System.Diagnostics.Debug.Assert(n>=1, "CMatrixRndOrthogonal: N<1!");
	140	a = new AP.Complex[n-1+1, n-1+1];
	141	for(i=0; i<=n-1; i++)
	142	{
	143	for(j=0; j<=n-1; j++)
	144	{
	145	if( i==j )
	146	{
	147	a[i,j] = 1;
	148	}
	149	else
	150	{
	151	a[i,j] = 0;
	152	}
	153	}
	154	}
	155	cmatrixrndorthogonalfromtheright(ref a, n, n);
	156	}
	157
	158
	159	/*************************************************************************
	160	Generation of random NxN complex matrix with given condition number C and
	161	norm2(A)=1
	162
	163	INPUT PARAMETERS:
	164	N - matrix size
	165	C - condition number (in 2-norm)
	166
	167	OUTPUT PARAMETERS:
	168	A - random matrix with norm2(A)=1 and cond(A)=C
	169
	170	-- ALGLIB routine --
	171	04.12.2009
	172	Bochkanov Sergey
	173	*************************************************************************/
	174	public static void cmatrixrndcond(int n,
	175	double c,
	176	ref AP.Complex[,] a)
	177	{
	178	int i = 0;
	179	int j = 0;
	180	double l1 = 0;
	181	double l2 = 0;
	182	hqrnd.hqrndstate state = new hqrnd.hqrndstate();
	183	AP.Complex v = 0;
	184
	185	System.Diagnostics.Debug.Assert(n>=1 & (double)(c)>=(double)(1), "CMatrixRndCond: N<1 or C<1!");
	186	a = new AP.Complex[n-1+1, n-1+1];
	187	if( n==1 )
	188	{
	189
	190	//
	191	// special case
	192	//
	193	hqrnd.hqrndrandomize(ref state);
	194	hqrnd.hqrndunit2(ref state, ref v.x, ref v.y);
	195	a[0,0] = v;
	196	return;
	197	}
	198	l1 = 0;
	199	l2 = Math.Log(1/c);
	200	for(i=0; i<=n-1; i++)
	201	{
	202	for(j=0; j<=n-1; j++)
	203	{
	204	a[i,j] = 0;
	205	}
	206	}
	207	a[0,0] = Math.Exp(l1);
	208	for(i=1; i<=n-2; i++)
	209	{
	210	a[i,i] = Math.Exp(AP.Math.RandomReal()*(l2-l1)+l1);
	211	}
	212	a[n-1,n-1] = Math.Exp(l2);
	213	cmatrixrndorthogonalfromtheleft(ref a, n, n);
	214	cmatrixrndorthogonalfromtheright(ref a, n, n);
	215	}
	216
	217
	218	/*************************************************************************
	219	Generation of random NxN symmetric matrix with given condition number and
	220	norm2(A)=1
	221
	222	INPUT PARAMETERS:
	223	N - matrix size
	224	C - condition number (in 2-norm)
	225
	226	OUTPUT PARAMETERS:
	227	A - random matrix with norm2(A)=1 and cond(A)=C
	228
	229	-- ALGLIB routine --
	230	04.12.2009
	231	Bochkanov Sergey
	232	*************************************************************************/
	233	public static void smatrixrndcond(int n,
	234	double c,
	235	ref double[,] a)
	236	{
	237	int i = 0;
	238	int j = 0;
	239	double l1 = 0;
	240	double l2 = 0;
	241
	242	System.Diagnostics.Debug.Assert(n>=1 & (double)(c)>=(double)(1), "SMatrixRndCond: N<1 or C<1!");
	243	a = new double[n-1+1, n-1+1];
	244	if( n==1 )
	245	{
	246
	247	//
	248	// special case
	249	//
	250	a[0,0] = 2*AP.Math.RandomInteger(2)-1;
	251	return;
	252	}
	253
	254	//
	255	// Prepare matrix
	256	//
	257	l1 = 0;
	258	l2 = Math.Log(1/c);
	259	for(i=0; i<=n-1; i++)
	260	{
	261	for(j=0; j<=n-1; j++)
	262	{
	263	a[i,j] = 0;
	264	}
	265	}
	266	a[0,0] = Math.Exp(l1);
	267	for(i=1; i<=n-2; i++)
	268	{
	269	a[i,i] = (2AP.Math.RandomInteger(2)-1)Math.Exp(AP.Math.RandomReal()*(l2-l1)+l1);
	270	}
	271	a[n-1,n-1] = Math.Exp(l2);
	272
	273	//
	274	// Multiply
	275	//
	276	smatrixrndmultiply(ref a, n);
	277	}
	278
	279
	280	/*************************************************************************
	281	Generation of random NxN symmetric positive definite matrix with given
	282	condition number and norm2(A)=1
	283
	284	INPUT PARAMETERS:
	285	N - matrix size
	286	C - condition number (in 2-norm)
	287
	288	OUTPUT PARAMETERS:
	289	A - random SPD matrix with norm2(A)=1 and cond(A)=C
	290
	291	-- ALGLIB routine --
	292	04.12.2009
	293	Bochkanov Sergey
	294	*************************************************************************/
	295	public static void spdmatrixrndcond(int n,
	296	double c,
	297	ref double[,] a)
	298	{
	299	int i = 0;
	300	int j = 0;
	301	double l1 = 0;
	302	double l2 = 0;
	303
	304
	305	//
	306	// Special cases
	307	//
	308	if( n<=0 \| (double)(c)<(double)(1) )
	309	{
	310	return;
	311	}
	312	a = new double[n-1+1, n-1+1];
	313	if( n==1 )
	314	{
	315	a[0,0] = 1;
	316	return;
	317	}
	318
	319	//
	320	// Prepare matrix
	321	//
	322	l1 = 0;
	323	l2 = Math.Log(1/c);
	324	for(i=0; i<=n-1; i++)
	325	{
	326	for(j=0; j<=n-1; j++)
	327	{
	328	a[i,j] = 0;
	329	}
	330	}
	331	a[0,0] = Math.Exp(l1);
	332	for(i=1; i<=n-2; i++)
	333	{
	334	a[i,i] = Math.Exp(AP.Math.RandomReal()*(l2-l1)+l1);
	335	}
	336	a[n-1,n-1] = Math.Exp(l2);
	337
	338	//
	339	// Multiply
	340	//
	341	smatrixrndmultiply(ref a, n);
	342	}
	343
	344
	345	/*************************************************************************
	346	Generation of random NxN Hermitian matrix with given condition number and
	347	norm2(A)=1
	348
	349	INPUT PARAMETERS:
	350	N - matrix size
	351	C - condition number (in 2-norm)
	352
	353	OUTPUT PARAMETERS:
	354	A - random matrix with norm2(A)=1 and cond(A)=C
	355
	356	-- ALGLIB routine --
	357	04.12.2009
	358	Bochkanov Sergey
	359	*************************************************************************/
	360	public static void hmatrixrndcond(int n,
	361	double c,
	362	ref AP.Complex[,] a)
	363	{
	364	int i = 0;
	365	int j = 0;
	366	double l1 = 0;
	367	double l2 = 0;
	368
	369	System.Diagnostics.Debug.Assert(n>=1 & (double)(c)>=(double)(1), "HMatrixRndCond: N<1 or C<1!");
	370	a = new AP.Complex[n-1+1, n-1+1];
	371	if( n==1 )
	372	{
	373
	374	//
	375	// special case
	376	//
	377	a[0,0] = 2*AP.Math.RandomInteger(2)-1;
	378	return;
	379	}
	380
	381	//
	382	// Prepare matrix
	383	//
	384	l1 = 0;
	385	l2 = Math.Log(1/c);
	386	for(i=0; i<=n-1; i++)
	387	{
	388	for(j=0; j<=n-1; j++)
	389	{
	390	a[i,j] = 0;
	391	}
	392	}
	393	a[0,0] = Math.Exp(l1);
	394	for(i=1; i<=n-2; i++)
	395	{
	396	a[i,i] = (2AP.Math.RandomInteger(2)-1)Math.Exp(AP.Math.RandomReal()*(l2-l1)+l1);
	397	}
	398	a[n-1,n-1] = Math.Exp(l2);
	399
	400	//
	401	// Multiply
	402	//
	403	hmatrixrndmultiply(ref a, n);
	404
	405	//
	406	// post-process to ensure that matrix diagonal is real
	407	//
	408	for(i=0; i<=n-1; i++)
	409	{
	410	a[i,i].y = 0;
	411	}
	412	}
	413
	414
	415	/*************************************************************************
	416	Generation of random NxN Hermitian positive definite matrix with given
	417	condition number and norm2(A)=1
	418
	419	INPUT PARAMETERS:
	420	N - matrix size
	421	C - condition number (in 2-norm)
	422
	423	OUTPUT PARAMETERS:
	424	A - random HPD matrix with norm2(A)=1 and cond(A)=C
	425
	426	-- ALGLIB routine --
	427	04.12.2009
	428	Bochkanov Sergey
	429	*************************************************************************/
	430	public static void hpdmatrixrndcond(int n,
	431	double c,
	432	ref AP.Complex[,] a)
	433	{
	434	int i = 0;
	435	int j = 0;
	436	double l1 = 0;
	437	double l2 = 0;
	438
	439
	440	//
	441	// Special cases
	442	//
	443	if( n<=0 \| (double)(c)<(double)(1) )
	444	{
	445	return;
	446	}
	447	a = new AP.Complex[n-1+1, n-1+1];
	448	if( n==1 )
	449	{
	450	a[0,0] = 1;
	451	return;
	452	}
	453
	454	//
	455	// Prepare matrix
	456	//
	457	l1 = 0;
	458	l2 = Math.Log(1/c);
	459	for(i=0; i<=n-1; i++)
	460	{
	461	for(j=0; j<=n-1; j++)
	462	{
	463	a[i,j] = 0;
	464	}
	465	}
	466	a[0,0] = Math.Exp(l1);
	467	for(i=1; i<=n-2; i++)
	468	{
	469	a[i,i] = Math.Exp(AP.Math.RandomReal()*(l2-l1)+l1);
	470	}
	471	a[n-1,n-1] = Math.Exp(l2);
	472
	473	//
	474	// Multiply
	475	//
	476	hmatrixrndmultiply(ref a, n);
	477
	478	//
	479	// post-process to ensure that matrix diagonal is real
	480	//
	481	for(i=0; i<=n-1; i++)
	482	{
	483	a[i,i].y = 0;
	484	}
	485	}
	486
	487
	488	/*************************************************************************
	489	Multiplication of MxN matrix by NxN random Haar distributed orthogonal matrix
	490
	491	INPUT PARAMETERS:
	492	A - matrix, array[0..M-1, 0..N-1]
	493	M, N- matrix size
	494
	495	OUTPUT PARAMETERS:
	496	A - A*Q, where Q is random NxN orthogonal matrix
	497
	498	-- ALGLIB routine --
	499	04.12.2009
	500	Bochkanov Sergey
	501	*************************************************************************/
	502	public static void rmatrixrndorthogonalfromtheright(ref double[,] a,
	503	int m,
	504	int n)
	505	{
	506	double tau = 0;
	507	double lambda = 0;
	508	int s = 0;
	509	int i = 0;
	510	double u1 = 0;
	511	double u2 = 0;
	512	double[] w = new double[0];
	513	double[] v = new double[0];
	514	hqrnd.hqrndstate state = new hqrnd.hqrndstate();
	515	int i_ = 0;
	516
	517	System.Diagnostics.Debug.Assert(n>=1 & m>=1, "RMatrixRndOrthogonalFromTheRight: N<1 or M<1!");
	518	if( n==1 )
	519	{
	520
	521	//
	522	// Special case
	523	//
	524	tau = 2*AP.Math.RandomInteger(2)-1;
	525	for(i=0; i<=m-1; i++)
	526	{
	527	a[i,0] = a[i,0]*tau;
	528	}
	529	return;
	530	}
	531
	532	//
	533	// General case.
	534	// First pass.
	535	//
	536	w = new double[m-1+1];
	537	v = new double[n+1];
	538	hqrnd.hqrndrandomize(ref state);
	539	for(s=2; s<=n; s++)
	540	{
	541
	542	//
	543	// Prepare random normal v
	544	//
	545	do
	546	{
	547	i = 1;
	548	while( i<=s )
	549	{
	550	hqrnd.hqrndnormal2(ref state, ref u1, ref u2);
	551	v[i] = u1;
	552	if( i+1<=s )
	553	{
	554	v[i+1] = u2;
	555	}
	556	i = i+2;
	557	}
	558	lambda = 0.0;
	559	for(i_=1; i_<=s;i_++)
	560	{
	561	lambda += v[i_]*v[i_];
	562	}
	563	}
	564	while( (double)(lambda)==(double)(0) );
	565
	566	//
	567	// Prepare and apply reflection
	568	//
	569	reflections.generatereflection(ref v, s, ref tau);
	570	v[1] = 1;
	571	reflections.applyreflectionfromtheright(ref a, tau, ref v, 0, m-1, n-s, n-1, ref w);
	572	}
	573
	574	//
	575	// Second pass.
	576	//
	577	for(i=0; i<=n-1; i++)
	578	{
	579	tau = 2*AP.Math.RandomInteger(2)-1;
	580	for(i_=0; i_<=m-1;i_++)
	581	{
	582	a[i_,i] = tau*a[i_,i];
	583	}
	584	}
	585	}
	586
	587
	588	/*************************************************************************
	589	Multiplication of MxN matrix by MxM random Haar distributed orthogonal matrix
	590
	591	INPUT PARAMETERS:
	592	A - matrix, array[0..M-1, 0..N-1]
	593	M, N- matrix size
	594
	595	OUTPUT PARAMETERS:
	596	A - Q*A, where Q is random MxM orthogonal matrix
	597
	598	-- ALGLIB routine --
	599	04.12.2009
	600	Bochkanov Sergey
	601	*************************************************************************/
	602	public static void rmatrixrndorthogonalfromtheleft(ref double[,] a,
	603	int m,
	604	int n)
	605	{
	606	double tau = 0;
	607	double lambda = 0;
	608	int s = 0;
	609	int i = 0;
	610	int j = 0;
	611	double u1 = 0;
	612	double u2 = 0;
	613	double[] w = new double[0];
	614	double[] v = new double[0];
	615	hqrnd.hqrndstate state = new hqrnd.hqrndstate();
	616	int i_ = 0;
	617
	618	System.Diagnostics.Debug.Assert(n>=1 & m>=1, "RMatrixRndOrthogonalFromTheRight: N<1 or M<1!");
	619	if( m==1 )
	620	{
	621
	622	//
	623	// special case
	624	//
	625	tau = 2*AP.Math.RandomInteger(2)-1;
	626	for(j=0; j<=n-1; j++)
	627	{
	628	a[0,j] = a[0,j]*tau;
	629	}
	630	return;
	631	}
	632
	633	//
	634	// General case.
	635	// First pass.
	636	//
	637	w = new double[n-1+1];
	638	v = new double[m+1];
	639	hqrnd.hqrndrandomize(ref state);
	640	for(s=2; s<=m; s++)
	641	{
	642
	643	//
	644	// Prepare random normal v
	645	//
	646	do
	647	{
	648	i = 1;
	649	while( i<=s )
	650	{
	651	hqrnd.hqrndnormal2(ref state, ref u1, ref u2);
	652	v[i] = u1;
	653	if( i+1<=s )
	654	{
	655	v[i+1] = u2;
	656	}
	657	i = i+2;
	658	}
	659	lambda = 0.0;
	660	for(i_=1; i_<=s;i_++)
	661	{
	662	lambda += v[i_]*v[i_];
	663	}
	664	}
	665	while( (double)(lambda)==(double)(0) );
	666
	667	//
	668	// Prepare and apply reflection
	669	//
	670	reflections.generatereflection(ref v, s, ref tau);
	671	v[1] = 1;
	672	reflections.applyreflectionfromtheleft(ref a, tau, ref v, m-s, m-1, 0, n-1, ref w);
	673	}
	674
	675	//
	676	// Second pass.
	677	//
	678	for(i=0; i<=m-1; i++)
	679	{
	680	tau = 2*AP.Math.RandomInteger(2)-1;
	681	for(i_=0; i_<=n-1;i_++)
	682	{
	683	a[i,i_] = tau*a[i,i_];
	684	}
	685	}
	686	}
	687
	688
	689	/*************************************************************************
	690	Multiplication of MxN complex matrix by NxN random Haar distributed
	691	complex orthogonal matrix
	692
	693	INPUT PARAMETERS:
	694	A - matrix, array[0..M-1, 0..N-1]
	695	M, N- matrix size
	696
	697	OUTPUT PARAMETERS:
	698	A - A*Q, where Q is random NxN orthogonal matrix
	699
	700	-- ALGLIB routine --
	701	04.12.2009
	702	Bochkanov Sergey
	703	*************************************************************************/
	704	public static void cmatrixrndorthogonalfromtheright(ref AP.Complex[,] a,
	705	int m,
	706	int n)
	707	{
	708	AP.Complex lambda = 0;
	709	AP.Complex tau = 0;
	710	int s = 0;
	711	int i = 0;
	712	AP.Complex[] w = new AP.Complex[0];
	713	AP.Complex[] v = new AP.Complex[0];
	714	hqrnd.hqrndstate state = new hqrnd.hqrndstate();
	715	int i_ = 0;
	716
	717	System.Diagnostics.Debug.Assert(n>=1 & m>=1, "CMatrixRndOrthogonalFromTheRight: N<1 or M<1!");
	718	if( n==1 )
	719	{
	720
	721	//
	722	// Special case
	723	//
	724	hqrnd.hqrndrandomize(ref state);
	725	hqrnd.hqrndunit2(ref state, ref tau.x, ref tau.y);
	726	for(i=0; i<=m-1; i++)
	727	{
	728	a[i,0] = a[i,0]*tau;
	729	}
	730	return;
	731	}
	732
	733	//
	734	// General case.
	735	// First pass.
	736	//
	737	w = new AP.Complex[m-1+1];
	738	v = new AP.Complex[n+1];
	739	hqrnd.hqrndrandomize(ref state);
	740	for(s=2; s<=n; s++)
	741	{
	742
	743	//
	744	// Prepare random normal v
	745	//
	746	do
	747	{
	748	for(i=1; i<=s; i++)
	749	{
	750	hqrnd.hqrndnormal2(ref state, ref tau.x, ref tau.y);
	751	v[i] = tau;
	752	}
	753	lambda = 0.0;
	754	for(i_=1; i_<=s;i_++)
	755	{
	756	lambda += v[i_]*AP.Math.Conj(v[i_]);
	757	}
	758	}
	759	while( lambda==0 );
	760
	761	//
	762	// Prepare and apply reflection
	763	//
	764	creflections.complexgeneratereflection(ref v, s, ref tau);
	765	v[1] = 1;
	766	creflections.complexapplyreflectionfromtheright(ref a, tau, ref v, 0, m-1, n-s, n-1, ref w);
	767	}
	768
	769	//
	770	// Second pass.
	771	//
	772	for(i=0; i<=n-1; i++)
	773	{
	774	hqrnd.hqrndunit2(ref state, ref tau.x, ref tau.y);
	775	for(i_=0; i_<=m-1;i_++)
	776	{
	777	a[i_,i] = tau*a[i_,i];
	778	}
	779	}
	780	}
	781
	782
	783	/*************************************************************************
	784	Multiplication of MxN complex matrix by MxM random Haar distributed
	785	complex orthogonal matrix
	786
	787	INPUT PARAMETERS:
	788	A - matrix, array[0..M-1, 0..N-1]
	789	M, N- matrix size
	790
	791	OUTPUT PARAMETERS:
	792	A - Q*A, where Q is random MxM orthogonal matrix
	793
	794	-- ALGLIB routine --
	795	04.12.2009
	796	Bochkanov Sergey
	797	*************************************************************************/
	798	public static void cmatrixrndorthogonalfromtheleft(ref AP.Complex[,] a,
	799	int m,
	800	int n)
	801	{
	802	AP.Complex tau = 0;
	803	AP.Complex lambda = 0;
	804	int s = 0;
	805	int i = 0;
	806	int j = 0;
	807	AP.Complex[] w = new AP.Complex[0];
	808	AP.Complex[] v = new AP.Complex[0];
	809	hqrnd.hqrndstate state = new hqrnd.hqrndstate();
	810	int i_ = 0;
	811
	812	System.Diagnostics.Debug.Assert(n>=1 & m>=1, "CMatrixRndOrthogonalFromTheRight: N<1 or M<1!");
	813	if( m==1 )
	814	{
	815
	816	//
	817	// special case
	818	//
	819	hqrnd.hqrndrandomize(ref state);
	820	hqrnd.hqrndunit2(ref state, ref tau.x, ref tau.y);
	821	for(j=0; j<=n-1; j++)
	822	{
	823	a[0,j] = a[0,j]*tau;
	824	}
	825	return;
	826	}
	827
	828	//
	829	// General case.
	830	// First pass.
	831	//
	832	w = new AP.Complex[n-1+1];
	833	v = new AP.Complex[m+1];
	834	hqrnd.hqrndrandomize(ref state);
	835	for(s=2; s<=m; s++)
	836	{
	837
	838	//
	839	// Prepare random normal v
	840	//
	841	do
	842	{
	843	for(i=1; i<=s; i++)
	844	{
	845	hqrnd.hqrndnormal2(ref state, ref tau.x, ref tau.y);
	846	v[i] = tau;
	847	}
	848	lambda = 0.0;
	849	for(i_=1; i_<=s;i_++)
	850	{
	851	lambda += v[i_]*AP.Math.Conj(v[i_]);
	852	}
	853	}
	854	while( lambda==0 );
	855
	856	//
	857	// Prepare and apply reflection
	858	//
	859	creflections.complexgeneratereflection(ref v, s, ref tau);
	860	v[1] = 1;
	861	creflections.complexapplyreflectionfromtheleft(ref a, tau, ref v, m-s, m-1, 0, n-1, ref w);
	862	}
	863
	864	//
	865	// Second pass.
	866	//
	867	for(i=0; i<=m-1; i++)
	868	{
	869	hqrnd.hqrndunit2(ref state, ref tau.x, ref tau.y);
	870	for(i_=0; i_<=n-1;i_++)
	871	{
	872	a[i,i_] = tau*a[i,i_];
	873	}
	874	}
	875	}
	876
	877
	878	/*************************************************************************
	879	Symmetric multiplication of NxN matrix by random Haar distributed
	880	orthogonal matrix
	881
	882	INPUT PARAMETERS:
	883	A - matrix, array[0..N-1, 0..N-1]
	884	N - matrix size
	885
	886	OUTPUT PARAMETERS:
	887	A - Q'AQ, where Q is random NxN orthogonal matrix
	888
	889	-- ALGLIB routine --
	890	04.12.2009
	891	Bochkanov Sergey
	892	*************************************************************************/
	893	public static void smatrixrndmultiply(ref double[,] a,
	894	int n)
	895	{
	896	double tau = 0;
	897	double lambda = 0;
	898	int s = 0;
	899	int i = 0;
	900	double u1 = 0;
	901	double u2 = 0;
	902	double[] w = new double[0];
	903	double[] v = new double[0];
	904	hqrnd.hqrndstate state = new hqrnd.hqrndstate();
	905	int i_ = 0;
	906
	907
	908	//
	909	// General case.
	910	//
	911	w = new double[n-1+1];
	912	v = new double[n+1];
	913	hqrnd.hqrndrandomize(ref state);
	914	for(s=2; s<=n; s++)
	915	{
	916
	917	//
	918	// Prepare random normal v
	919	//
	920	do
	921	{
	922	i = 1;
	923	while( i<=s )
	924	{
	925	hqrnd.hqrndnormal2(ref state, ref u1, ref u2);
	926	v[i] = u1;
	927	if( i+1<=s )
	928	{
	929	v[i+1] = u2;
	930	}
	931	i = i+2;
	932	}
	933	lambda = 0.0;
	934	for(i_=1; i_<=s;i_++)
	935	{
	936	lambda += v[i_]*v[i_];
	937	}
	938	}
	939	while( (double)(lambda)==(double)(0) );
	940
	941	//
	942	// Prepare and apply reflection
	943	//
	944	reflections.generatereflection(ref v, s, ref tau);
	945	v[1] = 1;
	946	reflections.applyreflectionfromtheright(ref a, tau, ref v, 0, n-1, n-s, n-1, ref w);
	947	reflections.applyreflectionfromtheleft(ref a, tau, ref v, n-s, n-1, 0, n-1, ref w);
	948	}
	949
	950	//
	951	// Second pass.
	952	//
	953	for(i=0; i<=n-1; i++)
	954	{
	955	tau = 2*AP.Math.RandomInteger(2)-1;
	956	for(i_=0; i_<=n-1;i_++)
	957	{
	958	a[i_,i] = tau*a[i_,i];
	959	}
	960	for(i_=0; i_<=n-1;i_++)
	961	{
	962	a[i,i_] = tau*a[i,i_];
	963	}
	964	}
	965	}
	966
	967
	968	/*************************************************************************
	969	Hermitian multiplication of NxN matrix by random Haar distributed
	970	complex orthogonal matrix
	971
	972	INPUT PARAMETERS:
	973	A - matrix, array[0..N-1, 0..N-1]
	974	N - matrix size
	975
	976	OUTPUT PARAMETERS:
	977	A - Q^HAQ, where Q is random NxN orthogonal matrix
	978
	979	-- ALGLIB routine --
	980	04.12.2009
	981	Bochkanov Sergey
	982	*************************************************************************/
	983	public static void hmatrixrndmultiply(ref AP.Complex[,] a,
	984	int n)
	985	{
	986	AP.Complex tau = 0;
	987	AP.Complex lambda = 0;
	988	int s = 0;
	989	int i = 0;
	990	AP.Complex[] w = new AP.Complex[0];
	991	AP.Complex[] v = new AP.Complex[0];
	992	hqrnd.hqrndstate state = new hqrnd.hqrndstate();
	993	int i_ = 0;
	994
	995
	996	//
	997	// General case.
	998	//
	999	w = new AP.Complex[n-1+1];
	1000	v = new AP.Complex[n+1];
	1001	hqrnd.hqrndrandomize(ref state);
	1002	for(s=2; s<=n; s++)
	1003	{
	1004
	1005	//
	1006	// Prepare random normal v
	1007	//
	1008	do
	1009	{
	1010	for(i=1; i<=s; i++)
	1011	{
	1012	hqrnd.hqrndnormal2(ref state, ref tau.x, ref tau.y);
	1013	v[i] = tau;
	1014	}
	1015	lambda = 0.0;
	1016	for(i_=1; i_<=s;i_++)
	1017	{
	1018	lambda += v[i_]*AP.Math.Conj(v[i_]);
	1019	}
	1020	}
	1021	while( lambda==0 );
	1022
	1023	//
	1024	// Prepare and apply reflection
	1025	//
	1026	creflections.complexgeneratereflection(ref v, s, ref tau);
	1027	v[1] = 1;
	1028	creflections.complexapplyreflectionfromtheright(ref a, tau, ref v, 0, n-1, n-s, n-1, ref w);
	1029	creflections.complexapplyreflectionfromtheleft(ref a, AP.Math.Conj(tau), ref v, n-s, n-1, 0, n-1, ref w);
	1030	}
	1031
	1032	//
	1033	// Second pass.
	1034	//
	1035	for(i=0; i<=n-1; i++)
	1036	{
	1037	hqrnd.hqrndunit2(ref state, ref tau.x, ref tau.y);
	1038	for(i_=0; i_<=n-1;i_++)
	1039	{
	1040	a[i_,i] = tau*a[i_,i];
	1041	}
	1042	tau = AP.Math.Conj(tau);
	1043	for(i_=0; i_<=n-1;i_++)
	1044	{
	1045	a[i,i_] = tau*a[i,i_];
	1046	}
	1047	}
	1048	}
	1049	}
	1050	}

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

Download in other formats:

Update cookies preferences