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] = (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 | }
|
---|