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

Line
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