Context Navigation

source: trunk/sources/ALGLIB/matgen.cs @ 2636

Visit:

Last change on this file since 2636 was 2563, checked in by gkronber, 15 years ago
Updated ALGLIB to latest version. #751 (Plugin for for data-modeling with ANN (integrated into CEDMA))
File size: 32.2 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
407	/*************************************************************************
408	Generation of random NxN Hermitian positive definite matrix with given
409	condition number and norm2(A)=1
410
411	INPUT PARAMETERS:
412	N - matrix size
413	C - condition number (in 2-norm)
414
415	OUTPUT PARAMETERS:
416	A - random HPD matrix with norm2(A)=1 and cond(A)=C
417
418	-- ALGLIB routine --
419	04.12.2009
420	Bochkanov Sergey
421	*************************************************************************/
422	public static void hpdmatrixrndcond(int n,
423	double c,
424	ref AP.Complex[,] a)
425	{
426	int i = 0;
427	int j = 0;
428	double l1 = 0;
429	double l2 = 0;
430
431
432	//
433	// Special cases
434	//
435	if( n<=0 \| (double)(c)<(double)(1) )
436	{
437	return;
438	}
439	a = new AP.Complex[n-1+1, n-1+1];
440	if( n==1 )
441	{
442	a[0,0] = 1;
443	return;
444	}
445
446	//
447	// Prepare matrix
448	//
449	l1 = 0;
450	l2 = Math.Log(1/c);
451	for(i=0; i<=n-1; i++)
452	{
453	for(j=0; j<=n-1; j++)
454	{
455	a[i,j] = 0;
456	}
457	}
458	a[0,0] = Math.Exp(l1);
459	for(i=1; i<=n-2; i++)
460	{
461	a[i,i] = Math.Exp(AP.Math.RandomReal()*(l2-l1)+l1);
462	}
463	a[n-1,n-1] = Math.Exp(l2);
464
465	//
466	// Multiply
467	//
468	hmatrixrndmultiply(ref a, n);
469	}
470
471
472	/*************************************************************************
473	Multiplication of MxN matrix by NxN random Haar distributed orthogonal matrix
474
475	INPUT PARAMETERS:
476	A - matrix, array[0..M-1, 0..N-1]
477	M, N- matrix size
478
479	OUTPUT PARAMETERS:
480	A - A*Q, where Q is random NxN orthogonal matrix
481
482	-- ALGLIB routine --
483	04.12.2009
484	Bochkanov Sergey
485	*************************************************************************/
486	public static void rmatrixrndorthogonalfromtheright(ref double[,] a,
487	int m,
488	int n)
489	{
490	double tau = 0;
491	double lambda = 0;
492	int s = 0;
493	int i = 0;
494	int j = 0;
495	double u1 = 0;
496	double u2 = 0;
497	double[] w = new double[0];
498	double[] v = new double[0];
499	double sm = 0;
500	hqrnd.hqrndstate state = new hqrnd.hqrndstate();
501	int i_ = 0;
502
503	System.Diagnostics.Debug.Assert(n>=1 & m>=1, "RMatrixRndOrthogonalFromTheRight: N<1 or M<1!");
504	if( n==1 )
505	{
506
507	//
508	// Special case
509	//
510	tau = 2*AP.Math.RandomInteger(2)-1;
511	for(i=0; i<=m-1; i++)
512	{
513	a[i,0] = a[i,0]*tau;
514	}
515	return;
516	}
517
518	//
519	// General case.
520	// First pass.
521	//
522	w = new double[m-1+1];
523	v = new double[n+1];
524	hqrnd.hqrndrandomize(ref state);
525	for(s=2; s<=n; s++)
526	{
527
528	//
529	// Prepare random normal v
530	//
531	do
532	{
533	i = 1;
534	while( i<=s )
535	{
536	hqrnd.hqrndnormal2(ref state, ref u1, ref u2);
537	v[i] = u1;
538	if( i+1<=s )
539	{
540	v[i+1] = u2;
541	}
542	i = i+2;
543	}
544	lambda = 0.0;
545	for(i_=1; i_<=s;i_++)
546	{
547	lambda += v[i_]*v[i_];
548	}
549	}
550	while( (double)(lambda)==(double)(0) );
551
552	//
553	// Prepare and apply reflection
554	//
555	reflections.generatereflection(ref v, s, ref tau);
556	v[1] = 1;
557	reflections.applyreflectionfromtheright(ref a, tau, ref v, 0, m-1, n-s, n-1, ref w);
558	}
559
560	//
561	// Second pass.
562	//
563	for(i=0; i<=n-1; i++)
564	{
565	tau = 2*AP.Math.RandomInteger(2)-1;
566	for(i_=0; i_<=m-1;i_++)
567	{
568	a[i_,i] = tau*a[i_,i];
569	}
570	}
571	}
572
573
574	/*************************************************************************
575	Multiplication of MxN matrix by MxM random Haar distributed orthogonal matrix
576
577	INPUT PARAMETERS:
578	A - matrix, array[0..M-1, 0..N-1]
579	M, N- matrix size
580
581	OUTPUT PARAMETERS:
582	A - Q*A, where Q is random MxM orthogonal matrix
583
584	-- ALGLIB routine --
585	04.12.2009
586	Bochkanov Sergey
587	*************************************************************************/
588	public static void rmatrixrndorthogonalfromtheleft(ref double[,] a,
589	int m,
590	int n)
591	{
592	double tau = 0;
593	double lambda = 0;
594	int s = 0;
595	int i = 0;
596	int j = 0;
597	double u1 = 0;
598	double u2 = 0;
599	double[] w = new double[0];
600	double[] v = new double[0];
601	double sm = 0;
602	hqrnd.hqrndstate state = new hqrnd.hqrndstate();
603	int i_ = 0;
604
605	System.Diagnostics.Debug.Assert(n>=1 & m>=1, "RMatrixRndOrthogonalFromTheRight: N<1 or M<1!");
606	if( m==1 )
607	{
608
609	//
610	// special case
611	//
612	tau = 2*AP.Math.RandomInteger(2)-1;
613	for(j=0; j<=n-1; j++)
614	{
615	a[0,j] = a[0,j]*tau;
616	}
617	return;
618	}
619
620	//
621	// General case.
622	// First pass.
623	//
624	w = new double[n-1+1];
625	v = new double[m+1];
626	hqrnd.hqrndrandomize(ref state);
627	for(s=2; s<=m; s++)
628	{
629
630	//
631	// Prepare random normal v
632	//
633	do
634	{
635	i = 1;
636	while( i<=s )
637	{
638	hqrnd.hqrndnormal2(ref state, ref u1, ref u2);
639	v[i] = u1;
640	if( i+1<=s )
641	{
642	v[i+1] = u2;
643	}
644	i = i+2;
645	}
646	lambda = 0.0;
647	for(i_=1; i_<=s;i_++)
648	{
649	lambda += v[i_]*v[i_];
650	}
651	}
652	while( (double)(lambda)==(double)(0) );
653
654	//
655	// Prepare and apply reflection
656	//
657	reflections.generatereflection(ref v, s, ref tau);
658	v[1] = 1;
659	reflections.applyreflectionfromtheleft(ref a, tau, ref v, m-s, m-1, 0, n-1, ref w);
660	}
661
662	//
663	// Second pass.
664	//
665	for(i=0; i<=m-1; i++)
666	{
667	tau = 2*AP.Math.RandomInteger(2)-1;
668	for(i_=0; i_<=n-1;i_++)
669	{
670	a[i,i_] = tau*a[i,i_];
671	}
672	}
673	}
674
675
676	/*************************************************************************
677	Multiplication of MxN complex matrix by NxN random Haar distributed
678	complex orthogonal matrix
679
680	INPUT PARAMETERS:
681	A - matrix, array[0..M-1, 0..N-1]
682	M, N- matrix size
683
684	OUTPUT PARAMETERS:
685	A - A*Q, where Q is random NxN orthogonal matrix
686
687	-- ALGLIB routine --
688	04.12.2009
689	Bochkanov Sergey
690	*************************************************************************/
691	public static void cmatrixrndorthogonalfromtheright(ref AP.Complex[,] a,
692	int m,
693	int n)
694	{
695	AP.Complex lambda = 0;
696	AP.Complex tau = 0;
697	int s = 0;
698	int i = 0;
699	int j = 0;
700	double u1 = 0;
701	double u2 = 0;
702	AP.Complex[] w = new AP.Complex[0];
703	AP.Complex[] v = new AP.Complex[0];
704	double sm = 0;
705	hqrnd.hqrndstate state = new hqrnd.hqrndstate();
706	int i_ = 0;
707
708	System.Diagnostics.Debug.Assert(n>=1 & m>=1, "CMatrixRndOrthogonalFromTheRight: N<1 or M<1!");
709	if( n==1 )
710	{
711
712	//
713	// Special case
714	//
715	hqrnd.hqrndrandomize(ref state);
716	hqrnd.hqrndunit2(ref state, ref tau.x, ref tau.y);
717	for(i=0; i<=m-1; i++)
718	{
719	a[i,0] = a[i,0]*tau;
720	}
721	return;
722	}
723
724	//
725	// General case.
726	// First pass.
727	//
728	w = new AP.Complex[m-1+1];
729	v = new AP.Complex[n+1];
730	hqrnd.hqrndrandomize(ref state);
731	for(s=2; s<=n; s++)
732	{
733
734	//
735	// Prepare random normal v
736	//
737	do
738	{
739	for(i=1; i<=s; i++)
740	{
741	hqrnd.hqrndnormal2(ref state, ref tau.x, ref tau.y);
742	v[i] = tau;
743	}
744	lambda = 0.0;
745	for(i_=1; i_<=s;i_++)
746	{
747	lambda += v[i_]*AP.Math.Conj(v[i_]);
748	}
749	}
750	while( lambda==0 );
751
752	//
753	// Prepare and apply reflection
754	//
755	creflections.complexgeneratereflection(ref v, s, ref tau);
756	v[1] = 1;
757	creflections.complexapplyreflectionfromtheright(ref a, tau, ref v, 0, m-1, n-s, n-1, ref w);
758	}
759
760	//
761	// Second pass.
762	//
763	for(i=0; i<=n-1; i++)
764	{
765	hqrnd.hqrndunit2(ref state, ref tau.x, ref tau.y);
766	for(i_=0; i_<=m-1;i_++)
767	{
768	a[i_,i] = tau*a[i_,i];
769	}
770	}
771	}
772
773
774	/*************************************************************************
775	Multiplication of MxN complex matrix by MxM random Haar distributed
776	complex orthogonal matrix
777
778	INPUT PARAMETERS:
779	A - matrix, array[0..M-1, 0..N-1]
780	M, N- matrix size
781
782	OUTPUT PARAMETERS:
783	A - Q*A, where Q is random MxM orthogonal matrix
784
785	-- ALGLIB routine --
786	04.12.2009
787	Bochkanov Sergey
788	*************************************************************************/
789	public static void cmatrixrndorthogonalfromtheleft(ref AP.Complex[,] a,
790	int m,
791	int n)
792	{
793	AP.Complex tau = 0;
794	AP.Complex lambda = 0;
795	int s = 0;
796	int i = 0;
797	int j = 0;
798	double u1 = 0;
799	double u2 = 0;
800	AP.Complex[] w = new AP.Complex[0];
801	AP.Complex[] v = new AP.Complex[0];
802	double sm = 0;
803	hqrnd.hqrndstate state = new hqrnd.hqrndstate();
804	int i_ = 0;
805
806	System.Diagnostics.Debug.Assert(n>=1 & m>=1, "CMatrixRndOrthogonalFromTheRight: N<1 or M<1!");
807	if( m==1 )
808	{
809
810	//
811	// special case
812	//
813	hqrnd.hqrndrandomize(ref state);
814	hqrnd.hqrndunit2(ref state, ref tau.x, ref tau.y);
815	for(j=0; j<=n-1; j++)
816	{
817	a[0,j] = a[0,j]*tau;
818	}
819	return;
820	}
821
822	//
823	// General case.
824	// First pass.
825	//
826	w = new AP.Complex[n-1+1];
827	v = new AP.Complex[m+1];
828	hqrnd.hqrndrandomize(ref state);
829	for(s=2; s<=m; s++)
830	{
831
832	//
833	// Prepare random normal v
834	//
835	do
836	{
837	for(i=1; i<=s; i++)
838	{
839	hqrnd.hqrndnormal2(ref state, ref tau.x, ref tau.y);
840	v[i] = tau;
841	}
842	lambda = 0.0;
843	for(i_=1; i_<=s;i_++)
844	{
845	lambda += v[i_]*AP.Math.Conj(v[i_]);
846	}
847	}
848	while( lambda==0 );
849
850	//
851	// Prepare and apply reflection
852	//
853	creflections.complexgeneratereflection(ref v, s, ref tau);
854	v[1] = 1;
855	creflections.complexapplyreflectionfromtheleft(ref a, tau, ref v, m-s, m-1, 0, n-1, ref w);
856	}
857
858	//
859	// Second pass.
860	//
861	for(i=0; i<=m-1; i++)
862	{
863	hqrnd.hqrndunit2(ref state, ref tau.x, ref tau.y);
864	for(i_=0; i_<=n-1;i_++)
865	{
866	a[i,i_] = tau*a[i,i_];
867	}
868	}
869	}
870
871
872	/*************************************************************************
873	Symmetric multiplication of NxN matrix by random Haar distributed
874	orthogonal matrix
875
876	INPUT PARAMETERS:
877	A - matrix, array[0..N-1, 0..N-1]
878	N - matrix size
879
880	OUTPUT PARAMETERS:
881	A - Q'AQ, where Q is random NxN orthogonal matrix
882
883	-- ALGLIB routine --
884	04.12.2009
885	Bochkanov Sergey
886	*************************************************************************/
887	public static void smatrixrndmultiply(ref double[,] a,
888	int n)
889	{
890	double tau = 0;
891	double lambda = 0;
892	int s = 0;
893	int i = 0;
894	int j = 0;
895	double u1 = 0;
896	double u2 = 0;
897	double[] w = new double[0];
898	double[] v = new double[0];
899	double sm = 0;
900	hqrnd.hqrndstate state = new hqrnd.hqrndstate();
901	int i_ = 0;
902
903
904	//
905	// General case.
906	//
907	w = new double[n-1+1];
908	v = new double[n+1];
909	hqrnd.hqrndrandomize(ref state);
910	for(s=2; s<=n; s++)
911	{
912
913	//
914	// Prepare random normal v
915	//
916	do
917	{
918	i = 1;
919	while( i<=s )
920	{
921	hqrnd.hqrndnormal2(ref state, ref u1, ref u2);
922	v[i] = u1;
923	if( i+1<=s )
924	{
925	v[i+1] = u2;
926	}
927	i = i+2;
928	}
929	lambda = 0.0;
930	for(i_=1; i_<=s;i_++)
931	{
932	lambda += v[i_]*v[i_];
933	}
934	}
935	while( (double)(lambda)==(double)(0) );
936
937	//
938	// Prepare and apply reflection
939	//
940	reflections.generatereflection(ref v, s, ref tau);
941	v[1] = 1;
942	reflections.applyreflectionfromtheright(ref a, tau, ref v, 0, n-1, n-s, n-1, ref w);
943	reflections.applyreflectionfromtheleft(ref a, tau, ref v, n-s, n-1, 0, n-1, ref w);
944	}
945
946	//
947	// Second pass.
948	//
949	for(i=0; i<=n-1; i++)
950	{
951	tau = 2*AP.Math.RandomInteger(2)-1;
952	for(i_=0; i_<=n-1;i_++)
953	{
954	a[i_,i] = tau*a[i_,i];
955	}
956	for(i_=0; i_<=n-1;i_++)
957	{
958	a[i,i_] = tau*a[i,i_];
959	}
960	}
961	}
962
963
964	/*************************************************************************
965	Hermitian multiplication of NxN matrix by random Haar distributed
966	complex orthogonal matrix
967
968	INPUT PARAMETERS:
969	A - matrix, array[0..N-1, 0..N-1]
970	N - matrix size
971
972	OUTPUT PARAMETERS:
973	A - Q^HAQ, where Q is random NxN orthogonal matrix
974
975	-- ALGLIB routine --
976	04.12.2009
977	Bochkanov Sergey
978	*************************************************************************/
979	public static void hmatrixrndmultiply(ref AP.Complex[,] a,
980	int n)
981	{
982	AP.Complex tau = 0;
983	AP.Complex lambda = 0;
984	int s = 0;
985	int i = 0;
986	int j = 0;
987	double u1 = 0;
988	double u2 = 0;
989	AP.Complex[] w = new AP.Complex[0];
990	AP.Complex[] v = new AP.Complex[0];
991	double sm = 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