Free cookie consent management tool by TermsFeed Policy Generator

source: branches/3.2/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.1.2.2591/ALGLIB-2.1.2.2591/matgen.cs @ 4539

Last change on this file since 4539 was 2645, checked in by mkommend, 15 years ago

extracted external libraries and adapted dependent plugins (ticket #837)

File size: 32.2 KB
Line 
1/*************************************************************************
2Copyright (c) 2007, Sergey Bochkanov (ALGLIB project).
3
4>>> SOURCE LICENSE >>>
5This program is free software; you can redistribute it and/or modify
6it under the terms of the GNU General Public License as published by
7the Free Software Foundation (www.fsf.org); either version 2 of the
8License, or (at your option) any later version.
9
10This program is distributed in the hope that it will be useful,
11but WITHOUT ANY WARRANTY; without even the implied warranty of
12MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13GNU General Public License for more details.
14
15A copy of the GNU General Public License is available at
16http://www.fsf.org/licensing/licenses
17
18>>> END OF LICENSE >>>
19*************************************************************************/
20
21using System;
22
23namespace 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] = (2*AP.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] = (2*AP.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'*A*Q, 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^H*A*Q, 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.